/*
    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.
*/

/*************************************************************************
 *
 * Defines the  Wavelet_Complex_Column class
 *
 ************************************************************************/
#include "header.h"

WL_ComplexColumn::WL_ComplexColumn():fat(TRUE)
{
//  fat = TRUE;
  length = 0;
  basis_spec = NULL;
  basis = NULL;
  representation = REALSPACE; //default

  embedding = SPARSE;
}


WL_ComplexColumn::WL_ComplexColumn(WL_Basis *b, char s[10])
{
  if(!strcmp(s,"fat"))
    fat = TRUE;
  else if (!strcmp(s,"slim"))
    fat = FALSE;
  else
    die("WL_Column(): Identity crisis - don't know if i'm slim or fat\n");

  basis = b;

  // We must be able to accept a null pointer!
  if(b != NULL) {
      basis_spec = b->basis_spec;

      // ipd: get rid of these - make everything  like PW_Column
      length = b->nbasis;
      // Grab the quantum numbers from the basis.
//      kvec = b->qnum->kvec;
    }

  // These are defaults.
  representation = REALSPACE;
  embedding = SPARSE;
}  


WL_ComplexColumn::WL_ComplexColumn(const WL_ComplexColumn &col,
				   char s[10])
{
  if(!strcmp(s,"fat"))
    fat = TRUE;
  else if (!strcmp(s,"slim"))
    fat = FALSE;
  else
    die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");

  length = col.length;
  basis_spec = col.basis_spec;
  basis = col.basis;
  representation = col.representation;
  embedding = col.embedding;

  //can we get rid of these?
//  kvec = col.kvec;
//  weight = col.weight;

  if(fat) {
    data.init(length);
    data = col.data;
  }
  else
    die("You can't copy a slim column, fatso!\n");
}

// destructor
WL_ComplexColumn::~WL_ComplexColumn()
{
  if(fat)
    data.free();
}

void
WL_ComplexColumn::init_wavefunction(WL_Basis &b)
{
  length = b.nbasis;
  basis_spec = b.basis_spec;
  basis = &b;
  representation = COEFFSPACE;
  embedding = SPARSE;
}


// We'll decide later if this should be different.
void
WL_ComplexColumn::init_scalarfield(WL_Basis &b, int space)
{
  // ipd: since *qnum is part of basis we want  scalarfields and wfs to
  // have different basis  - think how to fix it
  basis = &b; 
  basis_spec = b.basis_spec;
  representation = space;
  
  //ipd:replace this w/ an input flag (in basis/basis_spec)
  if(basis_spec->use_dense){
    embedding = DENSE;
    length = 8*b.nbasis;
  } else {
    embedding = SPARSE;
    length = b.nbasis;  
  }


  if(fat)
    data.init(length);
  else
    die("You cannot initialize a slim column to contain a scalar field\n");
}


void
WL_ComplexColumn::read(char *realfile, char *imagfile)
{
  FILE *fpin;
  int ii;
  
  double *databuf;

  databuf=(double *)mymalloc(length*sizeof(double),
                             "databuf", "WL_ComplexColumn::read(...)");
  
  fpin = dft_fopen(realfile,"r");
  
  // Store the data into databuf.
  if(dft_fread(databuf, sizeof(double), length,  fpin) != length)
    die("Error reading %d items for column initialization!",length);
  dft_fclose(fpin);
  
  // And store it in the real parts.
  for(ii=0; ii<length; ii++)
    data.d[ii].x = databuf[ii];

  // Now do the same thing for the complex parts.
  fpin = dft_fopen(imagfile,"r");
  
  // Store the data into databuf.
  if(dft_fread(databuf, sizeof(double), length,  fpin) != length)
    die("Error reading %d items for column initialization!",length);
  dft_fclose(fpin);

  // And store it in the real parts.
  for(ii=0; ii<length; ii++)
    data.d[ii].y = databuf[ii];
  myfree(databuf);
}


// Initializes a column with real data from an ascii file.
void
WL_ComplexColumn::read_text(char *filename)
{
  FILE *fpin;
  int ii;
  double databuf;

  fpin = dft_fopen(filename,"r");

  for(ii=0; ii<length; ii++)
    {
      fscanf(fpin,"%lf",&databuf);
      data.d[ii].x = databuf;
      data.d[ii].y = 0.0;
    }

  dft_fclose(fpin);
}


void
WL_ComplexColumn::print()
{
  int j;

  for (j = 0; j < length; j++)
    dft_log("%f\t%f\n", data.d[j].x, data.d[j].y);
}

void
WL_ComplexColumn::read(char *filename)
{
  data.read(filename);
}

void
WL_ComplexColumn::write(char *filename)
{
  data.write(filename);
}

