PSDWavefunctionStepper.C 3.51 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19 20 21
// PSDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "PSDWavefunctionStepper.h"
#include "Wavefunction.h"
#include "SlaterDet.h"
22
#include "Preconditioner.h"
23 24 25 26
#include <iostream>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
27
PSDWavefunctionStepper::PSDWavefunctionStepper(Wavefunction& wf,
28 29 30
  Preconditioner& prec, TimerMap& tmap) : WavefunctionStepper(wf,tmap),
  prec_(prec)
{}
Francois Gygi committed
31 32 33

////////////////////////////////////////////////////////////////////////////////
PSDWavefunctionStepper::~PSDWavefunctionStepper(void)
34
{}
35 36 37 38 39 40 41 42

////////////////////////////////////////////////////////////////////////////////
void PSDWavefunctionStepper::update(Wavefunction& dwf)
{
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
    {
Francois Gygi committed
43 44 45 46 47 48 49
      // compute A = V^T H V  and descent direction HV - VA
      tmap_["psd_residual"].start();
      if ( wf_.sd(ispin,ikp)->basis().real() )
      {
        // proxy real matrices c, cp
        DoubleMatrix c(wf_.sd(ispin,ikp)->c());
        DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
50

Francois Gygi committed
51
        DoubleMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
52

Francois Gygi committed
53 54 55 56
        // factor 2.0 in next line: G and -G
        a.gemm('t','n',2.0,c,cp,0.0);
        // rank-1 update correction
        a.ger(-1.0,c,0,cp,0);
57

Francois Gygi committed
58 59 60 61 62
        // cp = cp - c * a
        cp.gemm('n','n',-1.0,c,a,1.0);
      }
      else
      {
Francois Gygi committed
63 64
        ComplexMatrix& c = wf_.sd(ispin,ikp)->c();
        ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
Francois Gygi committed
65 66 67 68 69 70
        ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
        a.gemm('c','n',1.0,c,cp,0.0);
        // cp = cp - c * a
        cp.gemm('n','n',-1.0,c,a,1.0);
      }
      tmap_["psd_residual"].stop();
Francois Gygi committed
71 72
    }
  }
73

Francois Gygi committed
74
  // dwf.sd->c() now contains the descent direction (HV-VA)
75

76 77
  // update preconditioner
  prec_.update(wf_);
78

Francois Gygi committed
79 80 81 82 83
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
    {
      tmap_["psd_update_wf"].start();
Francois Gygi committed
84 85 86 87 88 89 90 91 92 93 94
      double* coeff = (double*) wf_.sd(ispin,ikp)->c().valptr();
      const double* dcoeff =
        (const double*) dwf.sd(ispin,ikp)->c().cvalptr();
      const int mloc = wf_.sd(ispin,ikp)->c().mloc();
      const int ngwl = wf_.sd(ispin,ikp)->basis().localsize();
      const int nloc = wf_.sd(ispin,ikp)->c().nloc();
      for ( int n = 0; n < nloc; n++ )
      {
        // note: double mloc length for complex<double> indices
        double* c = &coeff[2*mloc*n];
        const double* dc = &dcoeff[2*mloc*n];
Francois Gygi committed
95

Francois Gygi committed
96 97
        for ( int i = 0; i < ngwl; i++ )
        {
98
          const double fac = prec_.diag(ispin,ikp,n,i);
Francois Gygi committed
99 100 101 102
          const double delta_re = fac * dc[2*i];
          const double delta_im = fac * dc[2*i+1];
          c[2*i]   -= delta_re;
          c[2*i+1] -= delta_im;
103 104
        }
      }
Francois Gygi committed
105 106 107 108 109
      tmap_["psd_update_wf"].stop();

      tmap_["gram"].start();
      wf_.sd(ispin,ikp)->gram();
      tmap_["gram"].stop();
110 111 112
    }
  }
}