← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/detci/s3v.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/s3v.c
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <libciomr/libciomr.h>
4
 
#include <libqt/qt.h>
5
 
#include "structs.h"
6
 
#define EXTERN
7
 
#include "globals.h"
8
 
#include <pthread.h>
9
 
#include "tpool.h"
10
 
 
11
 
 
12
 
/* Global variables necessary for pthreads */
13
 
pthread_mutex_t inc_Ia_mutex;
14
 
struct stringwr *Ia_global;
15
 
int Ia_idx_global;
16
 
 
17
 
 
18
 
int form_ilist(struct stringwr *alplist, int Ja_list, int nas, int kl,
19
 
   int *L, int *R, double *Sgn);
20
 
int form_ilist_rotf(int *Cnt, int **Ridx, signed char **Sn, int **Ij,
21
 
   int nas, int kl, int *L, int *R, double *Sgn);
22
 
void s3_block_vdiag_pthread(void *threadarg);
23
 
void s3_block_v_pthread(void *threadarg);
24
 
 
25
 
 
26
 
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
27
 
 
28
 
 
29
 
/*
30
 
** S3_BLOCK_VDIAG()
31
 
**
32
 
** Calculate a block of the sigma3 vector in equation (9c) of
33
 
** Olsen, Roos, et al.  For diagonal blocks of sigma.
34
 
**
35
 
** currently assumes that (ij|ij)'s have not been halved
36
 
** Try to get the Olsen vector version working....again!!!!
37
 
*/
38
 
void s3_block_vdiag(struct stringwr *alplist, struct stringwr *betlist,
39
 
      double **C, double **S, double *tei, int nas, int nbs, int cnas,
40
 
      int Ib_list, int Ja_list, int Jb_list, int Ib_sym, int Jb_sym,
41
 
      double **Cprime, double *F, double *V, double *Sgn, int *L, int *R)
42
 
