← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/detci/oldcode/gencic.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/detci/oldcode/gencic.cc
 
1
/*
 
2
** GENCIC.CC
 
3
**
 
4
** This file contains code for the General CI (GENCI) module
 
5
**
 
6
** Matt Leininger
 
7
** Center for Computational Quantum Chemistry, UGA
 
8
** 1996
 
9
*/
 
10
 
 
11
 
 
12
#define EXTERN
 
13
 
 
14
/*** INCLUDES ***/
 
15
#include <cstdio>
 
16
#include <cstdlib>
 
17
#include <cmath>
 
18
#include <libciomr/libciomr.h>
 
19
#include "structs.h"
 
20
#include "globals.h"
 
21
#include "genci.h"
 
22
 
 
23
namespace psi { namespace detci {
 
24
 
 
25
/*** DEFINES ***/
 
26
#define NORM_TOL 1.0E-5
 
27
#define INDEX(i,j) ( (i>j) ? (ioff[(i)] + (j)): (ioff[(j)] + (i)) )
 
28
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
 
29
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
 
30
 
 
31
/*
 
32
** SCHMIDT_ADDOC(): Assume A is a orthogonal matrix.  
 
33
** This function Gram-Schmidt
 
34
** orthogonalizes a new vector v and adds it to matrix A.  A must contain
 
35
** a free row pointer for a new row.  Don't add orthogonalized v' if 
 
36
** norm(v') < NORM_TOL.
 
37
**
 
38
** David Sherrill, Feb 1994
 
39
** Matt Leininger, June 1996 Rewrote for out-of-core 
 
40
**
 
41
** Arguments:
 
42
**    buffer4 = buffer to store d or f vectors 
 
43
**    buffer5 = buffer to store b vectors 
 
44
**    buf_size = size of buffers
 
45
**    extra_buf = size of the last extra buffer of a vector
 
46
**    num_buf = number of buffers of size buf_size in a vector length N
 
47
**    f_index = 
 
48
**    N = dimension (length) of b, d, or f vectors 
 
49
**
 
50
** Returns: 1 if a vector is added to A, 0 otherwise
 
51
*/
 
52
 
 
53
int schmidt_addoc(double *buffer4, double *buffer5, int buf_size, int extra_buf,
 
54
                int num_buf, PSI_FPTR d_index, int N, int L,
 
55
                int b_file, int d_file)
 
56
{
 
57
   double *dotval, normval = 0.0, tval;
 
58
   int i, j, I;
 
59
   PSI_FPTR b_index = 0, jnk = 0, f_index;
 
60
 
 
61
   /* determine array of dot products b.f */
 
62
   f_index = d_index; 
 
63
   dotval = init_array(L); 
 
64
   for (i=0; i<num_buf; i++) { /* loop over num_buf's in b and f */
 
65
      wreadw(d_file, (char *) buffer4, sizeof(double)*buf_size, 
 
66
             d_index, &d_index); 
 
67
      b_index = i*buf_size*sizeof(double); 
 
68
      for (j=0; j<L; j++) { /* loop over b vectors */
 
69
         wreadw(b_file, (char *) buffer5, sizeof(double)*buf_size, 
 
70
                b_index, &jnk); 
 
71
         dot_arr(buffer4, buffer5, buf_size, &tval);
 
72
         dotval[j] += tval; 
 
73
         b_index += N*sizeof(double); 
 
74
         } /* end j loop */
 
75
      } /* end i loop */
 
76
   if (extra_buf != 0) {
 
77
   b_index = num_buf*buf_size*sizeof(double); 
 
78
   wreadw(d_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
79
          d_index, &d_index); 
 
80
   for (i=0; i<L; i++) {
 
81
      wreadw(b_file, (char *) buffer5, sizeof(double)*extra_buf, 
 
82
             b_index, &jnk);
 
83
      dot_arr(buffer4, buffer5, extra_buf, &tval);
 
84
      dotval[i] += tval; 
 
85
      b_index += N*sizeof(double);
 
86
      } /* end i loop */
 
87
     } /* end extra_buf */
 
88
 
 
89
   /* schmidt orthronormalize f to the set of b vectors */
 
90
   d_index -= N*sizeof(double); 
 
91
   for (i=0; i<num_buf; i++) {
 
92
      wreadw(d_file, (char *) buffer4, sizeof(double)*buf_size, 
 
93
             d_index, &d_index);
 
94
      b_index = i*buf_size*sizeof(double); 
 
95
      for (j=0; j<L; j++) { /* loop over b vectors */
 
96
         wreadw(b_file, (char *) buffer5, sizeof(double)*buf_size, 
 
97
                b_index, &jnk);
 
98
         b_index += N*sizeof(double);
 
99
         for (I=0; I<buf_size; I++) 
 
100
            buffer4[I] -= dotval[j]*buffer5[I]; 
 
101
         } /* end j loop */
 
102
      dot_arr(buffer4, buffer4, buf_size, &tval); 
 
103
      normval += tval; 
 
104
      wwritw(d_file, (char *) buffer4, sizeof(double)*buf_size, 
 
105
             f_index, &f_index); 
 
106
      } /* end i loop */
 
107
   if (extra_buf != 0) {
 
108
     b_index = num_buf*buf_size*sizeof(double);
 
109
     wreadw(d_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
110
            d_index, &d_index);
 
111
     for (i=0; i<L; i++) {
 
112
        wreadw(b_file, (char *) buffer5, sizeof(double)*extra_buf,
 
113
               b_index, &jnk);
 
114
        b_index += N*sizeof(double);
 
115
        for (I=0; I<extra_buf; I++)
 
116
           buffer4[I] -= dotval[i]*buffer5[I];
 
117
        } /* end i loop */
 
118
     dot_arr(buffer4, buffer4, extra_buf, &tval); 
 
119
     normval += tval; 
 
120
     wwritw(d_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
121
            f_index, &f_index); 
 
122
     } /* end extra_buf */
 
123
   normval = sqrt(normval); 
 
124
   /* fprintf(outfile,"\n normval = %15.10f\n", normval); */ 
 
125
 
 
126
  /* test norm of new b vector */
 
127
  free(dotval); 
 
128
  if (normval < NORM_TOL)
 
129
    return(0); 
 
130
  else {
 
131
      f_index -= N*sizeof(double); 
 
132
      b_index = N*L*sizeof(double); 
 
133
      for (i=0; i<num_buf; i++) {
 
134
         wreadw(d_file, (char *) buffer4, sizeof(double)*buf_size, 
 
135
                f_index, &f_index);
 
136
         for (I=0; I<buf_size; I++)  
 
137
            buffer4[I] /= normval; 
 
138
         wwritw(b_file, (char *) buffer4, sizeof(double)*buf_size, 
 
139
                b_index, &b_index); 
 
140
         } /* end i loop */
 
141
      if (extra_buf != 0) {
 
142
        wreadw(d_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
143
               f_index, &f_index);
 
144
        for (I=0; I<extra_buf; I++)  
 
145
           buffer4[I] /= normval; 
 
146
        wwritw(b_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
147
               b_index, &b_index); 
 
148
        } /* end extra_buf */
 
149
 
 
150
      return(1); 
 
151
      } /* end else */
 
152
 
 
153
}
 
154
 
 
155
 
 
156
/*
 
157
** NORMALIZE.C : Normalize a set of vectors
 
158
**
 
159
** Assume we're normalizing the ROWS
 
160
**
 
161
** David Sherrill, Feb 1994
 
162
** Matt Leininger, June 1996 modified for out-of-core vectors
 
163
*/
 
164
void v_normalize(double *A, PSI_FPTR index, int buf_size, 
 
165
               int extra_buf, int num_buf, int d_file)
 
166
{
 
167
   double normval = 0.0, tval;
 
168
   register int i, I;
 
169
   PSI_FPTR jnk = 0, index2; 
 
170
 
 
171
   index2 = index; 
 
172
  
 
173
   /* divide each row by the square root of its norm */
 
174
   for (i=0; i<num_buf; i++) {
 
175
      wreadw(d_file, (char *) A, sizeof(double)*buf_size, index, &index); 
 
176
      dot_arr(A, A, buf_size, &tval);
 
177
      normval += tval;
 
178
      } /* end i loop */
 
179
   if (extra_buf != 0) {
 
180
     wreadw(d_file, (char *) A, sizeof(double)*extra_buf, index, &index); 
 
181
     dot_arr(A, A, extra_buf, &tval);
 
182
     normval += tval;
 
183
     } /* end extra_buf */
 
184
   normval = sqrt(normval); 
 
185
 
 
186
   index = index2; 
 
187
   for (i=0; i<num_buf; i++) {
 
188
      wreadw(d_file, (char *) A, sizeof(double)*buf_size, index, &index);
 
189
      for (I=0; I<buf_size; I++) 
 
190
         A[I] /= normval;  
 
191
      wwritw(d_file, (char *) A, sizeof(double)*buf_size, index2, &index2); 
 
192
      } /* end i loop */
 
193
   if (extra_buf != 0) {
 
194
     wreadw(d_file, (char *) A, sizeof(double)*extra_buf, index, &index);
 
195
     for (I=0; I<extra_buf; I++)
 
196
        A[I] /= normval;
 
197
     wwritw(d_file, (char *) A, sizeof(double)*extra_buf, index2, &index2);  
 
198
     } /* end extra_buf */
 
199
 
 
200
}
 
201
 
 
202
/*
 
203
** V_SCHMIDT(): Assume A is a orthogonal matrix.  This function Gram-Schmidt
 
204
** orthogonalizes a set of given vectors.
 
205
** David Sherrill, Feb 1994
 
206
** Matt Leininger, June 1996 Rewrote for out-of-core 
 
207
**
 
208
** Arguments:
 
209
**    buffer4 = buffer to store b vectors 
 
210
**    buffer5 = buffer to store b vectors 
 
211
**    buf_size = size of buffers
 
212
**    extra_buf = size of the last extra buffer of a vector
 
213
**    num_buf = number of buffers of size buf_size in a vector length N
 
214
**    N = dimension (length) of b vectors 
 
215
**
 
216
*/
 
217
double *v_schmidt(double *buffer4, double *buffer5, int buf_size, int extra_buf,
 
218
                int num_buf, int N, int L, int b_file)
 
219
{
 
220
   double normval, tval, *dotval;
 
221
   int i, j, k, I, tridim;
 
222
   PSI_FPTR x_index = 0, jnk = 0, y_index, xwrit_index;
 
223
 
 
224
   /* determine array of dot products b.f */
 
225
   x_index = N*sizeof(double); 
 
226
   xwrit_index = x_index; 
 
227
   tridim = (L*(L+1))/2;
 
228
   dotval = init_array(tridim); 
 
229
 for (k=1; k<L; k++) { /* loop over b vectors */
 
230
   for (i=0; i<num_buf; i++) { /* loop over num_buf's in b vectors */
 
231
      wreadw(b_file, (char *) buffer4, sizeof(double)*buf_size, 
 
232
             x_index, &x_index); 
 
233
      y_index = i*buf_size*sizeof(double); 
 
234
      for (j=0; j<=k; j++) { /* loop over b vectors */
 
235
         wreadw(b_file, (char *) buffer5, sizeof(double)*buf_size, 
 
236
                y_index, &jnk); 
 
237
         dot_arr(buffer4, buffer5, buf_size, &tval);
 
238
         dotval[INDEX(k,j)] += tval; 
 
239
         y_index += N*sizeof(double); 
 
240
         } /* end j loop */
 
241
      } /* end i loop */
 
242
   if (extra_buf != 0) {
 
243
   y_index = num_buf*buf_size*sizeof(double); 
 
244
   wreadw(b_file, (char *) buffer4, sizeof(double)*extra_buf,
 
245
          x_index, &x_index); 
 
246
   for (i=0; i<=k; i++) {
 
247
      wreadw(b_file, (char *) buffer5, sizeof(double)*extra_buf, 
 
248
             y_index, &jnk);
 
249
      dot_arr(buffer4, buffer5, extra_buf, &tval);
 
250
      dotval[INDEX(i,k)] += tval; 
 
251
      y_index += N*sizeof(double);
 
252
      } /* end i loop */
 
253
     } /* end extra_buf */
 
254
 
 
255
   /* schmidt orthronormalize the set of b vectors */
 
256
   x_index -= N*sizeof(double); 
 
257
   normval = 0;
 
258
   for (i=0; i<num_buf; i++) {
 
259
      wreadw(b_file, (char *) buffer4, sizeof(double)*buf_size, 
 
260
             x_index, &x_index);
 
261
      y_index = i*buf_size*sizeof(double); 
 
262
      for (j=0; j<k; j++) { /* loop over b vectors */
 
263
         wreadw(b_file, (char *) buffer5, sizeof(double)*buf_size, 
 
264
                y_index, &jnk);
 
265
         y_index += N*sizeof(double);
 
266
         for (I=0; I<buf_size; I++) 
 
267
            buffer4[I] -= dotval[INDEX(j,k)]*buffer5[I]; 
 
268
         } /* end j loop */
 
269
      dot_arr(buffer4, buffer4, buf_size, &tval);
 
270
      normval += tval;
 
271
      wwritw(b_file, (char *) buffer4, sizeof(double)*buf_size, 
 
272
             xwrit_index, &xwrit_index); 
 
273
      } /* end i loop */
 
274
   if (extra_buf != 0) {
 
275
     y_index = num_buf*buf_size*sizeof(double);
 
276
     wreadw(b_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
277
            x_index, &x_index); 
 
278
     for (i=0; i<k; i++) {
 
279
        wreadw(b_file, (char *) buffer5, sizeof(double)*extra_buf, 
 
280
               y_index, &jnk);
 
281
        y_index += N*sizeof(double);
 
282
        for (I=0; I<extra_buf; I++)
 
283
           buffer4[I] -= dotval[INDEX(i,k)]*buffer5[I];
 
284
        } /* end i loop */
 
285
     dot_arr(buffer4, buffer4, extra_buf, &tval);
 
286
     normval += tval;
 
287
     wwritw(b_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
288
            xwrit_index, &xwrit_index); 
 
289
     } /* end extra_buf */
 
290
 
 
291
   /* Normalize vector */
 
292
   normval = sqrt(normval);
 
293
   x_index -= N*sizeof(double);
 
294
   for (i=0; i<num_buf; i++) {
 
295
      wreadw(b_file, (char *) buffer4, sizeof(double)*buf_size, 
 
296
             x_index, &jnk);
 
297
      for (I=0; I<buf_size; I++) 
 
298
         buffer4[I] /= normval;
 
299
      wwritw(b_file, (char *) buffer4, sizeof(double)*buf_size, 
 
300
             x_index, &x_index);
 
301
      } /* end i loop */
 
302
   if (extra_buf != 0) {
 
303
     wreadw(b_file, (char *) buffer4, sizeof(double)*extra_buf,
 
304
            x_index, &jnk);
 
305
     for (I=0; I<extra_buf; I++) 
 
306
         buffer4[I] /= normval;
 
307
      wwritw(b_file, (char *) buffer4, sizeof(double)*extra_buf, 
 
308
             x_index, &x_index);
 
309
     } /* end extra_buf */
 
310
   } /* end k loop */
 
311
 return(dotval);
 
312
}
 
313
 
 
314
void det2strings(BIGINT det, int *alp_code, int *alp_idx,
 
315
                 int *bet_code, int *bet_idx)
 
316
{
 
317
   int i;
 
318
 
 
319
   /* determine the CI block we're in */
 
320
   for (i=0; i<CIblks.num_blocks-1; i++) {
 
321
      if (CIblks.offset[i+1] > det) break;
 
322
      }
 
323
   *alp_code = CIblks.Ia_code[i];
 
324
   *bet_code = CIblks.Ib_code[i];
 
325
 
 
326
   *alp_idx = (det - CIblks.offset[i]) / CIblks.Ib_size[i];
 
327
   *bet_idx = (det - CIblks.offset[i]) % CIblks.Ib_size[i];
 
328
 
 
329
}
 
330
 
 
331
BIGINT strings2det(int alp_code, int alp_idx, int bet_code, int bet_idx)
 
332
{
 
333
   int blknum;
 
334
   BIGINT addr;
 
335
 
 
336
   blknum = CIblks.decode[alp_code][bet_code];
 
337
   if (blknum == -1) {
 
338
     fprintf(outfile, "CIvect::strings2det failed --- invalid block\n");
 
339
     exit(1);
 
340
   }
 
341
 
 
342
   addr = CIblks.offset[blknum];
 
343
   addr += alp_idx * CIblks.Ib_size[blknum] + bet_idx;
 
344
 
 
345
   return(addr);
 
346
}
 
347
 
 
348
 
 
349
/*
 
350
** unit_guess - uses a unit guess for the first b vector in the
 
351
**              davidson-liu algorithm
 
352
**
 
353
** 
 
354
*/
 
355
void unit_guess(int alp_code, int alp_idx, int bet_code, int bet_idx,
 
356
                int switch_buf3, double *buffer, int buf_size,
 
357
                int num_buf, int extra_buf, PSI_FPTR b_file,
 
358
                PSI_FPTR b_writ, int M, int N)
 
359
{
 
360
 int i;
 
361
 register int j;
 
362
  
 
363
 
 
364
   /* Form initial guess b vector */
 
365
   i = strings2det(CalcInfo.ref_alp_list, CalcInfo.ref_alp_rel,
 
366
                   CalcInfo.ref_bet_list, CalcInfo.ref_bet_rel);
 
367
   if (!switch_buf3) {
 
368
     buffer[i] = 1.0;
 
369
     wwritw(b_file, (char *) buffer, sizeof(double)*buf_size,
 
370
            b_writ, &b_writ);
 
371
     buffer[i] = 0.0;
 
372
     v_normalize(buffer, (b_writ-M*N*sizeof(double)), buf_size, 0,
 
373
                 1, b_file);
 
374
     }
 
375
   else if (switch_buf3) {
 
376
     for (j=0; j<num_buf; j++) {
 
377
        if ((i >= (j*buf_size)) && (i < ((j+1)*buf_size)))
 
378
          buffer[i-j*buf_size] = 1.0;
 
379
        wwritw(b_file, (char *) buffer, sizeof(double)*buf_size,
 
380
               b_writ, &b_writ);
 
381
        if ((i >= (j*buf_size)) && (i < ((j+1)*buf_size)))
 
382
          buffer[i-j*buf_size] = 0.0;
 
383
        }
 
384
     if (extra_buf != 0) {
 
385
       if ((i >= num_buf*buf_size) && (i < (extra_buf+num_buf*buf_size)))
 
386
         buffer[i-num_buf*buf_size] = 1.0;
 
387
       wwritw(b_file, (char *) buffer, sizeof(double)*extra_buf,
 
388
              b_writ, &b_writ);
 
389
       if ((i > num_buf*buf_size) && (i < (extra_buf+num_buf*buf_size)))
 
390
         buffer[i-num_buf*buf_size] = 0.0;
 
391
       }
 
392
     v_normalize(buffer, (b_writ-M*N*sizeof(double)), buf_size, extra_buf,
 
393
                 num_buf, b_file);
 
394
     }
 
395
 
 
396
   /* b_writ -= N*sizeof(double);
 
397
        wreadw(b_file, (char *) buffer, sizeof(double)*buf_size,
 
398
               b_writ, &b_writ);
 
399
        for (I=0; I<buf_size; I++)
 
400
           if (buffer[I] != 0.0) fprintf(outfile,"b[0][%d] = %lf\n\n",
 
401
               I, buffer[I]);
 
402
        printf("Done reading in b trial vector for check.\n"); */
 
403
 
 
404
}
 
405
 
 
406
 
 
407
/*
 
408
** max_element  - determines the maximum value in a list of diagonal
 
409
**               matrix elements.
 
410
** buffer       - array of matrix elements
 
411
** num_elements - number of elements in array buffer
 
412
** max          - pointer to max value
 
413
** max_num      - the element number to the maximum value
 
414
*/
 
415
void max_element(double *buffer, int num_elements, double *max, int *max_num)
 
416
{
 
417
 register int i;
 
418
 
 
419
 (*max) = buffer[0];
 
420
 (*max_num) = 0; 
 
421
 for (i=1; i<num_elements; i++) {
 
422
    if (buffer[i] > (*max)) {
 
423
      (*max) = buffer[i];
 
424
      (*max_num) = i;
 
425
      }
 
426
    } 
 
427
 
 
428
}
 
429
 
 
430
 
 
431
/*
 
432
** min_element  - determines the minimum value in a list of diagonal
 
433
**               matrix elements.
 
434
** buffer       - array of matrix elements
 
435
** num_elements - number of elements in array buffer
 
436
** min          - pointer to max value
 
437
** min_num      - the element number to the minimum value
 
438
*/
 
439
void min_element(double *buffer, int num_elements, double *min, int *min_num)
 
440
{
 
441
 register int i;
 
442
 
 
443
 (*min) = buffer[0];
 
444
 (*min_num) = 0;
 
445
 for (i=1; i<num_elements; i++) {
 
446
    if (buffer[i] < (*min)) {
 
447
      (*min) = buffer[i];
 
448
      (*min_num) = i;
 
449
      }
 
450
    }
 
451
 
 
452
}
 
453
 
 
454
 
 
455
/*
 
456
**  read_c() - reads in the CI vector for a restart of the calculation
 
457
**
 
458
*/
 
459
void read_c(int switch_buf3, double *buffer, int buf_size, int num_buf,
 
460
            int extra_buf, int b_file, PSI_FPTR b_writ, int c_file, 
 
461
            PSI_FPTR c_index)
 
462
{
 
463
 register int i;
 
464
 
 
465
 if (!switch_buf3) {
 
466
   wreadw(c_file, (char *) buffer, sizeof(double)*buf_size,
 
467
          c_index, &c_index);
 
468
   wwritw(b_file, (char *) buffer, sizeof(double)*buf_size,
 
469
          b_writ, &b_writ);
 
470
   }
 
471
 else if (switch_buf3) {
 
472
   for (i=0; i<num_buf; i++) {
 
473
      wreadw(c_file, (char *) buffer, sizeof(double)*buf_size,
 
474
          c_index, &c_index);
 
475
      wwritw(b_file, (char *) buffer, sizeof(double)*buf_size,
 
476
          b_writ, &b_writ);
 
477
      }
 
478
   if (extra_buf != 0) {
 
479
      wreadw(c_file, (char *) buffer, sizeof(double)*extra_buf,
 
480
          c_index, &c_index);
 
481
      wwritw(b_file, (char *) buffer, sizeof(double)*extra_buf,
 
482
          b_writ, &b_writ);
 
483
     }
 
484
   }
 
485
 
 
486
}
 
487
 
 
488
}} // namespace psi::detci
 
489