CalcByTPQ.c

Go to the documentation of this file.

/* HPhi  -  Quantum Lattice Model Simulator */
/* Copyright (C) 2015 The University of Tokyo */

/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */

/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
/* GNU General Public License for more details. */

/* You should have received a copy of the GNU General Public License */
/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
#include "FirstMultiply.h"
#include "Multiply.h"
#include "expec_energy_flct.h"
#include "expec_cisajs.h"
#include "expec_cisajscktaltdc.h"
#include "CalcByTPQ.h"
#include "FileIO.h"
#include "wrapperMPI.h"
#include "CalcTime.h"


int CalcByTPQ(
      const int NumAve,
      const int ExpecInterval,
      struct EDMainCalStruct *X
)
{
  char sdt[D_FileNameMax];
  char sdt_phys[D_FileNameMax];
  char sdt_norm[D_FileNameMax];
  char sdt_flct[D_FileNameMax];
  int rand_i, rand_max, iret;
  unsigned long int i_max;
  int step_iO=0;
  FILE *fp;
  double inv_temp, Ns;
  struct TimeKeepStruct tstruct;
  size_t byte_size;

  tstruct.tstart=time(NULL);
  
  rand_max = NumAve;
  step_spin = ExpecInterval;
  X->Bind.Def.St=0;
  fprintf(stdoutMPI, "%s", cLogTPQ_Start);
  for (rand_i = 0; rand_i<rand_max; rand_i++){
    sprintf(sdt_phys, cFileNameSSRand, rand_i);      
    sprintf(sdt_norm, cFileNameNormRand, rand_i);
    sprintf(sdt_flct, cFileNameFlctRand, rand_i);
    Ns = 1.0 * X->Bind.Def.NsiteMPI;
    fprintf(stdoutMPI, cLogTPQRand, rand_i+1, rand_max);
    iret=0;
    X->Bind.Def.irand=rand_i;

  //Make or Read initial vector
    if(X->Bind.Def.iReStart==RESTART_INOUT || X->Bind.Def.iReStart==RESTART_IN) {
      StartTimer(3600);
      TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cOutputVecStart, "a", rand_i, step_i);
      fprintf(stdoutMPI, "%s", cLogInputVecStart);
      sprintf(sdt, cFileNameInputVector, rand_i, myrank);
      childfopenALL(sdt, "rb", &fp);
      if(fp==NULL){
        fprintf(stdout, "A file of Inputvector does not exist.\n");
        fprintf(stdout, "Start to calculate in normal procedure.\n");
        iret=1;
      }
      byte_size = fread(&step_i, sizeof(step_i), 1, fp);
      byte_size = fread(&i_max, sizeof(long int), 1, fp);
      if(i_max != X->Bind.Check.idim_max){
        fprintf(stderr, "Error: A file of Inputvector is incorrect.\n");
        exitMPI(-1);
      }
      byte_size = fread(v0, sizeof(complex double), X->Bind.Check.idim_max+1, fp);
      TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cOutputVecFinish, "a", rand_i, step_i);
      fprintf(stdoutMPI, "%s", cLogInputVecFinish);
      fclose(fp);
      StopTimer(3600);
      X->Bind.Def.istep=step_i;
      StartTimer(3200);
      iret=expec_energy_flct(&(X->Bind));
      StopTimer(3200);
      if(iret != 0) return -1;

      step_iO=step_i-1;
      if (byte_size == 0) printf("byte_size: %d \n", (int)byte_size);
    }
    
    if(X->Bind.Def.iReStart==RESTART_NOT || X->Bind.Def.iReStart==RESTART_OUT || iret ==1) {
      StartTimer(3600);
      if (childfopenMPI(sdt_phys, "w", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%s", cLogSSRand);
      fclose(fp);
// for norm
      if (childfopenMPI(sdt_norm, "w", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%s", cLogNormRand);
      fclose(fp);
// for fluctuations
      if (childfopenMPI(sdt_flct, "w", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%s", cLogFlctRand);
      fclose(fp);

      StopTimer(3600);

      step_i = 0;

      StartTimer(3100);
      if(rand_i==0){
        TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cTPQStep, "w", rand_i, step_i);
      }
      else{
        TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cTPQStep, "a", rand_i, step_i);
      }
      FirstMultiply(rand_i, &(X->Bind));
      inv_temp = 0.0;
      StopTimer(3100);
      if (childfopenMPI(sdt_phys, "a", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%.16lf  %.16lf %.16lf %.16lf %.16lf %d\n", inv_temp, X->Bind.Phys.energy, X->Bind.Phys.var,
        X->Bind.Phys.doublon, X->Bind.Phys.num, step_i);
      fclose(fp);
      // for norm
      if (childfopenMPI(sdt_norm, "a", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%.16lf %.16lf %.16lf %d\n", inv_temp, global_1st_norm, global_1st_norm, step_i);
      fclose(fp);
      X->Bind.Def.istep = 0;
      StartTimer(3300);
      iret=expec_cisajs(&(X->Bind), v1);
      StopTimer(3300);
      if(iret !=0) return -1;

      StartTimer(3400);
      iret=expec_cisajscktaltdc(&(X->Bind), v1);
      StopTimer(3400);
      if(iret !=0) return -1;

      StartTimer(3200);
      iret=expec_energy_flct(&(X->Bind)); //v0 = H*v1
      StopTimer(3200);
      if(iret !=0) return -1;
      step_i += 1;
      inv_temp = (2.0 / Ns) / (LargeValue - X->Bind.Phys.energy / Ns);
      StartTimer(3600);
      if (childfopenMPI(sdt_phys, "a", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%.16lf  %.16lf %.16lf %.16lf %.16lf %d\n", inv_temp, X->Bind.Phys.energy, X->Bind.Phys.var,
              X->Bind.Phys.doublon, X->Bind.Phys.num, step_i);
      fclose(fp);
// for norm
      if (childfopenMPI(sdt_norm, "a", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%.16lf %.16lf %.16lf %d\n", inv_temp, global_norm, global_1st_norm, step_i);
      fclose(fp);
// for fluctuations
      if (childfopenMPI(sdt_flct, "a", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%.16lf %.16lf %.16lf %.16lf %.16lf %.16lf %.16lf %d\n", inv_temp,X->Bind.Phys.num,X->Bind.Phys.num2, X->Bind.Phys.doublon,X->Bind.Phys.doublon2, X->Bind.Phys.Sz,X->Bind.Phys.Sz2,step_i);
      fclose(fp);
//
      StopTimer(3600);
      step_i +=1;
      X->Bind.Def.istep = step_i;
      step_iO=0;
    }

    for (step_i = X->Bind.Def.istep; step_i<X->Bind.Def.Lanczos_max; step_i++){
      X->Bind.Def.istep=step_i;
      if(step_i %((X->Bind.Def.Lanczos_max-step_iO)/10)==0){
        fprintf(stdoutMPI, cLogTPQStep, step_i, X->Bind.Def.Lanczos_max);
      }
      X->Bind.Def.istep=step_i;
      StartTimer(3600);
      TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cTPQStep, "a", rand_i, step_i);
      StopTimer(3600);
      StartTimer(3500);
      Multiply(&(X->Bind));
      StopTimer(3500);

      StartTimer(3200);
      iret=expec_energy_flct(&(X->Bind));
      StopTimer(3200);
      if(iret !=0) return -1;

//
      inv_temp = (2.0*step_i / Ns) / (LargeValue - X->Bind.Phys.energy / Ns);

      StartTimer(3600);
      if(childfopenMPI(sdt_phys, "a", &fp)!=0){
        return FALSE;
      }
      fprintf(fp, "%.16lf  %.16lf %.16lf %.16lf %.16lf %d\n", inv_temp, X->Bind.Phys.energy, X->Bind.Phys.var, X->Bind.Phys.doublon, X->Bind.Phys.num ,step_i);
// for
      fclose(fp);

      if(childfopenMPI(sdt_norm, "a", &fp)!=0){
        return FALSE;
      }
      fprintf(fp, "%.16lf %.16lf %.16lf %d\n", inv_temp, global_norm, global_1st_norm, step_i);
      fclose(fp);

// for fluctuations
      if (childfopenMPI(sdt_flct, "a", &fp) != 0) {
        return -1;
      }
      fprintf(fp, "%.16lf %.16lf %.16lf %.16lf %.16lf %.16lf %.16lf %d\n", inv_temp,X->Bind.Phys.num,X->Bind.Phys.num2, X->Bind.Phys.doublon,X->Bind.Phys.doublon2, X->Bind.Phys.Sz,X->Bind.Phys.Sz2,step_i);
      fclose(fp);
//
      StopTimer(3600);


      if (step_i%step_spin == 0){
        StartTimer(3300);
        iret=expec_cisajs(&(X->Bind),v1);
        StopTimer(3300);
        if(iret !=0) return -1;

        StartTimer(3400);
        iret=expec_cisajscktaltdc(&(X->Bind), v1);
        StopTimer(3400);
        if(iret !=0) return -1;
      }
    }

    if(X->Bind.Def.iReStart== RESTART_OUT || X->Bind.Def.iReStart==RESTART_INOUT){
      TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cOutputVecStart, "a", rand_i, step_i);
      fprintf(stdoutMPI, "%s", cLogOutputVecStart);
      sprintf(sdt, cFileNameOutputVector, rand_i, myrank);
      if(childfopenALL(sdt, "wb", &fp)!=0){
        exitMPI(-1);
      }
      fwrite(&step_i, sizeof(step_i), 1, fp);
      fwrite(&X->Bind.Check.idim_max, sizeof(X->Bind.Check.idim_max),1,fp);
      fwrite(v1, sizeof(complex double),X->Bind.Check.idim_max+1, fp);
      fclose(fp);
      TimeKeeperWithRandAndStep(&(X->Bind), cFileNameTPQStep, cOutputVecFinish, "a", rand_i, step_i);
      fprintf(stdoutMPI, "%s", cLogOutputVecFinish);
    }
  }
  fprintf(stdoutMPI, "%s", cLogTPQ_End);

  tstruct.tend=time(NULL);
  fprintf(stdoutMPI, cLogTPQEnd, (int)(tstruct.tend-tstruct.tstart));
  return TRUE;
}