{
43
 
  struct stringwr *Ia;
44
 
  unsigned int Ia_ex;
45
 
  int ij, i, j, t, kl, I, J, RJ;
46
 
  double tval, VS, *CprimeI0, *CI0;
47
 
  int jlen, Jacnt, *Iaij, *orbsym, norbs, Ia_idx;
48
 
  unsigned int *Iaridx;
49
 
  signed char *Iasgn;
50
 
  double *Tptr;
51
 
  struct pthreads_s3diag **thread_info;
52
 
  int npthreads, rc, status;
53
 
  pthread_t *thread;
54
 
   
55
 
  npthreads = Parameters.nthreads-1;  /* subtract out the main thread */
56
 
 
57
 
  thread = (pthread_t *) malloc(sizeof(pthread_t)*Parameters.nthreads);
58
 
 
59
 
  thread_info = (struct pthreads_s3diag **)
60
 
                malloc(sizeof(struct pthreads_s3diag *) * nas);
61
 
 
62
 
  for (i=0; i<nas; i++) {
63
 
      thread_info[i] = (struct pthreads_s3diag *)
64
 
                       malloc(sizeof(struct pthreads_s3diag));
65
 
    }
66
 
  
67
 
  norbs = CalcInfo.num_ci_orbs;
68
 
  orbsym = CalcInfo.orbsym + CalcInfo.num_fzc_orbs;
69
 
 
70
 
  /* loop over i, j */
71
 
  for (i=0; i<norbs; i++) {
72
 
      for (j=0; j<=i; j++) {
73
 
          if ((orbsym[i] ^ orbsym[j] ^ Jb_sym ^ Ib_sym) != 0) continue;
74
 
          ij = ioff[i] + j;
75
 
          jlen = form_ilist(betlist, Jb_list, nbs, ij, L, R, Sgn);
76
 
           
77
 
          if (!jlen) continue;
78
 
          /*  fprintf(outfile,"S3_BLOCK_VDIAG: ij = %d\t jlen = %d\n", ij, jlen);
79
 
           */
80
 
          Tptr = tei + ioff[ij];
81
 
       
82
 
          /* gather operation */
83
 
          for (I=0; I<cnas; I++) {
84
 
              CprimeI0 = Cprime[I];
85
 
              CI0 = C[I];
86
 
              for (J=0; J<jlen; J++) {
87
 
                  tval = Sgn[J];
88
 
                  CprimeI0[J] = CI0[L[J]] * tval;
89
 
                }
90
 
            }
91
 
 
92
 
 
93
 
          /* loop over Ia */
94
 
          if (Parameters.pthreads) {
95
 
              detci_time.s3_mt_before_time = wall_time_new();
96
 
              tpool_queue_open(thread_pool);
97
 
              for (Ia=alplist, Ia_idx=0; Ia_idx<nas; Ia_idx++, Ia++) {
98
 
                  thread_info[Ia_idx]->nas = nas;
99
 
                  thread_info[Ia_idx]->jlen = jlen;
100
 
                  thread_info[Ia_idx]->ij = ij;
101
 
                  thread_info[Ia_idx]->Cprime = Cprime;    
102
 
                  thread_info[Ia_idx]->Ja_list = Ja_list;
103
 
                  thread_info[Ia_idx]->Tptr = Tptr;
104
 
                  thread_info[Ia_idx]->S = S;
105
 
                  thread_info[Ia_idx]->R = R;
106
 
                  thread_info[Ia_idx]->Ia_local = Ia;
107
 
                  thread_info[Ia_idx]->Ia_idx_local = Ia_idx;
108
 
                  tpool_add_work(thread_pool , s3_block_vdiag_pthread, (void *) thread_info[Ia_idx]); 
109
 
                }
110
 
              tpool_queue_close(thread_pool, 1);
111
 
              detci_time.s3_mt_after_time = wall_time_new();
112
 
              detci_time.s3_mt_total_time += detci_time.s3_mt_after_time - detci_time.s3_mt_before_time;
113
 
 
114
 
            }
115
 
          else {
116
 
              for (Ia=alplist, Ia_idx=0; Ia_idx<nas; Ia_idx++, Ia++) {
117
 
 
118
 
                  /* loop over excitations E^a_{kl} from |A(I_a)> */
119
 
                  Jacnt = Ia->cnt[Ja_list];
120
 
                  Iaridx = Ia->ridx[Ja_list];
121
 
                  Iasgn = Ia->sgn[Ja_list];
122
 
                  Iaij = Ia->ij[Ja_list];
123
 
 
124
 
                  zero_arr(V, jlen);
125
 
/*                  fprintf(outfile,"Ia = %x\t Ia_idx = %d\n", Ia, Ia_idx);
126
 
                    fflush(outfile);
127
 
 
128
 
                    if (Jacnt) {
129
 
                    fprintf(outfile,"S3_BLOCK_VDIAG: Jacnt = %d\n", Jacnt);
130
 
                    fflush(outfile);
131
 
                    }
132
 
*/
133
 
                  for (Ia_ex=0; Ia_ex < Jacnt && (kl = *Iaij++)<=ij; Ia_ex++) {
134
 
                      I = *Iaridx++;
135
 
                      tval = *Iasgn++;
136
 
                      if (ij == kl) tval *= 0.5;
137
 
                      VS = Tptr[kl] * tval;
138
 
                      CprimeI0 = Cprime[I];
139
 
           
140
 
#ifdef USE_BLAS
141
 
                      C_DAXPY(jlen, VS, CprimeI0, 1, V, 1);
142
 
#else
143
 
                      for (J=0; J<jlen; J++) {
144
 
                          V[J] += VS * CprimeI0[J];
145
 
                        }
146
 
#endif
147
 
           
148
 
                    }
149
 
 
150
 
         
151
 
                  /* scatter */
152
 
                  for (J=0; J<jlen; J++) {
153
 
                      RJ = R[J];
154
 
                      S[Ia_idx][RJ] += V[J];
155
 
                    }
156
 
 
157
 
                } /* end loop over Ia */
158
 
            }
159
 
       
160
 
        } /* end loop over j */
161
 
    } /* end loop over i */
162
 
 
163
 
  for (i=0; i<nas; i++) free(thread_info[i]);
164
 
  free(thread);
165
 
  
166
 
}              
167
 
 
168
 
 
169
 
/*
170
 
** S3_BLOCK_VDIAG_PTHREAD()
171
 
**
172
 
** Multithreaded component of the S3_BLOCK_VDIAG routine
173
 
**
174
 
** Calculate a block of the sigma3 vector in equation (9c) of
175
 
** Olsen, Roos, et al.  For diagonal blocks of sigma.
176
 
**
177
 
** currently assumes that (ij|ij)'s have not been halved
178
 
** Try to get the Olsen vector version working....again!!!!
179
 
*/
180
 
void s3_block_vdiag_pthread(void *threadarg)
181
 
{
182
 
 
183
 
  struct stringwr *Ia_local;
184
 
  unsigned int Ia_ex;
185
 
  int I, J, RJ, kl;
186
 
  double tval, VS, *CprimeI0, *CI0;
187
 
  int Jacnt, *Iaij;
188
 
  unsigned int *Iaridx;
189
 
  signed char *Iasgn;
190
 
  struct pthreads_s3diag *thread_info_i;
191
 
 
192
 
  int nas, jlen, ij, Ja_list;
193
 
  double **S, **Cprime, *Tptr, *V;
194
 
  int *R, Ia_idx_local;
195
 
  int loop_var = 0;
196
 
  int thread_id;
197
 
  
198
 
  thread_info_i = (struct pthreads_s3diag *) threadarg;
199
 
  nas = thread_info_i->nas;
200
 
  jlen = thread_info_i->jlen;
201
 
  ij = thread_info_i->ij;
202
 
  Ja_list = thread_info_i->Ja_list;
203
 
  S = thread_info_i->S;
204
 
  Cprime = thread_info_i->Cprime;
205
 
  Tptr = thread_info_i->Tptr;
206
 
  R = thread_info_i->R;
207
 
  Ia_idx_local = thread_info_i->Ia_idx_local;
208
 
  Ia_local = thread_info_i->Ia_local;
209
 
  
210
 
  V = init_array(jlen);
211
 
  /* loop over excitations E^a_{kl} from |A(I_a)> */
212
 
  Jacnt = Ia_local->cnt[Ja_list];
213
 
  Iaridx = Ia_local->ridx[Ja_list];
214
 
  Iasgn = Ia_local->sgn[Ja_list];
215
 
  Iaij = Ia_local->ij[Ja_list];
216
 
 
217
 
/*
218
 
  fprintf(outfile,"Ia_local = %x\t Ia_idx_local = %d\n", Ia_local, Ia_idx_local);
219
 
  fflush(outfile);
220
 
  
221
 
  fprintf(outfile,"Jacnt = %d\n", Jacnt);
222
 
*/
223
 
      
224
 
  zero_arr(V, jlen);
225
 
         
226
 
  for (Ia_ex=0; Ia_ex < Jacnt && (kl = *Iaij++)<=ij; Ia_ex++) {
227
 
      I = *Iaridx++;
228
 
      tval = *Iasgn++;
229
 
      if (ij == kl) tval *= 0.5;
230
 
      VS = Tptr[kl] * tval;
231
 
      CprimeI0 = Cprime[I];
232
 
           
233
 
#ifdef USE_BLAS
234
 
      C_DAXPY(jlen, VS, CprimeI0, 1, V, 1);
235
 
#else
236
 
      for (J=0; J<jlen; J++) {
237
 
          V[J] += VS * CprimeI0[J];
238
 
        }
239
 
#endif
240
 
    }
241
 
 
242
 
  /* scatter */
243
 
  for (J=0; J<jlen; J++) {
244
 
      RJ = R[J];
245
 
      S[Ia_idx_local][RJ] += V[J];
246
 
    }
247
 
 
248
 
  free(V);
249
 
  /* return 0; */
250
 
  
251
 
}              
252
 
 
253
 
 
254
 
/*
255
 
** S3_BLOCK_V()
256
 
**
257
 
** Calculate a block of the sigma3 vector in equation (9c) of
258
 
** Olsen, Roos, et al.  For non-diagonal blocks of s3
259
 
**
260
 
*/
261
 
void s3_block_v(struct stringwr *alplist, struct stringwr *betlist,
262
 
      double **C, double **S, double *tei, int nas, int nbs, int cnas,
263
 
      int Ib_list, int Ja_list, int Jb_list, int Ib_sym, int Jb_sym,
264
 
      double **Cprime, double *F, double *V, double *Sgn, int *L, int *R)
265
 
{
266
 
   struct stringwr *Ia;
267
 
   unsigned int Ia_ex;
268
 
   int ij, i, j, kl, ijkl, I, J, RJ;
269
 
   double tval, VS, *CprimeI0, *CI0;
270
 
   int jlen, Ia_idx, Jacnt, *Iaij, *orbsym, norbs;
271
 
   unsigned int *Iaridx;
272
 
   signed char *Iasgn;
273
 
   double *Tptr;
274
 
   struct pthreads_s3diag **thread_info;
275
 
   int t, npthreads, rc, status;
276
 
   pthread_t *thread;
277
 
   
278
 
   norbs = CalcInfo.num_ci_orbs;
279
 
   orbsym = CalcInfo.orbsym + CalcInfo.num_fzc_orbs;
280
 
 
281
 
   thread = (pthread_t *) malloc(sizeof(pthread_t)*Parameters.nthreads);
282
 
   thread_info = (struct pthreads_s3diag **)
283
 
                  malloc(sizeof(struct pthreads_s3diag *) * nas);
284
 
   for (i=0; i<nas; i++) {
285
 
       thread_info[i] = (struct pthreads_s3diag *)
286
 
                        malloc(sizeof(struct pthreads_s3diag));
287
 
     }
288
 
   
289
 
   /* loop over i, j */
290
 
   for (i=0; i<norbs; i++) {
291
 
     for (j=0; j<=i; j++) {
292
 
       if ((orbsym[i] ^ orbsym[j] ^ Jb_sym ^ Ib_sym) != 0) continue;
293
 
       ij = ioff[i] + j;
294
 
       jlen = form_ilist(betlist, Jb_list, nbs, ij, L, R, Sgn);
295
 
 
296
 
       if (!jlen) continue;
297
 
 
298
 
       Tptr = tei + ioff[ij];
299
 
 
300
 
       /* gather operation */
301
 
       for (I=0; I<cnas; I++) {
302
 
         CprimeI0 = Cprime[I];
303
 
         CI0 = C[I];
304
 
         for (J=0; J<jlen; J++) {
305
 
           tval = Sgn[J];
306
 
           CprimeI0[J] = CI0[L[J]] * tval;
307
 
         }
308
 
       }
309
 
 
310
 
 
311
 
       /* loop over Ia */
312
 
       if (Parameters.pthreads) {
313
 
           detci_time.s3_mt_before_time = wall_time_new();
314
 
           tpool_queue_open(thread_pool);
315
 
           for (Ia=alplist, Ia_idx=0; Ia_idx<nas; Ia_idx++, Ia++) {
316
 
               thread_info[Ia_idx]->nas = nas;
317
 
               thread_info[Ia_idx]->jlen = jlen;
318
 
               thread_info[Ia_idx]->ij = ij;
319
 
               thread_info[Ia_idx]->Cprime = Cprime;    
320
 
               thread_info[Ia_idx]->Ja_list = Ja_list;
321
 
               thread_info[Ia_idx]->Tptr = tei;
322
 
               thread_info[Ia_idx]->S = S;
323
 
               thread_info[Ia_idx]->R = R;
324
 
               thread_info[Ia_idx]->Ia_local = Ia;
325
 
               thread_info[Ia_idx]->Ia_idx_local = Ia_idx;
326
 
               tpool_add_work(thread_pool, s3_block_v_pthread, (void *) thread_info[Ia_idx]);
327
 
             }
328
 
           tpool_queue_close(thread_pool, 1);
329
 
           detci_time.s3_mt_after_time = wall_time_new();
330
 
           detci_time.s3_mt_total_time += detci_time.s3_mt_after_time - detci_time.s3_mt_before_time;
331
 
 
332
 
         }
333
 
       else {
334
 
           for (Ia=alplist, Ia_idx=0; Ia_idx<nas; Ia_idx++, Ia++) {
335
 
 
336
 
               /* loop over excitations E^a_{kl} from |A(I_a)> */
337
 
               Jacnt = Ia->cnt[Ja_list];
338
 
               Iaridx = Ia->ridx[Ja_list];
339
 
               Iasgn = Ia->sgn[Ja_list];
340
 
               Iaij = Ia->ij[Ja_list];
341
 
 
342
 
               zero_arr(V, jlen);
343
 
 
344
 
               for (Ia_ex=0; Ia_ex < Jacnt; Ia_ex++) {
345
 
                   kl = *Iaij++;
346
 
                   I = *Iaridx++;
347
 
                   tval = *Iasgn++;
348
 
                   ijkl = INDEX(ij,kl);
349
 
                   VS = tval * tei[ijkl];
350
 
                   CprimeI0 = Cprime[I];
351
 
 
352
 
#ifdef USE_BLAS
353
 
                   C_DAXPY(jlen, VS, CprimeI0, 1, V, 1);
354
 
#else
355
 
                   for (J=0; J<jlen; J++) {
356
 
                       V[J] += VS * CprimeI0[J];
357
 
                     }
358
 
#endif
359
 
 
360
 
                 }
361
 
 
362
 
 
363
 
               /* scatter */
364
 
               for (J=0; J<jlen; J++) {
365
 
                   RJ = R[J];
366
 
                   S[Ia_idx][RJ] += V[J];
367
 
                 }
368
 
 
369
 
             } /* end loop over Ia */
370
 
         }
371
 
       
372
 
     } /* end loop over j */
373
 
   } /* end loop over i */
374
 
  for (i=0; i<nas; i++) free(thread_info[i]);
375
 
  free(thread);
376
 
   
377
 
}
378
 
 
379
 
/*
380
 
** S3_BLOCK_V_PTHREAD()
381
 
**
382
 
** Multithreaded component of the S3_BLOCK_V routine
383
 
**
384
 
** Calculate a block of the sigma3 vector in equation (9c) of
385
 
** Olsen, Roos, et al.  For non-diagonal blocks of s3
386
 
**
387
 
*/
388
 
void s3_block_v_pthread(void *threadarg)
389
 
{
390
 
 
391
 
  struct stringwr *Ia_local;
392
 
  unsigned int Ia_ex;
393
 
  int I, J, RJ, kl, ijkl;
394
 
  double tval, VS, *CprimeI0, *CI0;
395
 
  int Jacnt, *Iaij;
396
 
  unsigned int *Iaridx;
397
 
  signed char *Iasgn;
398
 
  struct pthreads_s3diag *thread_info_i;
399
 
 
400
 
  int nas, jlen, ij, Ja_list;
401
 
  double **S, **Cprime, *tei, *V, *Tptr;
402
 
  int *R, Ia_idx_local;
403
 
  int loop_var = 0;
404
 
  int work_units = 0;
405
 
  int thread_id;
406
 
 
407
 
  thread_info_i = (struct pthreads_s3diag *) threadarg;
408
 
  nas = thread_info_i->nas;
409
 
  jlen = thread_info_i->jlen;
410
 
  ij = thread_info_i->ij;
411
 
  Ja_list = thread_info_i->Ja_list;
412
 
  S = thread_info_i->S;
413
 
  Cprime = thread_info_i->Cprime;
414
 
  tei = thread_info_i->Tptr;
415
 
  R = thread_info_i->R;
416
 
  thread_id = thread_info_i->thread_id;
417
 
  Ia_idx_local = thread_info_i->Ia_idx_local;
418
 
  Ia_local = thread_info_i->Ia_local;
419
 
  
420
 
  V = init_array(jlen);
421
 
 
422
 
  /* loop over excitations E^a_{kl} from |A(I_a)> */
423
 
  Jacnt = Ia_local->cnt[Ja_list];
424
 
  Iaridx = Ia_local->ridx[Ja_list];
425
 
  Iasgn = Ia_local->sgn[Ja_list];
426
 
  Iaij = Ia_local->ij[Ja_list];
427
 
 
428
 
  work_units++;
429
 
      
430
 
  zero_arr(V, jlen);
431
 
         
432
 
  for (Ia_ex=0; Ia_ex < Jacnt; Ia_ex++) {
433
 
      kl = *Iaij++;
434
 
      I = *Iaridx++;
435
 
      tval = *Iasgn++;
436
 
      ijkl = INDEX(ij,kl);
437
 
      VS = tval * tei[ijkl];
438
 
      CprimeI0 = Cprime[I];
439
 
 
440
 
#ifdef USE_BLAS
441
 
      C_DAXPY(jlen, VS, CprimeI0, 1, V, 1);
442
 
#else
443
 
      for (J=0; J<jlen; J++) {
444
 
          V[J] += VS * CprimeI0[J];
445
 
        }
446
 
#endif
447
 
    }
448
 
 
449
 
  /* scatter */
450
 
  for (J=0; J<jlen; J++) {
451
 
      RJ = R[J];
452
 
      S[Ia_idx_local][RJ] += V[J];
453
 
    }
454
 
 
455
 
  free(V);
456
 
  /* return 0; */
457
 
  
458
 
}
459
 
 
460
 
int form_ilist(struct stringwr *alplist, int Ja_list, int nas, int kl,
461
 
   int *L, int *R, double *Sgn)
462
 
{
463
 
 
464
 
  int inum=0, Ia_idx, Ia_ex, Iacnt, ij;
465
 
  int *Iaij;
466
 
  struct stringwr *Ia;
467
 
  unsigned int *Iaridx;
468
 
  signed char *Iasgn;
469
 
 
470
 
  /* loop over Ia */
471
 
  for (Ia=alplist, Ia_idx=0; Ia_idx < nas; Ia_idx++,Ia++) {
472
 
 
473
 
      /* loop over excitations E^a_{kl} from |A(I_a)> */
474
 
 
475
 
      Iacnt = Ia->cnt[Ja_list];
476
 
      if (!Iacnt) continue;
477
 
      Iaridx = Ia->ridx[Ja_list];
478
 
      Iasgn = Ia->sgn[Ja_list];
479
 
      Iaij = Ia->ij[Ja_list];
480
 
      Ia_ex=0;
481
 
      while (Ia_ex < Iacnt && (ij = *Iaij++)<kl) Ia_ex++;
482
 
      if (ij == kl) {
483
 
          *R++ = Ia_idx;
484
 
          *L++ = Iaridx[Ia_ex];
485
 
          *Sgn++ = (double) Iasgn[Ia_ex];
486
 
          inum++;
487
 
        }
488
 
    }  /* end loop over Ia */
489
 
/*  if(inum) {
490
 
      fprintf(outfile,"form_ilist: nas = %d\n", nas);
491
 
      fprintf(outfile,"form_ilist: jlen = %d\n", inum);
492
 
      fflush(outfile);
493
 
    }
494
 
*/
495
 
    return(inum);
496
 
}
497
 
 
498
 
 
499
 
/*
500
 
** S3_BLOCK_VDIAG_ROTF()
501
 
**
502
 
** Calculate a block of the sigma3 vector in equation (9c) of
503
 
** Olsen, Roos, et al.  For diagonal blocks of sigma.
504
 
**
505
 
** currently assumes that (ij|ij)'s have not been halved
506
 
** Try to get the Olsen vector version working....again!!!!
507
 
*/
508
 
 
509
 
void s3_block_vdiag_rotf(int *Cnt[2], int **Ij[2], int **Ridx[2],
510
 
      signed char **Sn[2], double **C, double **S, 
511
 
      double *tei, int nas, int nbs, int cnas,
512
 
      int Ib_list, int Ja_list, int Jb_list, int Ib_sym, int Jb_sym,
513
 
      double **Cprime, double *F, double *V, double *Sgn, int *L, int *R)
514
 
{
515
 
   int Ia_ex;
516
 
   int ij, i, j, kl, I, J, RJ;
517
 
   double tval, VS, *CprimeI0, *CI0;
518
 
   int jlen, Ia_idx, Jacnt, *Iaij, *orbsym, norbs;
519
 
   int *Iaridx;
520
 
   signed char *Iasgn;
521
 
   double *Tptr;
522
 
   
523
 
   norbs = CalcInfo.num_ci_orbs;
524
 
   orbsym = CalcInfo.orbsym + CalcInfo.num_fzc_orbs;
525
 
 
526
 
   /* loop over i, j */
527
 
   for (i=0; i<norbs; i++) {
528
 
     for (j=0; j<=i; j++) {
529
 
       if ((orbsym[i] ^ orbsym[j] ^ Jb_sym ^ Ib_sym) != 0) continue;
530
 
       ij = ioff[i] + j;
531
 
       jlen = form_ilist_rotf(Cnt[1], Ridx[1], Sn[1], Ij[1], 
532
 
          nbs, ij, L, R, Sgn);
533
 
       
534
 
       if (!jlen) continue;
535
 
       
536
 
       Tptr = tei + ioff[ij];
537
 
       
538
 
       /* gather operation */
539
 
       for (I=0; I<cnas; I++) {
540
 
         CprimeI0 = Cprime[I];
541
 
         CI0 = C[I];
542
 
         for (J=0; J<jlen; J++) {
543
 
           tval = Sgn[J];
544
 
           CprimeI0[J] = CI0[L[J]] * tval;
545
 
         }
546
 
       }
547
 
 
548
 
 
549
 
       /* loop over Ia */
550
 
       for (Ia_idx=0; Ia_idx<nas; Ia_idx++) {
551
 
 
552
 
         /* loop over excitations E^a_{kl} from |A(I_a)> */
553
 
         Jacnt = Cnt[0][Ia_idx];
554
 
         Iaridx = Ridx[0][Ia_idx];
555
 
         Iasgn = Sn[0][Ia_idx];
556
 
         Iaij = Ij[0][Ia_idx];
557
 
         
558
 
         zero_arr(V, jlen);
559
 
         
560
 
         /* rotf doesn't yet ensure kl's in order */
561
 
         for (Ia_ex=0; Ia_ex < Jacnt; Ia_ex++) {
562
 
           kl = *Iaij++;
563
 
           I = *Iaridx++;
564
 
           tval = *Iasgn++;
565
 
           if (kl > ij) continue;
566
 
           if (ij == kl) tval *= 0.5;
567
 
           VS = Tptr[kl] * tval;
568
 
           CprimeI0 = Cprime[I];
569
 
           
570
 
         #ifdef USE_BLAS
571
 
           C_DAXPY(jlen, VS, CprimeI0, 1, V, 1);
572
 
         #else
573
 
           for (J=0; J<jlen; J++) {
574
 
             V[J] += VS * CprimeI0[J];
575
 
           }
576
 
         #endif
577
 
         }
578
 
 
579
 
         
580
 
         /* scatter */
581
 
         for (J=0; J<jlen; J++) {
582
 
           RJ = R[J];
583
 
           S[Ia_idx][RJ] += V[J];
584
 
         }
585
 
 
586
 
       } /* end loop over Ia */
587
 
 
588
 
     } /* end loop over j */
589
 
   } /* end loop over i */
590
 
   
591
 
}              
592
 
 
593
 
 
594
 
/*
595
 
** S3_BLOCK_VROTF()
596
 
**
597
 
** Calculate a block of the sigma3 vector in equation (9c) of
598
 
** Olsen, Roos, et al.  For non-diagonal blocks of s3
599
 
**
600
 
*/
601
 
void s3_block_vrotf(int *Cnt[2], int **Ij[2], int **Ridx[2],
602
 
      signed char **Sn[2], double **C, double **S, 
603
 
      double *tei, int nas, int nbs, int cnas,
604
 
      int Ib_list, int Ja_list, int Jb_list, int Ib_sym, int Jb_sym,
605
 
      double **Cprime, double *F, double *V, double *Sgn, int *L, int *R)
606
 
{
607
 
   int Ia_ex;
608
 
   int ij, i, j, kl, ijkl, I, J, RJ;
609
 
   double tval, VS, *CprimeI0, *CI0;
610
 
   int jlen, Ia_idx, Jacnt, *Iaij, *orbsym, norbs;
611
 
   int *Iaridx;
612
 
   signed char *Iasgn;
613
 
   double *Tptr;
614
 
   
615
 
   norbs = CalcInfo.num_ci_orbs;
616
 
   orbsym = CalcInfo.orbsym + CalcInfo.num_fzc_orbs;
617
 
 
618
 
   /* loop over i, j */
619
 
   for (i=0; i<norbs; i++) {
620
 
     for (j=0; j<=i; j++) {
621
 
       if ((orbsym[i] ^ orbsym[j] ^ Jb_sym ^ Ib_sym) != 0) continue;
622
 
       ij = ioff[i] + j;
623
 
       jlen = form_ilist_rotf(Cnt[1], Ridx[1], Sn[1], Ij[1], 
624
 
          nbs, ij, L, R, Sgn);
625
 
       
626
 
       if (!jlen) continue;
627
 
       
628
 
       Tptr = tei + ioff[ij];
629
 
       
630
 
       /* gather operation */
631
 
       for (I=0; I<cnas; I++) {
632
 
         CprimeI0 = Cprime[I];
633
 
         CI0 = C[I];
634
 
         for (J=0; J<jlen; J++) {
635
 
           tval = Sgn[J];
636
 
           CprimeI0[J] = CI0[L[J]] * tval;
637
 
         }
638
 
       }
639
 
 
640
 
 
641
 
       /* loop over Ia */
642
 
       for (Ia_idx=0; Ia_idx<nas; Ia_idx++) {
643
 
 
644
 
         /* loop over excitations E^a_{kl} from |A(I_a)> */
645
 
         Jacnt = Cnt[0][Ia_idx];
646
 
         Iaridx = Ridx[0][Ia_idx];
647
 
         Iasgn = Sn[0][Ia_idx];
648
 
         Iaij = Ij[0][Ia_idx];
649
 
         
650
 
         zero_arr(V, jlen);
651
 
         
652
 
         for (Ia_ex=0; Ia_ex < Jacnt; Ia_ex++) {
653
 
           kl = *Iaij++;
654
 
           I = *Iaridx++;
655
 
           tval = *Iasgn++;
656
 
           ijkl = INDEX(ij,kl);
657
 
           VS = tval * tei[ijkl];
658
 
           CprimeI0 = Cprime[I];
659
 
           
660
 
         #ifdef USE_BLAS
661
 
           C_DAXPY(jlen, VS, CprimeI0, 1, V, 1);
662
 
         #else
663
 
           for (J=0; J<jlen; J++) {
664
 
             V[J] += VS * CprimeI0[J];
665
 
           }
666
 
         #endif
667
 
           
668
 
         }
669
 
 
670
 
         
671
 
         /* scatter */
672
 
         for (J=0; J<jlen; J++) {
673
 
           RJ = R[J];
674
 
           S[Ia_idx][RJ] += V[J];
675
 
         }
676
 
 
677
 
       } /* end loop over Ia */
678
 
 
679
 
     } /* end loop over j */
680
 
   } /* end loop over i */
681
 
   
682
 
}              
683
 
 
684
 
 
685
 
int form_ilist_rotf(int *Cnt, int **Ridx, signed char **Sn, int **Ij,
686
 
   int nas, int kl, int *L, int *R, double *Sgn)
687
 
{
688
 
 
689
 
   int inum=0, Ia_idx, Ia_ex, Iacnt, ij;
690
 
   int *Iaij;
691
 
   int *Iaridx;
692
 
   signed char *Iasgn;
693
 
 
694
 
   /* loop over Ia */
695
 
   for (Ia_idx=0; Ia_idx < nas; Ia_idx++) {
696
 
 
697
 
      /* loop over excitations E^a_{kl} from |A(I_a)> */
698
 
 
699
 
      Iacnt = Cnt[Ia_idx];
700
 
      if (!Iacnt) continue;
701
 
      Iaridx = Ridx[Ia_idx];
702
 
      Iasgn = Sn[Ia_idx];
703
 
      Iaij = Ij[Ia_idx];
704
 
      Ia_ex=0;
705
 
      for (Ia_ex=0; Ia_ex<Iacnt; Ia_ex++) {
706
 
         ij = *Iaij++;
707
 
         if (ij == kl) {
708
 
            *R++ = Ia_idx;
709
 
            *L++ = Iaridx[Ia_ex];
710
 
            *Sgn++ = (double) Iasgn[Ia_ex];
711
 
            inum++;
712
 
            }
713
 
         }
714
 
      }  /* end loop over Ia */
715
 
 
716
 
   return(inum);
717
 
}
718
 
 
719