/*
    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 <file>.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;
}