lapack_diag.cpp File Reference

#include "lapack_diag.hpp"
#include "matrixlapack.hpp"
#include "FileIO.hpp"

Go to the source code of this file.

Functions
int	lapack_diag (struct BindStruct *X)
	performing full diagonalization using lapack More...

Function Documentation

lapack_diag()

int lapack_diag

(

struct BindStruct *

X

)

performing full diagonalization using lapack

Parameters

[in,out]

X

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Returns

Parameters

[in,out]

X

Definition at line 35 of file lapack_diag.cpp.

References BindStruct::Check, childfopenMPI(), BindStruct::Def, PhysList::energy, CheckList::idim_max, DefineList::iNGPU, myrank, nproc, BindStruct::Phys, stdoutMPI, v0, v1, and ZHEEVall().

Referenced by CalcByFullDiag(), and CalcSpectrumByFullDiag().

  {

  FILE *fp;
  char sdt[D_FileNameMax] = "";
  long int i, j, i_max, xMsize;
#ifdef _SCALAPACK
  int rank, size, nprocs, nprow, npcol, myrow, mycol, ictxt;
  int i_negone=-1, i_zero=0, iam;
  long int mb, nb, mp, nq;
  int dims[2]={0,0};
#endif

  i_max = X->Check.idim_max;

  for (i = 0; i < i_max; i++) {
    for (j = 0; j < i_max; j++) {
      v0[i][j] = v0[i + 1][j];
    }
  }
  xMsize = i_max;
  if (X->Def.iNGPU == 0) {
#ifdef _SCALAPACK
    if(nproc >1) {
      fprintf(stdoutMPI, "Using SCALAPACK\n\n");
      MPI_Comm_rank(MPI_COMM_WORLD, &rank);
      MPI_Comm_size(MPI_COMM_WORLD, &size);
      MPI_Dims_create(size,2,dims);
      nprow=dims[0]; npcol=dims[1];

      blacs_pinfo_(&iam, &nprocs);
      blacs_get_(&i_negone, &i_zero, &ictxt);
      blacs_gridinit_(&ictxt, "R", &nprow, &npcol);
      blacs_gridinfo_(&ictxt, &nprow, &npcol, &myrow, &mycol);
      
      mb = GetBlockSize(xMsize, size);

      mp = numroc_(&xMsize, &mb, &myrow, &i_zero, &nprow);
      nq = numroc_(&xMsize, &mb, &mycol, &i_zero, &npcol);
      Z_vec = malloc(mp*nq*sizeof(std::complex<double>));
      //printf("xMsize %d\n", xMsize);
      /*if(rank == 0){
        for(i=0; i<xMsize; i++){
          for(j=0; j<xMsize; j++){
            printf("Ham %f %f ", real(Ham[i][j]), imag(Ham[i][j]));
          }
        }
      }*/
      diag_scalapack_cmp(xMsize, Ham, v0, Z_vec, descZ_vec);
      //printf("Z %f %f\n", real(Z_vec[0]), imag(Z_vec[1]));
    } else {
      ZHEEVall(xMsize, Ham, v0, v1);
    }
#else
    ZHEEVall(xMsize, v0, X->Phys.energy, v1);
#endif
  } else {
#ifdef _MAGMA
    if(myrank==0){
      if(diag_magma_cmp(xMsize, Ham, v0, v1, X->Def.iNGPU) != 0) {
        return -1;
      }
    }
#else
    fprintf(stdoutMPI, "Warning: MAGMA is not used in this calculation.");
    ZHEEVall(xMsize, v0, X->Phys.energy, v1);
#endif
  }
  for (i = 0; i < i_max; i++) {
    for (j = 0; j < i_max; j++) {
      v0[i + 1][j] = v1[i][j];
    }
  }

  strcpy(sdt, "Eigenvalue.dat");
  if (childfopenMPI(sdt, "w", &fp) != 0) {
    return -1;
  }
  for (i = 0; i < i_max; i++) {
    fprintf(fp, " %ld %.10lf \n", i, X->Phys.energy[i]);
  }
  fclose(fp);
  return 0;
}