HPhi++  3.1.0
Functions
CalcSpectrumByBiCG.cpp File Reference

File for givinvg functions of calculating spectrum by Lanczos. More...

#include "Common.hpp"
#include "FileIO.hpp"
#include "wrapperMPI.hpp"
#include "common/setmemory.hpp"
#include "mltply.hpp"
#include "CalcSpectrum.hpp"
#include "mltplyCommon.hpp"
#include <mpi.h>
#include <iostream>

Go to the source code of this file.

Functions

void ShiftedEq (int iter, int Nomega, int NdcSpectrum, int *lz_conv, std::complex< double > *alpha, std::complex< double > *beta, std::complex< double > *dcomega, std::complex< double > z_seed, std::complex< double > **pBiCG, std::complex< double > **res_proj, std::complex< double > **pi, std::complex< double > **dcSpectrum)
 
void SeedSwitch (int iter, int Nomega, int NdcSpectrum, int *lz_conv, int *iz_seed, std::complex< double > *z_seed, std::complex< double > *rho, std::complex< double > *dcomega, long int ndim, std::complex< double > **v2, std::complex< double > **v3, std::complex< double > **v4, std::complex< double > **v5, std::complex< double > **pi, std::complex< double > *alpha, std::complex< double > *beta, std::complex< double > **res_proj)
 Perform Seed Switch. More...
 
int CalcSpectrumByBiCG (struct EDMainCalStruct *X, std::complex< double > **v2, std::complex< double > **v4, int Nomega, int NdcSpectrum, std::complex< double > **dcSpectrum, std::complex< double > *dcomega, std::complex< double > **v1Org)
 A main function to calculate spectrum by BiCG method In this function, the \(K\omega\) library is used. The detailed procedure is written in the document of \(K\omega\). https://issp-center-dev.github.io/Komega/library/en/_build/html/komega_workflow_en.html#the-schematic-workflow-of-shifted-bicg-library. More...
 

Detailed Description

File for givinvg functions of calculating spectrum by Lanczos.

Author
Mitsuaki Kawamura (The University of Tokyo)

Definition in file CalcSpectrumByBiCG.cpp.

Function Documentation

◆ CalcSpectrumByBiCG()

int CalcSpectrumByBiCG ( struct EDMainCalStruct *  X,
std::complex< double > **  v2,
std::complex< double > **  v4,
int  Nomega,
int  NdcSpectrum,
std::complex< double > **  dcSpectrum,
std::complex< double > *  dcomega,
std::complex< double > **  v1Org 
)

A main function to calculate spectrum by BiCG method In this function, the \(K\omega\) library is used. The detailed procedure is written in the document of \(K\omega\). https://issp-center-dev.github.io/Komega/library/en/_build/html/komega_workflow_en.html#the-schematic-workflow-of-shifted-bicg-library.

Return values
0normally finished
-1error
Author
Mitsuaki Kawamura (The University of Tokyo)

  • Malloc vector for old residual vector ( \({\bf r}_{\rm old}\)) and old shadow residual vector ( \({\bf {\tilde r}}_{\rm old}\)).

  • Set initial result vector(+shadow result vector) Read residual vectors if restart

  • Input \(\alpha, \beta[iter]\), projected residual, or start from scratch

  • DO BiCG loop

    • \({\bf v}_{2}={\hat H}{\bf v}_{12}, {\bf v}_{4}={\hat H}{\bf v}_{14}\), where \({\bf v}_{12}, {\bf v}_{14}\) are old (shadow) residual vector.

    • Update projected result vector dcSpectrum.

  • END DO BiCG loop

  • Save \(\alpha, \beta[iter]\), projected residual

  • output vectors for recalculation
Parameters
[in,out]X
[in]v2[CheckList::idim_max] Right hand side vector, excited state.
[in,out]v4[CheckList::idim_max] Work space for residual vector \({\bf r}\)
[in]NomegaNumber of Frequencies
[in]NdcSpectrumNumber of left side operator
[out]dcSpectrum[Nomega] Spectrum
[in]dcomega[Nomega] Frequency
[in]v1OrgUnexcited vector

Definition at line 156 of file CalcSpectrumByBiCG.cpp.

