HΦ  3.2.0
Functions
CalcSpectrumByTPQ.c File Reference

Calculate spectrum function for the TPQ state.
Note: This method is trial and cannot be used in the release mode. More...

#include "CalcSpectrumByTPQ.h"
#include "Lanczos_EigenValue.h"
#include "FileIO.h"
#include "wrapperMPI.h"
#include "vec12.h"
#include "common/setmemory.h"
+ Include dependency graph for CalcSpectrumByTPQ.c:

Go to the source code of this file.

Functions

int ReadTPQData (struct EDMainCalStruct *X, double *ene, double *temp, double *specificHeat)
 Read TPQ data at "X->Bind.Large.itr" step in SS_rand file. More...
 
int GetCalcSpectrumTPQ (double complex dcomega, double dtemp, double dspecificheat, double ene, double *tmp_E, int Nsite, int idim_max, double complex *dc_tmpSpec)
 Calculate spectrum function from the TPQ state. More...
 
int CalcSpectrumByTPQ (struct EDMainCalStruct *X, double complex *tmp_v1, double dnorm, int Nomega, double complex *dcSpectrum, double complex *dcomega)
 A main function to calculate spectrum by TPQ (Note: This method is trial) More...
 

Detailed Description

Calculate spectrum function for the TPQ state.
Note: This method is trial and cannot be used in the release mode.

Version
1.2
Author
Kazuyoshi Yoshimi (The University of Tokyo)

Definition in file CalcSpectrumByTPQ.c.

Function Documentation

◆ CalcSpectrumByTPQ()

int CalcSpectrumByTPQ ( struct EDMainCalStruct *  X,
double complex *  tmp_v1,
double  dnorm,
int  Nomega,
double complex *  dcSpectrum,
double complex *  dcomega 
)

A main function to calculate spectrum by TPQ (Note: This method is trial)

Parameters
X[in,out] CalcStruct list for getting and pushing calculation information
tmp_v1[in] Normalized excited state.
dnorm[in] Norm of the excited state before normalization.
Nomega[in] Total number of frequencies.
dcSpectrum[out] Calculated spectrum.
dcomega[in] Target frequencies.
Return values
0normally finished
-1unnormally finished
Version
1.2
Author
Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 130 of file CalcSpectrumByTPQ.c.

References alpha, beta, EDMainCalStruct::Bind, c_CalcSpectrumFromTridiagonalEnd, c_CalcSpectrumFromTridiagonalStart, c_GetTridiagonalEnd, c_GetTridiagonalStart, c_InputSpectrumRecalcvecEnd, c_InputSpectrumRecalcvecStart, c_OutputSpectrumRecalcvecEnd, c_OutputSpectrumRecalcvecStart, DefineList::CDataFileHead, cFileNameOutputRestartVec, cFileNameTimeKeep, BindStruct::Check, childfopenALL(), D_FileNameMax, BindStruct::Def, exitMPI(), FALSE, GetCalcSpectrumTPQ(), CheckList::idim_max, DefineList::iFlgCalcSpec, Lanczos_GetTridiagonalMatrixComponents(), DefineList::Lanczos_max, DefineList::Lanczos_restart, myrank, DefineList::Nsite, DefineList::nvec, OutputTMComponents(), ReadTMComponents(), ReadTPQData(), stdoutMPI, TimeKeeper(), TRUE, v0, v1, and vec12().

Referenced by CalcSpectrum().