char *
WL_ComplexColumn::getrepresentation() const
{
  static char realspacename[] = "RealSpace";
  static char coeffspacename[] = "CoeffSpace";

  switch(representation){
  case REALSPACE:
    return realspacename;
  case COEFFSPACE:
    return coeffspacename;
  default:
    die("Unknown representation: %d\n", representation);
  }

  return (NULL);
}

char *
WL_ComplexColumn::getembedding() const
{
  static char sparsename[] = "Sparse";
  static char densename[] = "Dense";
  
  switch(embedding){
    case SPARSE:
      return sparsename;
    case DENSE:
      return densename;
    default:
      die("Unknown embedding: %d\n", embedding);
  }

  return (NULL);
}


// LINEAR ALGEBRA - just call the array algebra

void
WL_ComplexColumn::operator=(const WL_ComplexColumn &col)
{
  if(length != col.length)
    die("length != col.length in WL_ComplexColumn::operator=()\n");
  representation=col.representation;

  data = col.data;
}


void
WL_ComplexColumn::operator-=(const WL_ComplexColumn &col)
{
  if(length != col.length)
    die("length != col.length in WL_ComplexColumn::operator=()\n");
//  if(representation != col.representation)
//    die("representation != col.repr. in WL_ComplexColumn::operator-=()\n");

  data -= col.data;
}

void 
WL_ComplexColumn::operator+=(const WL_ComplexColumn &col)
{
  if(length != col.length)
    die("length != col.length in WL_ComplexColumn::operator=()\n");

  data += col.data;
}

void 
WL_ComplexColumn::operator*=(const real r)
{
  data*=r;
}

void 
WL_ComplexColumn::operator*=(const complex c)
{
  data *= c;
}

void
WL_ComplexColumn::operator*=(const WL_ComplexColumn &col)
{
  data *= col.data;
}

// NOTE:  Do we ever use this function?  It doesn't seem necessary.
//     -- mhe 3/8/02

// Realspace dot product: skip on ghosts
//
//  complex 
//  dot(const WL_ComplexColumn &col1, const WL_ComplexColumn &col2)
//  {

//  }

  
// the carrot operator is a little trickier this time - b/c of the grid
// structures
//ipd: this changes the input - i don't think we want it!!!
complex
WL_ComplexColumn::operator^(const WL_ComplexColumn &col)
{
  /*
  if(representation != col.representation)
    die("dude, u'r trying to do %s^%s\n",
        representation, col.representation);
  if(embedding!=col.embedding)
    die("trying %s^%s what r u thinking?!?\n",
        embedding, col.embedding);
  */
  switch (representation){
    case COEFFSPACE:
      return dot(data, col.data);
    case REALSPACE:{
      //ipd: this changes the input, so i'll die
      die("Are you sure you want to perform col^col REALSPACE??\n");
      Grid *grid = mkgridnew(basis_spec->gridspec,NULL);
      GetRe(grid,*this); zeroghostrec(grid); PutRe(grid,*this);
      GetIm(grid,*this); zeroghostrec(grid); PutIm(grid,*this);
      killgridnew(grid);
      return dot(data, col.data);
    }
    default:
      die("Unknown rerpesentation: %d\n",representation);
  }

  return 0.;
}


//ipd: would be nicer to do it in one loop - just change adjustdown/zeroghost
void
WL_ComplexColumn::randomize()
{
  Gridspec *gridspec;
    
    // choose the appropriate <sparse|dense> gridspec
    switch(embedding){
      case SPARSE:
        gridspec = basis_spec->gridspec;
        break;
      case DENSE:
        gridspec = basis_spec->gridspec_dense;
        break;
      default:
        die("randomize(): Unknown  input embedding %d\n",embedding);
    }

  Grid *grid = mkgridnew(gridspec,NULL);

  // Do the real part.
  for(int i=0; i<(grid->nxyzwholegrid); i++)
    grid->dat[i] = rand()/(RAND_MAX+1.0)-0.5;

  switch(representation)
    {
    case REALSPACE:
      // Realspace convention is:  ghost=parent
      adjustdown(grid);
      PutRe(grid,*this);
      break;
    case COEFFSPACE:
      // Coeffspace convention is:  ghost=0
      zeroghostrec(grid);
      PutRe(grid,*this);
      break;
    default:
      dft_log("\nInvalid representation in WL_ComplexColumn::randomize()!\n");
      break;
    }

  // Do the imaginary part.
  for(int i=0; i<(grid->nxyzwholegrid); i++)
    grid->dat[i] = rand()/(RAND_MAX+1.0)-0.5;

  switch(representation)
    {
    case REALSPACE:
      // Realspace convention is:  ghost=parent
      adjustdown(grid);
      PutIm(grid,*this);
      break;
    case COEFFSPACE:
      // Coeffspace convention is:  ghost=0
      zeroghostrec(grid);
      PutIm(grid,*this);
      break;
    default:
      dft_log("\nInvalid representation in WL_ComplexColumn::randomize()!\n");
      break;
    }

  killgridnew(grid);
}

