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

#include "header.h"

/* Returns (I-P)Y...the complement projector of P=O*C*Cdag where
 * Cdag is the Hermetian adjoint of C and C = Y*U(-1/2).
 * Thus the routine returns Pbar*Yin = (I-O*C*Cdag)*Yin. */
ColumnBundle
Pbar(const ColumnBundle &C,const ColumnBundle &Y)
{
  ColumnBundle result(Y);

  apply_Pbar(C, Y, result);

  return result;
}

/* Less memory intensive version of Pbar() above */
void
apply_Pbar(const ColumnBundle &C,
	   const ColumnBundle &Y,
	   ColumnBundle &PbarY)
{
#ifdef DFT_PROFILING
  timerOn(30); // turn on Pbar timer
#endif // DFT_PROFILING

  do_ColumnBundle_matrix_mult(C,-1.*(C^Y),PbarY,0);
  O(PbarY,PbarY);
  //PbarY.negate();
  PbarY += Y;

#ifdef DFT_PROFILING
  timerOff(30); // turn off Pbar timer
#endif // DFT_PROFILING
}


//
// Applies the local potential Vscloc to Y and accumulates
// the result into HspY.  It does this on the columns start_col
// to start_col+n_cols_todo-1 of Y.
//
// the last two arguments are necessary for reuse in case of threads
//
void
apply_Vsc(const ColumnBundle &Y, ColumnBundle &HspY,
          const RealSpaceScalarFieldColumn &V, 
          int start_col,int n_cols_todo)
{
  // Local part of self-consistent potential: Idag*Diag(Vsc)*I*Y
  //
  // HspY(:,i) += Idag*Diag(Vscloc)*I*Y(:,i) (i'th column of both sides)
  //
  int i;
  // Holds I of Ycol
  RealSpaceScalarFieldColumn IYcol(V);
  //ipd: since the columns of Y are slim you can't use the copy constructor
  //     in its current shape
  //CoeffSpaceWavefunctionColumn IdagVIYcol(*Y.col[start_col]);
  CoeffSpaceWavefunctionColumn IdagVIYcol(*Y.col[start_col],"fat");

  
  for (i = start_col; i < start_col+n_cols_todo; i++) {
    // Apply I to Ycol
    apply_I(*Y.col[i],IYcol);
    
    IYcol *= V;
    
    //ipd: this will generate 2 extra copies of a column
    apply_Idag(IYcol,IdagVIYcol);
    //ipd: avoid one of them by
    // apply_Idag(IYcol, IdagVIYcol);
    //ipd: save a column copy inside Idag by upgrading Idag and adding
    // Idag(IYcol, IYcol);
    // IYcol.map(IdagVIYcol);
      
    // Accumulate into the i'th column of HspY
    *(HspY.col[i]) += IdagVIYcol;
  }

}

#ifdef DFT_THREAD
//
// Threads run this routine.  It just applies Vsc to the columns
// this thread is working on by calling the above routine.
//
void *
apply_Vsc_thread(void *arg)
{
  // Decode the data passed to the thread
  dft_thread_data *data = (dft_thread_data *)arg;

  ColumnBundle *Y = (ColumnBundle *)data->p1;
  ColumnBundle *HspY = (ColumnBundle *)data->p2;
  RealSpaceScalarFieldColumn *V = (RealSpaceScalarFieldColumn *)data->p3;
  int start_col = data->start;
  int n_cols_todo = data->n;

  // Do the work
  apply_Vsc(*Y,*HspY,*V,start_col,n_cols_todo);

  // Free data-passing object and exit
  myfree(arg);
  return NULL;
}
#endif // DFT_THREAD


/*
 * Multiply ColumnBundle by single-particle Hamiltonian:
 *
 * Hsp = -0.5*L + Idag*Diag(Vscloc)*I
 *              + sum(ions,l,m,...) { Vnl*M*Vnl^ }
 *
 * i.e. Kinetic operator, local self-consistent potential, and non-local
 * potential respectively.
 *
 * Takes Vscloc and Everything is a hack for Vnl
 * It actaully calls the routine below it...
 */
ColumnBundle
Hsp(const ColumnBundle &Y, const RealSpaceScalarFieldColumn &Vscloc,
    Everything &e)
{
  /* Holds final result */
  ColumnBundle HspY(Y);

  /* Do the work!! */
  apply_Hsp(Y,HspY,Vscloc,e);

  return HspY;
}

//
// Actually does the work of applying Hsp onto Y.
//
void
apply_Hsp(const ColumnBundle &Y, ColumnBundle &HspY,
          const  RealSpaceScalarFieldColumn &Vscloc, Everything &e)
{
#ifdef DFT_PROFILING
  timerOn(26);   // Turn on apply_Hsp timer
#endif // DFT_PROFILING
  // Kinetic part: HspY = -0.5*L(Y);
  L(Y,HspY);

  HspY *= -0.5;
  
#ifdef DFT_PROFILING
  timerOn(42);  // Local Hsp application timer
#endif
  
  // Local part of self-consistent potential: Idag*Diag(Vsc)*I*Y 


#ifdef DFT_THREAD
  // Threads!  Distribute the work
  dft_call_threads(Y.my_ncols,
                   (void *)&Y,(void *)&HspY,(void *)&Vscloc,
                   NULL,NULL,
                   0,0,0,0,0,0,
                   apply_Vsc_thread);
#else
  // Serial mode:  do the work in one big call
  apply_Vsc(Y,HspY,Vscloc,0,Y.my_ncols);

//  (Y^Y).print(); exit(1);
#endif

#ifdef DFT_PROFILING
  timerOff(42);  // Local Hsp application timer
  timerOn(43); // Nonlocal timer
#endif

  // This is a slightly generalized version of a non-local pseudopotential
  // operator.  It applies Vnl to Y, and adds the result to HspY.
  //ipd: we have to decide about this guy - what do we need to pass
  // but whatever it is the bloch state could know about it
  accumulate_Vnl(Y,HspY,e);
  //(HspY^HspY).print(); exit(1);

#ifdef DFT_PROFILING
  timerOff(43);  // nonlocal timer
#endif

#ifdef DFT_PROFILING
  timerOff(26);   // Turn off apply_Hsp timer
#endif // DFT_PROFILING
}