← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cclambda/halftrans.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/bin/cclambda/halftrans.c
 
1
#include <libdpd/dpd.h>
 
2
#include <libqt/qt.h>
 
3
 
 
4
/* halftrans(): Routine to transform the last two indices of a dpdbuf4
 
5
** between the MO and SO bases.
 
6
** 
 
7
** dpdbuf4 *Buf1: Pointer to the MO dpdbuf4 (already initialized)
 
8
** dpdbuf4 *Buf2: Pointer to the SO dpdbuf4 (already initialized).
 
9
** double ***C:   Pointer to the transformation matrix (symmetry blocked, SO x MO)
 
10
** int nirreps:   The number of irreps in the point group
 
11
** int **mo_row:  A lookup array.  For a dpdbuf4 with MO indices (ij,ab),
 
12
**                given the irrep h of ij (= ab) and the irrep of orbital a, the
 
13
**                array returns the offset of the start of the set of b molecular
 
14
**                orbitals.
 
15
** int **so_row:  Like mo_row, but for a dpdbuf4 with the last two
 
16
**                indices in the SO basis.
 
17
** int *mospi:    The number of MO's per irrep.
 
18
** int *sospi:    The number of SO's per irrep.
 
19
** int type:      0 = MO --> SO; 1 = SO --> MO
 
20
** double alpha:  multiplicative factor for the transformation
 
21
** double beta:   multiplicative factor for the target
 
22
*/
 
23
 
 
24
void halftrans(dpdbuf4 *Buf1, int dpdnum1, dpdbuf4 *Buf2, int dpdnum2, double ***C, int nirreps, 
 
25
               int **mo_row, int **so_row, int *mospi, int *sospi, int type, double alpha, double beta)
 
26
{
 
27
  int h, Gc, Gd, cd, pq, ij;
 
28
  double **X;
 
29
 
 
30
  for(h=0; h < nirreps; h++) {
 
31
 
 
32
    dpd_set_default(dpdnum1);
 
33
    dpd_buf4_mat_irrep_init(Buf1, h);
 
34
 
 
35
    dpd_set_default(dpdnum2);
 
36
    dpd_buf4_mat_irrep_init(Buf2, h);
 
37
 
 
38
    if(type==0) { /* alpha * Buf1 --> beta * Buf2 */
 
39
      if(alpha != 0.0) { dpd_set_default(dpdnum1); dpd_buf4_mat_irrep_rd(Buf1, h); }
 
40
      if(beta != 0.0) { dpd_set_default(dpdnum2); dpd_buf4_mat_irrep_rd(Buf2, h); }
 
41
    }
 
42
    if(type==1) { /* alpha * Buf2 --> beta * Buf1 */
 
43
      if(alpha != 0.0) { dpd_set_default(dpdnum2); dpd_buf4_mat_irrep_rd(Buf2, h); }
 
44
      if(beta != 0.0) { dpd_set_default(dpdnum1); dpd_buf4_mat_irrep_rd(Buf1, h); }
 
45
    }
 
46
 
 
47
    for(Gc=0; Gc < nirreps; Gc++) {
 
48
      Gd = h^Gc;
 
49
 
 
50
      cd = mo_row[h][Gc];
 
51
      pq = so_row[h][Gc];
 
52
 
 
53
      if(mospi[Gc] && mospi[Gd] && sospi[Gc] && sospi[Gd]) {
 
54
 
 
55
        if(type == 0) {
 
56
          X = block_matrix(mospi[Gc],sospi[Gd]);
 
57
 
 
58
          for(ij=0; ij < Buf1->params->rowtot[h]; ij++) {
 
59
 
 
60
            C_DGEMM('n','t', mospi[Gc], sospi[Gd], mospi[Gd], 1.0,
 
61
                    &(Buf1->matrix[h][ij][cd]), mospi[Gd], &(C[Gd][0][0]), mospi[Gd],
 
62
                    0.0, &(X[0][0]), sospi[Gd]);
 
63
 
 
64
            C_DGEMM('n','n', sospi[Gc], sospi[Gd], mospi[Gc], alpha, 
 
65
                    &(C[Gc][0][0]), mospi[Gc], &(X[0][0]), sospi[Gd],
 
66
                    beta, &(Buf2->matrix[h][ij][pq]), sospi[Gd]);
 
67
          }
 
68
        }
 
69
        else {
 
70
          X = block_matrix(sospi[Gc],mospi[Gd]);
 
71
 
 
72
          for(ij=0; ij < Buf1->params->rowtot[h]; ij++) {
 
73
 
 
74
            C_DGEMM('n','n', sospi[Gc], mospi[Gd], sospi[Gd], 1.0,
 
75
                    &(Buf2->matrix[h][ij][pq]), sospi[Gd], &(C[Gd][0][0]), mospi[Gd],
 
76
                    0.0, &(X[0][0]), mospi[Gd]);
 
77
 
 
78
            C_DGEMM('t','n', mospi[Gc], mospi[Gd], sospi[Gc], alpha, 
 
79
                    &(C[Gc][0][0]), mospi[Gc], &(X[0][0]), mospi[Gd],
 
80
                    beta, &(Buf1->matrix[h][ij][cd]), mospi[Gd]);
 
81
 
 
82
          }
 
83
        }
 
84
 
 
85
        free_block(X);
 
86
      }
 
87
    }
 
88
 
 
89
    dpd_set_default(dpdnum1);
 
90
    if(type==1) dpd_buf4_mat_irrep_wrt(Buf1, h);
 
91
    dpd_buf4_mat_irrep_close(Buf1, h);
 
92
 
 
93
    dpd_set_default(dpdnum2);
 
94
    if(type==0) dpd_buf4_mat_irrep_wrt(Buf2, h);
 
95
    dpd_buf4_mat_irrep_close(Buf2, h);
 
96
 
 
97
  }
 
98
 
 
99
}