← Back to branch summary

~ubuntu-branches/ubuntu/trusty/psicode/trusty

« back to all changes in this revision

Viewing changes to src/lib/libqt/david.c

  • 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/lib/libqt/david.c
1
 
/*! 
2
 
    \file david.c
3
 
    \ingroup (QT)
4
 
*/
5
 
 
6
 
#include <stdio.h>
7
 
#include <stdlib.h>
8
 
#include <string.h>
9
 
#include <math.h>
10
 
#include <libciomr/libciomr.h>
11
 
#include <libqt/qt.h>
12
 
 
13
 
#define BIGNUM 1E100
14
 
#define MAXIT 1000
15
 
 
16
 
/*!
17
 
** david(): Computes the lowest few eigenvalues and eigenvectors of a
18
 
** symmetric matrix, A, using the Davidson-Liu algorithm.  The matrix
19
 
** must be small enough to fit entirely in core.  This algorithm is
20
 
** useful if one is interested in only a few roots of the matrix
21
 
** rather than the whole spectrum.
22
 
**
23
 
** NB: This implementation will keep up to eight guess vectors for each
24
 
** root desired before collapsing to one vector per root.  In
25
 
** addition, if smart_guess=1 (the default), guess vectors are
26
 
** constructed by diagonalization of a sub-matrix of A; otherwise,
27
 
** unit vectors are used.
28
 
**
29
 
** TDC, July-August 2002
30
 
**
31
 
**   \param A      = matrix to diagonalize
32
 
**   \param N      = dimension of A
33
 
**   \param M      = number of roots desired
34
 
**   \param eps    = eigenvalues
35
 
**   \param v      = eigenvectors
36
 
**   \param cutoff = tolerance for convergence of eigenvalues
37
 
**   \param print  = Boolean for printing additional information
38
 
**
39
 
** Returns: number of converged roots
40
 
** \ingroup (QT)
41
 
*/
42
 
 
43
 
int david(double **A, int N, int M, double *eps, double **v, 
44
 
          double cutoff, int print)
45
 
