/* DFT++ is a density functional package developed by the research group of Professor Tomas Arias Copyright 1996-2003 Sohrab Ismail-Beigi This file is part of DFT++. DFT++ 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 2 of the License, or (at your option) any later version. DFT++ 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 DFT++; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Please see the file CREDITS for a list of authors. For academic users, we request that publications using results obtained with this software reference "New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000). and, if using the wavelet basis, further reference "Multiresolution analysis of electronic structure: semicardinal and wavelet bases," T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999). and "Robust ab initio calculation of condensed matter: transparent convergence through semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003). For your convenience, preprints of the above articles may be obtained from http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively. */ /* * matrix.c -- implements matrix class */ /* $Id: ComplexMatrix.cpp,v 1.1.2.13 2003/05/29 18:53:56 ivan Exp $ */ #include "header.h" // // The init() function is the main initialization routine. // Sets up sizes and allocates memory for ComplexMatrix class. // All constructors call init() // void ComplexMatrix::init(int nrows,int ncols) { nr = nrows; nc = ncols; hermetian = 0; // If null constructor, set data to point to nothing if (nr*nc == 0) c = NULL; // Otherwise, get some memory else c = (scalar *)mymalloc(sizeof(scalar)*nrows*ncols, "c","ComplexMatrix.init()"); } // Free up memory for ComplexMatrix void ComplexMatrix::freemem(void) { myfree(c); c = NULL; } // constructor with sizes ComplexMatrix::ComplexMatrix(int nrows,int ncols) { init(nrows,ncols); } // Copy constructor ComplexMatrix::ComplexMatrix(const ComplexMatrix &m1) { // Setup copied object to have correct size and get memory for it init(m1.nr,m1.nc); hermetian = m1.hermetian; // Copy over m1 contents int i; for (i=0; i < nr*nc; i++) c[i] = m1.c[i]; } // destructor ComplexMatrix::~ComplexMatrix() { freemem(); } /* Assignment: nonstandard in that it returns void. To make it standard, * replace void -> ComplexMatrix and uncomment the return *this; */ void ComplexMatrix::operator=(const ComplexMatrix &m1) { int i; /* The sizes must agree */ if ( (nr != m1.nr) || (nc != m1.nc) ) die("In ComplexMatrix::operator=, sizes don't agree.\n\n"); hermetian = m1.hermetian; for (i=0; i < nr*nc; i++) c[i] = m1.c[i]; /* return *this; */ } /* Adding matrices, as a friend */ ComplexMatrix operator+(const ComplexMatrix &m1,const ComplexMatrix &m2) { int i; if ( (m1.nr != m2.nr) || (m1.nc != m2.nc) ) die("In ComplexMatrix::operator+, sizes don't agree.\n\n"); ComplexMatrix tm(m1); for (i=0; i < m1.nr*m1.nc; i++) tm.c[i] += m2.c[i]; return tm; } /* Subtracting matrices, as a friend */ ComplexMatrix operator-(const ComplexMatrix &m1,const ComplexMatrix &m2) { int i; if ( (m1.nr != m2.nr) || (m1.nc != m2.nc) ) die("In ComplexMatrix::operator-, sizes don't agree.\n\n"); ComplexMatrix tm(m1); for (i=0; i < m1.nr*m1.nc; i++) tm.c[i] -= m2.c[i]; return tm; } /* Multiplying matrices, as a friend */ ComplexMatrix operator*(const ComplexMatrix &m1, const ComplexMatrix &m2) { #ifdef DFT_PROFILING timerOn(28); // ComplexMatrix*ComplexMatrix timer // Flop count for multiply counterIncr(9,8.0*(double)m1.nr*(double)m2.nc*(double)m1.nc); counterIncr(10); // ComplexMatrix*ComplexMatrix counter #endif // check sizes if (m1.nc != m2.nr) die("In ComplexMatrix::operator*, m1.nc != m2.nr\n\n"); // allocate the result and zero it out ComplexMatrix mprod(m1.nr,m2.nc); mprod.zero_out(); // do the multiply and return matrix_matrix_block_matrix_mult(m1,m2,mprod); #ifdef DFT_PROFILING timerOff(28); #endif return mprod; } /* Scale a ComplexMatrix */ ComplexMatrix operator*(complex c,const ComplexMatrix &m) { int i; ComplexMatrix sm(m); for (i=0; i < m.nr*m.nc; i++) sm.c[i] *= c; return sm; } /* ComplexMatrix * scalar */ ComplexMatrix operator*(const ComplexMatrix &m, scalar c) { int i; ComplexMatrix sm(m); for (i=0; i < m.nr*m.nc; i++) sm.c[i] *= c; return sm; } ComplexMatrix operator*(real r,const ComplexMatrix &m) { int i; ComplexMatrix sm(m); for (i=0; i < m.nr*m.nc; i++) sm.c[i] *= r; return sm; } /* ComplexMatrix * scalar */ ComplexMatrix operator*(const ComplexMatrix &m, real r) { int i; ComplexMatrix sm(m); for (i=0; i < m.nr*m.nc; i++) sm.c[i] *= r; return sm; } /* diag_matrix*ComplexMatrix */ ComplexMatrix operator*(const diag_matrix &d, const ComplexMatrix &m) { int i,j; if (d.n != m.nr) die("In ComplexMatrix operator*: d.n != m.nr\n"); ComplexMatrix dm(d.n,m.nc); for (i=0; i < d.n; i++) for (j=0; j < m.nc; j++) dm(i,j) = d.c[i]*m(i,j); return dm; } /* ComplexMatrix*diag_matrix */ ComplexMatrix operator*(const ComplexMatrix &m,const diag_matrix &d) { int i,j; /*if (m.nc != d.n) die("In ComplexMatrix operator*: m.nr != d.n\n"); */ if (m.nc != d.n) die("m = %dx%d, d.n = %d\n%sIn ComplexMatrix operator*: m.nr != d.n\nmaddr=%u daddr=%u\n", m.nr,m.nc,d.n, &m, &d); ComplexMatrix md(m.nr,d.n); for (i=0; i < m.nr; i++) for (j=0; j < d.n; j++) md(i,j) = m(i,j)*d.c[j]; return md; } /* Scale a ComplexMatrix in place */ void ComplexMatrix::operator*=(scalar s) { int i; for (i=0; i < nc*nr; i++) c[i] *= s; } /* Accumulate to a ComplexMatrix in place */ void ComplexMatrix::operator+=(const ComplexMatrix &m) { if (nr!=m.nr || nc!=m.nc) die("void ComplexMatrix::operator+=(const ComplexMatrix &m) has size mismatch\n"); int i; for (i=0; i < nc*nr; i++) c[i] += m.c[i]; } /* decumulate to a ComplexMatrix in place */ void ComplexMatrix::operator-=(const ComplexMatrix &m) { if (nr!=m.nr || nc!=m.nc) die("void ComplexMatrix::operator-=(const ComplexMatrix &m) has size mismatch\n"); int i; for (i=0; i < nc*nr; i++) c[i] -= m.c[i]; } /* Diagonalization routine: uses the in-house jacobi.c routien or * uses the FORTRAN77 routines diagrs77_ and diagch77; * however I DO NOT declare them before using them * since then the compiler can't find the routines in the FORTRAN object * files since C++ naming conventions differ from the FORTRAN complier * naming options in the object files (using objdump -t .o to * see the naming of routines */ #if defined DFT_USE_JACOBI void Jacobi(ComplexMatrix &Ain,real *evals,ComplexMatrix &evecs); #elif defined DFT_USE_ESSL || defined DFT_USE_LAPACK /* * C interfaces to the LAPACK/ESSL diagonalization. */ int diagch77_C(real *eigs,complex *evecs,complex *A,int *n); int diagrs77_C(real *eigs,real *evecs,real *A,int *n); #else #error You must specify one of DFT_USE_LAPACK, DFT_USE_ESSL, or DFT_USE_JACOBI #endif void diagonalize_herm(real *eigs, ComplexMatrix &evecs, ComplexMatrix &a, int n) { #ifdef DFT_PROFILING timerOn(27); // Turn on diagonalize_herm timer counterIncr(11); // diagonalization counter #endif // DFT_PROFILING scalar z; if (evecs.nr != evecs.nc) die("evecs is not square in diagonalize_herm()\n\n"); if (a.nr != a.nc) die("a is not square in diagonalize_herm()\n\n"); if (a.nr != evecs.nr) die("Size mismatche of a and evecs in diagonalize_herm()\n\n"); if (a.nr != n) die("n != a.nr in diagonalize_herm()\n\n"); if (!a.hermetian) die("a is not hermetian in diagonalize_herm()\n"); // If we are not using LAPACK or ESSL, then just call the home-grown // Jacobi diagonalizer: slow but steady and does the trick #ifdef DFT_USE_JACOBI Jacobi(a,eigs,evecs); // Otherwise, we use LAPACK or ESSL, so we have to call the C-wrapper // function as well as fiddle with the order of storage in memory. #else int i,j; /* We have to flip the a as FORTRAN access arrays * in reverse order from C */ for (i=0; i < n; i++) for (j=i+1; j < n; j++){ z = a(i,j); a(i,j) = a(j,i); a(j,i) = z; } int status = 0; #if defined SCALAR_IS_COMPLEX status = diagch77_C(eigs,evecs.c,a.c,&n); #elif defined SCALAR_IS_REAL status = diagrs77_C(eigs,evecs.c,a.c,&n); #else #error scalar is neither real nor complex! #endif if (status != 0) die("Status!=0 from calling diag??77_C(): no memory for AP.\n"); /* We have to flip the eigenvector matrices as FORTRAN access arrays * in reverse order from C; we also flip a back to its original form. */ for (i=0; i < n; i++) for (j=i+1; j < n; j++) { z = evecs(i,j); evecs(i,j) = evecs(j,i); evecs(j,i) = z; z = a(i,j); a(i,j) = a(j,i); a(j,i) = z; } #endif // DFT_USE_JACOBI #ifdef DFT_PROFILING timerOff(27); // Turn off diagonalize_herm timer #endif // DFT_PROFILING } /* Return hermetian adjoint of a ComplexMatrix */ ComplexMatrix herm_adjoint(ComplexMatrix &a) { int i,j; ComplexMatrix adag(a.nc,a.nr); for (i=0; i < a.nc; i++) for (j=0; j < a.nr; j++) #if defined SCALAR_IS_COMPLEX adag(i,j) = conj(a(j,i)); #elif defined SCALAR_IS_REAL adag(i,j) = a(j,i); #else #error scalar is neither real nor complex! #endif return adag; } /* returns ^(-1/2) power of a hermetian ComplexMatrix U */ ComplexMatrix Uminusonehalf(ComplexMatrix &U,ComplexMatrix &W,real *u) { int i; if (!U.hermetian) die("In Uminusonehalf: U is not hermetian!\n"); if (U.nr != U.nc) die("In Uminusonehalf: U is not square\n"); if (W.nr != W.nc) die("In Uminusonehalf: W is not square\n"); if (U.nr != W.nr) die("In Uminusonehalf: U.n != W.n\n"); /* Diagonalize U: U = W*u*Wdag where W is unitary and u is diagonal */ diagonalize_herm(u,W,U,U.nr); /* do u -> u^(-1/2) and then multiply out to get U^(-1/2) */ diag_matrix umat(U.nr); ComplexMatrix Umhalf(U.nr,U.nc); for (i=0; i < U.nr; i++) { if (u[i] <= 0.0) die("In Uminusonehalf, eigenvalue %d is %e <= 0.0\n",i,u[i]); umat.c[i] = (real)1.0/sqrt((double)u[i]); } Umhalf = W*(umat*herm_adjoint(W)); Umhalf.hermetian = 1; return Umhalf; } /* write ComplexMatrix in binary format to file fname */ void ComplexMatrix::write(char *fname) { FILE *fp; /* dft_fopen() never returns NULL. that branch dies within the call, * and does the appropriate thing according to whether we are in a serial * or MPI context */ fp = dft_fopen(fname,"w"); dft_fwrite(c,sizeof(scalar),nr*nc,fp); dft_fclose(fp); } void ComplexMatrix::write(FILE *fp) { /* Make sure that if this is called from a non-io node on an MPI * system, then fp is NULL (ie. it was opened with dft_fopen() ) */ dft_fwrite(c,sizeof(scalar),nr*nc,fp); } /* read ComplexMatrix in binary format from file fname */ void ComplexMatrix::read(char *fname) { FILE *fp; fp = dft_fopen(fname,"r"); dft_fread(c,sizeof(scalar),nr*nc,fp); dft_fclose(fp); } void ComplexMatrix::read(FILE *fp) { /* Make sure that if this is called from a non-io node on an MPI * system, then fp is NULL (ie. it was opened with dft_fopen() ) */ dft_fread(c,sizeof(scalar),nr*nc,fp); } /* zero out the ComplexMatrix */ void ComplexMatrix::zero_out(void) { int i; for (i=0; i < nr*nc; i++) c[i] = (scalar)0.0; } /* print us */ void ComplexMatrix::print() { int i, j; for(i = 0; i < nr; i++){ for(j=0; j < nc; j++) dft_log("%5.2f+%5.2fi ", c[nc*i+j].x, c[nc*i+j].y); dft_log("\n"); } } /* print us in exp-format*/ void ComplexMatrix::printe() { int i, j; for(i = 0; i < nr; i++){ for(j=0; j < nc; j++) dft_log("%1.2e+%1.6ei ", c[nc*i+j].x, c[nc*i+j].y); dft_log("\n"); } } /* Allocate/free an array of matrices */ ComplexMatrix * alloc_Matrix_array_old(int nmats,int nrows,int ncols) { int i; ComplexMatrix *M; dft_log(" nmats= %d, nrows=%d ncols=%d\n", nmats, nrows, ncols); dft_log_flush(); M = (ComplexMatrix *)mymalloc(sizeof(ComplexMatrix)*nmats,"alloc_ComplexMatrix_array","M"); for (i=0; i < nmats; i++){ dft_log(" state= %d rows=%d cols=%d addr=%d\n", i,nrows,ncols, &M[i]); dft_log_flush(); M[i].init(nrows,ncols); } return M; } void free_Matrix_array_old(int nmats,ComplexMatrix *M) { int i; for (i=0; i < nmats; i++) M[i].freemem(); myfree(M); } /* Allocate/free an array of matrices */ ComplexMatrix ** alloc_Matrix_array(int nmats,int nrows,int ncols) { int i; ComplexMatrix **M; M = (ComplexMatrix **)mymalloc(sizeof(ComplexMatrix *)*nmats,"alloc_ComplexMatrix_array","M"); for (i=0; i < nmats; i++){ M[i] = new ComplexMatrix(nrows, ncols); } return M; } void free_Matrix_array(int nmats,ComplexMatrix **M) { int i; for (i=0; i < nmats; i++) delete M[i]; myfree(M); } /* Read/write an array of matrices */ void read_ComplexMatrix_array(char *fname,int nmatrices,ComplexMatrix *M) { FILE *filep; int i; filep = dft_fopen(fname,"r"); for (i=0; i < nmatrices; i++) dft_fread(M[i].c,sizeof(scalar),M[i].nr*M[i].nc,filep); dft_fclose(filep); } void write_ComplexMatrix_array(char *fname,int nmatrices,ComplexMatrix *M) { FILE *filep; int i; filep = dft_fopen(fname,"w"); for (i=0; i < nmatrices; i++) dft_fwrite(M[i].c,sizeof(scalar),M[i].nr*M[i].nc,filep); dft_fclose(filep); } /* * The ComplexMatrix calculates the linear operator Q(G) with a ComplexMatrix of * eigenvectors W and eigenvalues mu (these come from the overlap * operators where U = Wdag*mu*W). This is needed for the case of * variable fillings. * The definition of Q(G) is via: * * (Wdag*Q(G)*W)_{ij} = (Wdag*G*W)_{ij}/(sqrt(mu_i)+sqrt(mu_j)) * */ ComplexMatrix Q(const ComplexMatrix &G,ComplexMatrix &W,real *mu) { int i,j; if (G.nr != G.nc) die("Q(G) called with G not square\n"); if (W.nr != W.nc) die("Q(G) called with W not square\n"); if (G.nr != W.nr) die("Q(G) called with G.nr != W.nr\n"); /* Take G into the diagonal basis: i.e. W^*G*W, where ^ is herm. adjoint */ ComplexMatrix QG(G.nr,G.nc),Wdag(W.nr,W.nc); Wdag = herm_adjoint(W); QG = Wdag*G*W; /* Divide (i,j) entry by sqrt(mu_i)+sqrt(mu_j) */ for (i=0; i < QG.nr; i++) for (j=0; j < QG.nc; j++) { if (mu[i] < 0.0 || mu[j] < 0.0) die("Q(G) called with mu[%d]=%g mu[%d]=%g\n",i,mu[i],j,mu[j]); QG(i,j) = QG(i,j)/(sqrt(mu[i])+sqrt(mu[j])); } /* Multiply by W and W^ to bring back to standard basis */ QG = W*QG*Wdag; /* Done! */ QG.hermetian = G.hermetian; return QG; } /* Returns trace of ComplexMatrix */ scalar trace(const ComplexMatrix &m) { scalar tr; int i; if (m.nr != m.nc) die("trace(ComplexMatrix) called with non-square ComplexMatrix\n"); tr = 0.0; for (i=0; i < m.nr; i++) tr += m(i,i); return tr; } /* Returns diagonal of ComplexMatrix as a vector */ diag_matrix diag(const ComplexMatrix &m) { if (m.nr != m.nc) die("diag(ComplexMatrix) called with non-square ComplexMatrix\n"); diag_matrix d(m.nr); int i; for (i=0; i < m.nr; i++) d.c[i] = m(i,i); return d; } /* * The ComplexMatrix calculates the linear operator R(A) with a ComplexMatrix of * eigenvectors Z and eigenvalues beta (these come from the ComplexMatrix * B where B=Z*beta*Zdag and V=exp(iB) is the subspace rotation ComplexMatrix. * This is needed for the case of variable fillings. * The definition of R(A) is via: * * (Zdag*R(A)*Z)_{nm} = (Zdag*A*Z)_{nm}* * ( (exp(i*beta_m-i*beta_n)-1)/(beta_m-beta_n) ) * * The n=m case just gives i for the ratio. */ ComplexMatrix R(const ComplexMatrix &A,ComplexMatrix &Z,real *beta) { int i,j; if (A.nr != A.nc) die("R(A) called with G not square\n"); if (Z.nr != Z.nc) die("R(A) called with Z not square\n"); if (A.nr != Z.nr) die("R(A) called with A.nr != Z.nr\n"); /* Take A into the diagonal basis: i.e. Z^*A*Z, where ^ is herm. adjoint */ ComplexMatrix RA(A.nr,A.nc),Zdag(Z.nr,Z.nc); Zdag = herm_adjoint(Z); RA = Zdag*A*Z; /* Multiply (n,m) entry by (exp(i*beta_n-i*beta_m)-1)/(beta_n-beta_m) */ for (i=0; i < RA.nr; i++) for (j=0; j < RA.nc; j++) { complex ratio; real deltaji = beta[j]-beta[i]; /* Compute ratio=(exp(i*beta[n]-i*beta[m])-1)/(beta[n]-beta[m] */ if ( fabs(deltaji) < 1.0e-13*(fabs(beta[i])+fabs(beta[j])) ) { ratio.x = 0.0; ratio.y = 1.0; } else { complex z; z.x = 0.0; z.y = deltaji; ratio = (exp(z)-1.0)/deltaji; } /* Multiply RA by ratio */ RA(i,j) = RA(i,j)*ratio; } /* Multiply by Z and Z^ to bring back to standard basis */ RA = Z*RA*Zdag; /* Done! */ RA.hermetian = 0; return RA; } /* Does mout += s*min */ void scale_accumulate(scalar s,ComplexMatrix &min,ComplexMatrix &mout) { if (min.nr != mout.nr || min.nc != mout.nc) die("Incompatible sizes in scale_accumulate(scalar,ComplexMatrix,ComplexMatrix)\n"); for (int j=0; j < min.nr*min.nc; j++) mout.c[j] += s*min.c[j]; } /* Does mout[i] += s*min[i] */ void scale_accumulate(int nmat,scalar s,ComplexMatrix **min,ComplexMatrix **mout) { for (int i=0; i < nmat; i++) scale_accumulate(s,*min[i],*mout[i]); } /* Does mout = s1*m1 + s2*m2 */ void scaled_sum(scalar s1,ComplexMatrix &m1,scalar s2,ComplexMatrix &m2,ComplexMatrix &mout) { if (m1.nr!=m2.nr || m1.nr!=mout.nr || m1.nc!=m2.nc || m1.nc!=mout.nc) die("Incompatible ComplexMatrix sizes in scaled_sum(s1,m1,s2,m2,mout)\n"); for (int j=0; j < m1.nr*m1.nc; j++) mout.c[j] = s1*m1.c[j] + s2*m2.c[j]; } /* Does mout[i] = s1*m1[i] + s2*m2[i] */ void scaled_sum(int nmat,scalar s1,ComplexMatrix **m1, scalar s2,ComplexMatrix **m2,ComplexMatrix **mout) { for (int i=0; i < nmat; i++) scaled_sum(s1,*m1[i],s2,*m2[i],*mout[i]); } /* returns absolute square magnitude of sum of all ComplexMatrix elements */ real abs2(const ComplexMatrix &m) { real r(0); scalar z; for (int i=0; i < m.nr*m.nc; i++) { z = m.c[i]; r += z.x*z.x + z.y*z.y; } return r; } /* Same for an array of matrices */ real abs2(int nmat,ComplexMatrix **m) { int i; real r; r = 0.0; for (i=0; i < nmat; i++) r += abs2(*m[i]); return r; } /* * returns the "dot" product of two matrices...i.e. * the sum over all the products of their respective entries: * dot(A,B) = sum_{ij} { A(i,j)*B(i,j) } */ scalar dot(ComplexMatrix &A,ComplexMatrix &B) { if (A.nr!=B.nr || A.nc!=B.nc) die("In dot(ComplexMatrix,ComplexMatrix) the sizes don't match\n"); scalar d; int i; d = 0.0; for (i=0; i < A.nr*A.nc; i++) #if defined SCALAR_IS_COMPLEX d += conj(A.c[i])*B.c[i]; #elif defined SCALAR_IS_REAL #else d += A.c[i]*B.c[i]; #error scalar is neither real nor complex! #endif return d; } // Same for an array of matrices scalar dot(int nmat,ComplexMatrix **A,ComplexMatrix **B) { scalar d(0); for (int i=0; i < nmat; i++) d += dot(*A[i],*B[i]); return d; }