PSDAWavefunctionStepper.C 7 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 65 66 67 68 69 70 71
      }
      else
      {
        // not implemented in the complex case
        cout << "PSDA is not implemented for complex wave functions" << endl;
        assert(false);
      }
    }
  }
  tmap_["psda_residual"].stop();
Francois Gygi committed
72

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

Francois Gygi committed
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
  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
96
        {
97
          const double fac = prec_.diag(ispin,ikp,n,i);
Francois Gygi committed
98 99 100 101
          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
102
        }
Francois Gygi committed
103 104
      }
      tmap_["psda_prec"].stop();
Francois Gygi committed
105

Francois Gygi committed
106 107
      // dwf now contains the preconditioned descent
      // direction -K(HV-VA)
Francois Gygi committed
108

Francois Gygi committed
109 110 111 112 113 114 115
      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++ )
116
        {
Francois Gygi committed
117 118
          const double f = dc[i];
          const double delta_f = f - dc_last[i];
119

Francois Gygi committed
120 121
          // accumulate partial sums of a and b
          // a = delta_F * F
Francois Gygi committed
122

Francois Gygi committed
123 124 125 126 127 128 129 130 131 132
          a += f * delta_f;
          b += delta_f * delta_f;
        }
        // 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 )
        {
          for ( int n = 0; n < nloc; n++ )
133
          {
Francois Gygi committed
134 135 136 137 138 139 140
            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
141
          }
Francois Gygi committed
142
        }
143

Francois Gygi committed
144 145 146 147 148
        // 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];
149

Francois Gygi committed
150 151 152
        // compute theta = - a / b
        if ( b != 0.0 )
          theta = - a / b;
153

Francois Gygi committed
154 155
        if ( wf_.sdcontext()->onpe0() )
          cout << "  Anderson extrapolation: theta=" << theta;
156

Francois Gygi committed
157 158 159 160
        if ( theta < -1.0 )
        {
          theta = 0.0;
        }
Francois Gygi committed
161

Francois Gygi committed
162
        theta = min(2.0,theta);
Francois Gygi committed
163

Francois Gygi committed
164 165
        if ( wf_.sdcontext()->onpe0() )
          cout <<" (" << theta << ")" << endl;
Francois Gygi committed
166

Francois Gygi committed
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
        // 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
183
        }
Francois Gygi committed
184 185 186 187 188
      }
      else
      {
        // no extrapolation
        for ( int i = 0; i < 2*mloc*nloc; i++ )
Francois Gygi committed
189
        {
Francois Gygi committed
190 191 192
          // x_ = x_ + f_
          const double x = c[i];
          const double f = dc[i];
Francois Gygi committed
193

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

Francois Gygi committed
201 202 203 204
      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
205

Francois Gygi committed
206
      switch ( ortho )
Francois Gygi committed
207
      {
Francois Gygi committed
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
        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;
231
      }
Francois Gygi committed
232 233
    } // ikp
  } // ispin
234
  extrapolate_ = true;
235
}