lapack_diag.c File Reference

#include "lapack_diag.h"
#include "matrixlapack.h"
#include "FileIO.h"

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.c.

References cFileNameEigenvalue_Lanczos, childfopenMPI(), D_FileNameMax, Ham, L_vec, nproc, stdoutMPI, v0, X, 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++) {
      //printf("Ham %f %f ", creal(Ham[i+1][j+1]), cimag(Ham[i+1][j+1]));
      Ham[i][j] = Ham[i + 1][j + 1];
    }
  }
  xMsize = i_max;
  /*for(i=0; i<xMsize; i++){
    for(j=0; j<xMsize; j++){
      printf("Ham %f %f ", creal(Ham[i][j]), cimag(Ham[i][j]));
    }
  }*/
  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(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 ", creal(Ham[i][j]), cimag(Ham[i][j]));
          }
        }
      }*/
      diag_scalapack_cmp(xMsize, Ham, v0, Z_vec, descZ_vec);
      //printf("Z %f %f\n", creal(Z_vec[0]), cimag(Z_vec[1]));
    } else {
      ZHEEVall(xMsize, Ham, v0, L_vec);
    }
#else
    ZHEEVall(xMsize, Ham, v0, L_vec);
#endif
  } else {
#ifdef _MAGMA
    if(diag_magma_cmp(xMsize, Ham, v0, L_vec, X->Def.iNGPU) != 0) {
    return -1;
  }
#else
    fprintf(stdoutMPI, "Warning: MAGMA is not used in this calculation.");
    ZHEEVall(xMsize, Ham, v0, L_vec);
#endif
  }

  /*for (i = 0; i < i_max; i++) {
    for (j = 0; j < i_max; j++) {
      fprintf(stdoutMPI, "%f %f \n", creal(L_vec[i][j]), cimag(L_vec[i][j]));
    }
  }*/

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


  return 0;
}