matrixscalapack.cpp

Go to the documentation of this file.

/* HPhi  -  Quantum Lattice Model Simulator */
/* Copyright (C) 2015 The University of Tokyo */

/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */

/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
/* GNU General Public License for more details. */

/* You should have received a copy of the GNU General Public License */
/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
#include "matrixscalapack.hpp"
#ifdef _SCALAPACK
int use_scalapack = 0;

long int GetBlockSize(long int Msize, long int nproc) {
  long int block_size = 16;
  if(Msize*Msize/nproc > block_size*block_size)
    return block_size;
  return 1;
}

long int GetPArrayIndex(long int i, long int np, long int nb) {
  return (i/nb)%np;
}

long int GetLocalIndex(long int i, long int np, long int nb) {
  return (i/(np*nb))*nb + i%nb;
}

long int GetGlobalIndex(long int il, long int p, long int np, long int nb){
  return ((il/nb)*np+p)*nb + il%nb;
}

long int MatToRank(long int i, long int j, long int nprow, long int npcol, long int nb){
  long int iproc, jproc;
  iproc = GetPArrayIndex(i, nprow, nb);
  jproc = GetPArrayIndex(j, npcol, nb);
  return iproc+jproc*nprow;
}

long int GetMatRawInRank(long int lj, long int rank, long int npcol, long int nb){
  long int pcol = rank/npcol;
  return GetGlobalIndex(lj, pcol, npcol, nb);
}

long int *GetMatElementInRank(long int i, long int j, long int nprow, long int npcol, long int nb){
  long int *ij;
  ij = malloc(2*sizeof(int));
  ij[0] = GetLocalIndex(i, nprow, nb);
  ij[1] = GetLocalIndex(j, npcol, nb);
  return ij;
}

void DivMat(long int m, long int n, std::complex<double> Aorgmn, std::complex<double> *A, int *desca){
  long int mp = m+1, np = n+1;
  pzelset_(A, &mp, &np, desca, &Aorgmn);
}

void GetEigenVector(long int i, long int m, std::complex<double> *Z, int *descZ, std::complex<double> *vec) {
  std::complex<double> alpha;
  long int j, ip, jp;
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  ip = i+1;
  for(j=0; j<m; j++){
    jp = j+1;
    pzelget_("A", " ", &alpha, Z, &jp, &ip, descZ);
    if(rank==0) {
      vec[j] = alpha;
    }
  }
}

int diag_scalapack_cmp(long int xNsize, std::complex<double> **A, 
                       std::complex<double> *r, std::complex<double> *Z, int *descZ) {
  const int i_one=1, i_zero=0;
  const long int i_negone=-1;
  const double zero=0.0, one=1.0;
  long int m, n, mb, nb;
  int nprow, npcol;
  int myrow, mycol, info, lld;
  long int mp, nq;
  int ictxt;
  std::complex<double> *A_distr, *work, *rwork;
  double *W;
  int descA_distr[9];
  int rank, size, iam, nprocs;
  long int lwork, lrwork;
  int dims[2]={0,0};
  long int i, j, ip, jp;
  m=n=xNsize;

  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_((int *)&i_negone, &i_zero, &ictxt);
  blacs_gridinit_(&ictxt, "R", &nprow, &npcol);
  blacs_gridinfo_(&ictxt, &nprow, &npcol, &myrow, &mycol);
 
  mb = GetBlockSize(m, size);
  nb = GetBlockSize(n, size);

  mp = numroc_(&m, &mb, &myrow, &i_zero, &nprow);
  nq = numroc_(&n, &nb, &mycol, &i_zero, &npcol);
  W = malloc(n*sizeof(double));
  A_distr = malloc(mp*nq*sizeof(std::complex<double>));

  lld = (mp>0) ? mp : 1;
  descinit_(descA_distr, &m, &n, &mb, &nb, &i_zero, &i_zero, &ictxt, &lld, &info);
  descinit_(descZ, &m, &n, &mb, &nb, &i_zero, &i_zero, &ictxt, &lld, &info);

  for(i=0; i<m; i++){
    for(j=0; j<n; j++){
      DivMat(i, j, A[i][j], A_distr, descA_distr);
    }
  }

  std::complex<double> wkopt, rwkopt;
  pzheev_("V", "U", &n, A_distr, &i_one, &i_one, descA_distr, W, Z, &i_one, &i_one, descZ, &wkopt, &i_negone, &rwkopt, &i_negone, &info);

  lwork = (long int)wkopt;
  lrwork = (long int)rwkopt;
  work = malloc(lwork*sizeof(std::complex<double>));
  rwork = malloc(lrwork*sizeof(std::complex<double>));


  pzheev_("V", "U", &n, A_distr, &i_one, &i_one, descA_distr, W, Z, &i_one, &i_one, descZ, work, &lwork, rwork, &lrwork, &info);

  if(rank == 0){
    for(i=0; i<n; i++){
      r[i] = W[i];
    }
  }

  free(A_distr);
  free(work);
  free(rwork);
  free(W);

  use_scalapack = 1;

  return 0;
}

#endif