void
WL_ComplexColumn::setones()
{
  int ii;
  
  for(ii=0; ii<length; ii++){
    data.d[ii].x = 1.0;
    data.d[ii].y = 0.0;
  }

  if(representation == COEFFSPACE)
    apply_J(*this,*this);

  return;
}

void
WL_ComplexColumn::setmodones()
{
  int ii;
  real theta = rand(); // just a random phase
  complex one(cos(theta), sin(theta));

  for(ii=0; ii<length; ii++)
    data.d[ii] = one;

  //WL_ComplexColumn out(*this); // create a copy column
  if(this->representation == COEFFSPACE)
    apply_J(*this, *this);

  return;
}

// BASIS SPECIFIC


// Returns the real part of the data as a Grid structure.
void
GetRe(Grid *out, const WL_ComplexColumn &col)
{
#ifdef DFT_PROFILING
  timerOn(60);
  counterIncr(25);
#endif // DFT_PROFILING

  // Consistency checks.
  if(out == NULL)
    die("Passed NULL Grid to GetRe!\n");
  if(out->nxyzwholegrid != col.length)
    die("Grid structure doesn't match Column length in GetRe!%d %d\n",
        out->nxyzwholegrid, col.length);

  // In mkgridnew, the data is allocated as a single block (instead of 
  // level-by-level.  Thus, we should just need to copy the data as a vector.
  int ii;
  for(ii=0;ii<col.length;ii++)
    out->dat[ii] = col.data.d[ii].x;

#ifdef DFT_PROFILING
  timerOff(60); 
#endif // DFT_PROFILING

  return;
}

void
GetIm(Grid *out, const WL_ComplexColumn &col)
{
#ifdef DFT_PROFILING
  timerOn(60);
  counterIncr(25);
#endif // DFT_PROFILING

  // Consistency checks.
  if(out == NULL)
    die("Passed NULL Grid to GetIm!\n");
  if(out->nxyzwholegrid != col.length)
    die("Grid structure doesn't match Column length in GetIm!\n");

  // In mkgridnew, the data is allocated as a single block (instead of 
  // level-by-level.  Thus, we should just need to copy the data as a vector.
  int ii;
  for(ii=0;ii<col.length;ii++)
    out->dat[ii] = col.data.d[ii].y;

#ifdef DFT_PROFILING
  timerOff(60); 
#endif // DFT_PROFILING

  return;
}

void
PutRe(const Grid *in, WL_ComplexColumn &col)
{
#ifdef DFT_PROFILING
  timerOn(60);
  counterIncr(25);
#endif // DFT_PROFILING

  // Again, if we assume that the grid storage is allocated as a block, this is
  // pretty easy.

  int ii;
  for(ii=0;ii<col.length;ii++)
    col.data.d[ii].x = in->dat[ii];

#ifdef DFT_PROFILING
  timerOff(60); 
#endif // DFT_PROFILING

  return;
}

void 
PutIm(const Grid *in, WL_ComplexColumn &col)
{
#ifdef DFT_PROFILING
  timerOn(60);
  counterIncr(25);
#endif // DFT_PROFILING

  int ii;
  for(ii=0;ii<col.length;ii++)
    col.data.d[ii].y = in->dat[ii];

#ifdef DFT_PROFILING
  timerOff(60); 
#endif // DFT_PROFILING

  return;
}

// adds the square of the absolute value
void
add_scale_abs2(const complex &c, const  WL_ComplexColumn & in,
                    WL_ComplexColumn &out)
{
  add_scale_abs2(c, in.data, out.data);
}

WL_ComplexColumn
abs2(WL_ComplexColumn &in)
{
  int i;
  WL_ComplexColumn out(in);

  for(i=0; i<in.length; i++)
    out.data.d[i] = abs2(in.data.d[i]);

  return out;
}

WL_ComplexColumn
pointwise_mult(const WL_ComplexColumn &a, const WL_ComplexColumn &b)
{
  if(a.length != b.length)
    die("a.length != b.length in pointwise_mult!\n");

  WL_ComplexColumn ab(a);
  int i;

  for(i=0; i<a.length; i++)
    ab.data.d[i] *= b.data.d[i];

  return ab;
}

real
dot(const WL_ComplexColumn &d1, const WL_ComplexColumn &d2)
{
  // Copy code from ComplexArray::dot, but assume d.y = 0.
  int i;
  real dot12;

  if (d1.data.ndata != d2.data.ndata)
    die("d1.data.ndata != d2.data.ndata in real dot_WL_ComplexColumn()\n");


  dot12 = 0.0;
  for (i=0; i < d1.data.ndata; i++)
    dot12 += d1.data.d[i].x*d2.data.d[i].x;

  return dot12;
}