References EDMainCalStruct::Bind, DefineList::CDataFileHead, BindStruct::Check, childfopenALL(), childfopenMPI(), BindStruct::Def, eps_Lanczos, exitMPI(), fgetsMPI(), GetExcitedState(), CheckList::idim_max, DefineList::iFlgCalcSpec, DefineList::Lanczos_max, mltply(), myrank, NormMPI_dc(), SeedSwitch(), ShiftedEq(), stdoutMPI, TimeKeeper(), VecProdMPI(), and zclear().

Referenced by CalcSpectrum().

166 {
167  char sdt[D_FileNameMax];
168  char ctmp[256];
169  long int idim, i_max;
170  FILE *fp;
171  size_t byte_size;
172  int idcSpectrum;
173  std::complex<double> rho = 1.0, rho_old, z_seed, alpha_denom;
174  int iomega, iter_old;
175  int iter, iz_seed = 0, *lz_conv;
176  double resnorm, dtmp[4];
177  std::complex<double> **vL, **v12, **v14, **v3, **v5,
178  *alpha, *beta, **res_proj, **pi, **pBiCG;
179 
180  z_seed = dcomega[iz_seed];
181  fprintf(stdoutMPI, "##### Spectrum calculation with BiCG #####\n\n");
187  v12 = cd_2d_allocate(X->Bind.Check.idim_max + 1, 1);
188  v14 = cd_2d_allocate(X->Bind.Check.idim_max + 1, 1);
189  v3 = cd_2d_allocate(X->Bind.Check.idim_max + 1, 1);
190  v5 = cd_2d_allocate(X->Bind.Check.idim_max + 1, 1);
191  vL = cd_2d_allocate(X->Bind.Check.idim_max + 1, 1);
192  lz_conv = i_1d_allocate(Nomega);
197  if (X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC ||
198  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC) {
199  fprintf(stdoutMPI, " Start: Input vectors for recalculation.\n");
200  TimeKeeper(&(X->Bind), "%s_TimeKeeper.dat", "Input vectors for recalculation starts: %s", "a");
201 
202  sprintf(sdt, "%s_recalcvec_rank_%d.dat", X->Bind.Def.CDataFileHead, myrank);
203  if (childfopenALL(sdt, "rb", &fp) != 0) {
204  fprintf(stdoutMPI, "INFO: File for the restart is not found.\n");
205  fprintf(stdoutMPI, " Start from SCRATCH.\n");
206  zclear(X->Bind.Check.idim_max, &v2[1][0]);
207  GetExcitedState(&(X->Bind), 1, v2, v1Org, 0);
208 #pragma omp parallel for default(none) shared(v2,v4,v1Org,X) private(idim)
209  for (idim = 1; idim <= X->Bind.Check.idim_max; idim++)
210  v4[idim][0] = v2[idim][0];
211  }
212  else {
213  byte_size = fread(&iter_old, sizeof(int), 1, fp);
214  byte_size = fread(&i_max, sizeof(i_max), 1, fp);
215  if (i_max != X->Bind.Check.idim_max) {
216  fprintf(stderr, "Error: The size of the input vector is incorrect.\n");
217  printf("%s %ld %ld %d\n", sdt, i_max, X->Bind.Check.idim_max, iter_old);
218  exitMPI(-1);
219  }
220  byte_size = fread(&v2[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
221  byte_size = fread(&v3[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
222  byte_size = fread(&v4[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
223  byte_size = fread(&v5[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
224  fclose(fp);
225  fprintf(stdoutMPI, " End: Input vectors for recalculation.\n");
226  TimeKeeper(&(X->Bind), "%s_TimeKeeper.dat", "Input vectors for recalculation finishes: %s", "a");
227  if (byte_size == 0) printf("byte_size : %d\n", (int)byte_size);
228  }/*if (childfopenALL(sdt, "rb", &fp) == 0)*/
229  }/*if (X->Bind.Def.iFlgCalcSpec > RECALC_FROM_TMComponents)*/
230  else {
231  zclear(X->Bind.Check.idim_max, &v2[1][0]);
232  GetExcitedState(&(X->Bind), 1, v2, v1Org, 0);
233 #pragma omp parallel for default(none) shared(v2,v4,v1Org,X) private(idim)
234  for (idim = 1; idim <= X->Bind.Check.idim_max; idim++)
235  v4[idim][0] = v2[idim][0];
236  }
240  iter_old = 0;
241  fp = NULL;
242  if (X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents ||
243  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC ||
244  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC) {
245  sprintf(sdt, "%s_TMComponents.dat", X->Bind.Def.CDataFileHead);
246  if (childfopenALL(sdt, "rb", &fp) != 0) {
247  fprintf(stdoutMPI, "INFO: File for the restart is not found.\n");
248  fprintf(stdoutMPI, " Start from SCRATCH.\n");
249  }
250  else {
251  fgetsMPI(ctmp, sizeof(ctmp) / sizeof(char), fp);
252  sscanf(ctmp, "%d", &iter_old);
253  }
254  }
255  alpha = cd_1d_allocate(iter_old + X->Bind.Def.Lanczos_max + 1);
256  beta = cd_1d_allocate(iter_old + X->Bind.Def.Lanczos_max + 1);
257  res_proj = cd_2d_allocate(iter_old + X->Bind.Def.Lanczos_max + 1, NdcSpectrum);
258  pi = cd_2d_allocate(iter_old + X->Bind.Def.Lanczos_max + 1, Nomega);
259  pBiCG = cd_2d_allocate(Nomega, NdcSpectrum);
260  alpha[0] = std::complex<double>(1.0, 0.0);
261  beta[0] = std::complex<double>(0.0, 0.0);
262  for (iomega = 0; iomega < Nomega; iomega++) {
263  pi[0][iomega] = 1.0;
264  for (idcSpectrum = 0; idcSpectrum < NdcSpectrum; idcSpectrum++) {
265  pBiCG[iomega][idcSpectrum] = 0.0;
266  dcSpectrum[iomega][idcSpectrum] = 0.0;
267  }
268  }
269 
270  if (fp != NULL) {
271  if (X->Bind.Def.iFlgCalcSpec > RECALC_FROM_TMComponents)
272  X->Bind.Def.Lanczos_max = 0;
273  fgetsMPI(ctmp, sizeof(ctmp) / sizeof(char), fp);
274  sscanf(ctmp, "%lf %lf\n", &dtmp[0], &dtmp[1]);
275  z_seed = std::complex<double>(dtmp[0], dtmp[1]);
276 
277  iter = 1;
278  while (fgetsMPI(ctmp, sizeof(ctmp) / sizeof(char), fp) != NULL) {
279  sscanf(ctmp, "%lf %lf %lf %lf\n",
280  &dtmp[0], &dtmp[1], &dtmp[2], &dtmp[3]);
281  alpha[iter] = std::complex<double>(dtmp[0], dtmp[1]);
282  beta[iter] = std::complex<double>(dtmp[2], dtmp[3]);
283  for (idcSpectrum = 0; idcSpectrum < NdcSpectrum; idcSpectrum++)
284  sscanf(ctmp, "%lf %lf\n", &dtmp[0], &dtmp[1]);
285  res_proj[iter][idcSpectrum] = std::complex<double>(dtmp[0], dtmp[1]);
286  iter += 1;
287  }
288  fclose(fp);
289 
290  for (iter = 1; iter <= iter_old; iter++) {
291  ShiftedEq(iter, Nomega, NdcSpectrum, lz_conv, alpha, beta, dcomega,
292  z_seed, pBiCG, res_proj, pi, dcSpectrum);
293  }/*for (iter = 1; iter <= iter_old; iter++)*/
294 
295  rho = VecProdMPI(X->Bind.Check.idim_max, &v5[0][0], &v3[0][0]);
296 
297  SeedSwitch(iter, Nomega, NdcSpectrum, lz_conv, &iz_seed,
298  &z_seed, &rho, dcomega, X->Bind.Check.idim_max, v2, v3, v4, v5,
299  pi, alpha, beta, res_proj);
300 
301  resnorm = NormMPI_dc(X->Bind.Check.idim_max, &v2[0][0]);
302 
303  for (iomega = 0; iomega < Nomega; iomega++)
304  if (std::abs(resnorm / pi[iter][iomega]) < eps_Lanczos)
305  lz_conv[idcSpectrum] = 1;
306  }/*if (fp != NULL)*/
311  fprintf(stdoutMPI, " Start: Calculate tridiagonal matrix components.\n");
312  TimeKeeper(&(X->Bind), "%s_TimeKeeper.dat", "Calculating tridiagonal matrix components starts: %s", "a");
313  fprintf(stdoutMPI, "\n Iteration Seed Residual-2-Norm\n");
314  childfopenMPI("residual.dat", "w", &fp);
315 
316  for (iter = iter_old + 1; iter <= iter_old + X->Bind.Def.Lanczos_max; iter++) {
321  zclear(X->Bind.Check.idim_max, &v12[1][0]);
322  zclear(X->Bind.Check.idim_max, &v14[1][0]);
323  mltply(&X->Bind, 1, v12, v2);
324  mltply(&X->Bind, 1, v14, v4);
325 
326  for (idcSpectrum = 0; idcSpectrum < NdcSpectrum; idcSpectrum++) {
327  zclear(X->Bind.Check.idim_max, &vL[1][0]);
328  GetExcitedState(&(X->Bind), 1, vL, v1Org, idcSpectrum + 1);
329  res_proj[iter][idcSpectrum] = VecProdMPI(X->Bind.Check.idim_max, &vL[0][0], &v2[0][0]);
330  }
334  rho_old = rho;
335  rho = VecProdMPI(X->Bind.Check.idim_max, &v4[0][0], &v2[0][0]);
336  if (iter == 1)
337  beta[iter] = std::complex<double>(0.0, 0.0);
338  else
339  beta[iter] = rho / rho_old;
340 
341  for (idim = 1; idim <= X->Bind.Check.idim_max; idim++) {
342  v12[idim][0] = z_seed * v2[idim][0] - v12[idim][0];
343  v14[idim][0] = std::conj(z_seed) * v4[idim][0] - v14[idim][0];
344  }
345  alpha_denom = VecProdMPI(X->Bind.Check.idim_max, &v4[0][0], &v12[0][0])
346  - beta[iter] * rho / alpha[iter - 1];
347 
348  if (std::abs(alpha_denom) < 1.0e-50) {
349  printf("Error : The denominator of alpha is zero.\n");
350  exitMPI(-1);
351  }
352  else if (std::abs(rho) < 1.0e-50) {
353  printf("Error : rho is zero.\n");
354  exitMPI(-1);
355  }
356  alpha[iter] = rho / alpha_denom;
357  /*
358  Shifted equation
359  */
360  ShiftedEq(iter, Nomega, NdcSpectrum, lz_conv, alpha, beta, dcomega,
361  z_seed, pBiCG, res_proj, pi, dcSpectrum);
362  /*
363  Update residual
364  */
365  for (idim = 1; idim <= X->Bind.Check.idim_max; idim++) {
366  v12[idim][0] = (1.0 + alpha[iter] * beta[iter] / alpha[iter-1]) * v2[idim][0]
367  - alpha[iter] * v12[idim][0]
368  - alpha[iter] * beta[iter] / alpha[iter-1] * v3[idim][0];
369  v3[idim][0] = v2[idim][0];
370  v2[idim][0] = v12[idim][0];
371  v14[idim][0] = (1.0 + std::conj(alpha[iter] * beta[iter] / alpha[iter-1])) * v4[idim][0]
372  - std::conj(alpha[iter]) * v14[idim][0]
373  - std::conj(alpha[iter] * beta[iter] / alpha[iter-1]) * v5[idim][0];
374  v5[idim][0] = v4[idim][0];
375  v4[idim][0] = v14[idim][0];
376  }/*for (idim = 1; idim <= X->Bind.Check.idim_max; idim++)*/
377  /*
378  Seed Switching
379  */
380  SeedSwitch(iter, Nomega, NdcSpectrum, lz_conv, &iz_seed,
381  &z_seed, &rho, dcomega, X->Bind.Check.idim_max, v2, v3, v4, v5,
382  pi, alpha, beta, res_proj);
383  /*
384  Convergence check
385  */
386  resnorm = std::sqrt(NormMPI_dc(X->Bind.Check.idim_max, &v2[0][0]));
387 
388  for (iomega = 0; iomega < Nomega; iomega++)
389  if (std::abs(resnorm / pi[iter][iomega]) < eps_Lanczos)
390  lz_conv[idcSpectrum] = 1;
391 
392  fprintf(stdoutMPI, " %9d %9d %25.15e\n", iter, iz_seed, resnorm);
393  if (resnorm < eps_Lanczos) break;
394  }/*for (iter = 0; iter <= X->Bind.Def.Lanczos_max; iter++)*/
395 
396  if (iter >= iter_old + X->Bind.Def.Lanczos_max)
397  fprintf(stdoutMPI, "Remark : Not converged in iteration %d.", iter);
398  iter_old = iter;
403  fprintf(stdoutMPI, " End: Calculate tridiagonal matrix components.\n\n");
404  TimeKeeper(&(X->Bind), "%s_TimeKeeper.dat", "Calculating tridiagonal matrix components finishes: %s", "a");
408  if (X->Bind.Def.iFlgCalcSpec != RECALC_FROM_TMComponents) {
409  sprintf(sdt, "%s_TMComponents.dat", X->Bind.Def.CDataFileHead);
410  childfopenMPI(sdt, "w", &fp);
411  fprintf(fp, "%d \n", iter_old);
412  fprintf(fp, "%.10lf %.10lf\n", std::real(z_seed), std::imag(z_seed));
413  for (iter = 1; iter <= iter_old; iter++) {
414  fprintf(fp, "%25.16le %25.16le %25.16le %25.16le\n",
415  std::real(alpha[iter]), std::imag(alpha[iter]),
416  std::real(beta[iter]), std::imag(beta[iter]));
417  for (idcSpectrum = 0; idcSpectrum < NdcSpectrum; idcSpectrum++)
418  fprintf(fp, "%25.16le %25.16le\n",
419  std::real(res_proj[iter][idcSpectrum]),
420  std::imag(res_proj[iter][idcSpectrum]));
421  }/*for (iter = 0; iter < iter_old; iter++)*/
422  fclose(fp);
423  }
428  if (X->Bind.Def.iFlgCalcSpec == RECALC_OUTPUT_TMComponents_VEC ||
429  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC) {
430  fprintf(stdoutMPI, " Start: Output vectors for recalculation.\n");
431  TimeKeeper(&(X->Bind), "%s_TimeKeeper.dat", "Output vectors for recalculation starts: %s", "a");
432  sprintf(sdt, "%s_recalcvec_rank_%d.dat", X->Bind.Def.CDataFileHead, myrank);
433  if (childfopenALL(sdt, "wb", &fp) != 0) {
434  exitMPI(-1);
435  }
436  byte_size = fwrite(&iter, sizeof(iter), 1, fp);
437  byte_size = fwrite(&X->Bind.Check.idim_max, sizeof(X->Bind.Check.idim_max), 1, fp);
438  byte_size = fwrite(&v2[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
439  byte_size = fwrite(&v3[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
440  byte_size = fwrite(&v4[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
441  byte_size = fwrite(&v5[0][0], sizeof(std::complex<double>), X->Bind.Check.idim_max + 1, fp);
442  fclose(fp);
443 
444  fprintf(stdoutMPI, " End: Output vectors for recalculation.\n");
445  TimeKeeper(&(X->Bind), "%s_TimeKeeper.dat", "Output vectors for recalculation finishes: %s", "a");
446  }/*if (X->Bind.Def.iFlgCalcSpec > RECALC_FROM_TMComponents)*/
447 
448  free_cd_1d_allocate(alpha);
449  free_cd_1d_allocate(beta);
450  free_cd_2d_allocate(res_proj);
451  free_cd_2d_allocate(pi);
452  free_cd_2d_allocate(pBiCG);
453  free_cd_2d_allocate(v3);
454  free_cd_2d_allocate(v5);
455  free_cd_2d_allocate(v12);
456  free_cd_2d_allocate(v14);
457  free_cd_2d_allocate(vL);
458  return TRUE;
459 }/*int CalcSpectrumByBiCG*/
exitMPI
void exitMPI(int errorcode)
MPI Abortation wrapper.
Definition: wrapperMPI.cpp:86
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.hpp:395
GetExcitedState
int GetExcitedState(struct BindStruct *X, int nstate, std::complex< double > **tmp_v0, std::complex< double > **tmp_v1, int iEx)
Parent function to calculate the excited state.
Definition: CalcSpectrum.cpp:387
stdoutMPI
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()
Definition: global.cpp:75
childfopenALL
int childfopenALL(const char *_cPathChild, const char *_cmode, FILE **_fp)
All processes open file in output/ directory.
Definition: FileIO.cpp:50
NormMPI_dc
double NormMPI_dc(long int idim, std::complex< double > *_v1)
Compute norm of process-distributed vector .
Definition: wrapperMPI.cpp:321
mltply
int mltply(struct BindStruct *X, int nstate, std::complex< double > **tmp_v0, std::complex< double > **tmp_v1)
Parent function of multiplying the wavefunction by the Hamiltonian. . First, the calculation of diago...
Definition: mltply.cpp:56
VecProdMPI
std::complex< double > VecProdMPI(long int ndim, std::complex< double > *v1, std::complex< double > *v2)
Compute conjugate scaler product of process-distributed vector .
Definition: wrapperMPI.cpp:367
zclear
void zclear(long int n, std::complex< double > *x)
clear std::complex<double> array.
Definition: mltply.cpp:143
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.cpp:122
ShiftedEq
void ShiftedEq(int iter, int Nomega, int NdcSpectrum, int *lz_conv, std::complex< double > *alpha, std::complex< double > *beta, std::complex< double > *dcomega, std::complex< double > z_seed, std::complex< double > **pBiCG, std::complex< double > **res_proj, std::complex< double > **pi, std::complex< double > **dcSpectrum)
Definition: CalcSpectrumByBiCG.cpp:34
DefineList::Lanczos_max
int Lanczos_max
Maximum number of iterations.
Definition: struct.hpp:74
myrank
int myrank
Process ID, defined in InitializeMPI()
Definition: global.cpp:73
SeedSwitch
void SeedSwitch(int iter, int Nomega, int NdcSpectrum, int *lz_conv, int *iz_seed, std::complex< double > *z_seed, std::complex< double > *rho, std::complex< double > *dcomega, long int ndim, std::complex< double > **v2, std::complex< double > **v3, std::complex< double > **v4, std::complex< double > **v5, std::complex< double > **pi, std::complex< double > *alpha, std::complex< double > *beta, std::complex< double > **res_proj)
Perform Seed Switch.
Definition: CalcSpectrumByBiCG.cpp:73
TimeKeeper
int TimeKeeper(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType)
Functions for writing a time log.
Definition: log.cpp:42
DefineList::iFlgCalcSpec
int iFlgCalcSpec
Input parameter CalcSpec in teh CalcMod file.
Definition: struct.hpp:218
eps_Lanczos
double eps_Lanczos
Definition: global.cpp:65
BindStruct::Check
struct CheckList Check
Size of the Hilbert space.
Definition: struct.hpp:396
DefineList::CDataFileHead
char * CDataFileHead
Read from Calcmod in readdef.h. Header of output file such as Green&#39;s function.
Definition: struct.hpp:42
CheckList::idim_max
long int idim_max
The dimension of the Hilbert space of this process.
Definition: struct.hpp:305
EDMainCalStruct::Bind
struct BindStruct Bind
Binded struct.
Definition: struct.hpp:405
childfopenMPI
int childfopenMPI(const char *_cPathChild, const char *_cmode, FILE **_fp)
Only the root process open file in output/ directory.
Definition: FileIO.cpp:27

◆ SeedSwitch()

void SeedSwitch ( int  iter,
int  Nomega,
int  NdcSpectrum,
int *  lz_conv,
int *  iz_seed,
std::complex< double > *  z_seed,
std::complex< double > *  rho,
std::complex< double > *  dcomega,
long int  ndim,
std::complex< double > **  v2,
std::complex< double > **  v3,
std::complex< double > **  v4,
std::complex< double > **  v5,
std::complex< double > **  pi,
std::complex< double > *  alpha,
std::complex< double > *  beta,
std::complex< double > **  res_proj 
)

Perform Seed Switch.

Definition at line 73 of file CalcSpectrumByBiCG.cpp.

References exitMPI().

Referenced by CalcSpectrumByBiCG().

91  {
92  double pi_min;
93  std::complex<double> pi_seed;
94  int iz_seed0, iomega, jter, idcSpectrum;
95  long int idim;
96 
97  pi_min = std::abs(pi[iter][0]);
98  iz_seed0 = 0;
99  for (iomega = 0; iomega < Nomega; iomega++) {
100  if (lz_conv[iomega] == 0)
101  if (std::abs(pi[iter][iomega]) < pi_min) {
102  iz_seed0 = iomega;
103  pi_min = std::abs(pi[iter][iomega]);
104  }
105  }/*for (iomega = 0; iomega < Nomega; iomega++)*/
106 
107  if (std::abs(pi[iter][iz_seed0]) < 1.0e-50) {
108  printf("Error : pi at seed is 0.");
109  exitMPI(-1);
110  }
111 
112  if (iz_seed0 != *iz_seed) {
113 
114  *iz_seed = iz_seed0;
115  *z_seed = dcomega[iz_seed0];
116 
117  *rho /= std::pow(pi[iter-1][iz_seed0], 2);
118 
119  for (idim = 1; idim <= ndim; idim++) {
120  v2[idim][0] /= pi[iter][iz_seed0];
121  v4[idim][0] /= std::conj(pi[iter][iz_seed0]);
122  v3[idim][0] /= pi[iter-1][iz_seed0];
123  v5[idim][0] /= std::conj(pi[iter-1][iz_seed0]);
124  }
125  /*
126  For restarting
127  */
128  for (jter = 1; jter <= iter; jter++) {
129  alpha[jter] *= pi[jter - 1][iz_seed0] / pi[jter][iz_seed0];
130  if(jter != 1)
131  beta[jter] *= std::pow(pi[jter - 2][iz_seed0] / pi[jter - 1][iz_seed0], 2);
132  for (idcSpectrum = 0; idcSpectrum < NdcSpectrum; idcSpectrum++) {
133  res_proj[jter][idcSpectrum] /= pi[jter - 1][iz_seed0];
134  }
135  }
136 
137  for (jter = 1; jter <= iter; jter++) {
138  pi_seed = pi[jter][iz_seed0];
139  for (iomega = 0; iomega < Nomega; iomega++)
140  pi[jter][iomega] /= pi_seed;
141  }
142  }
143 
144 }/*void SeedSwitch*/
exitMPI
void exitMPI(int errorcode)
MPI Abortation wrapper.
Definition: wrapperMPI.cpp:86

◆ ShiftedEq()

void ShiftedEq ( int  iter,
int  Nomega,
int  NdcSpectrum,
int *  lz_conv,
std::complex< double > *  alpha,
std::complex< double > *  beta,
std::complex< double > *  dcomega,
std::complex< double >  z_seed,
std::complex< double > **  pBiCG,
std::complex< double > **  res_proj,
std::complex< double > **  pi,
std::complex< double > **  dcSpectrum 
)

Definition at line 34 of file CalcSpectrumByBiCG.cpp.

Referenced by CalcSpectrumByBiCG().

47  {
48  int iomega, idcSpectrum;
49  std::complex<double> pi_2;
50 
51  for (iomega = 0; iomega < Nomega; iomega++) {
52 
53  if (lz_conv[iomega] == 1) continue;
54 
55  if (iter == 1)
56  pi_2 = 1.0;
57  else
58  pi_2 = pi[iter - 2][iomega];
59 
60  pi[iter][iomega] = (1.0 + alpha[iter] * (dcomega[iomega] - z_seed)) * pi[iter - 1][iomega]
61  - alpha[iter] * beta[iter] / alpha[iter - 1] * (pi_2 - pi[iter - 1][iomega]);
62  for (idcSpectrum = 0; idcSpectrum < NdcSpectrum; idcSpectrum++) {
63  pBiCG[iomega][idcSpectrum] = res_proj[iter][idcSpectrum] / pi[iter - 1][iomega]
64  + std::pow(pi_2 / pi[iter - 1][iomega], 2) * beta[iter] * pBiCG[iomega][idcSpectrum];
65  dcSpectrum[iomega][idcSpectrum] +=
66  pi[iter-1][iomega] / pi[iter][iomega] * alpha[iter] * pBiCG[iomega][idcSpectrum];
67  }
68  }/*for (iomega = 0; iomega < Nomega; iomega++)*/
69 }