Lanczos_EigenValue.c File Reference

Calculate eigenvalues by the Lanczos method. More...

#include "Common.h"
#include "mfmemory.h"
#include "mltply.h"
#include "vec12.h"
#include "bisec.h"
#include "FileIO.h"
#include "matrixlapack.h"
#include "Lanczos_EigenValue.h"
#include "wrapperMPI.h"
#include "CalcTime.h"

Go to the source code of this file.

Functions
int	Lanczos_EigenValue (struct BindStruct *X)
	Main function for calculating eigen values by Lanczos method. The energy convergence is judged by the level of target energy determined by \( \verb|k_exct| \). . More...
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. More...
int	ReadInitialVector (struct BindStruct *X, double complex *_v0, double complex *_v1, unsigned long int *liLanczosStp_vec)
	Read initial vectors for the restart calculation. More...
int	OutputLanczosVector (struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1, unsigned long int liLanczosStp_vec)
	Output initial vectors for the restart calculation. More...
void	SetInitialVector (struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
	Set initial vector to start the calculation for Lanczos method. . More...
int	ReadTMComponents (struct BindStruct *X, double *_dnorm, unsigned long int *_i_max, int iFlg)
	Read tridiagonal matrix components obtained by the Lanczos method. . More...
int	OutputTMComponents (struct BindStruct *X, double *_alpha, double *_beta, double _dnorm, int liLanczosStp)
	Output tridiagonal matrix components obtained by the Lanczos method. More...

Detailed Description

Calculate eigenvalues by the Lanczos method.

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition in file Lanczos_EigenValue.c.

Function Documentation

Lanczos_EigenValue()

int Lanczos_EigenValue

(

struct BindStruct *

X

)

Main function for calculating eigen values by Lanczos method.
The energy convergence is judged by the level of target energy determined by \( \verb|k_exct| \).
.

Parameters

X	[in] Struct to give the information for calculating the eigen values.

Return values

-2	Fail to read the initial vectors or triangular matrix components.
-1	Fail to obtain the eigen values with in the \( \verb| Lanczos_max |\) step
0	Succeed to calculate the eigen values.

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 50 of file Lanczos_EigenValue.c.

References alpha, beta, cFileNameLanczosStep, cFileNameTimeKeep, childfopenMPI(), cLanczos_EigenValueFinish, cLanczos_EigenValueStart, cLanczos_EigenValueStep, cLogLanczos_EigenValueEnd, cLogLanczos_EigenValueNotConverged, cLogLanczos_EigenValueStart, D_FileNameMax, DSEVvalue(), eps_Lanczos, mfint, mltply(), OutputLanczosVector(), OutputTMComponents(), ReadInitialVector(), ReadTMComponents(), SetInitialVector(), StartTimer(), stdoutMPI, StopTimer(), SumMPI_dc(), SumMPI_li(), TimeKeeper(), TimeKeeperWithStep(), TRUE, v0, v1, vec12(), and X.

Referenced by CalcByLanczos().

                                             {

  fprintf(stdoutMPI, "%s", cLogLanczos_EigenValueStart);
  FILE *fp;
  char sdt[D_FileNameMax], sdt_2[D_FileNameMax];
  int stp;
  long int i, i_max;
  unsigned long int i_max_tmp;
  int k_exct, Target;
  int iconv = -1;
  double beta1, alpha1; //beta,alpha1 should be real
  double complex temp1, temp2;
  double complex cbeta1;
  double E[5], ebefor, E_target;

// for GC
  double dnorm=0.0;

  double **tmp_mat;
  double *tmp_E;
  int int_i, int_j, mfint[7];
  int iret=0;
  sprintf(sdt_2, cFileNameLanczosStep, X->Def.CDataFileHead);

  i_max = X->Check.idim_max;
  k_exct = X->Def.k_exct;
  unsigned long int liLanczosStp;
  liLanczosStp = X->Def.Lanczos_max;
  unsigned long int liLanczosStp_vec=0;

    if (X->Def.iReStart == RESTART_INOUT || X->Def.iReStart == RESTART_IN){
      if(ReadInitialVector( X, v0, v1, &liLanczosStp_vec)!=0){
        fprintf(stdoutMPI, "  Error: Fail to read initial vectors\n");
        return -2;
      }
      iret=ReadTMComponents(X, &dnorm, &liLanczosStp, 0);
      if(iret !=TRUE){
        fprintf(stdoutMPI, "  Error: Fail to read TMcomponents\n");
        return -2;
      }
      if(liLanczosStp_vec !=liLanczosStp){
        fprintf(stdoutMPI, "  Error: Input files for vector and TMcomponents are incoorect.\n");
        fprintf(stdoutMPI, "  Error: Input vector %ld th stps, TMcomponents  %ld th stps.\n", liLanczosStp_vec, liLanczosStp);
        return -2;
      }
      X->Def.Lanczos_restart=liLanczosStp;
      //Calculate EigenValue

      liLanczosStp = liLanczosStp+X->Def.Lanczos_max;
      alpha1=alpha[X->Def.Lanczos_restart];
      beta1=beta[X->Def.Lanczos_restart];
    }/*X->Def.iReStart == RESTART_INOUT || X->Def.iReStart == RESTART_IN*/
    else {
      SetInitialVector(X, v0, v1);
      //Eigenvalues by Lanczos method
      TimeKeeper(X, cFileNameTimeKeep, cLanczos_EigenValueStart, "a");
      StartTimer(4101);
      mltply(X, v0, v1);
      StopTimer(4101);
      stp = 1;
      TimeKeeperWithStep(X, cFileNameTimeKeep, cLanczos_EigenValueStep, "a", stp);

      alpha1 = creal(X->Large.prdct);// alpha = v^{\dag}*H*v

      alpha[1] = alpha1;
      cbeta1 = 0.0;

#pragma omp parallel for reduction(+:cbeta1) default(none) private(i) shared(v0, v1) firstprivate(i_max, alpha1)
      for (i = 1; i <= i_max; i++) {
        cbeta1 += conj(v0[i] - alpha1 * v1[i]) * (v0[i] - alpha1 * v1[i]);
      }
      cbeta1 = SumMPI_dc(cbeta1);
      beta1 = creal(cbeta1);
      beta1 = sqrt(beta1);
      beta[1] = beta1;
      ebefor = 0;
      liLanczosStp = X->Def.Lanczos_max;
      X->Def.Lanczos_restart =1;
    }/*else restart*/

  /*
   * Set Maximum number of loop to the dimention of the Wavefunction
   */
  i_max_tmp = SumMPI_li(i_max);
  if (i_max_tmp < liLanczosStp) {
    liLanczosStp = i_max_tmp;
  }
  if (i_max_tmp < X->Def.LanczosTarget) {
    liLanczosStp = i_max_tmp;
  }
  if (i_max_tmp == 1) {
    E[1] = alpha[1];
    StartTimer(4102);
    vec12(alpha, beta, stp, E, X);
    StopTimer(4102);
    X->Large.itr = stp;
    X->Phys.Target_energy = E[k_exct];
    iconv = 0;
    fprintf(stdoutMPI, "  LanczosStep  E[1] \n");
    fprintf(stdoutMPI, "  stp=%d %.10lf \n", stp, E[1]);
  }

  fprintf(stdoutMPI, "  LanczosStep  E[1] E[2] E[3] E[4] Target:E[%d] E_Max/Nsite\n", X->Def.LanczosTarget + 1);
  for (stp = X->Def.Lanczos_restart+1; stp <= liLanczosStp; stp++) {
#pragma omp parallel for default(none) private(i,temp1, temp2) shared(v0, v1) firstprivate(i_max, alpha1, beta1)
    for (i = 1; i <= i_max; i++) {
      temp1 = v1[i];
      temp2 = (v0[i] - alpha1 * v1[i]) / beta1;
      v0[i] = -beta1 * temp1;
      v1[i] = temp2;
    }

    StartTimer(4101);
    mltply(X, v0, v1);
    StopTimer(4101);
    TimeKeeperWithStep(X, cFileNameTimeKeep, cLanczos_EigenValueStep, "a", stp);
    alpha1 = creal(X->Large.prdct);
    alpha[stp] = alpha1;
    cbeta1 = 0.0;

#pragma omp parallel for reduction(+:cbeta1) default(none) private(i) shared(v0, v1) firstprivate(i_max, alpha1)
    for (i = 1; i <= i_max; i++) {
      cbeta1 += conj(v0[i] - alpha1 * v1[i]) * (v0[i] - alpha1 * v1[i]);
    }
    cbeta1 = SumMPI_dc(cbeta1);
    beta1 = creal(cbeta1);
    beta1 = sqrt(beta1);
    beta[stp] = beta1;

    Target = X->Def.LanczosTarget;

    if (stp == 2) {
      d_malloc2(tmp_mat, stp, stp);
      d_malloc1(tmp_E, stp + 1);

      for (int_i = 0; int_i < stp; int_i++) {
        for (int_j = 0; int_j < stp; int_j++) {
          tmp_mat[int_i][int_j] = 0.0;
        }
      }
      tmp_mat[0][0] = alpha[1];
      tmp_mat[0][1] = beta[1];
      tmp_mat[1][0] = beta[1];
      tmp_mat[1][1] = alpha[2];
      DSEVvalue(stp, tmp_mat, tmp_E);
      E[1] = tmp_E[0];
      E[2] = tmp_E[1];
      if (Target < 2) {
        E_target = tmp_E[Target];
        ebefor = E_target;
      }
      d_free1(tmp_E, stp + 1);
      d_free2(tmp_mat, stp, stp);

      childfopenMPI(sdt_2, "w", &fp);

      fprintf(fp, "LanczosStep  E[1] E[2] E[3] E[4] Target:E[%d] E_Max/Nsite\n", Target + 1);
      if (Target < 2) {
        fprintf(stdoutMPI, "  stp = %d %.10lf %.10lf xxxxxxxxxx xxxxxxxxx  %.10lf xxxxxxxxx \n", stp, E[1], E[2],
                E_target);
        fprintf(fp, "stp = %d %.10lf %.10lf xxxxxxxxxx xxxxxxxxx %.10lf xxxxxxxxx \n", stp, E[1], E[2], E_target);
      } else {
        fprintf(stdoutMPI, "  stp = %d %.10lf %.10lf xxxxxxxxxx xxxxxxxxx xxxxxxxxx xxxxxxxxx \n", stp, E[1], E[2]);
        fprintf(fp, "stp = %d %.10lf %.10lf xxxxxxxxxx xxxxxxxxx xxxxxxxxx xxxxxxxxx \n", stp, E[1], E[2]);
      }

      fclose(fp);
    }

    //if (stp > 2 && stp % 2 == 0) {
    if (stp > 2) {
      childfopenMPI(sdt_2, "a", &fp);

      d_malloc2(tmp_mat, stp, stp);
      d_malloc1(tmp_E, stp + 1);

      for (int_i = 0; int_i < stp; int_i++) {
        for (int_j = 0; int_j < stp; int_j++) {
          tmp_mat[int_i][int_j] = 0.0;
        }
      }
      tmp_mat[0][0] = alpha[1];
      tmp_mat[0][1] = beta[1];
      for (int_i = 1; int_i < stp - 1; int_i++) {
        tmp_mat[int_i][int_i] = alpha[int_i + 1];
        tmp_mat[int_i][int_i + 1] = beta[int_i + 1];
        tmp_mat[int_i][int_i - 1] = beta[int_i];
      }
      tmp_mat[int_i][int_i] = alpha[int_i + 1];
      tmp_mat[int_i][int_i - 1] = beta[int_i];
      StartTimer(4103);
      DSEVvalue(stp, tmp_mat, tmp_E);
      StopTimer(4103);
      E[1] = tmp_E[0];
      E[2] = tmp_E[1];
      E[3] = tmp_E[2];
      E[4] = tmp_E[3];
      E[0] = tmp_E[stp - 1];
      if (stp > Target) {
        E_target = tmp_E[Target];
      }
      d_free1(tmp_E, stp + 1);
      d_free2(tmp_mat, stp, stp);
      if (stp > Target) {
        fprintf(stdoutMPI, "  stp = %d %.10lf %.10lf %.10lf %.10lf %.10lf %.10lf\n", stp, E[1], E[2], E[3], E[4],
                E_target, E[0] / (double) X->Def.NsiteMPI);
        fprintf(fp, "stp=%d %.10lf %.10lf %.10lf %.10lf %.10lf %.10lf\n", stp, E[1], E[2], E[3], E[4], E_target,
                E[0] / (double) X->Def.NsiteMPI);
      } else {
        fprintf(stdoutMPI, "  stp = %d %.10lf %.10lf %.10lf %.10lf xxxxxxxxx %.10lf\n", stp, E[1], E[2], E[3], E[4],
                E[0] / (double) X->Def.NsiteMPI);
        fprintf(fp, "stp=%d %.10lf %.10lf %.10lf %.10lf xxxxxxxxx %.10lf\n", stp, E[1], E[2], E[3], E[4],
                E[0] / (double) X->Def.NsiteMPI);
      }
      fclose(fp);
      if (stp > Target) {
        if (fabs((E_target - ebefor) / E_target) < eps_Lanczos || fabs(beta[stp]) < pow(10.0, -14)) {
          /*
          if(X->Def.iReStart == RESTART_INOUT ||X->Def.iReStart == RESTART_OUT){
            break;
          }
           */
          d_malloc1(tmp_E, stp + 1);
          StartTimer(4102);
          vec12(alpha, beta, stp, tmp_E, X);
          StopTimer(4102);
          X->Large.itr = stp;
          X->Phys.Target_energy = E_target;
          X->Phys.Target_CG_energy = tmp_E[k_exct]; //for CG
          iconv = 0;
          d_free1(tmp_E, stp + 1);
          break;
        }
        ebefor = E_target;
      }

    }
  }
  if (X->Def.iReStart == RESTART_INOUT ||X->Def.iReStart == RESTART_OUT ){
    if(stp != X->Def.Lanczos_restart+2) { // 2 steps are needed to get the value: E[stp+2]-E[stp+1]
      OutputTMComponents(X, alpha, beta, dnorm, stp - 1);
      OutputLanczosVector(X, v0, v1, stp - 1);
    }
    //return 0;
  }

  sprintf(sdt, cFileNameTimeKeep, X->Def.CDataFileHead);
  if (iconv != 0) {
    sprintf(sdt, "%s", cLogLanczos_EigenValueNotConverged);
    fprintf(stdoutMPI, "  Lanczos Eigenvalue is not converged in this process.\n");
    return -1;
  }

  TimeKeeper(X, cFileNameTimeKeep, cLanczos_EigenValueFinish, "a");
  fprintf(stdoutMPI, "%s", cLogLanczos_EigenValueEnd);

  return 0;
}

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.

Parameters

X	[in] Struct to give the information to calculate triangular matrix components.
_alpha	[in,out] Triangular matrix components.
_beta	[in,out] Triangular matrix components.
tmp_v1	[in, out] A temporary vector to calculate triangular matrix components.
liLanczos_step	[in] The max iteration step.

Version

1.2

Returns

TRUE

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 321 of file Lanczos_EigenValue.c.

References alpha, beta, c_Lanczos_SpectrumStep, cFileNameLanczosStep, cFileNameTimeKeep, D_FileNameMax, mltply(), stdoutMPI, SumMPI_dc(), SumMPI_li(), TimeKeeperWithStep(), TRUE, v0, v1, and X.

Referenced by CalcSpectrumByTPQ().

   {

  char sdt[D_FileNameMax];
  int stp;
  long int i, i_max;
  i_max = X->Check.idim_max;

  unsigned long int i_max_tmp;
  double beta1, alpha1; //beta,alpha1 should be real
  double complex temp1, temp2;
  double complex cbeta1;

  sprintf(sdt, cFileNameLanczosStep, X->Def.CDataFileHead);

  /*
    Set Maximum number of loop to the dimension of the Wavefunction
  */
  i_max_tmp = SumMPI_li(i_max);
  if (i_max_tmp < *liLanczos_step || i_max_tmp < X->Def.LanczosTarget) {
    *liLanczos_step = i_max_tmp;
  }

  if (X->Def.Lanczos_restart == 0) { // initial procedure (not restart)
#pragma omp parallel for default(none) private(i) shared(v0, v1, tmp_v1) firstprivate(i_max)
    for (i = 1; i <= i_max; i++) {
      v0[i] = 0.0;
      v1[i] = tmp_v1[i];
    }
    stp = 0;
    mltply(X, v0, tmp_v1);
    TimeKeeperWithStep(X, cFileNameTimeKeep, c_Lanczos_SpectrumStep, "a", stp);
    alpha1 = creal(X->Large.prdct);// alpha = v^{\dag}*H*v
    _alpha[1] = alpha1;
    cbeta1 = 0.0;
    fprintf(stdoutMPI, "  Step / Step_max alpha beta \n");

#pragma omp parallel for reduction(+:cbeta1) default(none) private(i) shared(v0, v1) firstprivate(i_max, alpha1)
    for (i = 1; i <= i_max; i++) {
      cbeta1 += conj(v0[i] - alpha1 * v1[i]) * (v0[i] - alpha1 * v1[i]);
    }
    cbeta1 = SumMPI_dc(cbeta1);
    beta1 = creal(cbeta1);
    beta1 = sqrt(beta1);
    _beta[1] = beta1;
    X->Def.Lanczos_restart = 1;
  } else {  // restart case
    alpha1 = alpha[X->Def.Lanczos_restart];
    beta1 = beta[X->Def.Lanczos_restart];
  }

  for (stp = X->Def.Lanczos_restart + 1; stp <= *liLanczos_step; stp++) {

    if (fabs(_beta[stp - 1]) < pow(10.0, -14)) {
      *liLanczos_step = stp - 1;
      break;
    }

#pragma omp parallel for default(none) private(i, temp1, temp2) shared(v0, v1) firstprivate(i_max, alpha1, beta1)
    for (i = 1; i <= i_max; i++) {
      temp1 = v1[i];
      temp2 = (v0[i] - alpha1 * v1[i]) / beta1;
      v0[i] = -beta1 * temp1;
      v1[i] = temp2;
    }

    mltply(X, v0, v1);
    TimeKeeperWithStep(X, cFileNameTimeKeep, c_Lanczos_SpectrumStep, "a", stp);
    alpha1 = creal(X->Large.prdct);
    _alpha[stp] = alpha1;
    cbeta1 = 0.0;

#pragma omp parallel for reduction(+:cbeta1) default(none) private(i) shared(v0, v1) firstprivate(i_max, alpha1)
    for (i = 1; i <= i_max; i++) {
      cbeta1 += conj(v0[i] - alpha1 * v1[i]) * (v0[i] - alpha1 * v1[i]);
    }
    cbeta1 = SumMPI_dc(cbeta1);
    beta1 = creal(cbeta1);
    beta1 = sqrt(beta1);
    _beta[stp] = beta1;
    if(stp %10 == 0) {
      fprintf(stdoutMPI, "  stp = %d / %lu %.10lf  %.10lf \n", stp,  *liLanczos_step, alpha1, beta1);
    }
  }

  return TRUE;
}

OutputLanczosVector()

int OutputLanczosVector	(	struct BindStruct *	X,
		double complex *	tmp_v0,
		double complex *	tmp_v1,
		unsigned long int	liLanczosStp_vec
	)

Output initial vectors for the restart calculation.

Parameters

X	[in] Give the dimension for the vector _v0 and _v1.
tmp_v0	[in] The outputted vector for recalculation \( v_0 \).
tmp_v1	[in] The outputted vector for recalculation \( v_1 \).
liLanczosStp_vec	[in] The step for finishing the iteration.

Return values

-1	Fail to output the vector.
0	Succeed to output the vector.

Version

2.0

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 463 of file Lanczos_EigenValue.c.

References c_OutputSpectrumRecalcvecEnd, c_OutputSpectrumRecalcvecStart, cFileNameOutputRestartVec, cFileNameTimeKeep, childfopenALL(), D_FileNameMax, myrank, stdoutMPI, TimeKeeper(), and X.

Referenced by Lanczos_EigenValue().

                                                           {
  char sdt[D_FileNameMax];
  FILE *fp;

  fprintf(stdoutMPI, "    Start: Output vectors for recalculation.\n");
  TimeKeeper(X, cFileNameTimeKeep, c_OutputSpectrumRecalcvecStart, "a");

  sprintf(sdt, cFileNameOutputRestartVec, X->Def.CDataFileHead, myrank);
  if(childfopenALL(sdt, "wb", &fp)!=0){
    return -1;
  }
  fwrite(&liLanczosStp_vec, sizeof(liLanczosStp_vec),1,fp);
  fwrite(&X->Check.idim_max, sizeof(X->Check.idim_max),1,fp);
  fwrite(tmp_v0, sizeof(complex double),X->Check.idim_max+1, fp);
  fwrite(tmp_v1, sizeof(complex double),X->Check.idim_max+1, fp);
  fclose(fp);

  fprintf(stdoutMPI, "    End:   Output vectors for recalculation.\n");
  TimeKeeper(X, cFileNameTimeKeep, c_OutputSpectrumRecalcvecEnd, "a");
  return 0;
}

OutputTMComponents()

int OutputTMComponents	(	struct BindStruct *	X,
		double *	_alpha,
		double *	_beta,
		double	_dnorm,
		int	liLanczosStp
	)

Output tridiagonal matrix components obtained by the Lanczos method.

Parameters

X	[in] Give the input file name.
_alpha	[in] The array of tridiagonal matrix components.
_beta	[in] The array of tridiagonal matrix components.
_dnorm	[in] The norm.
liLanczosStp	[in] The iteration step.

Return values

FALSE	Fail to open the file for the output.
TRUE	Succeed to open the file for the output.

Version

1.2

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 691 of file Lanczos_EigenValue.c.

References cFileNameTridiagonalMatrixComponents, childfopenMPI(), D_FileNameMax, FALSE, TRUE, and X.

Referenced by CalcSpectrumByTPQ(), and Lanczos_EigenValue().

{
  char sdt[D_FileNameMax];
  unsigned long int i;
  FILE *fp;

  sprintf(sdt, cFileNameTridiagonalMatrixComponents, X->Def.CDataFileHead);
  if(childfopenMPI(sdt,"w", &fp)!=0){
    return FALSE;
  }
  fprintf(fp, "%d \n",liLanczosStp);
  fprintf(fp, "%.10lf \n",_dnorm);
  for( i = 1 ; i <= liLanczosStp; i++){
    fprintf(fp,"%.10lf %.10lf \n", _alpha[i], _beta[i]);
  }
  fclose(fp);
  return TRUE;
}

ReadInitialVector()

int ReadInitialVector	(	struct BindStruct *	X,
		double complex *	_v0,
		double complex *	_v1,
		unsigned long int *	liLanczosStp_vec
	)

Read initial vectors for the restart calculation.

Parameters

X	[in] Give the dimension for the vector _v0 and _v1.
_v0	[out] The inputted vector for recalculation \( v_0 \).
_v1	[out] The inputted vector for recalculation \( v_1 \).
liLanczosStp_vec	[in] The max iteration step.

Return values

-1	Fail to read the initial vector.
0	Succeed to read the initial vector.

Version

1.2

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 425 of file Lanczos_EigenValue.c.

References c_InputSpectrumRecalcvecEnd, c_InputSpectrumRecalcvecStart, cFileNameOutputRestartVec, cFileNameTimeKeep, childfopenALL(), D_FileNameMax, myrank, stdoutMPI, TimeKeeper(), and X.

Referenced by Lanczos_EigenValue().

{
  size_t byte_size;
  char sdt[D_FileNameMax];
  FILE *fp;
  unsigned long int i_max;

  fprintf(stdoutMPI, "  Start: Input vectors for recalculation.\n");
  TimeKeeper(X, cFileNameTimeKeep, c_InputSpectrumRecalcvecStart, "a");
  sprintf(sdt, cFileNameOutputRestartVec, X->Def.CDataFileHead, myrank);
  if (childfopenALL(sdt, "rb", &fp) != 0) {
    return -1;
  }
  byte_size = fread(liLanczosStp_vec, sizeof(*liLanczosStp_vec),1,fp);
  byte_size = fread(&i_max, sizeof(long int), 1, fp);
  if(i_max != X->Check.idim_max){
    fprintf(stderr, "Error: A size of Inputvector is incorrect.\n");
    return -1;
  }
  byte_size = fread(_v0, sizeof(complex double), X->Check.idim_max + 1, fp);
  byte_size = fread(_v1, sizeof(complex double), X->Check.idim_max + 1, fp);
  fclose(fp);
  fprintf(stdoutMPI, "  End:   Input vectors for recalculation.\n");
  TimeKeeper(X, cFileNameTimeKeep, c_InputSpectrumRecalcvecEnd, "a");
  if (byte_size == 0) printf("byte_size: %d \n", (int)byte_size);
  return 0;
}

ReadTMComponents()

int ReadTMComponents	(	struct BindStruct *	X,
		double *	_dnorm,
		unsigned long int *	_i_max,
		int	iFlg
	)

Read tridiagonal matrix components obtained by the Lanczos method.
.

Note

The arrays of tridiaonal components alpha and beta are global arrays.

Parameters

X	[in] Give the iteration number for the recalculation and the input file name.
_dnorm	[out] Get the norm.
_i_max	[in] The iteration step for the input data.
iFlg	[in] Flag for the recalculation.

Return values

FALSE	Fail to read the file.
TRUE	Succeed to read the file.

Version

1.2

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 621 of file Lanczos_EigenValue.c.

References alpha, beta, cFileNameTridiagonalMatrixComponents, childfopenMPI(), D_FileNameMax, FALSE, fgetsMPI(), TRUE, vec, and X.

Referenced by CalcSpectrumByTPQ(), and Lanczos_EigenValue().

 {
  char sdt[D_FileNameMax];
  char ctmp[256];

  unsigned long int idx, i, ivec;
  unsigned long int i_max;
  double dnorm;
  FILE *fp;
  idx=1;
  sprintf(sdt, cFileNameTridiagonalMatrixComponents, X->Def.CDataFileHead);
  if(childfopenMPI(sdt,"r", &fp)!=0){
    return FALSE;
  }

  fgetsMPI(ctmp, sizeof(ctmp)/sizeof(char), fp);
  sscanf(ctmp,"%ld \n", &i_max);
  if (X->Def.LanczosTarget > X->Def.nvec) {
    ivec = X->Def.LanczosTarget + 1;
  }
  else {
    ivec =X->Def.nvec + 1;
  }

  if(iFlg==0) {
    alpha = (double *) realloc(alpha, sizeof(double) * (i_max + X->Def.Lanczos_max + 1));
    beta = (double *) realloc(beta, sizeof(double) * (i_max + X->Def.Lanczos_max + 1));
    for (i = 0; i <  ivec; i++) {
      vec[i] = (complex double *) realloc(vec[i], (i_max + X->Def.Lanczos_max + 1) * sizeof(complex double));
    }
  }
  else if(iFlg==1){
    alpha=(double*)realloc(alpha, sizeof(double)*(i_max + 1));
    beta=(double*)realloc(beta, sizeof(double)*(i_max + 1));
    for (i = 0; i < ivec; i++) {
      vec[i] = (complex double *) realloc(vec[i], (i_max + X->Def.Lanczos_max + 1) * sizeof(complex double));
    }
  }
  else{
    fclose(fp);
    return FALSE;
  }
  fgetsMPI(ctmp, sizeof(ctmp)/sizeof(char), fp);
  sscanf(ctmp,"%lf \n", &dnorm);
  while(fgetsMPI(ctmp, sizeof(ctmp)/sizeof(char), fp) != NULL){
    sscanf(ctmp,"%lf %lf \n", &alpha[idx], &beta[idx]);
    idx++;
  }
  fclose(fp);
  *_dnorm=dnorm;
  *_i_max=i_max;
  return TRUE;
}

SetInitialVector()

void SetInitialVector	(	struct BindStruct *	X,
		double complex *	tmp_v0,
		double complex *	tmp_v1
	)

Set initial vector to start the calculation for Lanczos method.
.

Parameters

X	[in, out] Get the information of reading initisl vectors. Input: idim_max, iFlgMPI, k_exct, iInitialVecType. Output: Large.iv.
tmp_v0	[out] The initial vector whose components are zero.
tmp_v1	[out] The initial vector whose components are randomly given when initial_mode=1, otherwise, iv-th component is only given.

Definition at line 496 of file Lanczos_EigenValue.c.

References BcastMPI_li(), initial_mode, myrank, nproc, nthreads, stdoutMPI, SumMPI_dc(), SumMPI_li(), and X.

Referenced by Lanczos_EigenValue().

                                                                                            {
  int iproc;
  long int i, iv, i_max;
  unsigned long int i_max_tmp, sum_i_max;
  int mythread;

// for GC
  double dnorm;
  double complex cdnorm;
  long unsigned int u_long_i;
  dsfmt_t dsfmt;

  i_max = X->Check.idim_max;
  if (initial_mode == 0) {
    if(X->Def.iFlgMPI==0) {
      sum_i_max = SumMPI_li(X->Check.idim_max);
    }
    else{
      sum_i_max =X->Check.idim_max;
    }
    X->Large.iv = (sum_i_max / 2 + X->Def.initial_iv) % sum_i_max + 1;
    iv = X->Large.iv;
    fprintf(stdoutMPI, "  initial_mode=%d normal: iv = %ld i_max=%ld k_exct =%d \n\n", initial_mode, iv, i_max,
            X->Def.k_exct);
#pragma omp parallel for default(none) private(i) shared(tmp_v0, tmp_v1) firstprivate(i_max)
    for (i = 1; i <= i_max; i++) {
      tmp_v0[i] = 0.0;
      tmp_v1[i] = 0.0;
    }

    sum_i_max = 0;
    if(X->Def.iFlgMPI==0) {
      for (iproc = 0; iproc < nproc; iproc++) {
        i_max_tmp = BcastMPI_li(iproc, i_max);
        if (sum_i_max <= iv && iv < sum_i_max + i_max_tmp) {
          if (myrank == iproc) {
            tmp_v1[iv - sum_i_max + 1] = 1.0;
            if (X->Def.iInitialVecType == 0) {
              tmp_v1[iv - sum_i_max + 1] += 1.0 * I;
              tmp_v1[iv - sum_i_max + 1] /= sqrt(2.0);
            }
          }/*if (myrank == iproc)*/
        }/*if (sum_i_max <= iv && iv < sum_i_max + i_max_tmp)*/

        sum_i_max += i_max_tmp;

      }/*for (iproc = 0; iproc < nproc; iproc++)*/
    }
    else {
      tmp_v1[iv + 1] = 1.0;
      if (X->Def.iInitialVecType == 0) {
        tmp_v1[iv + 1] += 1.0 * I;
        tmp_v1[iv + 1] /= sqrt(2.0);
      }
    }
  }/*if(initial_mode == 0)*/
  else if (initial_mode == 1) {
    iv = X->Def.initial_iv;
    fprintf(stdoutMPI, "  initial_mode=%d (random): iv = %ld i_max=%ld k_exct =%d \n\n", initial_mode, iv, i_max,
            X->Def.k_exct);
#pragma omp parallel default(none) private(i, u_long_i, mythread, dsfmt) \
            shared(tmp_v0, tmp_v1, iv, X, nthreads, myrank) firstprivate(i_max)
    {

#pragma omp for
      for (i = 1; i <= i_max; i++) {
        tmp_v0[i] = 0.0;
      }
      /*
       Initialise MT
      */
#ifdef _OPENMP
      mythread = omp_get_thread_num();
#else
      mythread = 0;
#endif
      if(X->Def.iFlgMPI==0) {
        u_long_i = 123432 + labs(iv) + mythread + nthreads * myrank;
      }
      else{
        u_long_i = 123432 + labs(iv)+ mythread ;
      }
      dsfmt_init_gen_rand(&dsfmt, u_long_i);

      if (X->Def.iInitialVecType == 0) {
#pragma omp for
        for (i = 1; i <= i_max; i++)
          tmp_v1[i] = 2.0 * (dsfmt_genrand_close_open(&dsfmt) - 0.5) +
                      2.0 * (dsfmt_genrand_close_open(&dsfmt) - 0.5) * I;
      } else {
#pragma omp for
        for (i = 1; i <= i_max; i++)
          tmp_v1[i] = 2.0 * (dsfmt_genrand_close_open(&dsfmt) - 0.5);
      }

    }/*#pragma omp parallel*/

    cdnorm = 0.0;
#pragma omp parallel for default(none) private(i) shared(tmp_v1, i_max) reduction(+: cdnorm)
    for (i = 1; i <= i_max; i++) {
      cdnorm += conj(tmp_v1[i]) * tmp_v1[i];
    }
    if(X->Def.iFlgMPI==0) {
      cdnorm = SumMPI_dc(cdnorm);
    }
    dnorm = creal(cdnorm);
    dnorm = sqrt(dnorm);
#pragma omp parallel for default(none) private(i) shared(tmp_v1) firstprivate(i_max, dnorm)
    for (i = 1; i <= i_max; i++) {
      tmp_v1[i] = tmp_v1[i] / dnorm;
    }
  }/*else if(initial_mode==1)*/
}