138 {
139  char sdt[D_FileNameMax];
140  unsigned long int i, i_max;
141  FILE *fp;
142  int iret;
143  unsigned long int liLanczosStp = X->Bind.Def.Lanczos_max;
144  unsigned long int liLanczosStp_vec=0;
145  double dene, dtemp, dspecificHeat;
146  double *tmp_E;
147  double complex dctmp_Spectrum;
148  int stp;
149  size_t byte_size;
150 
151  //Read Ene, temp, C
152  if(ReadTPQData(X, &dene, &dtemp, &dspecificHeat)!=TRUE){
153  return FALSE;
154  }
155 
156  if(X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC ||
157  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC) {
158  fprintf(stdoutMPI, " Start: Input vectors for recalculation.\n");
159  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_InputSpectrumRecalcvecStart, "a");
160 
161  sprintf(sdt, cFileNameOutputRestartVec, X->Bind.Def.CDataFileHead, myrank);
162  if (childfopenALL(sdt, "rb", &fp) != 0) {
163  exitMPI(-1);
164  }
165  byte_size = fread(&liLanczosStp_vec, sizeof(liLanczosStp_vec),1,fp);
166  byte_size = fread(&i_max, sizeof(long int), 1, fp);
167  if(i_max != X->Bind.Check.idim_max){
168  fprintf(stderr, "Error: A size of Inputvector is incorrect.\n");
169  exitMPI(-1);
170  }
171  byte_size = fread(v0, sizeof(complex double), X->Bind.Check.idim_max + 1, fp);
172  byte_size = fread(v1, sizeof(complex double), X->Bind.Check.idim_max + 1, fp);
173  fclose(fp);
174  fprintf(stdoutMPI, " End: Input vectors for recalculation.\n");
175  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_InputSpectrumRecalcvecEnd, "a");
176  if (byte_size == 0) printf("byte_size: %d \n", (int)byte_size);
177  }
178 
179  //Read diagonal components
180  if(X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents ||
181  X->Bind.Def.iFlgCalcSpec ==RECALC_FROM_TMComponents_VEC||
182  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC)
183  {
184  int iFlgTMComp=0;
185  if(X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC ||
186  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC)
187  {
188  iFlgTMComp=0;
189  }
190  else{
191  iFlgTMComp=1;
192  }
193  iret=ReadTMComponents(&(X->Bind), &dnorm, &liLanczosStp, iFlgTMComp);
194  if(iret !=TRUE){
195  fprintf(stdoutMPI, " Error: Fail to read TMcomponents\n");
196  return FALSE;
197  }
198 
199  if(X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents){
200  X->Bind.Def.Lanczos_restart=0;
201  }
202  else if(X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC||
203  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC){
204  if(liLanczosStp_vec !=liLanczosStp){
205  fprintf(stdoutMPI, " Error: Input files for vector and TMcomponents are incoorect.\n");
206  fprintf(stdoutMPI, " Error: Input vector %ld th stps, TMcomponents %ld th stps.\n", liLanczosStp_vec, liLanczosStp);
207  return FALSE;
208  }
209  X->Bind.Def.Lanczos_restart=liLanczosStp;
210  liLanczosStp = liLanczosStp+X->Bind.Def.Lanczos_max;
211  }
212  }
213 
214  // calculate ai, bi
215  if (X->Bind.Def.iFlgCalcSpec == RECALC_NOT ||
216  X->Bind.Def.iFlgCalcSpec == RECALC_OUTPUT_TMComponents_VEC ||
217  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC ||
218  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC
219  )
220  {
221  fprintf(stdoutMPI, " Start: Calculate tridiagonal matrix components.\n");
222  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_GetTridiagonalStart, "a");
223  // Functions in Lanczos_EigenValue
224  iret = Lanczos_GetTridiagonalMatrixComponents(&(X->Bind), alpha, beta, tmp_v1, &(liLanczosStp));
225  if (iret != TRUE) {
226  //Error Message will be added.
227  return FALSE;
228  }
229  fprintf(stdoutMPI, " End: Calculate tridiagonal matrix components.\n\n");
230  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_GetTridiagonalEnd, "a");
231  OutputTMComponents(&(X->Bind), alpha,beta, dnorm, liLanczosStp);
232  }//X->Bind.Def.iFlgCalcSpec == RECALC_NOT || RECALC_FROM_TMComponents_VEC
233 
234  stp=liLanczosStp;
235  tmp_E = d_1d_allocate(stp + 1);
236  X->Bind.Def.nvec= stp;
237  vec12(alpha,beta,stp,tmp_E, &(X->Bind));
238  fprintf(stdoutMPI, " Start: Caclulate spectrum from tridiagonal matrix components.\n");
239  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumFromTridiagonalStart, "a");
240  for( i = 0 ; i < Nomega; i++) {
241  dctmp_Spectrum=0;
242  iret = GetCalcSpectrumTPQ(dcomega[i], dtemp, dspecificHeat, dene, tmp_E, X->Bind.Def.Nsite, stp, &dctmp_Spectrum);
243  if (iret != TRUE) {
244  //ReAlloc alpha, beta and Set alpha_start and beta_start in Lanczos_EigenValue
245  return FALSE;
246  }
247  dcSpectrum[i] = dnorm * dctmp_Spectrum;
248  }
249  fprintf(stdoutMPI, " End: Caclulate spectrum from tridiagonal matrix components.\n\n");
250  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumFromTridiagonalEnd, "a");
251 
252  free_d_1d_allocate(tmp_E);
253  //output vectors for recalculation
254  if(X->Bind.Def.iFlgCalcSpec==RECALC_OUTPUT_TMComponents_VEC ||
255  X->Bind.Def.iFlgCalcSpec==RECALC_INOUT_TMComponents_VEC){
256  fprintf(stdoutMPI, " Start: Output vectors for recalculation.\n");
257  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_OutputSpectrumRecalcvecStart, "a");
258 
259  sprintf(sdt, cFileNameOutputRestartVec, X->Bind.Def.CDataFileHead, myrank);
260  if(childfopenALL(sdt, "wb", &fp)!=0){
261  exitMPI(-1);
262  }
263  fwrite(&liLanczosStp, sizeof(liLanczosStp),1,fp);
264  fwrite(&X->Bind.Check.idim_max, sizeof(X->Bind.Check.idim_max),1,fp);
265  fwrite(v0, sizeof(complex double),X->Bind.Check.idim_max+1, fp);
266  fwrite(v1, sizeof(complex double),X->Bind.Check.idim_max+1, fp);
267  fclose(fp);
268 
269  fprintf(stdoutMPI, " End: Output vectors for recalculation.\n");
270  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_OutputSpectrumRecalcvecEnd, "a");
271  }
272 
273  return TRUE;
274 }
childfopenALL
int childfopenALL(const char *_cPathChild, const char *_cmode, FILE **_fp)
All processes open file in output/ directory.
Definition: FileIO.c:50
exitMPI
void exitMPI(int errorcode)
MPI Abortation wrapper.
Definition: wrapperMPI.c:86
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.h:410
CheckList::idim_max
unsigned long int idim_max
The dimension of the Hilbert space of this process.
Definition: struct.h:303
v1
double complex * v1
Definition: global.h:35
ReadTPQData
int ReadTPQData(struct EDMainCalStruct *X, double *ene, double *temp, double *specificHeat)
Read TPQ data at "X->Bind.Large.itr" step in SS_rand file.
Definition: CalcSpectrumByTPQ.c:40
D_FileNameMax
#define D_FileNameMax
Definition: global.h:23
Lanczos_GetTridiagonalMatrixComponents
int Lanczos_GetTridiagonalMatrixComponents(struct BindStruct *X, double *_alpha, double *_beta, double complex *tmp_v1, unsigned long int *liLanczos_step)
Calculate tridiagonal matrix components by Lanczos method.
Definition: Lanczos_EigenValue.c:324
cFileNameOutputRestartVec
const char * cFileNameOutputRestartVec
Definition: global.c:81
DefineList::nvec
unsigned int nvec
Read from Calcmod in readdef.h.
Definition: struct.h:46
TRUE
#define TRUE
Definition: global.h:26
c_OutputSpectrumRecalcvecEnd
const char * c_OutputSpectrumRecalcvecEnd
Definition: LogMessage.c:71
DefineList::Nsite
unsigned int Nsite
Number of sites in the INTRA process region.
Definition: struct.h:56
DefineList::Lanczos_restart
int Lanczos_restart
Number of iterations performed in the restart computation.
Definition: struct.h:75
c_InputSpectrumRecalcvecStart
const char * c_InputSpectrumRecalcvecStart
Definition: LogMessage.c:72
c_GetTridiagonalStart
const char * c_GetTridiagonalStart
Definition: LogMessage.c:66
c_CalcSpectrumFromTridiagonalEnd
const char * c_CalcSpectrumFromTridiagonalEnd
Definition: LogMessage.c:69
alpha
double * alpha
Definition: global.h:44
FALSE
#define FALSE
Definition: global.h:25
GetCalcSpectrumTPQ
int GetCalcSpectrumTPQ(double complex dcomega, double dtemp, double dspecificheat, double ene, double *tmp_E, int Nsite, int idim_max, double complex *dc_tmpSpec)
Calculate spectrum function from the TPQ state.
Definition: CalcSpectrumByTPQ.c:92
cFileNameTimeKeep
const char * cFileNameTimeKeep
Definition: global.c:23
c_InputSpectrumRecalcvecEnd
const char * c_InputSpectrumRecalcvecEnd
Definition: LogMessage.c:73
beta
double * beta
Definition: global.h:44
vec12
void vec12(double alpha[], double beta[], unsigned int ndim, double tmp_E[], struct BindStruct *X)
Diagonalize a tri-diagonal matrix and store eigenvectors into vec.
Definition: vec12.c:31
OutputTMComponents
int OutputTMComponents(struct BindStruct *X, double *_alpha, double *_beta, double _dnorm, int liLanczosStp)
Output tridiagonal matrix components obtained by the Lanczos method.
Definition: Lanczos_EigenValue.c:697
myrank
int myrank
Process ID, defined in InitializeMPI()
Definition: global.h:163
DefineList::iFlgCalcSpec
int iFlgCalcSpec
Input parameter CalcSpec in teh CalcMod file.
Definition: struct.h:216
c_CalcSpectrumFromTridiagonalStart
const char * c_CalcSpectrumFromTridiagonalStart
Definition: LogMessage.c:68
c_GetTridiagonalEnd
const char * c_GetTridiagonalEnd
Definition: LogMessage.c:67
ReadTMComponents
int ReadTMComponents(struct BindStruct *X, double *_dnorm, unsigned long int *_i_max, int iFlg)
Read tridiagonal matrix components obtained by the Lanczos method. .
Definition: Lanczos_EigenValue.c:624
v0
double complex * v0
Definition: global.h:34
c_OutputSpectrumRecalcvecStart
const char * c_OutputSpectrumRecalcvecStart
Definition: LogMessage.c:70
BindStruct::Check
struct CheckList Check
Size of the Hilbert space.
Definition: struct.h:411
DefineList::CDataFileHead
char * CDataFileHead
Read from Calcmod in readdef.h. Header of output file such as Green&#39;s function.
Definition: struct.h:42
EDMainCalStruct::Bind
struct BindStruct Bind
Binded struct.
Definition: struct.h:420
TimeKeeper
int TimeKeeper(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType)
Functions for writing a time log.
Definition: log.c:42
DefineList::Lanczos_max
unsigned int Lanczos_max
Maximum number of iterations.
Definition: struct.h:74
stdoutMPI
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()
Definition: global.h:165
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ GetCalcSpectrumTPQ()

