← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Tools/transform.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:
src/bin/cints/Tools/transform.cc
 
1
/*! \file transform.cc
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
    -------------------------------------------------------------------------------------------------
 
5
    Functions for performing four-index transformation from cartesian to pure angular momentum basis
 
6
    -------------------------------------------------------------------------------------------------*/
 
7
 
 
8
#include<cstdio>
 
9
#include<memory.h>
 
10
#include<libciomr/libciomr.h>
 
11
#include<libqt/qt.h>
 
12
#include<libint/libint.h>
 
13
 
 
14
#include"defines.h"
 
15
#define EXTERN
 
16
#include"global.h"
 
17
#include <stdexcept>
 
18
 
 
19
namespace psi { namespace CINTS {
 
20
 
 
21
#if !USE_MM
 
22
void transform_i(double *data, double *puream_data, double **c2p, int am_i, int nj, int nk, int nl)
 
23
{
 
24
  register int i,j,k,l;
 
25
  int step = nj*nk*nl;
 
26
  int offset;
 
27
  double *ii;
 
28
  int ni, nsph, isph;
 
29
  double tmp;
 
30
 
 
31
  ni = ioff[am_i + 1];
 
32
  nsph = 2*am_i + 1;
 
33
  
 
34
  for(j=0;j<nj;j++)
 
35
    for(k=0;k<nk;k++) {
 
36
      offset = (j*nk+k)*nl;
 
37
      for(l=0;l<nl;l++,offset++)
 
38
        for(isph=0;isph<nsph;isph++) {
 
39
          ii = data + offset; tmp = 0.0;
 
40
          for(i=0;i<ni;i++,ii+=step)
 
41
            tmp += c2p[isph][i]*(*ii);
 
42
          puream_data[step*isph+offset] = tmp;
 
43
        }
 
44
    }
 
45
 
 
46
  return;
 
47
}
 
48
          
 
49
 
 
50
void transform_j(double *data, double *puream_data, double **c2p, int am_j, int ni, int nk, int nl)
 
51
{
 
52
  register int i,j,k,l;
 
53
  int step = nk*nl;
 
54
  int offset,offset_sph;
 
55
  double *jj;
 
56
  int nj, nsph, jsph;
 
57
  double tmp;
 
58
 
 
59
  nj = ioff[am_j + 1];
 
60
  nsph = 2*am_j + 1;
 
61
  
 
62
  for(i=0;i<ni;i++)
 
63
    for(k=0;k<nk;k++,offset+=nl) {
 
64
      offset = (i*nj*nk + k)*nl;
 
65
      offset_sph = (i*nsph*nk + k)*nl;
 
66
      for(l=0;l<nl;l++,offset++,offset_sph++)
 
67
        for(jsph=0;jsph<nsph;jsph++) {
 
68
          jj = data + offset; tmp = 0.0;
 
69
          for(j=0;j<nj;j++,jj+=step)
 
70
            tmp += c2p[jsph][j]*(*jj);
 
71
          puream_data[step*jsph+offset_sph] = tmp;
 
72
        }
 
73
    }
 
74
 
 
75
  return;
 
76
}
 
77
 
 
78
 
 
79
 
 
80
void transform_k(double *data, double *puream_data, double **c2p, int am_k, int ni, int nj, int nl)
 
81
{
 
82
  register int i,j,k,l;
 
83
  int step = nl;
 
84
  int offset,offset_sph;
 
85
  double *kk;
 
86
  int nk, nsph, ksph;
 
87
  double tmp;
 
88
 
 
89
  nk = ioff[am_k + 1];
 
90
  nsph = 2*am_k + 1;
 
91
  
 
92
  for(i=0;i<ni;i++)
 
93
    for(j=0;j<nj;j++) {
 
94
      offset = (i*nj + j)*nk*nl;
 
95
      offset_sph = (i*nj + j)*nsph*nl;
 
96
      for(l=0;l<nl;l++,offset++,offset_sph++)
 
97
        for(ksph=0;ksph<nsph;ksph++) {
 
98
          kk = data + offset; tmp = 0.0;
 
99
          for(k=0;k<nk;k++,kk+=step)
 
100
            tmp += c2p[ksph][k]*(*kk);
 
101
          puream_data[step*ksph+offset_sph] = tmp;
 
102
        }
 
103
    }
 
104
 
 
105
  return;
 
106
}
 
107
 
 
108
 
 
109
 
 
110
void transform_l(double *data, double *puream_data, double **c2p, int am_l, int ni, int nj, int nk)
 
111
{
 
112
  register int i,j,k,l;
 
113
  int step = 1;
 
114
  int offset,offset_sph;
 
115
  double *ll;
 
116
  int nl, nsph, lsph;
 
117
  double tmp;
 
118
 
 
119
  nl = ioff[am_l + 1];
 
120
  nsph = 2*am_l + 1;
 
121
  
 
122
  for(i=0;i<ni;i++)
 
123
    for(j=0;j<nj;j++) {
 
124
      offset = (i*nj + j)*nk*nl;
 
125
      offset_sph = (i*nj + j)*nk*nsph;
 
126
      for(k=0;k<nk;k++,offset+=nl,offset_sph+=nsph)
 
127
        for(lsph=0;lsph<nsph;lsph++) {
 
128
          ll = data + offset; tmp = 0;
 
129
          for(l=0;l<nl;l++,ll+=step)
 
130
            tmp += c2p[lsph][l]*(*ll);
 
131
          puream_data[step*lsph+offset_sph] = tmp;
 
132
        }
 
133
    }
 
134
 
 
135
}
 
136
 
 
137
#elif USE_MM && !SPARSE_C2P
 
138
 
 
139
double *transform_ijkl(double *data, double *puream_data, double **cc2pp_ij, double **cc2pp_kl, int trans_ij,
 
140
                       int am_i, int am_j, int am_k, int am_l)
 
141
{
 
142
  int ij,kl,ijkl;
 
143
  int n_ij, n_kl, nsph_ij, nsph_kl;
 
144
  double tmp;
 
145
#if USE_BLAS
 
146
  int zero = 0;
 
147
  int one = 1;
 
148
  double zero_d = 0.0;
 
149
  double one_d = 1.0;
 
150
  char true_c = 't';
 
151
  char false_c = 'n';
 
152
#else
 
153
  double *tmp_ptr;
 
154
  static double *cij_ckl[CINTS_MAX_AM*CINTS_MAX_AM*(CINTS_MAX_AM+1)*(CINTS_MAX_AM+1)/4];
 
155
  static double *cij_pkl[CINTS_MAX_AM*(CINTS_MAX_AM+1)*(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)/2];
 
156
  static double *pij_pkl[(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)];
 
157
#endif
 
158
 
 
159
  n_ij = ioff[am_i+1]*ioff[am_j+1];
 
160
  n_kl = ioff[am_k+1]*ioff[am_l+1];
 
161
#if !USE_BLAS
 
162
  for(ij=0, tmp_ptr = data; ij<n_ij; ij++, tmp_ptr += n_kl)
 
163
    cij_ckl[ij] = tmp_ptr;
 
164
#endif
 
165
  
 
166
  nsph_kl = (2*am_k+1)*(2*am_l+1);
 
167
#if !USE_BLAS
 
168
  for(ij=0, tmp_ptr = puream_data; ij<n_ij; ij++, tmp_ptr += nsph_kl)
 
169
    cij_pkl[ij] = tmp_ptr;
 
170
  /* Insert any C-based matrix multiplication routine here */
 
171
  newmm_rking(&cij_ckl[0],0,cc2pp_kl,1,&cij_pkl[0],n_ij,n_kl,nsph_kl,1.0,0.0);
 
172
#else
 
173
/*  F_DGEMM(&true_c,&false_c,&nsph_kl,&n_ij,&n_kl,&one_d,&(cc2pp_kl[0][0]),&n_kl,
 
174
          data,&n_kl,&zero_d,puream_data,&nsph_kl);*/
 
175
    C_DGEMM('n', 't', n_ij, nsph_kl, n_kl, 1.0, data, n_kl, cc2pp_kl[0], n_kl,
 
176
            0.0, puream_data, nsph_kl);
 
177
#endif
 
178
 
 
179
  if (trans_ij) {
 
180
    nsph_ij = (2*am_i+1)*(2*am_j+1);
 
181
#if !USE_BLAS
 
182
    for(ij=0, tmp_ptr = data; ij<nsph_ij; ij++, tmp_ptr += nsph_kl)
 
183
      pij_pkl[ij] = tmp_ptr;
 
184
    /* Insert any C-based matrix multiplication routine here */
 
185
    newmm_rking(cc2pp_ij,0,&cij_pkl[0],0,&pij_pkl[0],nsph_ij,n_ij,nsph_kl,1.0,0.0);
 
186
#else
 
187
/*  F_DGEMM(&false_c,&false_c,&nsph_kl,&nsph_ij,&n_ij,&one_d,cc2pp_ij[0],&n_ij,
 
188
          puream_data,&nsph_kl,&zero_d,data,&nsph_kl);*/
 
189
    C_DGEMM('n', 'n', nsph_ij, nsph_kl, n_ij, 1.0, cc2pp_ij[0], n_ij,
 
190
            puream_data, nsph_kl, 0.0, data, nsph_kl);
 
191
 
 
192
#endif
 
193
 
 
194
    return data;
 
195
  }
 
196
  else
 
197
    return puream_data;
 
198
 
 
199
}
 
200
 
 
201
#elif USE_MM && SPARSE_C2P
 
202
 
 
203
double *transform_ijkl_sparse(double *data, double *puream_data,
 
204
                       mat_elem **pp2cc_sparse, mat_elem **cc2pp_sparse,
 
205
                       int *pp2cc_rowlength, int *cc2pp_rowlength, int trans_ij,
 
206
                       int am_i, int am_j, int am_k, int am_l)
 
207
{
 
208
  int ij,kl,ijkl,pq;
 
209
  int n_ij, n_kl, nsph_ij, nsph_kl;
 
210
  double tmp;
 
211
 
 
212
  double *tmp_ptr;
 
213
/*  static double *cij_ckl[CINTS_MAX_AM*CINTS_MAX_AM*(CINTS_MAX_AM+1)*(CINTS_MAX_AM+1)/4];
 
214
  static double *cij_pkl[CINTS_MAX_AM*(CINTS_MAX_AM+1)*(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)/2];
 
215
  static double *pij_pkl[(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)*(2*CINTS_MAX_AM-1)];*/
 
216
  double *cijckl_row, *cijpkl_row, *pijpkl_row;
 
217
 
 
218
  n_ij = ioff[am_i+1]*ioff[am_j+1];
 
219
  n_kl = ioff[am_k+1]*ioff[am_l+1];
 
220
  nsph_kl = (2*am_k+1)*(2*am_l+1);
 
221
  for(ij=0, cijckl_row = data, cijpkl_row = puream_data; ij<n_ij;
 
222
      ij++, cijckl_row += n_kl, cijpkl_row += nsph_kl) {
 
223
    memset(cijpkl_row,0,sizeof(double)*nsph_kl);
 
224
    for(pq=0; pq<n_kl; pq++) {
 
225
      tmp = cijckl_row[pq];
 
226
      for(kl=0; kl<cc2pp_rowlength[pq]; kl++) {
 
227
        cijpkl_row[cc2pp_sparse[pq][kl].column] += tmp*cc2pp_sparse[pq][kl].value;
 
228
      }
 
229
    }
 
230
  }
 
231
                         
 
232
  if (trans_ij) {
 
233
    nsph_ij = (2*am_i+1)*(2*am_j+1);
 
234
    for(ij=0, pijpkl_row = data; ij<nsph_ij; ij++, pijpkl_row += nsph_kl) {
 
235
      memset(pijpkl_row,0,sizeof(double)*nsph_kl);
 
236
      for(pq=0; pq<pp2cc_rowlength[ij]; pq++) {
 
237
        tmp = pp2cc_sparse[ij][pq].value;
 
238
        cijpkl_row = puream_data + nsph_kl*pp2cc_sparse[ij][pq].column;
 
239
        for(kl=0; kl<nsph_kl; kl++) {
 
240
          pijpkl_row[kl] += tmp*cijpkl_row[kl];
 
241
        }
 
242
      }
 
243
    }
 
244
    return data;
 
245
  }
 
246
 
 
247
  return puream_data;
 
248
 
 
249
}
 
250
#endif
 
251
};};