PSDAWavefunctionStepper.C 7.03 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
// PSDAWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////
27
PSDAWavefunctionStepper::PSDAWavefunctionStepper(Wavefunction& wf,
28
  Preconditioner& prec, TimerMap& tmap) : prec_(prec),
Francois Gygi committed
29
  WavefunctionStepper(wf,tmap), wf_last_(wf), dwf_last_(wf),
30
  extrapolate_(false)
31
{}
Francois Gygi committed
32 33 34

////////////////////////////////////////////////////////////////////////////////
PSDAWavefunctionStepper::~PSDAWavefunctionStepper(void)
35
{}
36 37 38 39

////////////////////////////////////////////////////////////////////////////////
void PSDAWavefunctionStepper::update(Wavefunction& dwf)
{
Francois Gygi committed
40
  tmap_["psda_residual"].start();
41 42 43 44
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
    {
Francois Gygi committed
45 46 47
      // compute A = V^T H V  and descent direction HV - VA

      if ( wf_.sd(ispin,ikp)->basis().real() )
48
      {
Francois Gygi committed
49 50 51 52 53 54 55 56 57 58 59 60 61
        // proxy real matrices c, cp
        DoubleMatrix c_proxy(wf_.sd(ispin,ikp)->c());
        DoubleMatrix cp_proxy(dwf.sd(ispin,ikp)->c());
        DoubleMatrix a(c_proxy.context(),c_proxy.n(),c_proxy.n(),
          c_proxy.nb(),c_proxy.nb());

        // factor 2.0 in next line: G and -G
        a.gemm('t','n',2.0,c_proxy,cp_proxy,0.0);
        // rank-1 update correction
        a.ger(-1.0,c_proxy,0,cp_proxy,0);

        // cp = cp - c * a
        cp_proxy.gemm('n','n',-1.0,c_proxy,a,1.0);
Francois Gygi committed
62 63 64
      }
      else
      {
65 66 67 68 69 70
        ComplexMatrix& c = wf_.sd(ispin,ikp)->c();
        ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
        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);
Francois Gygi committed
71 72 73 74
      }
    }
  }
  tmap_["psda_residual"].stop();
Francois Gygi committed
75

Francois Gygi committed
76
  // dwf.sd->c() now contains the descent direction (HV-VA) (residual)
77 78
  // update the preconditioner
  prec_.update(wf_);
Francois Gygi committed
79

Francois Gygi committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
    {
      tmap_["psda_prec"].start();
      // Apply preconditioner K and store -K(HV-VA) in dwf
      double* c = (double*) wf_.sd(ispin,ikp)->c().valptr();
      double* c_last = (double*) wf_last_.sd(ispin,ikp)->c().valptr();
      double* dc = (double*) dwf.sd(ispin,ikp)->c().valptr();
      double* dc_last = (double*) dwf_last_.sd(ispin,ikp)->c().valptr();
      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* dcn = &dc[2*mloc*n];
        for ( int i = 0; i < ngwl; i++ )
Francois Gygi committed
99
        {
100
          const double fac = prec_.diag(ispin,ikp,n,i);
Francois Gygi committed
101 102 103 104
          const double f0 = -fac * dcn[2*i];
          const double f1 = -fac * dcn[2*i+1];
          dcn[2*i] = f0;
          dcn[2*i+1] = f1;
Francois Gygi committed
105
        }
Francois Gygi committed
106 107
      }
      tmap_["psda_prec"].stop();
Francois Gygi committed
108

Francois Gygi committed
109 110
      // dwf now contains the preconditioned descent
      // direction -K(HV-VA)
Francois Gygi committed
111

