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

/*
 * Gabor Csanyi 8/1/2001
 *
 *
 * the rArray class holds a simple array of real numbers
 *
 */

#include "header.h"




RealArray::RealArray(int length):
  ndata(length)
{
}

void
RealArray::init(int length)
{
  ndata = length;
  if(ndata > 0){
    d = (real *)malloc(ndata*sizeof(real));
    //dft_log("Allocated RealArray\n");
  }
  else
    die("Can't initiliaze RealArray.init() with length < 0\n");
}


//////////////////////////////////////
//                                  //
// Operator functions               //
//                                  //
//////////////////////////////////////




/* Assignment:  nonstandard in that it returns void.  To make it standard,
 * replace void -> RealArray and uncomment the return *this; */
void
RealArray::operator=(const RealArray &rd)
{
  dft_log("RealArray::operator=\n");
  if (ndata != rd.ndata)
    die("Size mismatch in RealArray::operator=()\n");

  for (int i=0; i < ndata; i++)
    d[i] = rd.d[i];
  /* return *this; */
}

// Assignment of scalar:  all entries set to r
inline void
RealArray::operator=(real r)
{
  dft_log("RealArray::operator=(real r)\n");

  for (int i=0; i < ndata; i++)
    d[i] = r;

}




/* Add two RealArray */
RealArray
RealArray::operator+(const RealArray &rd)
{

  dft_log("RealArray::operator+\n");

  if (ndata != rd.ndata)
    die("Size mismatch in RealArray::operator+()\n");
  
  RealArray rd2(*this);

  for(int i = 0; i < ndata; i++)
    rd2.d[i] += rd.d[i]; 
  return rd2;
}


/* Accumulate sum of RealArray */
void
RealArray::operator+=(const RealArray &rd)
{

  dft_log("RealArray::operator+=\n");
  if (ndata != rd.ndata)
    die("Size mismatch in RealArray::operator+=()\n");

  for (int i=0; i < ndata; i++)
      d[i] += rd.d[i];
}


/* Subtract two RealArray */
RealArray
RealArray::operator-(const RealArray &rd)
{
  dft_log("RealArray::operator-\n");
  if (ndata != rd.ndata)
    die("Size mismatch in RealArray::operator-()\n");
  
  RealArray rd2(*this);

  for(int i = 0; i < ndata; i++)
    rd2.d[i] -= rd.d[i]; 
  return rd2;

}

/* Accumulate difference of column_bundles */
void
RealArray::operator-=(const RealArray &rd)
{
  dft_log("RealArray::operator-=\n");
  if (ndata != rd.ndata)
    die("Size mismatch in RealArray::operator-=()\n");

  for (int i=0; i < ndata; i++)
    d[i] -= rd.d[i];
    
}

/* Scale a RealArray rd by real r */
RealArray
operator*(real r,const RealArray &rd)
{
  dft_log("RealArray::operator*\n");
  
  RealArray rd2(rd);
  for(int i = 0; i < rd.ndata; i++)
    rd2.d[i] *= r;
  return rd2;
}

/* Scale a RealArray by real r */
RealArray
RealArray::operator*(real r)
{
  dft_log("RealArray::operator*\n");
  RealArray rd2(*this);
  for(int i = 0; i < ndata; i++)
    rd2.d[i] *= r;
  return rd2;
}

/* Scale a RealArray rd by complex number c */
ComplexArray
operator*(complex c,const RealArray &rd)
{
  dft_log("RealArray::operator*\n");
  
  ComplexArray cd(rd.ndata);
  for(int i =0; i < rd.ndata; i++){
    cd.d[i].x = c.x*rd.d[i];
    cd.d[i].y = c.y*rd.d[i];
  }
  return cd;
}

/* Scale a RealArray by a complex number c */
ComplexArray
RealArray::operator*(complex c)
{
  ComplexArray cd(ndata);
  for(int i = 0; i < ndata; i++){
    cd.d[i].x = d[i]*c.x;
    cd.d[i].y = d[i]*c.y;
  }
  return cd;
}


/* Scale a RealArray real r in place */
void
RealArray::operator*=(real r)
{
  dft_log("RealArray::operator*=\n");

  for (int i=0; i < ndata; i++)
    d[i] *= r;
}


/////////////////////////////////////////////
//                                         //
//  MEMBER FUNCTIONS:                      //
//                                         //
/////////////////////////////////////////////
// zeroes out the whole array
void
RealArray::zero_out(void)
{
  for(int i=0;i < ndata; i++)
      d[i] = (real)0.0;
}

/* Negates all the entries in the array */
void
RealArray::negate(void)
{
  dft_log("RealArray::negate()\n");

  for(int i=0;i < ndata; i++)
    d[i] = -d[i];

}


/* Randomizes w/ uniform distribution in [-.5, .5]*/
void
RealArray::randomize(void)
{
  //dft_log("ComplexArray::negate()\n");

  for(int i=0;i < ndata; i++)
    d[i] = rand()/(RAND_MAX+1.) -.5;

}




// Sum of absolute squares of all elements in rd
real abs2(const RealArray &rd)
{
  real rdrd(0);
  real r;

  for(int i = 0; i < rd.ndata; i++){
    r = rd.d[i];
    rdrd += r*r;
  }

  return rdrd;
}


/* Take "dot-product" of two RealArray:  sum the diagonals of rd1^rd2 */
real
dot(const RealArray &rd1,const RealArray &rd2)
{
  dft_log("RealArray::dot()\n");
  int i;
  real dot12;

  if (rd1.ndata != rd2.ndata)
    die("rd1.ndata != rd2.ndata in dot_RealArray()\n");


  dot12 = 0.0;
  for (i=0; i < rd1.ndata; i++)
    dot12 += rd1.d[i]*rd2.d[i];

  return dot12;
}


/* Does rd2 += r * rd1 */
void
scale_accumulate(real r,const RealArray &rd1, RealArray &rd2)
{
  dft_log("RealArray::scale_accumulate()\n");
  for (int i=0; i < rd1.ndata; i++)
      rd2.d[i] += r * rd1.d[i];
}

/* Does rd3 = r1*rd1 + r2*rd2 */
void
scaled_sum(real r1,RealArray &rd1,
	   real r2,RealArray &rd2,
	   RealArray &rd3)
{
  dft_log("RealArray::scaled_sum()\n");
  if (rd1.ndata != rd2.ndata || rd1.ndata != rd3.ndata)
    die("Incomaptible sizes in scaled_sum(r1,rd1,r2,rd2,rd3)");

  for (int i=0; i < rd1.ndata; i++)
      rd3.d[i] = r1*rd1.d[i] + r2*rd2.d[i];
}