fleur_elemental.cpp 4.15 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
/*
   Copyright (c) 2014, Daniel Wortmann
   All rights reserved.

   This file provides an interface from FLEUR to Elemental 
*/
#include "elemental.hpp"
using namespace std;
using namespace elem;

// Typedef our real or complex types to 'C' for convenience
#ifdef CPP_INVERSION
typedef double C;
#else
typedef Complex<double> C;
#endif

class global_data {
 public:
  Grid *g;
  mpi::Comm mpi_comm;
  int matrix_dimension;
  DistMatrix<C> *H_mat,*S_mat;
  DistMatrix<C,STAR,VC> eigenvectors;
  DistMatrix<double,VC,STAR> eigenvalues;

  
};

//Global variables
  

static global_data *gd;


DistMatrix<C>* fleur_matrix(int n,C* buffer)
{
  // Create the distributed matrix
  DistMatrix<C,STAR,VC> *mat;
  // The Matrix should be a n x n matrix in 1-D cyclic distribution
  mat= new DistMatrix<C,STAR,VC>(n,n,*(gd->g));
  C* localbuffer=mat->Buffer(); // this is the local buffer of the matrix
  const int buffersize1=mat->LocalWidth();
  const int buffersize2=mat->LocalHeight();
  int localindex=0;
  int fleur_index=0; 
  // Now copy all the data into local buffer to initialize the matrix
  // Loop over columns of local data
  for (int i=mpi::CommRank(gd->mpi_comm);i<n;i+=mpi::CommSize(gd->mpi_comm))
    {
      for (int j=0;j<=i;j++)
	{
	  localbuffer[localindex++]=buffer[fleur_index++];
	}
      // further Off-diagonal elements are set to zero !Probably not needed
      for (int j=0;j<n-i-1;j++)
	{
	  localbuffer[localindex++]=0.0;
	}
    }
  DistMatrix<C> *mat2=new DistMatrix<C>(*mat);
  delete mat;
  //*mat2=*mat;
  return mat2;
}
extern "C"{ 
void fl_el_initialize(int n, C* hbuf, C* sbuf, int mpi_used_comm)
// Set the two matrices
{
  // Initialize the Library
  int argc=0; char** argv;
  Initialize( argc, argv );
  
  //Store the matrix dimension& the mpi_communicator
  gd = new global_data;
  gd->mpi_comm=MPI_Comm_f2c(mpi_used_comm);
  gd->matrix_dimension=n;
  
  
  // First we need a mpi-grid
  gd->g= new Grid(gd->mpi_comm);
  // Store the Matrices
  gd->H_mat=fleur_matrix(n,hbuf);
  gd->S_mat=fleur_matrix(n,sbuf);
  
}

void fl_el_diagonalize(int no_of_eigenpairs)
// Diagonalize the Matrix and return the number of local eigenvalues
{
  
  /* this is for the development version
  // The subset determines the no of eigenvalues found
  HermitianEigSubset<double> subset;
  subset.indexSubset=true;
  subset.lowerIndex=0;
  subset.upperIndex=no_of_global_eigenpairs;
  
  // Space for eigenvalues
  DistMatrix<double> eigenval(g);
  DistMatrix<C> eigenvec(g);
  //default sorting
  const SortType sort = static_cast<SortType>(0);

  //call diagonalization
  HermitianGenDefEig( AXBX, LOWER, H_mat, S_mat, eigenval, eigenvec, sort, subset );
  */
  DistMatrix<double, VR, STAR> eigenval(*(gd->g));
  DistMatrix<C> evec(*(gd->g));
  HermitianGenDefiniteEigType eigtype=AXBX;
  UpperOrLower uplo=UPPER;
  if (mpi::CommRank(gd->mpi_comm)==0) {
    cout<<"H/S-matrix of size "<<gd->matrix_dimension<<endl;
  }
  Display(*(gd->H_mat));
  Display(*(gd->S_mat));
  HermitianGenDefiniteEig(eigtype, uplo, *(gd->H_mat), *(gd->S_mat), eigenval, evec,0,no_of_eigenpairs);
  //redistribute matrices
  //eigenvalues are of type DistMatrix(C,STAR,VC);
  gd->eigenvalues=eigenval;
  gd->eigenvectors=evec;

  no_of_eigenpairs=gd->eigenvectors.LocalWidth();
    }
  

void fl_el_eigenvalues(int neig, double* eig){
  //Return the eigenvalues

  double* buf=gd->eigenvalues.Buffer();
  if (neig > gd->eigenvalues.LocalWidth()*gd->eigenvalues.LocalHeight())
    {
      cerr<<"Error in dimensions in fleur_elemental\n";
    }
  for (int i=0; i<neig;i++){
    eig[i]=buf[i];
  }
}

void fl_el_eigenvectors(int neig, double* eig, C* eigvec){
  //Return all the local eigenvectors&eigenvalues
  double* eigbuf=gd->eigenvalues.Buffer();
  Display(gd->eigenvalues);
  Display(gd->eigenvectors);
  C* eigbuff=gd->eigenvectors.Buffer();
  int local_index=0;
  for (int i=0; i<neig;i++)
    {
      //Copy eigenvalue
      int pe=mpi::CommRank(gd->mpi_comm);
      int in=i*mpi::CommSize(gd->mpi_comm)+pe;
      cout<< "PE:"<<pe<<":"<<i<<"->"<<in<<endl;
      eig[i]=eigbuf[i];
      //Copy eigenvector
      for (int j=0;j<gd->matrix_dimension; j++){
	eigvec[local_index]=eigbuff[local_index];
        local_index++;
      }
    }
}
}