3
\brief Enter brief description of file here
7
** Simultaneous Expansion Method for the Iterative Solution of
8
** Several of the Lowest Eigenvalues and Corresponding Eivenvectors of
9
** Large Real-Symmetric Matrices
11
** Algorithm due to Bowen Liu
12
** IBM Research Laboratory
14
** Implemented for Schaefer Group by David Sherrill
15
** Center for Computational Quantum Chemistry, UGA
17
** In-core version for now!
20
** Updated 12/7/95 for testing within rasci code
21
** Updated 11/21/97 for least squares extrapolation and debugging of sem
28
#include <libciomr/libciomr.h>
34
namespace psi { namespace detci {
36
#define MAX_B_ROWS 200
37
#define MIN_F_DENOM 1.0E-3
40
/*** do a little test routine
44
double *evals, **evecs ;
49
ffile(&outfile, "output.dat", 0) ;
52
A = init_matrix(50,50) ;
53
evals = init_array(4) ;
54
evecs = init_matrix(4, 50) ;
55
for (i=0; i<50; i++) {
56
for (j=0; j<=i; j++) {
57
if (i!=j) {A[i][j] = 1.0; A[j][i] = 1.0; }
58
else if (i<5) A[i][j] = 1.0 + 0.1 * (double) i ;
59
else A[i][j] = 2.0 * (double) i + 1.0;
62
sem(A, 50, 4, evecs, evals, 1.0E-10, 6, &used);
65
fprintf(outfile, "Ok, the eigenvectors are sideways!\n");
66
eivout(evecs, evals, 4, 50, outfile);
67
fprintf(outfile, "\nused %d expansion vectors\n", used);
76
** sem(): Use Liu's Simultaneous Expansion Method to find the lowest
77
** few eigenvalues of a real matrix.
80
** A = matrix to find eigenvalues of
82
** M = number of eigenvalues to solve for
83
** L = number of initial vectors in subspace
84
** evecs = matrix for eigenvectors
85
** evals = array for eigenvalues
86
** b = set of subspace vectors (dimensions maxiter x N)
87
** space for rows i < L should not yet be allocated!!
88
** conv_e = convergence tolerance. The lowest energy eigenvalue must
89
** be converged to within this range. It is interesting
90
** that higher roots _may_ converge faster.
91
** conv_rms = the required tolerance for convergence of the CI correction
93
** maxiter = max number of iterations allowed
94
** offst = offset to add to eigenvalues in printing (e.g. enuc)
95
** vu = pointer to int to hold how many expansion vectors used
96
** outfile = output file
100
void sem_test(double **A, int N, int M, int L, double **evecs, double *evals,
101
double **b, double conv_e, double conv_rms, int maxiter, double offst,
102
int *vu, int maxnvect, FILE *outfile)
104
double *tmp_vec, **tmp_mat ;
112
double *lambda, **alpha, **f ;
113
double *converged_root;
114
double **m_lambda, ***m_alpha;
115
double tval, *dvecnorm;
116
int converged=0, iter=1;
117
int iter2=0; /* iterations since last collapse */
119
int lse_do=0, last_lse_collapse_num=-Parameters.lse_collapse, collapse_num=0;
120
double lse_tolerance=5e-3;
121
double **sigma_overlap, ***Mmatrix;
124
/* check parameters */
125
if (evecs == NULL || evals == NULL) {
126
printf("(sem): passed uncallocated pointers for evecs or evals\n") ;
132
A[I][I] -= CalcInfo.efzc;
135
/* make space for temp vector */
136
tmp_vec = init_array(N);
137
lastroot = init_array(N);
138
converged_root = init_array(M);
139
dvecnorm = init_array(M);
141
/* allocate other arrays with ~fixed dimensions during iteration */
142
d = init_matrix(M, N); /* M and N are both fixed */
143
f = init_matrix(M, N);
144
G = init_matrix(maxnvect, maxnvect);
145
tmp_mat = init_matrix(maxnvect, N);
146
jnk = init_matrix(maxnvect, N);
147
lambda = init_array(maxnvect);
148
alpha = init_matrix(maxnvect, maxnvect);
149
sigma_overlap = init_matrix(maxnvect, maxnvect);
150
Mmatrix = (double ***) malloc (sizeof(double **) * M);
152
Mmatrix[i] = init_matrix(maxnvect, maxnvect);
153
m_lambda = init_matrix(M, maxnvect);
154
m_alpha = (double ***) malloc (sizeof(double **) * maxnvect);
155
for (i=0; i<maxnvect; i++) {
156
m_alpha[i] = init_matrix(maxnvect, maxnvect);
158
Lvec = init_int_array(maxnvect);
162
while (!converged && iter <= maxiter) {
166
mmult(b, 0, A, 0, tmp_mat, 0, L, N, N, 0); /* tmp = B * A */
167
mmult(tmp_mat, 0, b, 1, G, 0, L, N, L, 0); /* G = tmp * B(T) */
169
/* solve the L x L eigenvalue problem G a = lambda a for M roots */
170
sq_rsp(L, L, G, lambda, 1, alpha, 1E-14);
172
if (N<100 && Parameters.print_lvl >=3) {
173
fprintf(outfile,"\n b matrix\n");
174
print_mat(b,L,N,outfile);
175
fprintf(outfile,"\n sigma matrix\n");
176
print_mat(tmp_mat,L,N,outfile);
177
fprintf(outfile,"\n G matrix (%d)\n", iter-1);
178
print_mat(G,L,L,outfile);
179
fprintf(outfile,"\n Eigenvectors and eigenvalues of G matrix (%d)\n", iter-1);
180
eivout(alpha, lambda, L, L, outfile);
184
if (Parameters.lse && (maxnvect-L <= M*Parameters.collapse_size) && L>2 &&
185
(lse_tolerance > fabs(lambda[0]-lastroot[0])) && iter>=3 &&
186
((collapse_num-last_lse_collapse_num)>= Parameters.lse_collapse))
189
/* Form sigma_overlap matrix */
190
zero_mat(sigma_overlap,maxnvect,maxnvect);
191
mmult(b, 0, A, 0, tmp_mat, 0, L, N, N, 0);
192
mmult(tmp_mat, 0, tmp_mat, 1, sigma_overlap, 0, L, N, L, 0);
194
/* Form Mij matrix */
195
for (k=0; k<M; k++) {
196
zero_mat(Mmatrix[k],maxnvect,maxnvect);
197
for (i=0; i<L; i++) {
198
for (j=i; j<L; j++) {
199
Mmatrix[k][i][j] = Mmatrix[k][j][i] = sigma_overlap[i][j]
200
-2.0 * lambda[k] * G[i][j];
201
if (i==j) Mmatrix[k][i][i] += pow(lambda[k],2.0);
204
} /* end loop over k (nroots) */
206
if (Parameters.print_lvl > 2) {
207
fprintf(outfile, "\nsigma_overlap matrix (%2d) = \n", iter-1);
208
print_mat(sigma_overlap, L, L, outfile);
210
for (k=0; k<M; k++) {
211
fprintf(outfile, "\nM matrix (%2d) for root %d = \n", iter, k);
212
print_mat(Mmatrix[k], L, L, outfile);
213
fprintf(outfile, "\n");
217
/* solve the L x L eigenvalue problem M a = lambda a for M roots */
218
for (k=0; k<M; k++) {
219
sq_rsp(L, L, Mmatrix[k], m_lambda[k], 1, m_alpha[k], 1.0E-14);
220
if (Parameters.print_lvl > 2) {
221
fprintf(outfile, "\n M eigenvectors and eigenvalues root %d:\n",k);
222
eivout(m_alpha[k], m_lambda[k], L, L, outfile);
226
} /* end if lse_do */
228
if ((Parameters.collapse_size>0) && (iter2-Parameters.collapse_size+1 > 0)
229
&& (Lvec[iter2-Parameters.collapse_size+1]+M*Parameters.collapse_size
230
> maxnvect) && iter!=maxiter) {
233
if (lse_do) last_lse_collapse_num = collapse_num;
234
/* copy ci vector into d matrix */
240
d[k][I] += m_alpha[k][i][0] * b[i][I];
246
d[k][I] += alpha[i][k] * b[i][I];
248
/* copy d matrix to end of b matrix */
251
b[maxnvect-1-i][I] = d[i][I];
253
/* reorder b matrix pointers */
254
for (i=0; i<L; i++) jnk[i] = b[i];
255
for (i=0; i<L; i++) b[i] = jnk[maxnvect-1-i];
260
/* zero out old parts of b matrix */
261
for (i=L; i<maxnvect; i++) zero_arr(b[i], N);
263
/* reform G matrix */
264
mmult(b, 0, A, 0, tmp_mat, 0, L, N, N, 0); /* tmp = B * A */
265
mmult(tmp_mat, 0, b, 1, G, 0, L, N, L, 0); /* G = tmp * B(T) */
267
/* solve the L x L eigenvalue problem G a = lambda a for M roots */
268
sq_rsp(L, L, G, lambda, 1, alpha, 1E-14);
270
if (N<100 && Parameters.print_lvl >= 3) {
271
fprintf(outfile," Reformed G matrix (%d)\n",iter-1);
272
print_mat(G,L,L,outfile);
273
fprintf(outfile,"\n");
276
if (lse_do) fprintf(outfile," Least Squares Extrapolation\n");
277
fprintf(outfile," Collapse Davidson subspace to %d vectors\n", L);
280
/* form the d part of the correction vector */
282
for (k=0; k<M; k++) {
283
for (i=0; i<L; i++) {
284
mmult(A,0,&(b[i]),1,&(tmp_vec),1,N,N,1,0); /* tmp=A*b[i] */
285
for (I=0; I<N; I++) {
286
d[k][I] += alpha[i][k] * (tmp_vec[I] - lambda[k] * b[i][I]);
291
if (N<100 && Parameters.print_lvl >= 3) {
292
fprintf(outfile," D vectors for iter (%d)\n",iter-1);
293
print_mat(d,M,N,outfile);
296
/* check for convergence */
298
for (i=0; i<M; i++) {
299
dot_arr(d[i], d[i], N, &tval);
302
if (dvecnorm[i] <= conv_rms && fabs(lambda[i] - lastroot[i]) <= conv_e) converged_root[i] = 1;
304
converged_root[i] = 0;
307
fprintf(outfile, "Iter %3d Root %d = %13.9lf",
308
iter-1, i+1, (lambda[i] + offst));
309
fprintf(outfile, " Delta_E %10.3E Delta_C %10.3E %c\n",
310
lambda[i] - lastroot[i], tval, converged_root[i] ? 'c' : ' ');
315
fprintf(outfile, "\n");
319
if (converged || iter == maxiter) {
320
for (i=0; i<M; i++) {
321
evals[i] = lambda[i];
322
for (j=0; j<L; j++) {
325
evecs[i][I] += tval * b[j][I];
331
for (i=0; i<M; i++) lastroot[i] = lambda[i];
334
/* form the correction vector and normalize */
335
for (k=0; k<M; k++) {
336
for (I=0; I<N; I++) {
337
tval = lambda[k] - A[I][I];
338
/* make sure denom isn't 0. If so, make some arbitrary val */
339
/* It might be interesting to figure the best way to do this */
341
/* previous way to do it
342
if (fabs(tval) < MIN_F_DENOM) tval = 0.1;
343
f[k][I] = d[k][I] / tval;
346
/* the way GUGA does it */
347
if (fabs(tval) < 1.0E-8) f[k][I] = 0.0;
348
else f[k][I] = d[k][I] / tval;
353
/* Schmidt orthog and append f's to b */
355
if (converged_root[i] == 0)
356
if (schmidt_add(b, L, N, f[i])) L++;
357
fprintf(outfile," Number of b vectors = %d\n", L);
360
fprintf(outfile, "(test_sem): L(%2d) > maxnvect(%2d)!",L,maxnvect);
361
fprintf(outfile, " Aborting!\n");
365
/* Again Schmidt orthog b's (minimize numerical error) */
366
/* Doesn't this mess up the sigma vectors slightly */
367
schmidt(b, L, N, outfile);
378
free_matrix(b, maxnvect);
379
free_matrix(G, maxnvect);
380
free_matrix(tmp_mat, maxnvect);
381
free_matrix(alpha, maxnvect);
382
free_matrix(sigma_overlap, maxnvect);
383
for (i=0; i<M; i++) free_matrix(Mmatrix[i], maxnvect);
386
}} // namespace psi::detci
added
removed
src/bin/detci/sem_test.cc
