~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/lib/libdpd/dot14.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*! \file
 
2
    \ingroup DPD
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <libqt/qt.h>
 
7
#include "dpd.h"
 
8
 
 
9
extern "C" {
 
10
 
 
11
/* the dot functions have not been tested for transposed cases */
 
12
 
 
13
int dpd_dot14(dpdfile2 *T, dpdbuf4 *I, dpdfile2 *Z,
 
14
    int transt, int transz, double alpha, double beta)
 
15
{
 
16
  int h, Gp, Gq, Gr, Gs, GT, GI, GZ, Tblock, Zblock;
 
17
  int p, q, r, s;
 
18
  int P, Q, R, S;
 
19
  int row, col;
 
20
  int nirreps;
 
21
  double **X;
 
22
  double value;
 
23
 
 
24
  nirreps = T->params->nirreps;
 
25
  GT = T->my_irrep;
 
26
  GI = I->file.my_irrep;
 
27
  GZ = Z->my_irrep;
 
28
 
 
29
  /* Get the two-index quantities from disk */
 
30
  dpd_file2_mat_init(T);
 
31
  dpd_file2_mat_rd(T);
 
32
  dpd_file2_scm(Z, beta);
 
33
  dpd_file2_mat_init(Z);
 
34
  dpd_file2_mat_rd(Z);
 
35
 
 
36
#ifdef DPD_TIMER
 
37
  timer_on("dot14");
 
38
#endif
 
39
 
 
40
  /* loop over irreps of buffer I; h = Gpq */
 
41
  for(h=0; h < nirreps; h++) {
 
42
 
 
43
    dpd_buf4_mat_irrep_init(I, h);
 
44
    dpd_buf4_mat_irrep_rd(I, h);
 
45
 
 
46
    /* Loop over irreps of the target */
 
47
    for(Gq=0; Gq < nirreps; Gq++) {
 
48
      /*Gr = Gq;  Gp = Gs = h^Gq; */
 
49
      Gp = h^Gq; Gr = Gq^GZ; Gs = h^Gq^GT;
 
50
      if (!transt) Tblock = Gp; else Tblock = Gs;
 
51
      if (!transz) Zblock = Gq; else Zblock = Gr;
 
52
 
 
53
      /* Allocate space for the X buffer */
 
54
      if(T->params->ppi[Gp] && T->params->qpi[Gs])
 
55
        X = dpd_block_matrix(T->params->ppi[Gp],T->params->qpi[Gs]);
 
56
 
 
57
      /* Loop over orbitals of the target */
 
58
      for(q=0; q < Z->params->ppi[Gq]; q++) {
 
59
        Q = Z->params->poff[Gq] + q;
 
60
        for(r=0; r < Z->params->qpi[Gr]; r++) {
 
61
          R = Z->params->qoff[Gr] + r;
 
62
 
 
63
          /* Loop over orbitals of the two-index term */
 
64
          for(p=0; p < T->params->ppi[Gp]; p++) {
 
65
            P = T->params->poff[Gp] + p;
 
66
            for(s=0; s < T->params->qpi[Gs]; s++) {
 
67
              S = T->params->qoff[Gs] + s;
 
68
 
 
69
              /* Calculate row and column indices in I */
 
70
              if(!transt && !transz) {
 
71
                row = I->params->rowidx[P][Q];
 
72
                col = I->params->colidx[R][S];
 
73
              }
 
74
              else if(transt && !transz) {
 
75
                row = I->params->rowidx[S][Q];
 
76
                col = I->params->colidx[R][P];
 
77
              }
 
78
              else if(!transt && transz) {
 
79
                row = I->params->rowidx[P][R];
 
80
                col = I->params->colidx[Q][S];
 
81
              }
 
82
              else if(transt && transz) {
 
83
                row = I->params->rowidx[S][R];
 
84
                col = I->params->colidx[Q][P];
 
85
              }
 
86
 
 
87
              /* Build the X buffer */
 
88
              X[p][s] = I->matrix[h][row][col]; 
 
89
 
 
90
            }
 
91
          }
 
92
 
 
93
          value = dot_block(T->matrix[Tblock], X, T->params->ppi[Gp],
 
94
              T->params->qpi[Gs], alpha); 
 
95
 
 
96
          Z->matrix[Zblock][q][r] += value;
 
97
        }
 
98
      }
 
99
      if(T->params->ppi[Gp] && T->params->qpi[Gs])
 
100
        dpd_free_block(X,T->params->ppi[Gp],T->params->qpi[Gs]);
 
101
    }
 
102
    dpd_buf4_mat_irrep_close(I, h);
 
103
  }
 
104
 
 
105
#ifdef DPD_TIMER
 
106
  timer_off("dot14");
 
107
#endif
 
108
 
 
109
  /* Close the two-index quantities */
 
110
  dpd_file2_mat_close(T);
 
111
  dpd_file2_mat_wrt(Z);
 
112
  dpd_file2_mat_close(Z);
 
113
 
 
114
  return 0;
 
115
}
 
116
 
 
117
} /* extern "C" */