Francois Gygi committed
112 113 114 115 116 117 118
      tmap_["psda_update_wf"].start();
      // Anderson extrapolation
      if ( extrapolate_ )
      {
        double theta = 0.0;
        double a = 0.0, b = 0.0;
        for ( int i = 0; i < 2*mloc*nloc; i++ )
119
        {
Francois Gygi committed
120 121
          const double f = dc[i];
          const double delta_f = f - dc_last[i];
122

Francois Gygi committed
123 124
          // accumulate partial sums of a and b
          // a = delta_F * F
Francois Gygi committed
125

Francois Gygi committed
126 127 128
          a += f * delta_f;
          b += delta_f * delta_f;
        }
129 130

        if ( wf_.sd(ispin,ikp)->basis().real() )
Francois Gygi committed
131
        {
132 133 134 135 136
          // correct for double counting of asum and bsum on first row
          // factor 2.0: G and -G
          a *= 2.0;
          b *= 2.0;
          if ( wf_.sdcontext()->myrow() == 0 )
137
          {
138 139 140 141 142 143 144 145 146 147
            for ( int n = 0; n < nloc; n++ )
            {
              const int i = 2*mloc*n;
              const double f0 = dc[i];
              const double f1 = dc[i+1];
              const double delta_f0 = f0 - dc_last[i];
              const double delta_f1 = f1 - dc_last[i+1];
              a -= f0 * delta_f0 + f1 * delta_f1;
              b -= delta_f0 * delta_f0 + delta_f1 * delta_f1;
            }
Francois Gygi committed
148
          }
Francois Gygi committed
149
        }
150

Francois Gygi committed
151 152 153 154 155
        // a and b contain the partial sums of a and b
        double tmpvec[2] = { a, b };
        wf_.sdcontext()->dsum(2,1,&tmpvec[0],2);
        a = tmpvec[0];
        b = tmpvec[1];
156

Francois Gygi committed
157 158 159
        // compute theta = - a / b
        if ( b != 0.0 )
          theta = - a / b;
160

Francois Gygi committed
161 162 163 164
        if ( theta < -1.0 )
        {
          theta = 0.0;
        }
Francois Gygi committed
165

Francois Gygi committed
166
        theta = min(2.0,theta);
Francois Gygi committed
167

Francois Gygi committed
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
        // extrapolation
        for ( int i = 0; i < 2*mloc*nloc; i++ )
        {
          // x_bar = x_ + theta * ( x_ - xlast_ ) (store in x_)
          const double x = c[i];
          const double xlast = c_last[i];
          const double xbar = x + theta * ( x - xlast );

          // f_bar = f + theta * ( f - flast ) (store in f)
          const double f = dc[i];
          const double flast = dc_last[i];
          const double fbar = f + theta * ( f - flast );

          c[i] = xbar + fbar;
          c_last[i] = x;
          dc_last[i] = f;
Francois Gygi committed
184
        }
Francois Gygi committed
185 186 187 188 189
      }
      else
      {
        // no extrapolation
        for ( int i = 0; i < 2*mloc*nloc; i++ )
Francois Gygi committed
190
        {
Francois Gygi committed
191 192 193
          // x_ = x_ + f_
          const double x = c[i];
          const double f = dc[i];
Francois Gygi committed
194

Francois Gygi committed
195 196 197
          c[i] = x + f;
          c_last[i] = x;
          dc_last[i] = f;
198
        }
Francois Gygi committed
199 200
      }
      tmap_["psda_update_wf"].stop();
Francois Gygi committed
201

Francois Gygi committed
202 203 204 205
      enum ortho_type { GRAM, LOWDIN, ORTHO_ALIGN, RICCATI };
      //const ortho_type ortho = GRAM;
      //const ortho_type ortho = LOWDIN;
      const ortho_type ortho = ORTHO_ALIGN;
Francois Gygi committed
206

Francois Gygi committed
207
      switch ( ortho )
Francois Gygi committed
208
      {
Francois Gygi committed
209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231
        case GRAM:
          tmap_["gram"].start();
          wf_.sd(ispin,ikp)->gram();
          tmap_["gram"].stop();
          break;

        case LOWDIN:
          tmap_["lowdin"].start();
          wf_.sd(ispin,ikp)->lowdin();
          tmap_["lowdin"].stop();
          break;

        case ORTHO_ALIGN:
          tmap_["ortho_align"].start();
          wf_.sd(ispin,ikp)->ortho_align(*wf_last_.sd(ispin,ikp));
          tmap_["ortho_align"].stop();
          break;

        case RICCATI:
          tmap_["riccati"].start();
          wf_.sd(ispin,ikp)->riccati(*wf_last_.sd(ispin,ikp));
          tmap_["riccati"].stop();
          break;
232
      }
Francois Gygi committed
233 234
    } // ikp
  } // ispin
235
  extrapolate_ = true;
236
}