5
#include <libciomr/libciomr.h>
6
#include <libpsio/psio.h>
7
#include <libdpd/dpd.h>
14
** DIIS: Direct inversion in the iterative subspace routine to
15
** accelerate convergence of the CCSD amplitude equations.
17
** Substantially improved efficiency of this routine:
18
** (1) Keeping at most two error vectors in core at once.
19
** (2) Limiting direct product (overlap) calculation to unique pairs.
20
** (3) Using LAPACK's linear equation solver DGESV instead of flin.
22
** These improvements have been applied only to RHF cases so far.
25
** Modified for CCRESPONSE, TDC, 5/03
28
void diis(int iter, char *pert, char *cart, int irrep, double omega)
30
int nvector=8; /* Number of error vectors to keep */
32
int row, col, word, p, q;
36
dpdfile2 T1, T1a, T1b;
37
dpdbuf4 T2, T2a, T2b, T2c;
38
psio_address start, end, next;
40
double **B, *C, **vector;
41
double product, determinant, maximum;
44
nirreps = moinfo.nirreps;
46
if(params.ref == 0) { /** RHF **/
47
/* Compute the length of a single error vector */
48
dpd_file2_init(&T1, CC_MISC, irrep, 0, 1, "XXX");
49
dpd_buf4_init(&T2, CC_MISC, irrep, 0, 5, 0, 5, 0, "XXX");
50
for(h=0; h < nirreps; h++) {
51
vector_length += T1.params->rowtot[h] * T1.params->coltot[h^irrep];
52
vector_length += T2.params->rowtot[h] * T2.params->coltot[h^irrep];
57
/* Set the diis cycle value */
58
diis_cycle = (iter-1) % nvector;
60
/* Build the current error vector and dump it to disk */
61
error = dpd_block_matrix(1,vector_length);
64
sprintf(lbl, "New X_%s_%1s_IA (%5.3f)", pert, cart, omega);
65
dpd_file2_init(&T1a, CC_OEI, irrep, 0, 1, lbl);
66
dpd_file2_mat_init(&T1a);
67
dpd_file2_mat_rd(&T1a);
68
sprintf(lbl, "X_%s_%1s_IA (%5.3f)", pert, cart, omega);
69
dpd_file2_init(&T1b, CC_OEI, irrep, 0, 1, lbl);
70
dpd_file2_mat_init(&T1b);
71
dpd_file2_mat_rd(&T1b);
72
for(h=0; h < nirreps; h++)
73
for(row=0; row < T1a.params->rowtot[h]; row++)
74
for(col=0; col < T1a.params->coltot[h^irrep]; col++)
75
error[0][word++] = T1a.matrix[h][row][col] - T1b.matrix[h][row][col];
76
dpd_file2_mat_close(&T1a);
77
dpd_file2_close(&T1a);
78
dpd_file2_mat_close(&T1b);
79
dpd_file2_close(&T1b);
81
sprintf(lbl, "New X_%s_%1s_IjAb (%5.3f)", pert, cart, omega);
82
dpd_buf4_init(&T2a, CC_LR, irrep, 0, 5, 0, 5, 0, lbl);
83
sprintf(lbl, "X_%s_%1s_IjAb (%5.3f)", pert, cart, omega);
84
dpd_buf4_init(&T2b, CC_LR, irrep, 0, 5, 0, 5, 0, lbl);
85
for(h=0; h < nirreps; h++) {
86
dpd_buf4_mat_irrep_init(&T2a, h);
87
dpd_buf4_mat_irrep_rd(&T2a, h);
88
dpd_buf4_mat_irrep_init(&T2b, h);
89
dpd_buf4_mat_irrep_rd(&T2b, h);
90
for(row=0; row < T2a.params->rowtot[h]; row++)
91
for(col=0; col < T2a.params->coltot[h^irrep]; col++)
92
error[0][word++] = T2a.matrix[h][row][col] - T2b.matrix[h][row][col];
93
dpd_buf4_mat_irrep_close(&T2a, h);
94
dpd_buf4_mat_irrep_close(&T2b, h);
99
start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
100
sprintf(lbl, "DIIS %s %1s Error Vectors", pert, cart);
101
psio_write(CC_DIIS_ERR, lbl , (char *) error[0],
102
vector_length*sizeof(double), start, &end);
104
/* Store the current amplitude vector on disk */
107
sprintf(lbl, "New X_%s_%1s_IA (%5.3f)", pert, cart, omega);
108
dpd_file2_init(&T1a, CC_OEI, irrep, 0, 1, lbl);
109
dpd_file2_mat_init(&T1a);
110
dpd_file2_mat_rd(&T1a);
111
for(h=0; h < nirreps; h++)
112
for(row=0; row < T1a.params->rowtot[h]; row++)
113
for(col=0; col < T1a.params->coltot[h^irrep]; col++)
114
error[0][word++] = T1a.matrix[h][row][col];
115
dpd_file2_mat_close(&T1a);
116
dpd_file2_close(&T1a);
118
sprintf(lbl, "New X_%s_%1s_IjAb (%5.3f)", pert, cart, omega);
119
dpd_buf4_init(&T2a, CC_LR, irrep, 0, 5, 0, 5, 0, lbl);
120
for(h=0; h < nirreps; h++) {
121
dpd_buf4_mat_irrep_init(&T2a, h);
122
dpd_buf4_mat_irrep_rd(&T2a, h);
123
for(row=0; row < T2a.params->rowtot[h]; row++)
124
for(col=0; col < T2a.params->coltot[h^irrep]; col++)
125
error[0][word++] = T2a.matrix[h][row][col];
126
dpd_buf4_mat_irrep_close(&T2a, h);
128
dpd_buf4_close(&T2a);
130
start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
131
sprintf(lbl, "DIIS %s %1s Amplitude Vectors", pert, cart);
132
psio_write(CC_DIIS_AMP, lbl , (char *) error[0],
133
vector_length*sizeof(double), start, &end);
135
/* If we haven't run through enough iterations, set the correct dimensions
136
for the extrapolation */
137
if(!(iter >= (nvector))) {
138
if(iter < 2) { /* Leave if we can't extrapolate at all */
139
dpd_free_block(error, 1, vector_length);
145
/* Build B matrix of error vector products */
146
vector = dpd_block_matrix(2, vector_length);
147
B = block_matrix(nvector+1,nvector+1);
148
for(p=0; p < nvector; p++) {
150
start = psio_get_address(PSIO_ZERO, p*vector_length*sizeof(double));
152
sprintf(lbl, "DIIS %s %1s Error Vectors", pert, cart);
153
psio_read(CC_DIIS_ERR, lbl, (char *) vector[0],
154
vector_length*sizeof(double), start, &end);
156
dot_arr(vector[0], vector[0], vector_length, &product);
160
for(q=0; q < p; q++) {
162
start = psio_get_address(PSIO_ZERO, q*vector_length*sizeof(double));
164
sprintf(lbl, "DIIS %s %1s Error Vectors", pert, cart);
165
psio_read(CC_DIIS_ERR, lbl, (char *) vector[1],
166
vector_length*sizeof(double), start, &end);
168
dot_arr(vector[1], vector[0], vector_length, &product);
170
B[p][q] = B[q][p] = product;
173
dpd_free_block(vector, 2, vector_length);
175
for(p=0; p < nvector; p++) {
180
B[nvector][nvector] = 0;
182
/* Find the maximum value in B and scale all its elements */
183
maximum = fabs(B[0][0]);
184
for(p=0; p < nvector; p++)
185
for(q=0; q < nvector; q++)
186
if(fabs(B[p][q]) > maximum) maximum = fabs(B[p][q]);
188
for(p=0; p < nvector; p++)
189
for(q=0; q < nvector; q++)
192
/* Build the constant vector */
193
C = init_array(nvector+1);
196
/* Solve the linear equations */
197
ipiv = init_int_array(nvector+1);
199
errcod = C_DGESV(nvector+1, 1, &(B[0][0]), nvector+1, &(ipiv[0]), &(C[0]), nvector+1);
201
fprintf(outfile, "\nError in DGESV return in diis.\n");
202
exit(PSI_RETURN_FAILURE);
205
/* Build a new amplitude vector from the old ones */
206
vector = dpd_block_matrix(1, vector_length);
207
for(p=0; p < vector_length; p++) error[0][p] = 0.0;
208
for(p=0; p < nvector; p++) {
210
start = psio_get_address(PSIO_ZERO, p*vector_length*sizeof(double));
212
sprintf(lbl, "DIIS %s %1s Amplitude Vectors", pert, cart);
213
psio_read(CC_DIIS_AMP, lbl, (char *) vector[0],
214
vector_length*sizeof(double), start, &end);
216
for(q=0; q < vector_length; q++)
217
error[0][q] += C[p] * vector[0][q];
220
dpd_free_block(vector, 1, vector_length);
222
/* Now place these elements into the DPD amplitude arrays */
224
sprintf(lbl, "New X_%s_%1s_IA (%5.3f)", pert, cart, omega);
225
dpd_file2_init(&T1a, CC_OEI, irrep, 0, 1, lbl);
226
dpd_file2_mat_init(&T1a);
227
for(h=0; h < nirreps; h++)
228
for(row=0; row < T1a.params->rowtot[h]; row++)
229
for(col=0; col < T1a.params->coltot[h^irrep]; col++)
230
T1a.matrix[h][row][col] = error[0][word++];
231
dpd_file2_mat_wrt(&T1a);
232
dpd_file2_mat_close(&T1a);
233
dpd_file2_close(&T1a);
235
sprintf(lbl, "New X_%s_%1s_IjAb (%5.3f)", pert, cart, omega);
236
dpd_buf4_init(&T2a, CC_LR, irrep, 0, 5, 0, 5, 0, lbl);
237
for(h=0; h < nirreps; h++) {
238
dpd_buf4_mat_irrep_init(&T2a, h);
239
for(row=0; row < T2a.params->rowtot[h]; row++)
240
for(col=0; col < T2a.params->coltot[h^irrep]; col++)
241
T2a.matrix[h][row][col] = error[0][word++];
242
dpd_buf4_mat_irrep_wrt(&T2a, h);
243
dpd_buf4_mat_irrep_close(&T2a, h);
245
dpd_buf4_close(&T2a);
247
/* Release memory and return */
248
/* free_matrix(vector, nvector); */
252
dpd_free_block(error, 1, vector_length);
added
removed
src/bin/ccresponse/diis.c
