← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/psicode/intrepid

« back to all changes in this revision

Viewing changes to src/lib/libr12/r_vrr_build.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libr12/r_vrr_build.c
 
1
#include <math.h>
 
2
#include <stdio.h>
 
3
#include <stdlib.h>
 
4
#include <libint/libint.h>
 
5
#include "libr12.h"
 
6
 
 
7
extern void punt(char *);
 
8
static int hash(int a[2][3], int b[2]);
 
9
 
 
10
REALTYPE *r_vrr_build_xxxx(int am_in[2], prim_data *Data, REALTYPE *vp, const REALTYPE *i0, const REALTYPE *i1, REALTYPE *i2,
 
11
    const REALTYPE *i3, const REALTYPE *i4, const REALTYPE *i5)
 
12
{
 
13
  int i, j, k, l;
 
14
  int a;
 
15
  int flag = 0;
 
16
  int am[2][3];
 
17
  int t1, t2, t3, t4, t2max;
 
18
  int xyz;
 
19
  int la, lc;
 
20
  REALTYPE PA[3], U1[3], loo4zn, loo2z, loo2p, r12int;
 
21
  static int io[] = {0,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153};
 
22
 
 
23
  la = am_in[0];
 
24
  lc = am_in[1];
 
25
  loo4zn = Data->oo2z*Data->oo2n;
 
26
  loo2p = Data->oo2p;
 
27
  if (la == 0) { /*--- Decrement on C ---*/
 
28
    a = 1;
 
29
    t2max = io[lc];
 
30
    PA[0] = Data->U[2][0];
 
31
    PA[1] = Data->U[2][1];
 
32
    PA[2] = Data->U[2][2];
 
33
    U1[0] = Data->U[2][0]*Data->oo2z + Data->U[3][0]*Data->oo2n;
 
34
    U1[1] = Data->U[2][1]*Data->oo2z + Data->U[3][1]*Data->oo2n;
 
35
    U1[2] = Data->U[2][2]*Data->oo2z + Data->U[3][2]*Data->oo2n;
 
36
    loo2z = Data->oo2n;
 
37
  }
 
38
  else {         /*--- Decrement on A ---*/
 
39
    a = 0;
 
40
    t2max = io[la]*io[lc+1];
 
41
    PA[0] = Data->U[0][0];
 
42
    PA[1] = Data->U[0][1];
 
43
    PA[2] = Data->U[0][2];
 
44
    U1[0] = Data->U[0][0]*Data->oo2n + Data->U[1][0]*Data->oo2z;
 
45
    U1[1] = Data->U[0][1]*Data->oo2n + Data->U[1][1]*Data->oo2z;
 
46
    U1[2] = Data->U[0][2]*Data->oo2n + Data->U[1][2]*Data->oo2z;
 
47
    loo2z = Data->oo2z;
 
48
  }
 
49
 
 
50
  t2 = t3 = t4 = 0;
 
51
 
 
52
  for(i = 0; i <= la; i++){
 
53
    am[0][0] = la - i;
 
54
    for(j = 0; j <= i; j++){
 
55
      am[0][1] = i - j;
 
56
      am[0][2] = j;
 
57
 
 
58
      for(k = 0; k <= lc; k++){
 
59
        am[1][0] = lc - k;
 
60
        for(l = 0; l <= k; l++){
 
61
          am[1][1] = k - l;
 
62
          am[1][2] = l;
 
63
          
 
64
          if(am[a][2]) xyz = 2;
 
65
          if(am[a][1]) xyz = 1;
 
66
          if(am[a][0]) xyz = 0;
 
67
          
 
68
          if (t2 == t2max) {
 
69
            /*--- reset indices (read Justin Fermann's thesis, pp 36-41 ---*/
 
70
            /*--- (a-1,0|c0) ---*/
 
71
            am[a][xyz] = am[a][xyz] - 1;
 
72
            am_in[a] = am_in[a] - 1;
 
73
            t2 = hash(am,am_in);
 
74
            
 
75
            /*--- (a-2,0|c0) ---*/
 
76
            if (am_in[a]) {
 
77
              am[a][xyz] = am[a][xyz] - 1;
 
78
              am_in[a] = am_in[a] - 1;
 
79
              t3 = hash(am,am_in);
 
80
              am[a][xyz] = am[a][xyz] + 1;
 
81
              am_in[a] = am_in[a] + 1;
 
82
            }
 
83
            
 
84
            /*--- (a-1,0|c-1,0) ---*/
 
85
            if (am_in[a^1]) {
 
86
              am[a^1][0] = am[a^1][0] - 1;
 
87
              am_in[a^1] = am_in[a^1] - 1;
 
88
              t4 = hash(am,am_in);
 
89
              am[a^1][0] = am[a^1][0] + 1;
 
90
              am_in[a^1] = am_in[a^1] + 1;
 
91
            }
 
92
 
 
93
            am[a][xyz] = am[a][xyz] + 1;
 
94
            am_in[a] = am_in[a] + 1;
 
95
          }
 
96
          r12int = PA[xyz]*i0[t2] - U1[xyz]*i3[t2] + loo2p*(*i2);
 
97
          t2++;
 
98
          if(am[a][xyz] > 1){
 
99
            r12int += (am[a][xyz]-1)*(loo2z*i1[t3] - loo4zn*i4[t3]);
 
100
            t3++;
 
101
          }
 
102
          if(am[a^1][xyz] > 0){
 
103
            r12int -= am[a^1][xyz]*loo4zn*i5[t4];
 
104
            t4++;
 
105
          }
 
106
          *vp = r12int;
 
107
          vp++;
 
108
          i2++;
 
109
        }
 
110
      }
 
111
    }
 
112
  }
 
113
  
 
114
  return vp;
 
115
}
 
116
 
 
117
 
 
118
int hash(a, b)
 
119
  int a[2][3];
 
120
  int b[2];
 
121
{
 
122
  int c[2] = {0,0};
 
123
  int i;
 
124
  static int io[] = {0,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153};
 
125
 
 
126
  if(b[0]){
 
127
    i=b[0]-a[0][0];
 
128
    c[0]=i+io[i]-a[0][1];
 
129
    }
 
130
  if(b[1]){
 
131
    i=b[1]-a[1][0];
 
132
    c[1]=i+io[i]-a[1][1];
 
133
    }
 
134
 
 
135
  return c[0]*io[b[1]+1]+c[1];
 
136
}