int GetCalcSpectrumTPQ ( double complex  dcomega,
double  dtemp,
double  dspecificheat,
double  ene,
double *  tmp_E,
int  Nsite,
int  idim_max,
double complex *  dc_tmpSpec 
)

Calculate spectrum function from the TPQ state.

Parameters
dcomega[in] Target frequencies.
dtemp[in] Temperature corresponding to the target TPQ state.
dspecificheat[in] Specific heat.
ene[in] Energy for the target TPQ state.
tmp_E[in] Energies included in the excited TPQ state obtained by the continued fraction expansions.
Nsite[in] Total number of sites.
idim_max[in] Dimension of the Hilbert space.
dc_tmpSpec[out] Calculated spectrum.
Return values
FALSEfail to calculate spectrum.
TRUEsucsceed to calculate spectrum.

Definition at line 92 of file CalcSpectrumByTPQ.c.

References FALSE, TRUE, and vec.

Referenced by CalcSpectrumByTPQ().

94 {
95  int l;
96  double tmp_dcSpec;
97  double factor, pre_factor;
98  pre_factor=2.0*dtemp*dtemp*dspecificheat;
99  factor=M_PI*pre_factor;
100  factor=1.0/sqrt(factor);
101  tmp_dcSpec=0;
102 
103  if(cimag(dcomega)>0) {
104  for (l = 1; l <= idim_max; l++) {
105  //TODO: Check omega is real ?
106  //fprintf(stdoutMPI, "Debug: %lf, %lf\n", creal(dcomega) - tmp_E[l] + ene, pre_factor);
107  tmp_dcSpec += (double)(vec[l][1] * conj(vec[l][1])) * exp(-pow((creal(dcomega) - tmp_E[l] + ene),2)/(pre_factor));
108  }
109  }
110  else{
111  fprintf(stderr, " an imarginary part of omega must be positive.\n");
112  return FALSE;
113  }
114  tmp_dcSpec *=factor;
115  *dc_tmpSpec = tmp_dcSpec;
116  return TRUE;
117 }
vec
double complex ** vec
Definition: global.h:45
TRUE
#define TRUE
Definition: global.h:26
FALSE
#define FALSE
Definition: global.h:25
+ Here is the caller graph for this function:

