← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/detcas/diis.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/detcas/diis.cc
 
1
/*! \file
 
2
    \ingroup DETCAS
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
/*
 
6
** DIIS.C
 
7
** 
 
8
** Routines for Direct Inversion of the Iterative Subspace Interpolation
 
9
** of orbital rotation angles, P. Pulay, Chem. Phys. Lett. 73, 393 (1980).
 
10
**
 
11
** C. David Sherrill
 
12
** University of California, Berkeley
 
13
** May 1998
 
14
*/
 
15
 
 
16
#include <cstdlib>
 
17
#include <cstdio>
 
18
#include <cmath>
 
19
#include <libciomr/libciomr.h>
 
20
#include "globaldefs.h"
 
21
#include "globals.h"
 
22
 
 
23
namespace psi { namespace detcas {
 
24
 
 
25
#define DIIS_MIN_DET 1.0E-16
 
26
 
 
27
/*
 
28
** diis()
 
29
**
 
30
** This top-level routine manages all the DIIS stuff.  
 
31
**
 
32
** Parameters:
 
33
**   veclen = length of the vectors to extrapolate
 
34
**   vec    = new vector to add
 
35
**   errvec = new error vector to add
 
36
**
 
37
** Returns:
 
38
**   1 if DIIS step taken, otherwise 0
 
39
*/
 
40
int diis(int veclen, double *vec, double *errvec)
 
41
{
 
42
  int i, j, k;
 
43
  int num_vecs, new_num_vecs, offset, diis_iter, do_diis;
 
44
  double **vecs, **errvecs, **bmat, *bvec, tval, det, scale_factor;
 
45
  FILE *fp;
 
46
 
 
47
  
 
48
  /* add the vector and error vector to subspace */
 
49
  ffileb_noexit(&fp,"diis.dat",2);
 
50
  if (fp != NULL) {
 
51
 
 
52
    if (fread(&num_vecs, sizeof(int), 1, fp) != 1) {
 
53
      fprintf(outfile, "(diis): Error reading number of diis vectors.\n");
 
54
      return(0);
 
55
    }
 
56
 
 
57
    if (fread(&diis_iter, sizeof(int), 1, fp) != 1) {
 
58
      fprintf(outfile, "(diis): Error reading diis iteration number.\n");
 
59
      return(0);
 
60
    }
 
61
 
 
62
    vecs = block_matrix(num_vecs+1, veclen);
 
63
    errvecs = block_matrix(num_vecs+1, veclen);
 
64
 
 
65
    for (i=0; i<num_vecs; i++) {
 
66
 
 
67
      if (fread(vecs[i], sizeof(double), veclen, fp) != veclen) { 
 
68
        fprintf(outfile, "(diis): Error reading diis vector %d\n", i);
 
69
        return(0);
 
70
      }
 
71
 
 
72
      if (fread(errvecs[i], sizeof(double), veclen, fp) != veclen) { 
 
73
        fprintf(outfile, "(diis): Error reading diis error vector %d\n", i);
 
74
        return(0);
 
75
      }
 
76
 
 
77
    }
 
78
 
 
79
    fclose(fp);
 
80
  }
 
81
 
 
82
  else {
 
83
    num_vecs = 0;
 
84
    diis_iter = 0;
 
85
    vecs = block_matrix(num_vecs+1, veclen);
 
86
    errvecs = block_matrix(num_vecs+1, veclen);
 
87
  } 
 
88
 
 
89
  /* will we do a diis this time? */
 
90
  diis_iter++;
 
91
  num_vecs++;
 
92
  if ((diis_iter % Params.diis_freq) || (num_vecs < Params.diis_min_vecs)) 
 
93
    do_diis = 0; 
 
94
  else
 
95
    do_diis = 1;
 
96
 
 
97
  for (i=0; i<veclen; i++) vecs[num_vecs-1][i] = vec[i];
 
98
  for (i=0; i<veclen; i++) errvecs[num_vecs-1][i] = errvec[i];
 
99
 
 
100
  offset = 0;
 
101
  if (num_vecs > Params.diis_max_vecs) 
 
102
    offset = num_vecs - Params.diis_max_vecs;
 
103
 
 
104
  if (Params.print_lvl > 2) 
 
105
    fprintf(outfile, "Diis: iter %2d, vecs %d, do_diis %d, offset %d\n", 
 
106
            diis_iter, num_vecs, do_diis, offset);
 
107
 
 
108
  new_num_vecs = num_vecs - offset;
 
109
 
 
110
  /* write out the diis info */
 
111
  ffileb_noexit(&fp,"diis.dat",0);
 
112
  if (fp == NULL) {
 
113
    fprintf(outfile, "(diis): Error opening diis.dat\n");
 
114
    return(0);
 
115
  } 
 
116
 
 
117
  if (fwrite(&new_num_vecs, sizeof(int), 1, fp) != 1) {
 
118
    fprintf(outfile, "(diis): Error writing number of diis vectors.\n");
 
119
    return(0);
 
120
  }
 
121
 
 
122
  if (fwrite(&diis_iter, sizeof(int), 1, fp) != 1) {
 
123
    fprintf(outfile, "(diis): Error writing diis iteration number.\n");
 
124
    return(0);
 
125
  }
 
126
 
 
127
  for (i=offset; i<num_vecs; i++) {
 
128
    if (fwrite(vecs[i], sizeof(double), veclen, fp) != veclen) {
 
129
      fprintf(outfile, "(diis): Error writing vector %d.\n", i);
 
130
    }
 
131
    if (fwrite(errvecs[i], sizeof(double), veclen, fp) != veclen) {
 
132
      fprintf(outfile, "(diis): Error writing error vector %d.\n", i);
 
133
    }
 
134
  }
 
135
 
 
136
  fclose(fp);
 
137
 
 
138
 
 
139
  /* don't take a diis step if it's not time */
 
140
  if (!do_diis) {
 
141
    free_block(vecs);
 
142
    free_block(errvecs);
 
143
    return(0); 
 
144
  }
 
145
 
 
146
  /* form diis matrix, solve equations */
 
147
  if (Params.print_lvl) 
 
148
    fprintf(outfile, "Attempting a DIIS step\n");
 
149
 
 
150
  bmat = block_matrix(new_num_vecs+1, new_num_vecs+1);
 
151
  bvec = init_array(new_num_vecs+1);
 
152
 
 
153
  for (i=0; i<new_num_vecs; i++) {
 
154
    for (j=0; j<=i; j++) {
 
155
      tval = 0.0;
 
156
      for (k=0; k<veclen; k++) {
 
157
        tval += errvecs[i+offset][k] * errvecs[j+offset][k];
 
158
      }
 
159
      bmat[i][j] = tval;
 
160
      bmat[j][i] = tval;
 
161
    }
 
162
  }
 
163
 
 
164
  for (i=0; i<new_num_vecs; i++) {
 
165
    bmat[i][new_num_vecs] = -1.0;
 
166
    bmat[new_num_vecs][i] = -1.0;
 
167
  }
 
168
 
 
169
  bmat[new_num_vecs][new_num_vecs] = 0.0;
 
170
  bvec[new_num_vecs] = -1.0;
 
171
 
 
172
  if (Params.print_lvl > 2) {
 
173
    fprintf(outfile, "DIIS B Matrix:\n");
 
174
    print_mat(bmat, new_num_vecs+1, new_num_vecs+1, outfile);
 
175
  }
 
176
 
 
177
  /* scale the B matrix */
 
178
  scale_factor = 1.0 / bmat[0][0];
 
179
  for (i=0; i<new_num_vecs; i++) {
 
180
    for (j=0; j<new_num_vecs; j++) {
 
181
      bmat[i][j] = bmat[i][j] * scale_factor;
 
182
    }
 
183
  }
 
184
 
 
185
  if (Params.print_lvl > 2) {
 
186
    fprintf(outfile, "DIIS B Matrix:\n");
 
187
    print_mat(bmat, new_num_vecs+1, new_num_vecs+1, outfile);
 
188
  }
 
189
 
 
190
  /* now solve the linear equations */ 
 
191
  flin(bmat, bvec, new_num_vecs+1, 1, &det); 
 
192
 
 
193
  if (fabs(det) < DIIS_MIN_DET) {
 
194
    fprintf(outfile, "Warning: diis matrix near-singular\n");
 
195
    fprintf(outfile, "Determinant is %6.3E\n", det);
 
196
  }
 
197
 
 
198
  if (Params.print_lvl > 3) {
 
199
    fprintf(outfile, "\nCoefficients of DIIS extrapolant:\n");
 
200
    for (i=0; i<new_num_vecs; i++) {
 
201
      fprintf(outfile, "%12.6lf\n", bvec[i]);
 
202
    }
 
203
    fprintf(outfile, "Lambda:\n");
 
204
    fprintf(outfile, "%12.6lf\n", bvec[i]);
 
205
  }
 
206
 
 
207
  /* get extrapolated vector */
 
208
  zero_arr(CalcInfo.theta_cur, veclen);
 
209
  for (i=0; i<new_num_vecs; i++) {
 
210
    for (j=0; j<veclen; j++) {
 
211
      CalcInfo.theta_cur[j] += bvec[i] * vecs[i+offset][j];
 
212
    }
 
213
  }
 
214
 
 
215
  free_block(vecs);
 
216
  free_block(errvecs);
 
217
  free_block(bmat); 
 
218
  free(bvec);
 
219
 
 
220
  return(1);
 
221
}
 
222
 
 
223
}} // end namespace psi::detcas
 
224