← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/lib/libr12/emit_hrr_t_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/emit_hrr_t_build.c
 
1
#include <math.h>
 
2
#include <stdio.h>
 
3
#include "build_libr12.h"
 
4
 
 
5
extern FILE *outfile, *hrr_header;
 
6
extern Libr12Params_t Params;
 
7
 
 
8
extern void punt(char *);
 
9
static int hash(int a[2][3], int b[2]);
 
10
 
 
11
int emit_hrr_t_build()
 
12
{
 
13
  int new_am = Params.new_am;
 
14
  int max_class_size = Params.max_class_size;
 
15
 
 
16
  FILE *code;
 
17
  int p,q,r,s;
 
18
  int ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz;
 
19
  int t0, t1, t2, t3, t4;
 
20
  int i,j,nj,i_i0,i_i1;
 
21
  int k,l,nl,k_i0,k_i1;
 
22
  int cp1dm1_num,cdm1_num;
 
23
  int a, b;
 
24
  int flag;
 
25
  int am_in[2];
 
26
  int am[2][3];
 
27
  int xyz;
 
28
  int class_size;
 
29
  int split;
 
30
  int la, lb;
 
31
  int ld, lc, ld_max;
 
32
  int curr_count,curr_subfunction;
 
33
  int num_subfunctions, subbatch_length;
 
34
  int f;
 
35
  static int io[] = {1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153};
 
36
  static const char am_letter[] = "0pdfghiklmnoqrtuvwxyz";
 
37
  char code_name[21];
 
38
  char function_name[18];
 
39
  char **subfunction_name;
 
40
  
 
41
  for(lc=0;lc<=new_am;lc++) {
 
42
    ld_max = (lc+1)/2;
 
43
    for(ld=1;ld<=ld_max;ld++) {
 
44
 
 
45
      /*-----------------------
 
46
        HRR on centers C and D
 
47
       -----------------------*/
 
48
 
 
49
      am_in[0] = lc-ld;
 
50
      am_in[1] = ld;
 
51
 
 
52
      class_size = ((am_in[0]+1)*(am_in[0]+2)*(am_in[1]+1)*(am_in[1]+2))/4;
 
53
      nl = (am_in[1]*(am_in[1]+1))/2;
 
54
      cp1dm1_num = (am_in[0]+2)*(am_in[0]+3)*nl/2;
 
55
      cdm1_num = (am_in[0]+1)*(am_in[0]+2)*nl/2;
 
56
 
 
57
      /* Decide if the routine has to be split into several routines producing "subbatches" */
 
58
      if (class_size > max_class_size) {
 
59
        split = 1;
 
60
        num_subfunctions = ceil((double)class_size/max_class_size);
 
61
        subbatch_length = 1 + class_size/num_subfunctions;
 
62
      }
 
63
      else {
 
64
        split = 0;
 
65
      }
 
66
      
 
67
      sprintf(function_name,"t2hrr3_build_%c%c",am_letter[am_in[0]],am_letter[am_in[1]]);
 
68
      if (split) {
 
69
        subfunction_name = (char **) malloc (num_subfunctions*sizeof(char *));
 
70
        for(i=0;i<num_subfunctions;i++) {
 
71
          subfunction_name[i] = (char *) malloc(22*sizeof(char));
 
72
          sprintf(subfunction_name[i],"_%s_%d",
 
73
                  function_name,i);
 
74
        }
 
75
      }
 
76
      sprintf(code_name,"%s.cc",function_name);
 
77
      code = fopen(code_name,"w");
 
78
 
 
79
      /* include the function into the hrr_header.h */
 
80
      fprintf(hrr_header,"void %s(REALTYPE *, REALTYPE *, REALTYPE *, const REALTYPE *, const REALTYPE *, ",function_name);
 
81
      fprintf(hrr_header,"const REALTYPE *, const REALTYPE *, const REALTYPE *, int, int);\n");
 
82
 
 
83
      fprintf(code,"  /* This machine-generated function computes a quartet of [r12,T2]|%c%c) integrals */\n\n",
 
84
              am_letter[am_in[0]],am_letter[am_in[1]]);
 
85
      fprintf(code,"#include<libint/libint.h>\n\n");
 
86
      if (split) {
 
87
        for(i=0;i<num_subfunctions;i++) {
 
88
          fprintf(code,"REALTYPE *%s(REALTYPE *, REALTYPE *, REALTYPE *, const REALTYPE *, const REALTYPE *, ",
 
89
                  subfunction_name[i]);
 
90
          fprintf(code,"const REALTYPE *, const REALTYPE *, const REALTYPE *, int, int);\n");
 
91
        }
 
92
        fprintf(code,"\n");
 
93
      }
 
94
      fprintf(code,"void %s(REALTYPE *CD, REALTYPE *AC, REALTYPE *vp, const REALTYPE *I0, const REALTYPE *I1, ",function_name);
 
95
      fprintf(code,"const REALTYPE *I2, const REALTYPE *I3, const REALTYPE *I4, int la, int lb)\n{\n");
 
96
      if (split == 1) {
 
97
        curr_subfunction = 0;
 
98
        curr_count = 0;
 
99
      }
 
100
      else {
 
101
        fprintf(code,"  const REALTYPE CD0 = CD[0];\n");
 
102
        fprintf(code,"  const REALTYPE CD1 = CD[1];\n");
 
103
        fprintf(code,"  const REALTYPE CD2 = CD[2];\n");
 
104
        fprintf(code,"  const REALTYPE AC0 = AC[0];\n");
 
105
        fprintf(code,"  const REALTYPE AC1 = AC[1];\n");
 
106
        fprintf(code,"  const REALTYPE AC2 = AC[2];\n");
 
107
        fprintf(code,"  static int io[] = { 1");
 
108
        for(i=1;i<=new_am+1;i++)
 
109
          fprintf(code,", %d",(i+1)*(i+2)/2);
 
110
        fprintf(code,"};\n");
 
111
        fprintf(code,"  int bcdm1_num = io[lb]*%d;\n\n",cdm1_num);
 
112
      }
 
113
      fprintf(code,"  int pa, qa, b;\n");
 
114
 
 
115
      fprintf(code,"  for(pa=0;pa<=la;pa++)\n");
 
116
      fprintf(code,"    for(qa=0;qa<=pa;qa++)\n");
 
117
      fprintf(code,"      for(b=0;b<((lb+1)*(lb+2)/2);b++) {\n");
 
118
 
 
119
      if (split == 1) {
 
120
        for(f=0;f<num_subfunctions;f++)
 
121
          fprintf(code,"    vp = %s(CD, AC, vp, I0, I1, I2, I3, I4, pa, lb);\n",
 
122
                subfunction_name[f]);
 
123
        fprintf(code,"        I0 += %d;\n",cp1dm1_num);
 
124
        fprintf(code,"        I1 += %d;\n",cdm1_num);
 
125
        fprintf(code,"        I2 += %d;\n",cp1dm1_num);
 
126
        fprintf(code,"        I3 += %d;\n",cdm1_num);
 
127
        fprintf(code,"        I4 += %d;\n",cdm1_num);
 
128
        fprintf(code,"  }\n}\n\n");
 
129
 
 
130
        fprintf(code,"REALTYPE *%s(REALTYPE *CD, REALTYPE *AC, REALTYPE *vp, const REALTYPE *I0, const REALTYPE *I1, ",
 
131
              subfunction_name[0]);
 
132
        fprintf(code,"const REALTYPE *I2, const REALTYPE *I3, const REALTYPE *I4, int pa, int lb)\n{\n");
 
133
        fprintf(code,"  const REALTYPE CD0 = CD[0];\n");
 
134
        fprintf(code,"  const REALTYPE CD1 = CD[1];\n");
 
135
        fprintf(code,"  const REALTYPE CD2 = CD[2];\n");
 
136
        fprintf(code,"  const REALTYPE AC0 = AC[0];\n");
 
137
        fprintf(code,"  const REALTYPE AC1 = AC[1];\n");
 
138
        fprintf(code,"  const REALTYPE AC2 = AC[2];\n");
 
139
        fprintf(code,"  static int io[] = { 1");
 
140
        for(i=1;i<=new_am+1;i++)
 
141
          fprintf(code,", %d",(i+1)*(i+2)/2);
 
142
        fprintf(code,"};\n");
 
143
        fprintf(code,"  int bcdm1_num = io[lb]*%d;\n\n",cdm1_num);
 
144
      }
 
145
      
 
146
      for(p = 0; p <= am_in[0]; p++){
 
147
        am[0][0] = am_in[0] - p;
 
148
        for(q = 0; q <= p; q++){
 
149
          am[0][1] = p - q;
 
150
          am[0][2] = q;
 
151
          
 
152
          for(r = 0; r <= am_in[1]; r++){
 
153
            am[1][0] = am_in[1] - r;
 
154
            for(s = 0; s <= r; s++){
 
155
              am[1][1] = r - s;
 
156
              am[1][2] = s;
 
157
 
 
158
              if (am[1][0]) /* build along x */
 
159
                xyz = 0;
 
160
              else if (am[1][1]) /* build along y */
 
161
                xyz = 1;
 
162
              else /*build along z */
 
163
                xyz = 2;
 
164
 
 
165
              am[0][xyz] += 1;
 
166
              am_in[0] += 1;
 
167
              am[1][xyz] -= 1;
 
168
              am_in[1] -= 1;
 
169
              t0 = hash(am,am_in);
 
170
              am[0][xyz] -= 1;
 
171
              am_in[0] -= 1;
 
172
              t1 = hash(am,am_in);
 
173
              am[1][xyz] += 1;
 
174
              am_in[1] += 1;
 
175
              t2 = t0;
 
176
              t3 = t1;
 
177
              t4 = t1;
 
178
              
 
179
              fprintf(code,"        *(vp++) = I0[%d] + CD%d*I1[%d] + I2[%d] - AC%d*I4[%d] - ",t0,xyz,t1,t2,xyz,t4);
 
180
              if (xyz == 0)
 
181
                fprintf(code,"I3[%d];\n",t3);
 
182
              else
 
183
                fprintf(code,"I3[(pa+%d)*bcdm1_num+%d];\n",xyz,t3);
 
184
 
 
185
              curr_count++;
 
186
              if (curr_count == subbatch_length && split == 1) {
 
187
                curr_count = 0;
 
188
                curr_subfunction++;
 
189
                fprintf(code,"  return vp;\n}\n\n");
 
190
                fprintf(code,"REALTYPE *%s(REALTYPE *CD, REALTYPE *AC, REALTYPE *vp, const REALTYPE *I0, const REALTYPE *I1, ",
 
191
                        subfunction_name[curr_subfunction]);
 
192
                fprintf(code,"const REALTYPE *I2, const REALTYPE *I3, const REALTYPE *I4, int pa, int lb)\n{\n");
 
193
                fprintf(code,"  const REALTYPE CD0 = CD[0];\n");
 
194
                fprintf(code,"  const REALTYPE CD1 = CD[1];\n");
 
195
                fprintf(code,"  const REALTYPE CD2 = CD[2];\n");
 
196
                fprintf(code,"  const REALTYPE AC0 = AC[0];\n");
 
197
                fprintf(code,"  const REALTYPE AC1 = AC[1];\n");
 
198
                fprintf(code,"  const REALTYPE AC2 = AC[2];\n");
 
199
                fprintf(code,"  static int io[] = { 1");
 
200
                for(i=1;i<=new_am+1;i++)
 
201
                  fprintf(code,", %d",(i+1)*(i+2)/2);
 
202
                fprintf(code,"};\n");
 
203
                fprintf(code,"  int bcdm1_num = io[lb]*%d;\n\n",cdm1_num);
 
204
              }
 
205
            }
 
206
          }
 
207
        }
 
208
      }
 
209
      if (split == 0) {
 
210
        fprintf(code,"        I0 += %d;\n",cp1dm1_num);
 
211
        fprintf(code,"        I1 += %d;\n",cdm1_num);
 
212
        fprintf(code,"        I2 += %d;\n",cp1dm1_num);
 
213
        fprintf(code,"        I3 += %d;\n",cdm1_num);
 
214
        fprintf(code,"        I4 += %d;\n",cdm1_num);
 
215
        fprintf(code,"  }\n}\n");
 
216
      }
 
217
      else {
 
218
        fprintf(code,"  return vp;\n}\n");
 
219
      }
 
220
      fclose(code);
 
221
      printf("Done with %s\n",code_name);
 
222
 
 
223
      
 
224
      
 
225
      /*-----------------------
 
226
        HRR on centers A and B
 
227
       -----------------------*/
 
228
 
 
229
      la = lc-ld;  lb = ld;
 
230
      am_in[0] = la;
 
231
      am_in[1] = lb;
 
232
 
 
233
      class_size = ((am_in[0]+1)*(am_in[0]+2)*(am_in[1]+1)*(am_in[1]+2))/4;
 
234
 
 
235
      /*--- Whether to split has been decided already ---*/
 
236
      
 
237
      sprintf(function_name,"t1hrr1_build_%c%c",am_letter[am_in[0]],am_letter[am_in[1]]);
 
238
      if (split) {
 
239
        subfunction_name = (char **) malloc (num_subfunctions*sizeof(char *));
 
240
        for(i=0;i<num_subfunctions;i++) {
 
241
          subfunction_name[i] = (char *) malloc(22*sizeof(char));
 
242
          sprintf(subfunction_name[i],"_%s_%d",
 
243
                  function_name,i);
 
244
        }
 
245
      }
 
246
      sprintf(code_name,"t1hrr1_build_%c%c.cc",am_letter[am_in[0]],am_letter[am_in[1]]);
 
247
      code = fopen(code_name,"w");
 
248
 
 
249
      /* include the function into the hrr_header.h */
 
250
      fprintf(hrr_header,"void %s(REALTYPE *, REALTYPE *, REALTYPE *, const REALTYPE *, const REALTYPE *, ",function_name);
 
251
      fprintf(hrr_header,"const REALTYPE *, const REALTYPE *, const REALTYPE *, int, int);\n");
 
252
      
 
253
      fprintf(code,"  /* This machine-generated function computes a quartet of (%c%c|[r12,T1] integrals */\n\n",
 
254
              am_letter[am_in[0]],am_letter[am_in[1]]);
 
255
      fprintf(code,"#include<libint/libint.h>\n\n");
 
256
      
 
257
      if (split) {
 
258
        for(i=0;i<num_subfunctions;i++) {
 
259
          fprintf(code,"REALTYPE *%s(REALTYPE *, REALTYPE *, REALTYPE *, const REALTYPE *, const REALTYPE *, ",
 
260
                  subfunction_name[i]);
 
261
          fprintf(code,"const REALTYPE *, const REALTYPE *, const REALTYPE *, int, int);\n");
 
262
        }
 
263
        fprintf(code,"\n");
 
264
      }
 
265
      fprintf(code,"void %s(REALTYPE *AB, REALTYPE *AC, REALTYPE *vp, const REALTYPE *I0, const REALTYPE *I1, ",
 
266
              function_name);
 
267
      fprintf(code,"const REALTYPE *I2, const REALTYPE *I3, const REALTYPE *I4, int lc, int ld)\n{\n");
 
268
      if (split == 1) {
 
269
        curr_subfunction = 0;
 
270
        curr_count = 0;
 
271
      }
 
272
      else {
 
273
        fprintf(code,"  int cd, cd_num, c_num, cp1_num, d_num;\n");
 
274
        fprintf(code,"  int pc, qc, d, ind_c, ind_cp1d;\n");
 
275
        fprintf(code,"  const REALTYPE *i0, *i1, *i2, *i3, *i4;\n");
 
276
        fprintf(code,"  const REALTYPE AB0 = AB[0];\n");
 
277
        fprintf(code,"  const REALTYPE AB1 = AB[1];\n");
 
278
        fprintf(code,"  const REALTYPE AB2 = AB[2];\n");
 
279
        fprintf(code,"  const REALTYPE AC0 = AC[0];\n");
 
280
        fprintf(code,"  const REALTYPE AC1 = AC[1];\n");
 
281
        fprintf(code,"  const REALTYPE AC2 = AC[2];\n");
 
282
        fprintf(code,"  static int io[] = { 1");
 
283
        for(i=1;i<=new_am+1;i++)
 
284
          fprintf(code,", %d",(i+1)*(i+2)/2);
 
285
        fprintf(code,"};\n\n");
 
286
        fprintf(code,"  c_num = io[lc];\n");
 
287
        fprintf(code,"  cp1_num = io[lc+1];\n");
 
288
        fprintf(code,"  d_num = io[ld];\n\n");
 
289
      }
 
290
 
 
291
      if (split == 1) {
 
292
        for(f=0;f<num_subfunctions;f++)
 
293
          fprintf(code,"  vp = %s(AB, AC, vp, I0, I1, I2, I3, I4, lc, ld);\n",
 
294
                subfunction_name[f]);
 
295
        fprintf(code,"}\n\n");
 
296
 
 
297
        fprintf(code,"REALTYPE *%s(REALTYPE *AB, REALTYPE *AC, REALTYPE *vp, const REALTYPE *I0, const REALTYPE *I1, ",
 
298
                subfunction_name[0]);
 
299
        fprintf(code,"const REALTYPE *I2, const REALTYPE *I3, const REALTYPE *I4, int lc, int ld)\n{\n");
 
300
        fprintf(code,"  int cd, cd_num, c_num, cp1_num, d_num;\n");
 
301
        fprintf(code,"  int pc, qc, d, ind_c, ind_cp1d;\n");
 
302
        fprintf(code,"  const REALTYPE *i0, *i1, *i2, *i3, *i4;\n");
 
303
        fprintf(code,"  const REALTYPE AB0 = AB[0];\n");
 
304
        fprintf(code,"  const REALTYPE AB1 = AB[1];\n");
 
305
        fprintf(code,"  const REALTYPE AB2 = AB[2];\n");
 
306
        fprintf(code,"  const REALTYPE AC0 = AC[0];\n");
 
307
        fprintf(code,"  const REALTYPE AC1 = AC[1];\n");
 
308
        fprintf(code,"  const REALTYPE AC2 = AC[2];\n");
 
309
        fprintf(code,"  static int io[] = { 1");
 
310
        for(i=1;i<=new_am+1;i++)
 
311
          fprintf(code,", %d",(i+1)*(i+2)/2);
 
312
        fprintf(code,"};\n\n");
 
313
        fprintf(code,"  c_num = io[lc];\n");
 
314
        fprintf(code,"  cp1_num = io[lc+1];\n");
 
315
        fprintf(code,"  d_num = io[ld];\n\n");
 
316
      }
 
317
      
 
318
      nj = (lb*(lb+1))/2;
 
319
 
 
320
      for(p = 0; p <= am_in[0]; p++){
 
321
        am[0][0] = am_in[0] - p;
 
322
        for(q = 0; q <= p; q++){
 
323
          am[0][1] = p - q;
 
324
          am[0][2] = q;
 
325
          
 
326
          for(r = 0; r <= am_in[1]; r++){
 
327
            am[1][0] = am_in[1] - r;
 
328
            for(s = 0; s <= r; s++){
 
329
              am[1][1] = r - s;
 
330
              am[1][2] = s;
 
331
 
 
332
              if (am[1][0]) /* build along x */
 
333
                xyz = 0;
 
334
              else if (am[1][1]) /* build along y */
 
335
                xyz = 1;
 
336
              else /* build along z */
 
337
                xyz = 2;
 
338
 
 
339
              am[0][xyz] += 1;
 
340
              am_in[0] += 1;
 
341
              am[1][xyz] -= 1;
 
342
              am_in[1] -= 1;
 
343
              t0 = hash(am,am_in);
 
344
              am[0][xyz] -= 1;
 
345
              am_in[0] -= 1;
 
346
              t1 = hash(am,am_in);
 
347
              am[1][xyz] += 1;
 
348
              am_in[1] += 1;
 
349
              t2 = t0;
 
350
              t3 = t1;
 
351
              t4 = t1;
 
352
              
 
353
              if (t0) {
 
354
                fprintf(code,"  i0 = I0 + %d*(c_num*d_num);\n",t0);
 
355
                fprintf(code,"  i2 = I2 + %d*(c_num*d_num);\n",t2);
 
356
              }
 
357
              else {
 
358
                fprintf(code,"  i0 = I0;\n");
 
359
                fprintf(code,"  i2 = I2;\n");
 
360
              }
 
361
              if (t1) {
 
362
                fprintf(code,"  i1 = I1 + %d*(c_num*d_num);\n",t1);
 
363
                fprintf(code,"  i4 = I4 + %d*(c_num*d_num);\n",t4);
 
364
                fprintf(code,"  i3 = I3 + %d*(cp1_num*d_num);\n",t3);
 
365
              }
 
366
              else {
 
367
                fprintf(code,"  i1 = I1;\n");
 
368
                fprintf(code,"  i3 = I3;\n");
 
369
                fprintf(code,"  i4 = I4;\n");
 
370
              }
 
371
 
 
372
              if (xyz == 0) {
 
373
                fprintf(code,"  for(cd=0;cd<c_num*d_num;cd++)\n");
 
374
                fprintf(code,"    *(vp++) = *(i0++) + AB0*(*(i1++)) + *(i2++) + AC0*(*(i4++)) - *(i3++);\n\n");
 
375
              }
 
376
              else {
 
377
                fprintf(code,"  ind_c = %d;\n",xyz);
 
378
                fprintf(code,"  for(pc=0;pc<=lc;pc++) {\n");
 
379
                fprintf(code,"    for(qc=0;qc<=pc;qc++) {\n");
 
380
                fprintf(code,"      ind_cp1d = (ind_c + pc)*d_num;\n");
 
381
                fprintf(code,"      for(d=0;d<d_num;d++) {\n");
 
382
                fprintf(code,"        *(vp++) = *(i0++) + AB%d*(*(i1++)) + *(i2++) + AC%d*(*(i4++)) - i3[ind_cp1d];\n",
 
383
                        xyz,xyz);
 
384
                fprintf(code,"        ind_cp1d++;\n");
 
385
                fprintf(code,"      }\n");
 
386
                fprintf(code,"      ind_c++;\n");
 
387
                fprintf(code,"    }\n");
 
388
                fprintf(code,"  }\n\n");
 
389
              }
 
390
 
 
391
              curr_count++;
 
392
              if (curr_count == subbatch_length && split == 1) {
 
393
                curr_count = 0;
 
394
                curr_subfunction++;
 
395
                fprintf(code,"  return vp;\n}\n\n");
 
396
                fprintf(code,"REALTYPE *%s(REALTYPE *AB, REALTYPE *AC, REALTYPE *vp, const REALTYPE *I0, const REALTYPE *I1, ",
 
397
                        subfunction_name[curr_subfunction]);
 
398
                fprintf(code,"const REALTYPE *I2, const REALTYPE *I3, const REALTYPE *I4, int lc, int ld)\n{\n");
 
399
                fprintf(code,"  int cd, cd_num, c_num, cp1_num, d_num;\n");
 
400
                fprintf(code,"  int pc, qc, d, ind_c, ind_cp1d;\n");
 
401
                fprintf(code,"  const REALTYPE *i0, *i1, *i2, *i3, *i4;\n");
 
402
                fprintf(code,"  const REALTYPE AB0 = AB[0];\n");
 
403
                fprintf(code,"  const REALTYPE AB1 = AB[1];\n");
 
404
                fprintf(code,"  const REALTYPE AB2 = AB[2];\n");
 
405
                fprintf(code,"  const REALTYPE AC0 = AC[0];\n");
 
406
                fprintf(code,"  const REALTYPE AC1 = AC[1];\n");
 
407
                fprintf(code,"  const REALTYPE AC2 = AC[2];\n");
 
408
                fprintf(code,"  static int io[] = { 1");
 
409
                for(i=1;i<=new_am+1;i++)
 
410
                  fprintf(code,", %d",(i+1)*(i+2)/2);
 
411
                fprintf(code,"};\n\n");
 
412
                fprintf(code,"  c_num = io[lc];\n");
 
413
                fprintf(code,"  cp1_num = io[lc+1];\n");
 
414
                fprintf(code,"  d_num = io[ld];\n\n");
 
415
              }
 
416
            }
 
417
          }
 
418
        }
 
419
      }
 
420
      if (split == 1) {
 
421
        fprintf(code,"  return vp;\n");
 
422
      }
 
423
      fprintf(code,"}\n");
 
424
      fclose(code);
 
425
      printf("Done with %s\n",code_name);
 
426
    }
 
427
  }
 
428
}
 
429
 
 
430
 
 
431
/*----------------------------------------------------------------------------------
 
432
  hash(a,b) returns a number of the (a[0] a[1]) type product within a doublet.
 
433
  a contains x y and z exponents of functions on centers A and B, and b contains
 
434
  their angular momenta
 
435
 ----------------------------------------------------------------------------------*/
 
436
 
 
437
int hash(a, b)
 
438
  int a[2][3];
 
439
  int b[2];
 
440
{
 
441
  int c[2] = {0,0};
 
442
  int i;
 
443
  static int io[] = {0,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153};
 
444
 
 
445
  if(b[0]){
 
446
    i=b[0]-a[0][0];
 
447
    c[0]=i+io[i]-a[0][1];
 
448
    }
 
449
  if(b[1]){
 
450
    i=b[1]-a[1][0];
 
451
    c[1]=i+io[i]-a[1][1];
 
452
    }
 
453
 
 
454
  return c[0]*io[b[1]+1]+c[1];
 
455
}