7
#include <libciomr/libciomr.h>
11
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
12
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
18
void IWL::sort_buffer(IWL *Inbuf, IWL *Outbuf,
19
double *ints, int fpq, int lpq, int *ioff, int *ioff2,
20
int nbfso, int elbert, int intermediate, int no_pq_perm,
21
int qdim, int add, int printflg, FILE *outfile)
24
Value *valptr; /* array of integral values */
25
Label *lblptr; /* array of integral labels */
26
int idx; /* index for curr integral (0..ints_per_buf) */
27
int lastbuf; /* last buffer flag */
28
int p, q, qmax, qmin, r, rmin, rmax, s, smin, smax, pq, rs;
29
long int pqrs, offset;
30
int first_p, first_q, first_pq, last_p, last_q;
34
fprintf(outfile, "\nsortbuf for pq=%d to %d\n", fpq, lpq);
37
if (no_pq_perm && !intermediate) {
38
fprintf(outfile,"(sortbuf): illegal parameter combination.\n");
39
fprintf(stderr, "(sortbuf): illegal parameter combination.\n");
42
nbstri = nbfso * (nbfso + 1) / 2;
44
/* figure out ranges on things */
45
/* I believe this section works fine, even with different ioff arrays */
47
while (fpq >= ioff[i] && i < BIGNUM) i++;
49
fprintf(outfile, "(sortbuf): parameter error\n") ;
52
first_p = i-1 ; first_q = fpq - ioff[i-1];
53
first_pq = ioff[first_p] + first_q;
54
if (first_pq != fpq) {
55
fprintf(outfile, "(sortbuf): fpq != first_pq.\n");
56
fprintf(stderr, "(sortbuf): fpq != first_pq.\n");
60
if (elbert) offset = ioff2[first_pq] + first_pq ;
61
else offset = ioff[first_pq] ;
66
while (lpq >= ioff[i] && i < BIGNUM) i++ ;
68
fprintf(outfile, "(sortbuf): parameter error\n") ;
71
last_p = i-1 ; last_q = lpq - ioff[i-1] ;
74
lblptr = Inbuf->labels_;
75
valptr = Inbuf->values_;
77
/* read a buffer at a time until we're done */
81
lastbuf = Inbuf->lastbuf_;
82
for (idx=4*Inbuf->idx_; Inbuf->idx_<Inbuf->inbuf_; Inbuf->idx_++) {
83
p = (int) lblptr[idx++];
84
q = (int) lblptr[idx++];
85
r = (int) lblptr[idx++];
86
s = (int) lblptr[idx++];
88
/* if (no_pq_perm) ioff is the appropriate offset array for the left
89
indices (ioff[p] = nvirt * p for MP2); ioff2 is then the usual
90
ioff offset array, used for the right indices */
93
rs = ioff2[MAX0(r,s)] + MIN0(r,s);
96
pq = ioff[MAX0(p,q)] + MIN0(p,q);
97
rs = ioff[MAX0(r,s)] + MIN0(r,s);
102
pqrs = ioff2[pq] + rs;
104
pqrs = ioff[MAX0(pq,rs)];
109
pqrs = (pq - first_pq);
114
if (printflg && ints[pqrs-offset] != 0.0)
115
fprintf(outfile, "Adding %10.6f to el %d %d %d %d = %10.6lf\n",
116
valptr[Inbuf->idx_], p, q, r, s, ints[pqrs-offset]);
118
if (add) ints[pqrs-offset] += valptr[Inbuf->idx_];
119
else ints[pqrs-offset] += valptr[Inbuf->idx_];
122
fprintf(outfile, "<%d %d %d %d | %d %d [%ld] = %10.6f\n",
123
p, q, r, s, pq, rs, pqrs, ints[pqrs-offset]) ;
127
/* now write them out again, in order */
128
lblptr = Outbuf->labels_;
129
valptr = Outbuf->values_;
133
for (p=first_p; p<=last_p; p++) {
134
qmax = (p==last_p) ? last_q : p ;
135
qmin = (p==first_p) ? first_q : 0 ;
137
qmax = (p==last_p) ? last_q : (qdim-1);
138
qmin = (p==first_p) ? first_q: 0;
140
for (q=qmin; q<=qmax; q++) {
141
pq = ioff[p] + q ; /* This should be fine even with MP2 */
144
rmin = (elbert) ? p : 0 ;
145
rmax = (elbert) ? nbfso : p+1 ;
147
else { /* This should be fine with MP2, also */
152
for (r=rmin; r<rmax; r++) {
157
smin = (p==r) ? q : 0 ;
160
smax = (p==r) ? (q+1) : (r+1) ;
164
else { /* This should be fine with MP2, also */
169
for (s=smin; s < smax; s++) {
170
if(no_pq_perm) rs = ioff2[r] + s;
171
else rs = ioff[r] + s;
173
/* Again, this should be fine with MP2 */
174
if (elbert) pqrs = ioff2[pq] + rs ;
175
else if (intermediate) {
176
pqrs = (pq - first_pq);
180
else pqrs = ioff[pq] + rs ;
182
if (fabs(ints[pqrs-offset]) > Outbuf->cutoff_) {
183
idx = 4*Outbuf->idx_;
188
valptr[Outbuf->idx_] = ints[pqrs-offset];
190
fprintf(outfile, ">%d %d %d %d | %d %d [%ld] = %10.6f\n",
191
p, q, r, s, pq, rs, pqrs, ints[pqrs-offset]) ;
194
if (Outbuf->idx_ == Outbuf->ints_per_buf_) {
195
Outbuf->lastbuf_ = 0;
196
Outbuf->inbuf_ = Outbuf->idx_;
212
** Function reads a file of two-electron integrals into
213
** core and writes them back out again in canonical order. Used in
214
** Yoshimine sorts where we have a file containing all rs for a few
215
** values of pq, but the ints are not in canonical order. At the
216
** very least, we need to sort to make sure that all (pq|rs) for a
217
** given pq are grouped together, since the transformation program
218
** wants to work with all rs for a single pq value at one time.
219
** We may or may not use the restriction pq >= rs (not used if
220
** intermediate = 1, which is how this routine is always called
223
** One interesting issue here is that the intermediate array ('ints')
224
** must be big enough to hold the integrals in the current buffer, but
225
** we don't generally want it to be much larger than necessary! Thus
226
** we calculate an 'offset' which is the canonical index of the first
227
** integral in the buffer, and we use this so that the first integral
228
** in the buffer is stored in ints[0]. What's different for the
229
** Elbert ordering is that we have all rs for a given pq but rs >= pq!
230
** If we want our canonical indices to be consecutive WITHIN THE
231
** CURRENT BUFFER, we MUST use upper triangle ordering rather than
232
** lower triangle! That's what ioff2 is used for. Obviously the
233
** offset must also be calculated with ioff2 for the Elbert order.
234
** Formula for ioff2: ioff2[0] = 0; ioff2[i] = ioff2[i-1] + n - i;
235
** Note that this is not the case when this routine is used for MP2
236
** sorts. The definitions of the ioff arrays can be confusing, so
237
** care should be taken when using this routine.
239
** \param Inbuf = IWL buffer for input
240
** \param Outbuf = IWL buffer for output
241
** \param ints = array to hold integrals in
242
** \param fpq = first pq for this tape
243
** \param lpq = last pq for this tape
244
** \param ioff = offset array for the left indices
245
** \param ioff2 = offset array for Elbert sorts or for the right
246
** indices when no_pq_perm=1
247
** \param nbfso = number of basis functions in SO's
248
** \param lastsort = 1 if this is the last intape, 0 otherwise
249
** \param elbert = integrals obey rs >= pq. Use ioff2 to get offset.
250
** \param intermediate= 1 if sorting a intermediate in the transformation
251
** which is indexed as X[ij][kl] where ij runs from
252
** fpq to lpq and kl runs from 0 to nbstri
253
** \param no_pq_perm = don't use permutational symmetry to swap p and q
254
** (appropriate for MP2 where one is occ and one is virt)
255
** \param qdim = dimensions for the q index...nvirt for MP2
256
** \param add = add contributions to the same integral during sort
257
** \param printflg = 1 for printing, 0 otherwise
258
** \param outfile = output file pointer
262
** Revised 6/27/96 by CDS for new IWL format
263
** N.B. Now need to iwl_flush the output buffer...not done in here!!
266
void sortbuf(struct iwlbuf *Inbuf, struct iwlbuf *Outbuf,
267
double *ints, int fpq, int lpq, int *ioff, int *ioff2,
268
int nbfso, int elbert, int intermediate, int no_pq_perm,
269
int qdim, int add, int printflg, FILE *outfile)
272
Value *valptr; /* array of integral values */
273
Label *lblptr; /* array of integral labels */
274
int idx; /* index for curr integral (0..ints_per_buf) */
275
int lastbuf; /* last buffer flag */
276
int p, q, qmax, qmin, r, rmin, rmax, s, smin, smax, pq, rs;
277
long int pqrs, offset;
278
int first_p, first_q, first_pq, last_p, last_q;
282
fprintf(outfile, "\nsortbuf for pq=%d to %d\n", fpq, lpq);
285
if (no_pq_perm && !intermediate) {
286
fprintf(outfile,"(sortbuf): illegal parameter combination.\n");
287
fprintf(stderr, "(sortbuf): illegal parameter combination.\n");
290
nbstri = nbfso * (nbfso + 1) / 2;
292
/* figure out ranges on things */
293
/* I believe this section works fine, even with different ioff arrays */
295
while (fpq >= ioff[i] && i < BIGNUM) i++;
297
fprintf(outfile, "(sortbuf): parameter error\n") ;
300
first_p = i-1 ; first_q = fpq - ioff[i-1];
301
first_pq = ioff[first_p] + first_q;
302
if (first_pq != fpq) {
303
fprintf(outfile, "(sortbuf): fpq != first_pq.\n");
304
fprintf(stderr, "(sortbuf): fpq != first_pq.\n");
308
if (elbert) offset = ioff2[first_pq] + first_pq ;
309
else offset = ioff[first_pq] ;
314
while (lpq >= ioff[i] && i < BIGNUM) i++ ;
316
fprintf(outfile, "(sortbuf): parameter error\n") ;
319
last_p = i-1 ; last_q = lpq - ioff[i-1] ;
322
lblptr = Inbuf->labels;
323
valptr = Inbuf->values;
325
/* read a buffer at a time until we're done */
328
iwl_buf_fetch(Inbuf);
329
lastbuf = Inbuf->lastbuf;
330
for (idx=4*Inbuf->idx; Inbuf->idx<Inbuf->inbuf; Inbuf->idx++) {
331
p = (int) lblptr[idx++];
332
q = (int) lblptr[idx++];
333
r = (int) lblptr[idx++];
334
s = (int) lblptr[idx++];
336
/* if (no_pq_perm) ioff is the appropriate offset array for the left
337
indices (ioff[p] = nvirt * p for MP2); ioff2 is then the usual
338
ioff offset array, used for the right indices */
341
rs = ioff2[MAX0(r,s)] + MIN0(r,s);
344
pq = ioff[MAX0(p,q)] + MIN0(p,q);
345
rs = ioff[MAX0(r,s)] + MIN0(r,s);
350
pqrs = ioff2[pq] + rs;
352
pqrs = ioff[MAX0(pq,rs)];
357
pqrs = (pq - first_pq);
362
if (printflg && ints[pqrs-offset] != 0.0)
363
fprintf(outfile, "Adding %10.6lf to el %d %d %d %d = %10.6lf\n",
364
valptr[Inbuf->idx], p, q, r, s, ints[pqrs-offset]);
366
if (add) ints[pqrs-offset] += valptr[Inbuf->idx];
367
else ints[pqrs-offset] += valptr[Inbuf->idx];
370
fprintf(outfile, "<%d %d %d %d | %d %d [%ld] = %10.6lf\n",
371
p, q, r, s, pq, rs, pqrs, ints[pqrs-offset]) ;
375
/* now write them out again, in order */
376
lblptr = Outbuf->labels;
377
valptr = Outbuf->values;
381
for (p=first_p; p<=last_p; p++) {
382
qmax = (p==last_p) ? last_q : p ;
383
qmin = (p==first_p) ? first_q : 0 ;
385
qmax = (p==last_p) ? last_q : (qdim-1);
386
qmin = (p==first_p) ? first_q: 0;
388
for (q=qmin; q<=qmax; q++) {
389
pq = ioff[p] + q ; /* This should be fine even with MP2 */
392
rmin = (elbert) ? p : 0 ;
393
rmax = (elbert) ? nbfso : p+1 ;
395
else { /* This should be fine with MP2, also */
400
for (r=rmin; r<rmax; r++) {
405
smin = (p==r) ? q : 0 ;
408
smax = (p==r) ? (q+1) : (r+1) ;
412
else { /* This should be fine with MP2, also */
417
for (s=smin; s < smax; s++) {
418
if(no_pq_perm) rs = ioff2[r] + s;
419
else rs = ioff[r] + s;
421
/* Again, this should be fine with MP2 */
422
if (elbert) pqrs = ioff2[pq] + rs ;
423
else if (intermediate) {
424
pqrs = (pq - first_pq);
428
else pqrs = ioff[pq] + rs ;
430
if (fabs(ints[pqrs-offset]) > Outbuf->cutoff) {
436
valptr[Outbuf->idx] = ints[pqrs-offset];
438
fprintf(outfile, ">%d %d %d %d | %d %d [%ld] = %10.6lf\n",
439
p, q, r, s, pq, rs, pqrs, ints[pqrs-offset]) ;
442
if (Outbuf->idx == Outbuf->ints_per_buf) {
444
Outbuf->inbuf = Outbuf->idx;
b'\\ No newline at end of file'
added
removed
src/lib/libiwl/sortbuf.cc
