← Back to branch summary

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

« back to all changes in this revision

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