{
46
 
  int i, j, k, L, I;
47
 
  double minimum;
48
 
  int min_pos, numf, iter, *conv, converged, maxdim, skip_check;
49
 
  int *small2big, init_dim;
50
 
  int smart_guess = 1;
51
 
  double *Adiag, **b, **bnew, **sigma, **G;
52
 
  double *lambda, **alpha, **f, *lambda_old;
53
 
  double norm, denom, diff;
54
 
 
55
 
  maxdim = 8 * M;
56
 
 
57
 
  b = block_matrix(maxdim, N);  /* current set of guess vectors,
58
 
                                   stored by row */
59
 
  bnew = block_matrix(M, N); /* guess vectors formed from old vectors,
60
 
                                stored by row*/
61
 
  sigma = block_matrix(N, maxdim); /* sigma vectors, stored by column */
62
 
  G = block_matrix(maxdim, maxdim); /* Davidson mini-Hamitonian */
63
 
  f = block_matrix(maxdim, N); /* residual eigenvectors, stored by row */
64
 
  alpha = block_matrix(maxdim, maxdim); /* eigenvectors of G */
65
 
  lambda = init_array(maxdim); /* eigenvalues of G */
66
 
  lambda_old = init_array(maxdim); /* approximate roots from previous
67
 
                                      iteration */
68
 
 
69
 
  if(smart_guess) { /* Use eigenvectors of a sub-matrix as initial guesses */
70
 
 
71
 
    if(N > 7*M) init_dim = 7*M;
72
 
    else init_dim = M;
73
 
    Adiag = init_array(N);
74
 
    small2big = init_int_array(7*M);
75
 
    for(i=0; i < N; i++) { Adiag[i] = A[i][i]; }
76
 
    for(i=0; i < init_dim; i++) {
77
 
      minimum = Adiag[0];
78
 
      min_pos = 0;
79
 
      for(j=1; j < N; j++)
80
 
        if(Adiag[j] < minimum) {        
81
 
          minimum = Adiag[j]; 
82
 
          min_pos = j; 
83
 
          small2big[i] = j; 
84
 
        }
85
 
 
86
 
      Adiag[min_pos] = BIGNUM;
87
 
      lambda_old[i] = minimum;
88
 
    }
89
 
    for(i=0; i < init_dim; i++) {
90
 
      for(j=0; j < init_dim; j++)
91
 
        G[i][j] = A[small2big[i]][small2big[j]];
92
 
    }
93
 
 
94
 
    sq_rsp(init_dim, init_dim, G, lambda, 1, alpha, 1e-12);
95
 
 
96
 
    for(i=0; i < init_dim; i++) {
97
 
      for(j=0; j < init_dim; j++)
98
 
        b[i][small2big[j]] = alpha[j][i];
99
 
    }
100
 
 
101
 
    free(Adiag);
102
 
    free(small2big);
103
 
  }
104
 
  else { /* Use unit vectors as initial guesses */
105
 
    Adiag = init_array(N);
106
 
    for(i=0; i < N; i++) { Adiag[i] = A[i][i]; }
107
 
    for(i=0; i < M; i++) {
108
 
      minimum = Adiag[0];
109
 
      min_pos = 0;
110
 
      for(j=1; j < N; j++)
111
 
        if(Adiag[j] < minimum) { minimum = Adiag[j]; min_pos = j; }
112
 
 
113
 
      b[i][min_pos] = 1.0; 
114
 
      Adiag[min_pos] = BIGNUM; 
115
 
      lambda_old[i] = minimum;
116
 
    }
117
 
    free(Adiag);
118
 
  }
119
 
 
120
 
  L = init_dim;
121
 
  iter =0;
122
 
  converged = 0;
123
 
  conv = init_int_array(M); /* boolean array for convergence of each
124
 
                               root */
125
 
  while(converged < M && iter < MAXIT) {
126
 
 
127
 
    skip_check = 0;
128
 
    if(print) printf("\niter = %d\n", iter); 
129
 
 
130
 
    /* form mini-matrix */
131
 
    C_DGEMM('n','t', N, L, N, 1.0, &(A[0][0]), N, &(b[0][0]), N,
132
 
            0.0, &(sigma[0][0]), maxdim);
133
 
    C_DGEMM('n','n', L, L, N, 1.0, &(b[0][0]), N,
134
 
            &(sigma[0][0]), maxdim, 0.0, &(G[0][0]), maxdim);
135
 
 
136
 
    /* diagonalize mini-matrix */
137
 
    sq_rsp(L, L, G, lambda, 1, alpha, 1e-12);
138
 
 
139
 
    /* form preconditioned residue vectors */
140
 
    for(k=0; k < M; k++)
141
 
      for(I=0; I < N; I++) {
142
 
        f[k][I] = 0.0;
143
 
        for(i=0; i < L; i++) {
144
 
          f[k][I] += alpha[i][k] * (sigma[I][i] - lambda[k] * b[i][I]);
145
 
        }
146
 
        denom = lambda[k] - A[I][I];
147
 
        if(fabs(denom) > 1e-6) f[k][I] /= denom;
148
 
        else f[k][I] = 0.0;
149
 
      }
150
 
 
151
 
    /* normalize each residual */
152
 
    for(k=0; k < M; k++) {
153
 
      norm = 0.0;
154
 
      for(I=0; I < N; I++) {
155
 
        norm += f[k][I] * f[k][I];
156
 
      }
157
 
      norm = sqrt(norm);
158
 
      for(I=0; I < N; I++) {
159
 
        if(norm > 1e-6) f[k][I] /= norm;
160
 
        else f[k][I] = 0.0;
161
 
      }
162
 
    }
163
 
 
164
 
    /* schmidt orthogonalize the f[k] against the set of b[i] and add
165
 
       new vectors */
166
 
    for(k=0,numf=0; k < M; k++)
167
 
      if(schmidt_add(b, L, N, f[k])) { L++; numf++; }
168
 
 
169
 
    /* If L is close to maxdim, collapse to one guess per root */
170
 
    if(maxdim - L < M) {
171
 
      if(print) {
172
 
        printf("Subspace too large: maxdim = %d, L = %d\n", maxdim, L);
173
 
        printf("Collapsing eigenvectors.\n");
174
 
      }
175
 
      for(i=0; i < M; i++) {
176
 
        memset((void *) bnew[i], 0, N*sizeof(double));
177
 
        for(j=0; j < L; j++) {
178
 
          for(k=0; k < N; k++) {
179
 
            bnew[i][k] += alpha[j][i] * b[j][k];
180
 
          }
181
 
        }
182
 
      }
183
 
 
184
 
      /* copy new vectors into place */
185
 
      for(i=0; i < M; i++) 
186
 
        for(k=0; k < N; k++)
187
 
          b[i][k] = bnew[i][k];
188
 
 
189
 
      skip_check = 1;
190
 
 
191
 
      L = M;
192
 
    }
193
 
 
194
 
    /* check convergence on all roots */
195
 
    if(!skip_check) {
196
 
      converged = 0;
197
 
      zero_int_array(conv, M);
198
 
      if(print) {
199
 
        printf("Root      Eigenvalue       Delta  Converged?\n");
200
 
        printf("---- -------------------- ------- ----------\n");
201
 
      }
202
 
      for(k=0; k < M; k++) {
203
 
        diff = fabs(lambda[k] - lambda_old[k]);
204
 
        if(diff < cutoff) {
205
 
          conv[k] = 1;
206
 
          converged++;
207
 
        }
208
 
        lambda_old[k] = lambda[k];
209
 
        if(print) {
210
 
          printf("%3d  %20.14f %4.3e    %1s\n", k, lambda[k], diff,
211
 
                 conv[k] == 1 ? "Y" : "N");
212
 
        }
213
 
      }
214
 
    }
215
 
 
216
 
    iter++;
217
 
  }
218
 
 
219
 
  /* generate final eigenvalues and eigenvectors */
220
 
  if(converged == M) {
221
 
    for(i=0; i < M; i++) {
222
 
      eps[i] = lambda[i];
223
 
      for(j=0; j < L; j++) {
224
 
        for(I=0; I < N; I++) {
225
 
          v[I][i] += alpha[j][i] * b[j][I];
226
 
        }
227
 
      }
228
 
    }
229
 
    if(print) printf("Davidson algorithm converged in %d iterations.\n", iter);
230
 
  }
231
 
 
232
 
  free(conv);
233
 
  free_block(b);
234
 
  free_block(bnew);
235
 
  free_block(sigma);
236
 
  free_block(G);
237
 
  free_block(f);
238
 
  free_block(alpha);
239
 
  free(lambda);
240
 
  free(lambda_old);
241
 
 
242
 
  return converged;
243
 
}