5
/* dpd_trace42_13(): Take a "trace" of the specified indices of a buf4 and put the
6
** result into a file2.
8
** In this version of the trace42 functions, the summation indices are 1 and 3:
10
** B(q,s) = alpha * \Sum_{p==r} A(pq,rs) + beta * B(q,s)
13
** dpdbuf4 *A: Pointer to the input four-index buffer.
14
** dpdfile2 *B: Pointer to the target two-index file.
15
** int transb: Boolean to indicate whether B should be transposed to match the
16
** target index ordering of A.
17
** double alpha: Prefactor for A.
18
** double beta: Prefactor for B.
21
int dpd_trace42_13(dpdbuf4 *A, dpdfile2 *B, int transb, double alpha, double beta)
23
int h, Gp, Gq, Gr, Gs;
29
nirreps = A->params->nirreps;
31
dpd_file2_scm(B, beta);
32
dpd_file2_mat_init(B);
39
/* Read all of A into core */
40
for(h=0; h < nirreps; h++) {
41
dpd_buf4_mat_irrep_init(A, h);
42
dpd_buf4_mat_irrep_rd(A, h);
45
for(h=0; h < nirreps; h++) {
47
for(Gp=0; Gp < nirreps; Gp++) {
52
/* Loop over target indices */
53
for(q=0; q < A->params->qpi[Gq]; q++) {
54
Q = A->params->qoff[Gq] + q;
56
for(s=0; s < A->params->spi[Gs]; s++) {
57
S = A->params->soff[Gs] + s;
59
/* Loop over elements for which P==R */
60
for(p=0; p < A->params->ppi[Gp]; p++) {
61
P = A->params->poff[Gp] + p;
64
pq = A->params->rowidx[P][Q];
65
rs = A->params->colidx[R][S];
67
if(!transb) B->matrix[Gq][q][s] += alpha * A->matrix[h][pq][rs];
68
else B->matrix[Gq][s][q] += alpha * A->matrix[h][pq][rs];
78
for(h=0; h < nirreps; h++) {
79
dpd_buf4_mat_irrep_close(A, h);
86
/* Close the two-index quantities */
88
dpd_file2_mat_close(B);