◆ ReadTPQData()

int ReadTPQData ( struct EDMainCalStruct *  X,
double *  ene,
double *  temp,
double *  specificHeat 
)

Read TPQ data at "X->Bind.Large.itr" step in SS_rand file.

Parameters
[in]XCalcStruct list for getting and pushing calculation information
[out]eneenergy
[out]temptemperature
[out]specificHeatspecific heat
Return values
TRUEsucceed to read data
FALSEfail to read data

Definition at line 40 of file CalcSpectrumByTPQ.c.

Referenced by CalcSpectrumByTPQ().

45  {
46  FILE *fp;
47  char sdt[D_FileNameMax];
48  char ctmp2[256];
49  double dinv_temp;
50  double dene, dHvar, dn, ddoublon;
51  int istp;
52  sprintf(sdt, cFileNameSSRand, X->Bind.Def.irand);
53  childfopenMPI(sdt, "r", &fp);
54  if(fp==NULL){
55  fprintf(stderr, " Error: SS_rand%d.dat does not exist.\n", X->Bind.Def.irand);
56  fclose(fp);
57  return FALSE;
58  }
59  fgetsMPI(ctmp2, 256, fp);
60  while(fgetsMPI(ctmp2, 256, fp) != NULL) {
61  sscanf(ctmp2, "%lf %lf %lf %lf %lf %d\n",
62  &dinv_temp,
63  &dene,
64  &dHvar,
65  &dn,
66  &ddoublon,
67  &istp
68  );
69  if(istp==X->Bind.Large.itr) break;
70  }
71  fclose(fp);
72 
73  *ene = dene;
74  *temp = 1.0/dinv_temp;
75  *specificHeat = (dHvar-dene*dene)*(dinv_temp*dinv_temp);
76 
77  return TRUE;
78 }
DefineList::irand
int irand
Input keyword TargetTPQRand ???
Definition: struct.h:79
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.h:410
LargeList::itr
int itr
Iteration number.
Definition: struct.h:315
D_FileNameMax
#define D_FileNameMax
Definition: global.h:23
BindStruct::Large
struct LargeList Large
Variables for Matrix-Vector product.
Definition: struct.h:412
childfopenMPI
int childfopenMPI(const char *_cPathChild, const char *_cmode, FILE **_fp)
Only the root process open file in output/ directory.
Definition: FileIO.c:27
TRUE
#define TRUE
Definition: global.h:26
FALSE
#define FALSE
Definition: global.h:25
cFileNameSSRand
const char * cFileNameSSRand
Definition: global.c:56
fgetsMPI
char * fgetsMPI(char *InputString, int maxcount, FILE *fp)
MPI file I/O (get a line, fgets) wrapper. Only the root node (myrank = 0) reads and broadcast string...
Definition: wrapperMPI.c:122
EDMainCalStruct::Bind
struct BindStruct Bind
Binded struct.
Definition: struct.h:420
+ Here is the caller graph for this function: