← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/detci/s2v.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/s2v.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
 
extern unsigned char ***Occs;
12
 
 
13
 
extern void b2brepl(unsigned char **occs, int *Jcnt, int **Jij, int **Joij, 
14
 
      int **Jridx, signed char **Jsgn, struct olsen_graph *Graph,
15
 
      int Ilist, int Jlist, int len);
16
 
void s2_block_vfci_pthread(void *threadarg);
17
 
void s2_block_vras_pthread(void *threadarg);
18
 
 
19
 
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
20
 
 
21
 
/*
22
 
** S2_BLOCK_VFCI()
23
 
** 
24
 
** Calculate the sigma_2 vector as described by
25
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
26
 
**
27
 
** This sigma1 routine is for Full CI's only, assumes (ij|ij)'s have not
28
 
** been halved, and attempts to follow Olsen's vectorized algorithm more
29
 
** closely than previous versions, using sparsity of F.
30
 
** 
31
 
** David Sherrill, 18 April 1996
32
 
** Based on many previous versions by David Sherrill 1994-5
33
 
*/
34
 
void s2_block_vfci(struct stringwr **alplist, struct stringwr **betlist, 
35
 
      double **C, double **S, double *oei, double *tei, double *F,
36
 
      int nlists, int nas, int nbs, int Ia_list, int Ja_list, 
37
 
      int Ja_list_nas)
38
 
{
39
 
   struct stringwr *Ia, *Ka;
40
 
   unsigned int Ia_idx, Ib_idx, Ka_idx, Ja_idx;
41
 
   unsigned int Iacnt, Kacnt, Ka_list, Ia_ex, Ka_ex;
42
 
   unsigned int *Iaridx, *Karidx;
43
 
   int *Iaij, *Kaij;
44
 
   signed char *Iasgn, *Kasgn;
45
 
   int ij,kl,ijkl;
46
 
   double Ka_sgn, Ja_sgn;
47
 
   double tval;
48
 
   double *Sptr, *Cptr;
49
 
 
50
 
 
51
 
   /* loop over I_a */
52
 
   for (Ia=alplist[Ia_list], Ia_idx=0; Ia_idx < nas; Ia_idx++, Ia++) {
53
 
 
54
 
      Sptr = S[Ia_idx];
55
 
      zero_arr(F, Ja_list_nas);
56
 
 
57
 
      /* loop over excitations E^a_{kl} from |A(I_a)> */
58
 
      for (Ka_list=0; Ka_list < nlists; Ka_list++) {
59
 
            Iacnt = Ia->cnt[Ka_list];
60
 
            Iaridx = Ia->ridx[Ka_list];
61
 
            Iasgn = Ia->sgn[Ka_list];
62
 
            Iaij = Ia->ij[Ka_list];
63
 
            for (Ia_ex=0; Ia_ex < Iacnt; Ia_ex++) {
64
 
               kl = *Iaij++;
65
 
               Ka_idx = *Iaridx++;
66
 
               Ka_sgn = (double) *Iasgn++;
67
 
 
68
 
               /* A(K_a) = sgn(kl) * E^a_{kl} |A(I_a)> */
69
 
               Ka = alplist[Ka_list] + Ka_idx;
70
 
               if (Ka_list == Ja_list) F[Ka_idx] += Ka_sgn * oei[kl];
71
 
 
72
 
               /* loop over excitations E^a_{ij} from |A(K_a)> */
73
 
               /* Ja_list pre-determined because of C blocking */
74
 
               Kacnt = Ka->cnt[Ja_list];
75
 
               Karidx = Ka->ridx[Ja_list];
76
 
               Kasgn = Ka->sgn[Ja_list];
77
 
               Kaij = Ka->ij[Ja_list];
78
 
               for (Ka_ex=0; Ka_ex < Kacnt; Ka_ex++) {
79
 
                  Ja_idx = *Karidx++;
80
 
                  Ja_sgn = (double) *Kasgn++;
81
 
                  ij = *Kaij++;
82
 
                  ijkl = INDEX(ij,kl);
83
 
                  F[Ja_idx] += 0.5 * Ka_sgn * Ja_sgn * tei[ijkl] ;
84
 
                  }
85
 
               } /* end loop over Ia excitations */
86
 
         } /* end loop over Ka_list */
87
 
 
88
 
      
89
 
      /*
90
 
      for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
91
 
         tval = 0.0;
92
 
         for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
93
 
            tval += C[Ja_idx][Ib_idx] * F[Ja_idx];
94
 
            }
95
 
         S[Ia_idx][Ib_idx] += tval;
96
 
         }
97
 
      */
98
 
 
99
 
      for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
100
 
         if ((tval=F[Ja_idx]) == 0.0) continue;
101
 
         Cptr = C[Ja_idx];
102
 
 
103
 
      #ifdef USE_BLAS
104
 
         C_DAXPY(nbs, tval, Cptr, 1, Sptr, 1);
105
 
      #else   
106
 
         for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
107
 
            Sptr[Ib_idx] += tval * Cptr[Ib_idx];
108
 
            }
109
 
      #endif
110
 
         }
111
 
 
112
 
      } /* end loop over Ia */
113
 
 
114
 
}
115
 
 
116
 
 
117
 
/*
118
 
** S2_BLOCK_VFCI_THREAD()
119
 
** 
120
 
** Calculate the sigma_2 vector as described by
121
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
122
 
**
123
 
** This sigma1 routine is for Full CI's only, assumes (ij|ij)'s have not
124
 
** been halved, and attempts to follow Olsen's vectorized algorithm more
125
 
** closely than previous versions, using sparsity of F.
126
 
** 
127
 
** David Sherrill, 18 April 1996
128
 
** Based on many previous versions by David Sherrill 1994-5
129
 
*/
130
 
void s2_block_vfci_thread(struct stringwr **alplist, struct stringwr **betlist, 
131
 
      double **C, double **S, double *oei, double *tei, double *F,
132
 
      int nlists, int nas, int nbs, int Ia_list, int Ja_list, 
133
 
      int Ja_list_nas)
134
 
{
135
 
  struct stringwr *Ia, *Ka;
136
 
  unsigned int Ia_idx, Ib_idx, Ka_idx, Ja_idx;
137
 
  unsigned int Iacnt, Kacnt, Ka_list, Ia_ex, Ka_ex;
138
 
  unsigned int *Iaridx, *Karidx;
139
 
  int *Iaij, *Kaij;
140
 
  signed char *Iasgn, *Kasgn;
141
 
  int ij,kl,ijkl;
142
 
  double Ka_sgn, Ja_sgn;
143
 
  double tval;
144
 
  double *Sptr, *Cptr;
145
 
  struct pthreads_s2vfci **thread_info;
146
 
  int i;
147
 
   
148
 
  thread_info = (struct pthreads_s2vfci **)
149
 
                malloc(sizeof(struct pthreads_s2vfci *) * nas);
150
 
  for (i=0; i<nas; i++) {
151
 
      thread_info[i] = (struct pthreads_s2vfci *)
152
 
                       malloc(sizeof(struct pthreads_s2vfci));
153
 
    }
154
 
 
155
 
  tpool_queue_open(thread_pool);
156
 
 
157
 
  detci_time.s2_mt_before_time = wall_time_new();
158
 
  /* loop over I_a */
159
 
  for (Ia=alplist[Ia_list], Ia_idx=0; Ia_idx < nas; Ia_idx++, Ia++) {
160
 
      thread_info[Ia_idx]->alplist=alplist;
161
 
      thread_info[Ia_idx]->betlist=betlist;
162
 
      thread_info[Ia_idx]->C=C;
163
 
      thread_info[Ia_idx]->S=S;
164
 
      thread_info[Ia_idx]->oei=oei;
165
 
      thread_info[Ia_idx]->tei=tei;
166
 
      thread_info[Ia_idx]->nlists=nlists;
167
 
      thread_info[Ia_idx]->nas=nas;
168
 
      thread_info[Ia_idx]->nbs=nbs;
169
 
      thread_info[Ia_idx]->Ia_list=Ia_list;
170
 
      thread_info[Ia_idx]->Ja_list=Ja_list;
171
 
      thread_info[Ia_idx]->Ja_list_nas=Ja_list_nas;
172
 
      thread_info[Ia_idx]->Ia=Ia;
173
 
      thread_info[Ia_idx]->Ia_idx=Ia_idx;
174
 
      tpool_add_work(thread_pool, s2_block_vfci_pthread, (void *) thread_info[Ia_idx]);
175
 
    } /* end loop over Ia */
176
 
  tpool_queue_close(thread_pool, 1);
177
 
 
178
 
  detci_time.s2_mt_after_time = wall_time_new();
179
 
  detci_time.s2_mt_total_time += detci_time.s2_mt_after_time - detci_time.s2_mt_before_time;
180
 
 
181
 
  for (i=0; i<nas; i++) free(thread_info[i]);
182
 
}
183
 
 
184
 
 
185
 
/*
186
 
** S2_BLOCK_VFCI_PTHREAD()
187
 
** 
188
 
** Calculate the sigma_2 vector as described by
189
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
190
 
**
191
 
** This sigma1 routine is for Full CI's only, assumes (ij|ij)'s have not
192
 
** been halved, and attempts to follow Olsen's vectorized algorithm more
193
 
** closely than previous versions, using sparsity of F.
194
 
** 
195
 
** David Sherrill, 18 April 1996
196
 
** Based on many previous versions by David Sherrill 1994-5
197
 
*/
198
 
void s2_block_vfci_pthread(void *threadarg)
199
 
{
200
 
  struct stringwr *Ia, *Ka, **alplist, **betlist;
201
 
  unsigned int Ia_idx, Ib_idx, Ka_idx, Ja_idx;
202
 
  unsigned int Iacnt, Kacnt, Ka_list, Ia_ex, Ka_ex;
203
 
  unsigned int *Iaridx, *Karidx;
204
 
  int *Iaij, *Kaij;
205
 
  signed char *Iasgn, *Kasgn;
206
 
  int ij,kl,ijkl;
207
 
  double Ka_sgn, Ja_sgn;
208
 
  double tval;
209
 
  double *Sptr, *Cptr, *oei, *tei;
210
 
  double *F, **C, **S;
211
 
  struct pthreads_s2vfci *thread_info;
212
 
  int nlists, nas, nbs, Ia_list, Ja_list, Ja_list_nas;
213
 
  
214
 
  thread_info = (struct pthreads_s2vfci *) threadarg;
215
 
  alplist = thread_info->alplist;
216
 
  betlist = thread_info->betlist;
217
 
  C = thread_info->C;
218
 
  S = thread_info->S;
219
 
  oei = thread_info->oei;
220
 
  tei = thread_info->tei;
221
 
  nlists = thread_info->nlists;
222
 
  nas = thread_info->nas;
223
 
  nbs = thread_info->nbs;
224
 
  Ia_list = thread_info->Ia_list;
225
 
  Ja_list = thread_info->Ja_list;
226
 
  Ja_list_nas = thread_info->Ja_list_nas;
227
 
  Ia = thread_info->Ia;
228
 
  Ia_idx = thread_info->Ia_idx;
229
 
 
230
 
  F = init_array(Ja_list_nas);
231
 
  Sptr = S[Ia_idx];
232
 
  zero_arr(F, Ja_list_nas);
233
 
 
234
 
  /* loop over excitations E^a_{kl} from |A(I_a)> */
235
 
  for (Ka_list=0; Ka_list < nlists; Ka_list++) {
236
 
      Iacnt = Ia->cnt[Ka_list];
237
 
      Iaridx = Ia->ridx[Ka_list];
238
 
      Iasgn = Ia->sgn[Ka_list];
239
 
      Iaij = Ia->ij[Ka_list];
240
 
      for (Ia_ex=0; Ia_ex < Iacnt; Ia_ex++) {
241
 
          kl = *Iaij++;
242
 
          Ka_idx = *Iaridx++;
243
 
          Ka_sgn = (double) *Iasgn++;
244
 
 
245
 
          /* A(K_a) = sgn(kl) * E^a_{kl} |A(I_a)> */
246
 
          Ka = alplist[Ka_list] + Ka_idx;
247
 
          if (Ka_list == Ja_list) F[Ka_idx] += Ka_sgn * oei[kl];
248
 
 
249
 
          /* loop over excitations E^a_{ij} from |A(K_a)> */
250
 
          /* Ja_list pre-determined because of C blocking */
251
 
          Kacnt = Ka->cnt[Ja_list];
252
 
          Karidx = Ka->ridx[Ja_list];
253
 
          Kasgn = Ka->sgn[Ja_list];
254
 
          Kaij = Ka->ij[Ja_list];
255
 
          for (Ka_ex=0; Ka_ex < Kacnt; Ka_ex++) {
256
 
              Ja_idx = *Karidx++;
257
 
              Ja_sgn = (double) *Kasgn++;
258
 
              ij = *Kaij++;
259
 
              ijkl = INDEX(ij,kl);
260
 
              F[Ja_idx] += 0.5 * Ka_sgn * Ja_sgn * tei[ijkl] ;
261
 
            }
262
 
        } /* end loop over Ia excitations */
263
 
    } /* end loop over Ka_list */
264
 
 
265
 
      
266
 
  /*
267
 
    for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
268
 
    tval = 0.0;
269
 
    for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
270
 
    tval += C[Ja_idx][Ib_idx] * F[Ja_idx];
271
 
    }
272
 
    S[Ia_idx][Ib_idx] += tval;
273
 
    }
274
 
  */
275
 
 
276
 
  for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
277
 
      if ((tval=F[Ja_idx]) == 0.0) continue;
278
 
      Cptr = C[Ja_idx];
279
 
 
280
 
#ifdef USE_BLAS
281
 
      C_DAXPY(nbs, tval, Cptr, 1, Sptr, 1);
282
 
#else   
283
 
      for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
284
 
          Sptr[Ib_idx] += tval * Cptr[Ib_idx];
285
 
        }
286
 
#endif
287
 
    }
288
 
  free(F);
289
 
}
290
 
 
291
 
/*
292
 
** S2_BLOCK_VRAS()
293
 
** 
294
 
** Calculate the sigma_2 vector as described by
295
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
296
 
**
297
 
** This sigma2 routine is for RAS CI's.
298
 
** currently assumes that (ij|ij)'s have not been halved!! 
299
 
** 
300
 
** David Sherrill, 10 May 1996
301
 
** Based on previous code by David Sherrill, 1994-5
302
 
**
303
 
** Updated 3/27/94 to include g matrix for RAS
304
 
** Modified 4/8/94 to make C and s one-dimensional
305
 
** Modified 4/10/94 to make FCI-only (for now) and use new string structs
306
 
** Modified 6/21/95 for use in new RAS program
307
 
** Obtained 7/22/95 from s1 routine by changing a's to b's and vice versa
308
 
** Modified 5/10/96 for more vectorized approach
309
 
*/
310
 
void s2_block_vras(struct stringwr **alplist, struct stringwr **betlist, 
311
 
      double **C, double **S, double *oei, double *tei, double *F,
312
 
      int nlists, int nas, int nbs, int Ia_list, int Ja_list, 
313
 
      int Ja_list_nas)
314
 
{
315
 
   struct stringwr *Ia, *Ka;
316
 
   unsigned int Ia_idx, Ib_idx, Ka_idx, Ja_idx;
317
 
   unsigned int Iacnt, Kacnt, Ka_list, Ia_ex, Ka_ex;
318
 
   unsigned int *Iaridx, *Karidx;
319
 
   int nirreps, *Iaij, *Kaij, *Iaoij, *Kaoij;
320
 
   signed char *Iasgn, *Kasgn;
321
 
   int ij,kl,ijkl,oij,okl;
322
 
   double Ka_sgn, Ja_sgn;
323
 
   double tval;
324
 
   double *Sptr, *Cptr;
325
 
 
326
 
   nirreps = CalcInfo.nirreps;
327
 
 
328
 
   /* loop over I_a */
329
 
   for (Ia=alplist[Ia_list], Ia_idx=0; Ia_idx < nas; Ia_idx++, Ia++) {
330
 
 
331
 
      Sptr = S[Ia_idx];
332
 
      zero_arr(F, Ja_list_nas);
333
 
 
334
 
      /* loop over excitations E^a_{kl} from |A(I_a)> */
335
 
      for (Ka_list=0; Ka_list < nlists; Ka_list++) {
336
 
            Iacnt = Ia->cnt[Ka_list];
337
 
            Iaridx = Ia->ridx[Ka_list];
338
 
            Iasgn = Ia->sgn[Ka_list];
339
 
            Iaij = Ia->ij[Ka_list];
340
 
            Iaoij = Ia->oij[Ka_list];
341
 
            for (Ia_ex=0; Ia_ex < Iacnt; Ia_ex++) {
342
 
               kl = *Iaij++;
343
 
               okl = *Iaoij++;
344
 
               Ka_idx = *Iaridx++;
345
 
               Ka_sgn = (double) *Iasgn++;
346
 
 
347
 
               /* A(K_a) = sgn(kl) * E^a_{kl} |A(I_a)> */
348
 
               Ka = alplist[Ka_list] + Ka_idx;
349
 
               /* note okl on next line, not kl */
350
 
               if (Ka_list == Ja_list) F[Ka_idx] += Ka_sgn * oei[okl];
351
 
 
352
 
               /* loop over excitations E^a_{ij} from |A(K_a)> */
353
 
               /* Ja_list pre-determined because of C blocking */
354
 
               Kacnt = Ka->cnt[Ja_list];
355
 
               Karidx = Ka->ridx[Ja_list];
356
 
               Kasgn = Ka->sgn[Ja_list];
357
 
               Kaij = Ka->ij[Ja_list];
358
 
               Kaoij = Ka->oij[Ja_list];
359
 
               for (Ka_ex=0; Ka_ex < Kacnt; Ka_ex++) {
360
 
                  Ja_idx = *Karidx++;
361
 
                  Ja_sgn = (double) *Kasgn++;
362
 
                  ij = *Kaij++;
363
 
                  oij = *Kaoij++;
364
 
                  ijkl = INDEX(ij,kl);
365
 
                  if (oij > okl) 
366
 
                     F[Ja_idx] += Ka_sgn * Ja_sgn * tei[ijkl] ;
367
 
                  else if (oij == okl) 
368
 
                     F[Ja_idx] += 0.5 * Ka_sgn * Ja_sgn * tei[ijkl] ;
369
 
                  }
370
 
               } /* end loop over Ia excitations */
371
 
         } /* end loop over Ka_list */
372
 
      
373
 
      /*
374
 
      for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
375
 
         tval = 0.0;
376
 
         for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
377
 
            tval += C[Ja_idx][Ib_idx] * F[Ja_idx];
378
 
            }
379
 
         S[Ia_idx][Ib_idx] += tval;
380
 
         }
381
 
      */
382
 
 
383
 
      for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
384
 
         if ((tval=F[Ja_idx]) == 0.0) continue;
385
 
         Cptr = C[Ja_idx];
386
 
      #ifdef USE_BLAS
387
 
         C_DAXPY(nbs, tval, Cptr, 1, Sptr, 1);
388
 
      #else
389
 
         for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
390
 
            Sptr[Ib_idx] += tval * Cptr[Ib_idx];
391
 
            }
392
 
      #endif
393
 
         }
394
 
 
395
 
      } /* end loop over Ia */
396
 
 
397
 
}
398
 
 
399
 
 
400
 
/*
401
 
** S2_BLOCK_VRAS_THREAD()
402
 
** 
403
 
** Calculate the sigma_2 vector as described by
404
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
405
 
**
406
 
** This sigma2 routine is for RAS CI's.
407
 
** currently assumes that (ij|ij)'s have not been halved!! 
408
 
** 
409
 
** David Sherrill, 10 May 1996
410
 
** Based on previous code by David Sherrill, 1994-5
411
 
**
412
 
** Updated 3/27/94 to include g matrix for RAS
413
 
** Modified 4/8/94 to make C and s one-dimensional
414
 
** Modified 4/10/94 to make FCI-only (for now) and use new string structs
415
 
** Modified 6/21/95 for use in new RAS program
416
 
** Obtained 7/22/95 from s1 routine by changing a's to b's and vice versa
417
 
** Modified 5/10/96 for more vectorized approach
418
 
*/
419
 
void s2_block_vras_thread(struct stringwr **alplist, struct stringwr **betlist, 
420
 
      double **C, double **S, double *oei, double *tei, double *F,
421
 
      int nlists, int nas, int nbs, int Ia_list, int Ja_list, 
422
 
      int Ja_list_nas)
423
 
{
424
 
  struct stringwr *Ia;
425
 
  unsigned int Ia_idx;
426
 
  struct pthreads_s2vfci **thread_info;
427
 
  int i;
428
 
 
429
 
  thread_info = (struct pthreads_s2vfci **)
430
 
                malloc(sizeof(struct pthreads_s2vfci *) * nas);
431
 
  for (i=0; i<nas; i++) {
432
 
      thread_info[i] = (struct pthreads_s2vfci *)
433
 
                       malloc(sizeof(struct pthreads_s2vfci));
434
 
    }
435
 
 
436
 
  tpool_queue_open(thread_pool);
437
 
  detci_time.s2_mt_before_time = wall_time_new();
438
 
 
439
 
 
440
 
  /* loop over I_a */
441
 
  for (Ia=alplist[Ia_list], Ia_idx=0; Ia_idx < nas; Ia_idx++, Ia++) {
442
 
      thread_info[Ia_idx]->alplist=alplist;
443
 
      thread_info[Ia_idx]->betlist=betlist;
444
 
      thread_info[Ia_idx]->C=C;
445
 
      thread_info[Ia_idx]->S=S;
446
 
      thread_info[Ia_idx]->oei=oei;
447
 
      thread_info[Ia_idx]->tei=tei;
448
 
      thread_info[Ia_idx]->nlists=nlists;
449
 
      thread_info[Ia_idx]->nas=nas;
450
 
      thread_info[Ia_idx]->nbs=nbs;
451
 
      thread_info[Ia_idx]->Ia_list=Ia_list;
452
 
      thread_info[Ia_idx]->Ja_list=Ja_list;
453
 
      thread_info[Ia_idx]->Ja_list_nas=Ja_list_nas;
454
 
      thread_info[Ia_idx]->Ia=Ia;
455
 
      thread_info[Ia_idx]->Ia_idx=Ia_idx;
456
 
      tpool_add_work(thread_pool, s2_block_vras_pthread, (void *) thread_info[Ia_idx]);
457
 
    } /* end loop over Ia */
458
 
 
459
 
  tpool_queue_close(thread_pool, 1);
460
 
  detci_time.s2_mt_after_time = wall_time_new();
461
 
  detci_time.s2_mt_total_time += detci_time.s2_mt_after_time - detci_time.s2_mt_before_time;
462
 
 
463
 
 
464
 
  for (i=0; i<nas; i++) free(thread_info[i]);
465
 
   
466
 
}
467
 
 
468
 
 
469
 
/*
470
 
** S2_BLOCK_VRAS_PTHREAD()
471
 
** 
472
 
** Calculate the sigma_2 vector as described by
473
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
474
 
**
475
 
** This sigma2 routine is for RAS CI's.
476
 
** currently assumes that (ij|ij)'s have not been halved!! 
477
 
** 
478
 
** David Sherrill, 10 May 1996
479
 
** Based on previous code by David Sherrill, 1994-5
480
 
**
481
 
** Updated 3/27/94 to include g matrix for RAS
482
 
** Modified 4/8/94 to make C and s one-dimensional
483
 
** Modified 4/10/94 to make FCI-only (for now) and use new string structs
484
 
** Modified 6/21/95 for use in new RAS program
485
 
** Obtained 7/22/95 from s1 routine by changing a's to b's and vice versa
486
 
** Modified 5/10/96 for more vectorized approach
487
 
*/
488
 
void s2_block_vras_pthread(void *threadarg)
489
 
{
490
 
  struct stringwr *Ia, *Ka, **alplist, **betlist;
491
 
  unsigned int Ia_idx, Ib_idx, Ka_idx, Ja_idx;
492
 
  unsigned int Iacnt, Kacnt, Ka_list, Ia_ex, Ka_ex;
493
 
  unsigned int *Iaridx, *Karidx;
494
 
  int nirreps, *Iaij, *Kaij, *Iaoij, *Kaoij;
495
 
  signed char *Iasgn, *Kasgn;
496
 
  int ij,kl,ijkl,oij,okl;
497
 
  double Ka_sgn, Ja_sgn;
498
 
  double tval;
499
 
  double *Sptr, *Cptr, *oei, *tei, **C, **S, *F;
500
 
  struct pthreads_s2vfci *thread_info;
501
 
  int nlists, nas, nbs, Ia_list, Ja_list, Ja_list_nas;
502
 
 
503
 
  thread_info = (struct pthreads_s2vfci *) threadarg;
504
 
  alplist = thread_info->alplist;
505
 
  betlist = thread_info->betlist;
506
 
  C = thread_info->C;
507
 
  S = thread_info->S;
508
 
  oei = thread_info->oei;
509
 
  tei = thread_info->tei;
510
 
  nlists = thread_info->nlists;
511
 
  nas = thread_info->nas;
512
 
  nbs = thread_info->nbs;
513
 
  Ia_list = thread_info->Ia_list;
514
 
  Ja_list = thread_info->Ja_list;
515
 
  Ja_list_nas = thread_info->Ja_list_nas;
516
 
  Ia = thread_info->Ia;
517
 
  Ia_idx = thread_info->Ia_idx;
518
 
 
519
 
  F = init_array(Ja_list_nas);
520
 
  Sptr = S[Ia_idx];
521
 
  zero_arr(F, Ja_list_nas);
522
 
 
523
 
  nirreps = CalcInfo.nirreps;
524
 
 
525
 
  /* loop over excitations E^a_{kl} from |A(I_a)> */
526
 
  for (Ka_list=0; Ka_list < nlists; Ka_list++) {
527
 
      Iacnt = Ia->cnt[Ka_list];
528
 
      Iaridx = Ia->ridx[Ka_list];
529
 
      Iasgn = Ia->sgn[Ka_list];
530
 
      Iaij = Ia->ij[Ka_list];
531
 
      Iaoij = Ia->oij[Ka_list];
532
 
      for (Ia_ex=0; Ia_ex < Iacnt; Ia_ex++) {
533
 
          kl = *Iaij++;
534
 
          okl = *Iaoij++;
535
 
          Ka_idx = *Iaridx++;
536
 
          Ka_sgn = (double) *Iasgn++;
537
 
 
538
 
          /* A(K_a) = sgn(kl) * E^a_{kl} |A(I_a)> */
539
 
          Ka = alplist[Ka_list] + Ka_idx;
540
 
          /* note okl on next line, not kl */
541
 
          if (Ka_list == Ja_list) F[Ka_idx] += Ka_sgn * oei[okl];
542
 
 
543
 
          /* loop over excitations E^a_{ij} from |A(K_a)> */
544
 
          /* Ja_list pre-determined because of C blocking */
545
 
          Kacnt = Ka->cnt[Ja_list];
546
 
          Karidx = Ka->ridx[Ja_list];
547
 
          Kasgn = Ka->sgn[Ja_list];
548
 
          Kaij = Ka->ij[Ja_list];
549
 
          Kaoij = Ka->oij[Ja_list];
550
 
          for (Ka_ex=0; Ka_ex < Kacnt; Ka_ex++) {
551
 
              Ja_idx = *Karidx++;
552
 
              Ja_sgn = (double) *Kasgn++;
553
 
              ij = *Kaij++;
554
 
              oij = *Kaoij++;
555
 
              ijkl = INDEX(ij,kl);
556
 
              if (oij > okl) 
557
 
                  F[Ja_idx] += Ka_sgn * Ja_sgn * tei[ijkl] ;
558
 
              else if (oij == okl) 
559
 
                  F[Ja_idx] += 0.5 * Ka_sgn * Ja_sgn * tei[ijkl] ;
560
 
            }
561
 
        } /* end loop over Ia excitations */
562
 
    } /* end loop over Ka_list */
563
 
      
564
 
  /*
565
 
    for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
566
 
    tval = 0.0;
567
 
    for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
568
 
    tval += C[Ja_idx][Ib_idx] * F[Ja_idx];
569
 
    }
570
 
    S[Ia_idx][Ib_idx] += tval;
571
 
    }
572
 
  */
573
 
 
574
 
  for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
575
 
      if ((tval=F[Ja_idx]) == 0.0) continue;
576
 
      Cptr = C[Ja_idx];
577
 
#ifdef USE_BLAS
578
 
      C_DAXPY(nbs, tval, Cptr, 1, Sptr, 1);
579
 
#else
580
 
      for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
581
 
          Sptr[Ib_idx] += tval * Cptr[Ib_idx];
582
 
        }
583
 
#endif
584
 
    }
585
 
  free(F);
586
 
}
587
 
 
588
 
 
589
 
/*
590
 
** S2_BLOCK_VRAS_ROTF()
591
 
** 
592
 
** s2_block_vras_rotf(): Calculate the sigma_2 vector as described by
593
 
** equation (20) of RAS Paper (Olsen, Roos, Jorgensen, Aa. Jensen JCP 1988)
594
 
**
595
 
** String replacements on-the-fly version
596
 
** currently assumes that (ij|ij)'s have not been halved!! 
597
 
** 
598
 
** This sigma2 routine is for RAS CI's.
599
 
** 
600
 
** David Sherrill, 13 May 1996
601
 
** Based on previous code by David Sherrill, 1994-5
602
 
**
603
 
** Updated 3/27/94 to include g matrix for RAS
604
 
** Modified 4/8/94 to make C and s one-dimensional
605
 
** Modified 4/10/94 to make FCI-only (for now) and use new string structs
606
 
** Modified 6/21/95 for use in new RAS program
607
 
** Obtained 7/22/95 from s1 routine by changing a's to b's and vice versa
608
 
** Modified 5/13/96 for new sparse-F vectorized version
609
 
**
610
 
*/
611
 
void s2_block_vras_rotf(int *Cnt[2], int **Ij[2], int **Oij[2],
612
 
      int **Ridx[2], signed char **Sgn[2], unsigned char **Toccs,
613
 
      double **C, double **S,
614
 
      double *oei, double *tei, double *F, int nlists, int nas, int nbs,
615
 
      int Ia_list, int Ja_list, int Ja_list_nas)
616
 
{
617
 
   int Ia_idx, Ib_idx, Ka_idx, Ja_idx;
618
 
   int Iacnt, Kacnt, Ka_list, Ia_ex, Ka_ex;
619
 
   int *Iaridx, *Karidx;
620
 
   int nirreps, *Iaij, *Kaij, *Iaoij, *Kaoij;
621
 
   signed char *Iasgn, *Kasgn;
622
 
   int i,ij,kl,ijkl,oij,okl;
623
 
   double Ka_sgn, Ja_sgn;
624
 
   double tval, *Cptr, *Sptr;
625
 
 
626
 
   nirreps = CalcInfo.nirreps;
627
 
 
628
 
   for (Ka_list=0; Ka_list < nlists; Ka_list++) {
629
 
      b2brepl(Occs[Ia_list], Cnt[0], Ij[0], Oij[0], Ridx[0],
630
 
         Sgn[0], BetaG, Ia_list, Ka_list, nas);
631
 
 
632
 
      /* loop over I_a */
633
 
      for (Ia_idx=0; Ia_idx < nas; Ia_idx++) {
634
 
 
635
 
         if ((Iacnt = Cnt[0][Ia_idx]) < 0) continue;
636
 
         Sptr = S[Ia_idx];
637
 
         zero_arr(F, Ja_list_nas);
638
 
 
639
 
         /* loop over excitations E^a_{kl} from |A(I_a)> */
640
 
         Iaridx = Ridx[0][Ia_idx];
641
 
         Iasgn = Sgn[0][Ia_idx];
642
 
         Iaij = Ij[0][Ia_idx];
643
 
         Iaoij = Oij[0][Ia_idx];
644
 
 
645
 
         for (i=0; i<Iacnt; i++) 
646
 
            Toccs[i] = Occs[Ka_list][Iaridx[i]];
647
 
 
648
 
         b2brepl(Toccs, Cnt[1], Ij[1], Oij[1], Ridx[1], Sgn[1],
649
 
            AlphaG, Ka_list, Ja_list, Iacnt);
650
 
 
651
 
         for (Ia_ex=0; Ia_ex < Iacnt; Ia_ex++) {
652
 
            kl = *Iaij++;
653
 
            okl = *Iaoij++;
654
 
            Ka_idx = *Iaridx++;
655
 
            Ka_sgn = (double) *Iasgn++;
656
 
 
657
 
            /* A(K_a) = sgn(kl) * E^a_{kl} |A(I_a)> */
658
 
            /* note okl on next line, not kl */
659
 
            if (Ka_list == Ja_list) F[Ka_idx] += Ka_sgn * oei[okl];
660
 
 
661
 
            /* loop over excitations E^a_{ij} from |A(K_a)> */
662
 
            /* Ja_list pre-determined because of C blocking */
663
 
            Kacnt = Cnt[1][Ia_ex];
664
 
            Karidx = Ridx[1][Ia_ex];
665
 
            Kasgn = Sgn[1][Ia_ex];
666
 
            Kaij = Ij[1][Ia_ex];
667
 
            Kaoij = Oij[1][Ia_ex];
668
 
            for (Ka_ex=0; Ka_ex < Kacnt; Ka_ex++) {
669
 
               Ja_idx = *Karidx++;
670
 
               Ja_sgn = (double) *Kasgn++;
671
 
               ij = *Kaij++;
672
 
               oij = *Kaoij++;
673
 
               ijkl = INDEX(ij,kl);
674
 
               if (oij > okl) 
675
 
                  F[Ja_idx] += Ka_sgn * Ja_sgn * tei[ijkl] ;
676
 
               else if (oij == okl) 
677
 
                  F[Ja_idx] += 0.5 * Ka_sgn * Ja_sgn * tei[ijkl] ;
678
 
               }
679
 
            } /* end loop over Ia excitations */
680
 
 
681
 
      /*
682
 
      for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
683
 
         tval = 0.0;
684
 
         for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
685
 
            tval += C[Ja_idx][Ib_idx] * F[Ja_idx];
686
 
            }
687
 
         S[Ia_idx][Ib_idx] += tval;
688
 
         }
689
 
      */
690
 
 
691
 
      for (Ja_idx=0; Ja_idx < Ja_list_nas; Ja_idx++) {
692
 
         if ((tval=F[Ja_idx]) == 0.0) continue;
693
 
         Cptr = C[Ja_idx];
694
 
         for (Ib_idx=0; Ib_idx < nbs; Ib_idx++) {
695
 
            Sptr[Ib_idx] += tval * Cptr[Ib_idx];
696
 
            }
697
 
         }
698
 
 
699
 
      } /* end loop over Ia */
700
 
   } /* end loop over Ka_list */
701
 
 
702
 
}
703