← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccsort/build_abcd_packed.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/ccsort/build_abcd_packed.c
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
#include <math.h>
 
4
#include <libciomr/libciomr.h>
 
5
#include <libpsio/psio.h>
 
6
#include <libiwl/iwl.h>
 
7
#include <ccfiles.h>
 
8
#include <psifiles.h>
 
9
#include <libdpd/dpd.h>
 
10
#include "MOInfo.h"
 
11
#include "Params.h"
 
12
#define EXTERN
 
13
#include "globals.h"
 
14
 
 
15
void idx_permute_multipass(dpdfile4 *File, int this_bucket,
 
16
                           int **bucket_map, int **bucket_offset,
 
17
                           int p, int q, int r, int s,
 
18
                           int perm_pr, int perm_qs, int perm_prqs,
 
19
                           double value, FILE *outfile);
 
20
 
 
21
int build_abcd_packed(int inputfile, double tolerance, int keep)
 
22
{
 
23
  struct iwlbuf InBuf;
 
24
  int lastbuf;
 
25
  long int memoryb, memoryd, core_left, row_length;
 
26
  int h, nirreps, n, row, nump, numq, nbuckets;
 
27
  int **bucket_map, **bucket_offset, **bucket_rowdim;
 
28
  long int **bucket_size;
 
29
  Value *valptr;
 
30
  Label *lblptr;
 
31
  int idx, p, q, r, s;
 
32
  int ab, cd, c, d, CD, DC;
 
33
  double value, abcd, abdc;
 
34
  struct iwlbuf *SortBuf;
 
35
  psio_address next;
 
36
  dpdfile4 B;
 
37
  dpdbuf4 B_s, B_a;
 
38
  double **B_diag;
 
39
  int rows_per_bucket, rows_left, row_start, nvirt;
 
40
  int m, Gc, C, cc;
 
41
 
 
42
  dpd_file4_init_nocache(&B, CC_BINTS, 0, 8, 5, "B(+) <ab|cd>");  /* junk target */
 
43
  dpd_buf4_init(&B_s, CC_BINTS, 0, 8, 8, 8, 8, 0, "B(+) <ab|cd> + <ab|dc>");
 
44
  dpd_buf4_scm(&B_s, 0.0);
 
45
  dpd_buf4_init(&B_a, CC_BINTS, 0, 9, 9, 9, 9, 0, "B(-) <ab|cd> - <ab|dc>");
 
46
  dpd_buf4_scm(&B_a, 0.0);
 
47
 
 
48
  nirreps = B.params->nirreps;
 
49
 
 
50
  fndcor(&memoryb, infile, outfile);
 
51
  memoryd = memoryb/sizeof(double);
 
52
 
 
53
  /* It's annoying that I have to compute this here */
 
54
  for(h=0,nump=0,numq=0; h < B.params->nirreps; h++) {
 
55
    nump += B.params->ppi[h];
 
56
    numq += B.params->qpi[h];
 
57
  }
 
58
  bucket_map = init_int_matrix(nump,numq);
 
59
  for(p=0; p < nump; p++)
 
60
    for(q=0; q < numq; q++)
 
61
      bucket_map[p][q] = -1;  /* initialize with junk */
 
62
 
 
63
  /* Room for one bucket to begin with */
 
64
  bucket_offset = (int **) malloc(sizeof(int *));
 
65
  bucket_offset[0] = init_int_array(nirreps);
 
66
  bucket_rowdim = (int **) malloc(sizeof(int *));
 
67
  bucket_rowdim[0] = init_int_array(nirreps);
 
68
  bucket_size = (long int **) malloc(sizeof(long int *));
 
69
  bucket_size[0] = init_long_int_array(nirreps);
 
70
    
 
71
  /* Figure out how many passes we need and where each p,q goes */
 
72
  for(h=0,core_left=memoryd,nbuckets=1; h < nirreps; h++) {
 
73
 
 
74
    row_length = (long int) B.params->coltot[h^(B.my_irrep)];
 
75
               
 
76
    for(row=0; row < B.params->rowtot[h]; row++) {
 
77
 
 
78
      if((core_left - row_length) >= 0) {
 
79
        core_left -= row_length;
 
80
        bucket_rowdim[nbuckets-1][h]++;
 
81
        bucket_size[nbuckets-1][h] += row_length;
 
82
      }
 
83
      else {
 
84
        nbuckets++;
 
85
        core_left = memoryd - row_length;
 
86
 
 
87
        /* Make room for another bucket */
 
88
        bucket_offset = (int **) realloc((void *) bucket_offset,
 
89
                                         nbuckets * sizeof(int *));
 
90
        bucket_offset[nbuckets-1] = init_int_array(nirreps);
 
91
        bucket_offset[nbuckets-1][h] = row;
 
92
 
 
93
        bucket_rowdim = (int **) realloc((void *) bucket_rowdim,
 
94
                                         nbuckets * sizeof(int *));
 
95
        bucket_rowdim[nbuckets-1] = init_int_array(nirreps);
 
96
        bucket_rowdim[nbuckets-1][h] = 1;
 
97
 
 
98
        bucket_size = (long int **) realloc((void *) bucket_size,
 
99
                                            nbuckets * sizeof(long int *));
 
100
        bucket_size[nbuckets-1] = init_long_int_array(nirreps);
 
101
        bucket_size[nbuckets-1][h] = row_length;
 
102
      }
 
103
 
 
104
      p = B.params->roworb[h][row][0];
 
105
      q = B.params->roworb[h][row][1];
 
106
      bucket_map[p][q] = nbuckets - 1;
 
107
    }
 
108
  }
 
109
 
 
110
  fprintf(outfile, "\tSorting File: %s nbuckets = %d\n", B.label, nbuckets);
 
111
  fflush(outfile);
 
112
 
 
113
  next = PSIO_ZERO;
 
114
  for(n=0; n < nbuckets; n++) { /* nbuckets = number of passes */
 
115
 
 
116
    /* Prepare target matrix */
 
117
    for(h=0; h < nirreps; h++) {
 
118
      B.matrix[h] = block_matrix(bucket_rowdim[n][h], B.params->coltot[h]);
 
119
    }
 
120
 
 
121
    iwl_buf_init(&InBuf, inputfile, tolerance, 1, 1);
 
122
 
 
123
    lblptr = InBuf.labels;
 
124
    valptr = InBuf.values;
 
125
    lastbuf = InBuf.lastbuf;
 
126
 
 
127
    for(idx=4*InBuf.idx; InBuf.idx < InBuf.inbuf; InBuf.idx++) {
 
128
      p = (int) lblptr[idx++];
 
129
      q = (int) lblptr[idx++];
 
130
      r = (int) lblptr[idx++];
 
131
      s = (int) lblptr[idx++];
 
132
 
 
133
      value = (double) valptr[InBuf.idx];
 
134
 
 
135
      idx_permute_multipass(&B,n,bucket_map,bucket_offset,
 
136
                            p,q,r,s,1,1,1,value,outfile);
 
137
 
 
138
    } /* end loop through current buffer */
 
139
 
 
140
    /* Now run through the rest of the buffers in the file */
 
141
    while (!lastbuf) {
 
142
      iwl_buf_fetch(&InBuf);
 
143
      lastbuf = InBuf.lastbuf;
 
144
 
 
145
      for (idx=4*InBuf.idx; InBuf.idx < InBuf.inbuf; InBuf.idx++) {
 
146
        p = (int) lblptr[idx++];
 
147
        q = (int) lblptr[idx++];
 
148
        r = (int) lblptr[idx++];
 
149
        s = (int) lblptr[idx++];
 
150
 
 
151
        value = (double) valptr[InBuf.idx];
 
152
 
 
153
        idx_permute_multipass(&B,n,bucket_map,bucket_offset,
 
154
                              p,q,r,s,1,1,1,value,outfile);
 
155
      
 
156
      } /* end loop through current buffer */
 
157
    } /* end loop over reading buffers */
 
158
 
 
159
    iwl_buf_close(&InBuf, 1); /* close buffer for next pass */
 
160
 
 
161
    for(h=0; h < nirreps;h++) {
 
162
      if(bucket_size[n][h]) {
 
163
 
 
164
/*      psio_write(B.filenum, B.label, (char *) B.matrix[h][0], */
 
165
/*                 bucket_size[n][h]*((long int) sizeof(double)), next, &next); */
 
166
 
 
167
        dpd_buf4_mat_irrep_row_init(&B_s, h);
 
168
        dpd_buf4_mat_irrep_row_init(&B_a, h);
 
169
        for(ab=0; ab < bucket_rowdim[n][h]; ab++) {
 
170
          for(cd=0; cd < B_s.params->coltot[h]; cd++) {
 
171
            c = B_s.params->colorb[h][cd][0];
 
172
            d = B_s.params->colorb[h][cd][1];
 
173
            CD = B.params->colidx[c][d];
 
174
            DC = B.params->colidx[d][c];
 
175
            abcd = B.matrix[h][ab][CD];
 
176
            abdc = B.matrix[h][ab][DC];
 
177
            B_s.matrix[h][0][cd] = abcd + abdc;
 
178
            B_a.matrix[h][0][cd] = abcd - abdc;
 
179
          }
 
180
          dpd_buf4_mat_irrep_row_wrt(&B_s, h, ab+bucket_offset[n][h]);
 
181
          dpd_buf4_mat_irrep_row_wrt(&B_a, h, ab+bucket_offset[n][h]);
 
182
        }
 
183
        dpd_buf4_mat_irrep_row_close(&B_s, h);
 
184
        dpd_buf4_mat_irrep_row_close(&B_a, h);
 
185
      }
 
186
      free_block(B.matrix[h]);
 
187
    }
 
188
 
 
189
  } /* end loop over buckets/passes */
 
190
 
 
191
  /* Get rid of the input integral file */
 
192
  psio_open(inputfile, PSIO_OPEN_OLD);
 
193
  psio_close(inputfile, keep);
 
194
 
 
195
  free_int_matrix(bucket_map,nump);
 
196
 
 
197
  for(n=0; n < nbuckets; n++) {
 
198
    free(bucket_offset[n]);
 
199
    free(bucket_rowdim[n]);
 
200
    free(bucket_size[n]);
 
201
  }
 
202
  free(bucket_offset);
 
203
  free(bucket_rowdim);
 
204
  free(bucket_size);
 
205
 
 
206
  dpd_file4_close(&B);
 
207
  dpd_buf4_close(&B_s);
 
208
  dpd_buf4_close(&B_a);
 
209
 
 
210
  /* Generate <ab|cc> components of B(+) */
 
211
  for(h=0,nvirt=0; h < moinfo.nirreps; h++) nvirt += moinfo.virtpi[h];
 
212
  dpd_buf4_init(&B_s, CC_BINTS, 0, 8, 8, 8, 8, 0, "B(+) <ab|cd> + <ab|dc>");
 
213
 
 
214
  rows_per_bucket = dpd_memfree()/(B_s.params->coltot[0] + nvirt);
 
215
  if(rows_per_bucket > B_s.params->rowtot[0]) rows_per_bucket = B_s.params->rowtot[0];
 
216
  nbuckets = ceil((double) B_s.params->rowtot[0]/(double) rows_per_bucket);
 
217
  rows_left = B_s.params->rowtot[0] % rows_per_bucket;
 
218
 
 
219
  dpd_buf4_mat_irrep_init_block(&B_s, 0, rows_per_bucket);
 
220
  B_diag = dpd_block_matrix(rows_per_bucket, nvirt);
 
221
 
 
222
  next = PSIO_ZERO;
 
223
  for(m=0; m < (rows_left ? nbuckets-1:nbuckets); m++) {
 
224
    row_start = m * rows_per_bucket;
 
225
    dpd_buf4_mat_irrep_rd_block(&B_s, 0, row_start, rows_per_bucket);
 
226
    for(ab=0; ab < rows_per_bucket; ab++)
 
227
      for(Gc=0; Gc < moinfo.nirreps; Gc++)
 
228
        for(C=0; C < moinfo.virtpi[Gc]; C++) {
 
229
          c = moinfo.vir_off[Gc] + C;
 
230
          cc = B_s.params->colidx[c][c];
 
231
          B_diag[ab][c] = B_s.matrix[0][ab][cc];
 
232
        }
 
233
    psio_write(CC_BINTS, "B(+) <ab|cc>", (char *) B_diag[0], rows_per_bucket*nvirt*sizeof(double), next, &next);
 
234
  }
 
235
  if(rows_left) {
 
236
    row_start = m * rows_per_bucket;
 
237
    dpd_buf4_mat_irrep_rd_block(&B_s, 0, row_start, rows_left);
 
238
    for(ab=0; ab < rows_left; ab++)
 
239
      for(Gc=0; Gc < moinfo.nirreps; Gc++)
 
240
        for(C=0; C < moinfo.virtpi[Gc]; C++) {
 
241
          c = moinfo.vir_off[Gc] + C;
 
242
          cc = B_s.params->colidx[c][c];
 
243
          B_diag[ab][c] = B_s.matrix[0][ab][cc];
 
244
        }
 
245
    psio_write(CC_BINTS, "B(+) <ab|cc>", (char *) B_diag[0], rows_left*nvirt*sizeof(double), next, &next);
 
246
  }
 
247
  dpd_free_block(B_diag, rows_per_bucket, nvirt);
 
248
  dpd_buf4_mat_irrep_close_block(&B_s, 0, rows_per_bucket);
 
249
  dpd_buf4_close(&B_s);
 
250
 
 
251
  return 0;
 
252
}