Psi4
tensors.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2017 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This program is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License along
22  * with this program; if not, write to the Free Software Foundation, Inc.,
23  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24  *
25  * @END LICENSE
26  */
27 
28 #ifndef _dfocc_tensors_h_
29 #define _dfocc_tensors_h_
30 
31 #include "psi4/libpsio/psio.hpp"
32 
33 #define index2(i,j) ((i>j) ? ((i*(i+1)/2)+j) : ((j*(j+1)/2)+i))
34 #define index4(i,j,k,l) index2(index2(i,j),index2(k,l))
35 #define MIN0(a,b) (((a)<(b)) ? (a) : (b))
36 #define MAX0(a,b) (((a)>(b)) ? (a) : (b))
37 #define idx_asym(i,j) ((i>j) ? ((i*(i-1)/2)+j) : ((j*(j-1)/2)+i))
38 
39 
40 using namespace psi;
41 using namespace std;
42 
43 namespace psi{ namespace dfoccwave{
44 
45 class Tensor1d;
46 class Tensor2d;
47 class Tensor3d;
48 class Tensor1i;
49 class Tensor2i;
50 class Tensor3i;
51 
52 typedef std::shared_ptr<Tensor1d> SharedTensor1d;
53 typedef std::shared_ptr<Tensor2d> SharedTensor2d;
54 typedef std::shared_ptr<Tensor3d> SharedTensor3d;
55 typedef std::shared_ptr<Tensor1i> SharedTensor1i;
56 typedef std::shared_ptr<Tensor2i> SharedTensor2i;
57 typedef std::shared_ptr<Tensor3i> SharedTensor3i;
58 
59 class Tensor1d
60 {
61 
62  private:
63  double *A1d_;
64  int dim1_;
65  string name_; // Name of the array
66 
67  public:
68  Tensor1d(int d1);
69  Tensor1d(string name, int d1);
70  Tensor1d(); //default constructer
71  ~Tensor1d(); //destructer
72 
73  void init(string name, int d1);
74  void init(int d1);
75  void memalloc();
76  void zero();
77  void print();
78  void print(std::string OutFileRMR);
79  void release();
80  void set(int i, double value);
81  void set(double *vec);
82  void set(const SharedTensor1d &vec);
83  void add(const SharedTensor1d& Adum);
84  void add(int i, double value);// add value to ith element of the vector
85  void subtract(const SharedTensor1d& Adum);
86  void subtract(int i, double value);
87  double get(int i);
88  void to_shared_vector(SharedVector A);
89  // rms: rms of A1d_
90  double rms();
91  // rms: rms of (A1d_ - Atemp)
92  double rms(const SharedTensor1d& Atemp);
93  // dot: return result of A1d_' * y
94  double dot(const SharedTensor1d &y);
95  // gemv: A1d_ = alpha * A * b + beta, where A is a general matrix
96  // gemv: C(m) = \sum_{n} A(m,n) b(n)
97  void gemv(bool transa, const SharedTensor2d& a, const SharedTensor1d& b, double alpha, double beta);
98  // gemv: C(m) = \sum_{n} A(m,n) b(n)
99  void gemv(bool transa, int m, int n, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
100  void gemv(bool transa, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
101  void gemv(bool transa, int m, int n, const SharedTensor2d& a, const SharedTensor2d& b, int start_a, int start_b, double alpha, double beta);
102  void gemv(bool transa, int m, int n, const SharedTensor2d& a, const SharedTensor2d& b, int start_a, int start_b, int start_c, double alpha, double beta);
103  // gbmv: This function may NOT working correctly!!!!
104  void gbmv(bool transa, const SharedTensor2d& a, const SharedTensor1d& b, double alpha, double beta);
105  // xay: return result of A1d_' * A * y
106  double xay(const SharedTensor2d &a, const SharedTensor1d &y);
107  // axpy: Y <-- a * X + Y
108  void axpy(const SharedTensor1d &a, double alpha);
109  void scale(double a);
110  void copy(double *x);
111  void copy(const SharedTensor1d &x);
112  // row_vector: set A1d to nth row of A, dim1_ = A->dim2
113  void row_vector(SharedTensor2d &A, int n);
114  // column_vector: set A1d to nth column of A, dim1_ = A->dim1
115  void column_vector(SharedTensor2d &A, int n);
116  int dim1() const { return dim1_; }
117  // dirprd: A1d_[i] = a[i] * b[i]
118  void dirprd(SharedTensor1d &a, SharedTensor1d &b);
119  // symm_packed A(p>=q) = A(p,q) * (2 - \delta_{pq})
120  void symm_packed(const SharedTensor2d &A);
121  // ltm: lower triangular matrix: A(p>=q) = A(p,q)
122  void ltm(const SharedTensor2d &A);
123 
124  friend class Tensor2d;
125  friend class Tensor3d;
126 };
127 
128 
129 class Tensor2d
130 {
131 
132  private:
133  double **A2d_;
134  int dim1_, dim2_, d1_, d2_, d3_, d4_;
135  int **row_idx_, **col_idx_;
136  int *row2d1_, *row2d2_, *col2d1_, *col2d2_;
137  string name_; // Name of the array
138 
139  public:
140  Tensor2d(int d1,int d2);
141  Tensor2d(string name, int d1,int d2);
142  Tensor2d(psi::PSIO* psio, unsigned int fileno, string name, int d1,int d2);
143  Tensor2d(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, string name, int d1,int d2);
144  Tensor2d(psi::PSIO& psio, unsigned int fileno, string name, int d1,int d2);
145  Tensor2d(string name, int d1, int d2, int d3, int d4);
146  Tensor2d(string name, int d1, int d2, int d3);
147  Tensor2d(); //default constructer
148  ~Tensor2d(); //destructer
149 
150  void init(string name, int d1,int d2);
151  void init(int d1,int d2);
152  void memalloc();
153  void zero();
154  void zero_diagonal();
155  void print();
156  void print(std::string OutFileRMR);
157  void release();
158  void set(int i, int j, double value);
159  void set(double **A);
160  void set(double *A);
161  void set(SharedTensor2d &A);
162  void set(SharedMatrix A);
163  void set2(SharedMatrix A);
164  void set(SharedTensor1d &A);
165  // A2d_[n][ij] = A(i,j)
166  void set_row(const SharedTensor2d &A, int n);
167  // A2d_[ij][n] = A(i,j)
168  void set_column(const SharedTensor2d &A, int n);
169  double get(int i, int j);
170  // A2d_[ij] = A(n, ij)
171  void get_row(const SharedTensor2d &A, int n);
172  // A2d_[ij] = A(ij, n)
173  void get_column(const SharedTensor2d &A, int n);
174  // A2d = alpha * Adum
175  void add(const SharedTensor2d &a);
176  void add(double **a);
177  void add(double alpha, const SharedTensor2d &a);
178  void add(int i, int j, double value);
179  // A2d_[n][ij] += A(i,j)
180  void add2row(const SharedTensor2d &A, int n);
181  // A2d_[ij][n] += A(i,j)
182  void add2col(const SharedTensor2d &A, int n);
183  void subtract(const SharedTensor2d &a);
184  void subtract(int i, int j, double value);
185  // axpy: Y <-- a * X + Y
186  void axpy(double **a, double alpha);
187  void axpy(const SharedTensor2d &a, double alpha);
188  void axpy(ULI length, int inc_a, const SharedTensor2d &a, int inc_2d, double alpha);
189  void axpy(ULI length, int start_a, int inc_a, const SharedTensor2d &A, int start_2d, int inc_2d, double alpha);
190  double **transpose2();
191  SharedTensor2d transpose();
192  void trans(const SharedTensor2d &A);
193  void trans(double **A);
194  void copy(double **a);
195  void copy(const SharedTensor2d &Adum);
196  void copy(ULI length, const SharedTensor2d &A, int inc_a, int inc_2d);
197  void copy(const SharedTensor2d &A, int start);
198  // partial copy
199  void pcopy(const SharedTensor2d &A, int dim_copy, int dim_skip);
200  void pcopy(const SharedTensor2d &A, int dim_copy, int dim_skip, int start);
201  double get_max_element();
202  // diagonalize: diagonalize via rsp
203  void diagonalize(const SharedTensor2d &eigvectors, const SharedTensor1d &eigvalues, double cutoff);
204  // cdsyev: diagonalize via lapack
205  void cdsyev(char jobz, char uplo, const SharedTensor2d& eigvectors, const SharedTensor1d& eigvalues);
206  // davidson: diagonalize via davidson algorithm
207  void davidson(int n_eigval, const SharedTensor2d& eigvectors, const SharedTensor1d& eigvalues, double cutoff, int print);
208  // cdgesv: solve a linear equation via lapack
209  void cdgesv(const SharedTensor1d& Xvec);
210  void cdgesv(double* Xvec);
211  void cdgesv(const SharedTensor1d& Xvec, int errcod);
212  void cdgesv(double* Xvec, int errcod);
213  // lineq_flin: solve a linear equation via FLIN
214  void lineq_flin(const SharedTensor1d& Xvec, double *det);
215  void lineq_flin(double* Xvec, double *det);
216  // pople: solve a linear equation via Pople's algorithm
217  void lineq_pople(const SharedTensor1d& Xvec, int num_vecs, double cutoff);
218  void lineq_pople(double* Xvec, int num_vecs, double cutoff);
219 
220  // gemm: matrix multiplication C = A * B
221  void gemm(bool transa, bool transb, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
222  // contract: general contraction C(m,n) = \sum_{k} A(m,k) * B(k,n)
223  void contract(bool transa, bool transb, int m, int n, int k, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
224  void contract(bool transa, bool transb, int m, int n, int k, const SharedTensor2d& a, const SharedTensor2d& b, int start_a, int start_b, double alpha, double beta);
225  void contract(bool transa, bool transb, int m, int n, int k, const SharedTensor2d& a, const SharedTensor2d& b,
226  int start_a, int start_b, int start_c, double alpha, double beta);
227  // contract323: C[Q](m,n) = \sum_{k} A[Q](m,k) * B(k,n). Note: contract332 should be called with beta=1.0
228  void contract323(bool transa, bool transb, int m, int n, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
229  // contract233: C[Q](m,n) = \sum_{k} A(m,k) * B[Q](k,n)
230  void contract233(bool transa, bool transb, int m, int n, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
231  // contract332: C(m,n) = \sum_{k} A[Q](m,k) * B[Q](k,n)
232  void contract332(bool transa, bool transb, int k, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
233  // contract424: Z(pq,rs) = \sum_{o} X(pq,ro) * Y(o,s): where target_x = 4, target_y = 1, and Z = A2d_
234  void contract424(int target_x, int target_y, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
235  // contract442: C(p,q) \sum_{rst} A(pr,st) B(qr,st) , where row/col pair indices are related to B and sorted B.
236  void contract442(int target_a, int target_b, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
237  // gemv: C(mn) = \sum_{k} A(mn,k) * b(k)
238  void gemv(bool transa, const SharedTensor2d& a, const SharedTensor1d& b, double alpha, double beta);
239  // gemv: C(mn) = \sum_{kl} A(mn,kl) * b(kl)
240  void gemv(bool transa, const SharedTensor2d& a, const SharedTensor2d& b, double alpha, double beta);
241 
242  // level_shift: A[i][i] = A[i][i] - value
243  void level_shift(double value);
244  // outer_product: A = x * y'
245  void outer_product(const SharedTensor1d &x, const SharedTensor1d &y);
246  void scale(double a);
247  // scale_row: scales mth row with a
248  void scale_row(int m, double a);
249  // scale_column: scales nth column with a
250  void scale_column(int n, double a);
251  // identity: A = I
252  void identity();
253  double trace();
254  double norm();
255  // transform: A = L' * B * L
256  void transform(const SharedTensor2d& a, const SharedTensor2d& transformer);
257  // back_transform: A = L * B * L'
258  void back_transform(const SharedTensor2d& a, const SharedTensor2d& transformer);
259  // back_transform: A = alpha* L * B * L' + (beta*A)
260  void back_transform(const SharedTensor2d& a, const SharedTensor2d& transformer, double alpha, double beta);
261  // pseudo_transform: A = L * B * L
262  void pseudo_transform(const SharedTensor2d& a, const SharedTensor2d& transformer);
263  // triple_gemm: A2d_ = a * b * c
264  void triple_gemm(const SharedTensor2d& a, const SharedTensor2d& b, const SharedTensor2d& c);
265  // vector_dot: value = Tr(A' * B)
266  double vector_dot(const SharedTensor2d &rhs);
267  double vector_dot(double **rhs);
268  double **to_block_matrix();
269  double *to_lower_triangle();
270  void to_shared_matrix(SharedMatrix A);
271  void to_matrix(SharedMatrix A);
272  void to_pointer(double *A);
273  // mgs: orthogonalize with a Modified Gram-Schmid algorithm
274  void mgs();
275  // gs: orthogonalize with a Classical Gram-Schmid algorithm
276  void gs();
277  // row_vector: return nth row as a vector
278  double *row_vector(int n);
279  // column_vector: return nth column as a vector
280  double *column_vector(int n);
281  int dim1() const { return dim1_; }
282  int dim2() const { return dim2_; }
283 
284  void write(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
285  void write(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, psio_address start, psio_address *end);
286  void write(psi::PSIO* const psio, unsigned int fileno);
287  void write(psi::PSIO* psio, unsigned int fileno, psio_address start, psio_address *end);
288  void write(psi::PSIO& psio, unsigned int fileno);
289  void write(psi::PSIO& psio, unsigned int fileno, psio_address start, psio_address *end);
290  void write(std::shared_ptr<psi::PSIO> psio, const string& filename, unsigned int fileno);
291  void write(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, bool three_index, bool symm);
292  void write(std::shared_ptr<psi::PSIO> psio, const string& filename, unsigned int fileno, bool three_index, bool symm);
293  void write_symm(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
294  void write_anti_symm(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
295 
296  void read(psi::PSIO* psio, unsigned int fileno);
297  void read(psi::PSIO* psio, unsigned int fileno, psio_address start, psio_address *end);
298  void read(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
299  void read(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, psio_address start, psio_address *end);
300  void read(psi::PSIO& psio, unsigned int fileno);
301  void read(psi::PSIO& psio, unsigned int fileno, psio_address start, psio_address *end);
302  void read(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, bool three_index, bool symm);
303  void read_symm(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
304  void read_anti_symm(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
305 
306  bool read(PSIO* psio, int itap, const char *label, int dim);
307  bool read(std::shared_ptr<psi::PSIO> psio, int itap, const char *label, int dim);
308  void save(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
309  void save(psi::PSIO* const psio, unsigned int fileno);
310  void save(psi::PSIO& psio, unsigned int fileno);
311  void load(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, string name, int d1,int d2);
312  void load(psi::PSIO* const psio, unsigned int fileno, string name, int d1,int d2);
313  void load(psi::PSIO& psio, unsigned int fileno, string name, int d1,int d2);
314 
315  void mywrite(const string& filename);
316  void mywrite(int fileno);
317  void mywrite(int fileno, bool append);
318 
319  void myread(const string& filename);
320  void myread(int fileno);
321  void myread(int fileno, bool append);
322  void myread(int fileno, ULI start);
323 
324  // sort (for example 1432 sort): A2d_(ps,rq) = A(pq,rs)
325  // A2d_ = alpha*A + beta*A2d_
326  void sort(int sort_type, const SharedTensor2d &A, double alpha, double beta);
327  // A2d_[p][qr] = sort(A[p][qr])
328  void sort3a(int sort_type, int d1, int d2, int d3, const SharedTensor2d &A, double alpha, double beta);
329  // A2d_[pq][r] = sort(A[pq][r])
330  void sort3b(int sort_type, int d1, int d2, int d3, const SharedTensor2d &A, double alpha, double beta);
331  // apply_denom: T(ij,ab) /= D(ij,ab)
332  void apply_denom(int frzc, int occ, const SharedTensor2d &fock);
333  // apply_denom_os: T(Ij,Ab) /= D(Ij,Ab)
334  void apply_denom_os(int frzc, int occA, int occB, const SharedTensor2d &fockA, const SharedTensor2d &fockB);
335  // apply_denom_chem: T(ia,jb) /= D(ij,ab)
336  void apply_denom_chem(int frzc, int occ, const SharedTensor2d &fock);
337  void reg_denom(int frzc, int occ, const SharedTensor2d &fock, double reg);
338  void reg_denom_os(int frzc, int occA, int occB, const SharedTensor2d &fockA, const SharedTensor2d &fockB, double reg);
339  void reg_denom_chem(int frzc, int occ, const SharedTensor2d &fock, double reg);
340  // dirprd: A2d_[i][j] = a[i][j] * b[i][j]
341  void dirprd(const SharedTensor2d &a, const SharedTensor2d &b);
342  // dirprd123: A2d_[Q][ij] = a[Q] * b[i][j]
343  void dirprd123(const SharedTensor1d &a, const SharedTensor2d &b, double alpha, double beta);
344  void dirprd123(bool transb, const SharedTensor1d &a, const SharedTensor2d &b, double alpha, double beta);
345  // dirprd112: A2d_[i][j] = a[i] * b[j]
346  void dirprd112(const SharedTensor1d &a, const SharedTensor1d &b);
347  // dirprd112: A2d_[i][j] = alpha *a[i] * b[j] + beta * A2d_[i][j]
348  void dirprd112(const SharedTensor1d &a, const SharedTensor1d &b, double alpha, double beta);
349  // dirprd224: A2d_[ij][kl] = a[i][j] * b[k][l]
350  void dirprd224(const SharedTensor2d &a, const SharedTensor2d &b);
351  void dirprd224(const SharedTensor2d &a, const SharedTensor2d &b, double alpha, double beta);
352  double* to_vector(const SharedTensor2i &pair_idx);
353  double* to_vector();
354  double rms();
355  double rms(const SharedTensor2d& a);
356 
357  void set_act_oo(int aocc, const SharedTensor2d &a);
358  void set_act_oo(int frzc, int aocc, const SharedTensor2d &a);
359  void set_act_vv(int occ, int avir, const SharedTensor2d &A);
360  void set_act_vv(const SharedTensor2d &A);
361  void set_oo(const SharedTensor2d &a);
362  void set_ov(const SharedTensor2d &A);
363  void set_vo(const SharedTensor2d &A);
364  void set_vv(int occ, const SharedTensor2d &A);
365  void set_act_ov(int frzc, const SharedTensor2d &A);
366  void set_act_vo(int frzc, const SharedTensor2d &A);
367 
368  void add_oo(const SharedTensor2d &A, double alpha, double beta);
369  void add_vv(int occ, const SharedTensor2d &A, double alpha, double beta);
370  void add_ov(const SharedTensor2d &A, double alpha, double beta);
371  void add_vo(const SharedTensor2d &A, double alpha, double beta);
372  void add_aocc_fc(const SharedTensor2d &A, double alpha, double beta);
373  void add_fc_aocc(const SharedTensor2d &A, double alpha, double beta);
374 
375  void add3_oo(const SharedTensor2d &A, double alpha, double beta);
376  void set3_oo(const SharedTensor2d &A);
377  void set3_ov(const SharedTensor2d &A);
378  void set3_vo(const SharedTensor2d &A);
379  void set3_vv(const SharedTensor2d &A, int occ);
380 
381  void set3_act_ov(int frzc, int aocc, int avir, int vir, const SharedTensor2d &a);
382  void set3_act_oo(int frzc, const SharedTensor2d &A);
383  void set3_act_vv(const SharedTensor2d &A);
384 
385  void swap_3index_col(const SharedTensor2d &A);
386 
387  void form_oo(const SharedTensor2d &A);
388  void form_act_oo(int frzc, const SharedTensor2d &A);
389  void form_vv(int occ, const SharedTensor2d &A);
390  void form_act_vv(int occ, const SharedTensor2d &A);
391  void form_vo(const SharedTensor2d &A);
392  void form_vo(int occ, const SharedTensor2d &A);
393  void form_act_vo(int frzc, const SharedTensor2d &A);
394  void form_act_vo(int frzc, int occ, const SharedTensor2d &A);
395  void form_ov(const SharedTensor2d &A);
396  void form_ov(int occ, const SharedTensor2d &A);
397  void form_act_ov(int frzc, const SharedTensor2d &A);
398  void form_act_ov(int frzc, int occ, const SharedTensor2d &A);
399  void form_ooAB(const SharedTensor2d &A);
400 
401  void form_b_ij(int frzc, const SharedTensor2d &A);
402  void form_b_ia(int frzc, const SharedTensor2d &A);
403  void form_b_ab(const SharedTensor2d &A);
404  // form_b_kl: k is active occupied, and l is frozen core
405  void form_b_kl(const SharedTensor2d &A);
406  // form_b_ki: k is active occupied, and i is all occupied
407  void form_b_ki(const SharedTensor2d &A);
408  // form_b_li: l is frozen core, and i is all occupied
409  void form_b_li(const SharedTensor2d &A);
410  // form_b_il: l is frozen core, and i is all occupied
411  void form_b_il(const SharedTensor2d &A);
412  // form_b_ka: k is active occupied, and a is all virtual
413  void form_b_ka(const SharedTensor2d &A);
414  // form_b_la: l is frozen core, and a is all virtual
415  void form_b_la(const SharedTensor2d &A);
416 
417  // B_pq = 1/2 (A_pq + A_qp)
418  void symmetrize();
419  void symmetrize(const SharedTensor2d &A);
420  // B(Q,pq) = 1/2 [ A(Q,pq) + A(Q,qp) ]
421  void symmetrize3(const SharedTensor2d &A);
422  // A(Q, p>=q) = A(Q,pq) * (2 -\delta_{pq})
423  void symm_packed(const SharedTensor2d &A);
424  // A(Q, p>=q) = A(Q,pq)
425  void ltm(const SharedTensor2d &A);
426  // A(p,qr) = A(p,q>=r)
427  void expand23(int d1, int d2, int d3, const SharedTensor2d &A);
428 
429  // (+)A(p>=q, r>=s) = 1/2 [A(pq,rs) + A(qp,rs)]
430  void symm4(const SharedTensor2d &a);
431  // (-)A(p>=q, r>=s) = 1/2 [A(pq,rs) - A(qp,rs)]
432  void antisymm4(const SharedTensor2d &a);
433  // (+)A(p>=q, r>=s) = 1/2 [A(pq,rs) + A(pq,sr)]
434  void symm_col4(const SharedTensor2d &a);
435  // (-)A(p>=q, r>=s) = 1/2 [A(pq,rs) - A(pq,sr)]
436  void antisymm_col4(const SharedTensor2d &a);
437  // (+)At(p>=q, r>=s) = 1/2 [A(pq,rs) + A(qp,rs)] * (2 -\delta_{pq})
438  void symm_row_packed4(const SharedTensor2d &a);
439  // (+)At(p>=q, r>=s) = 1/2 [A(pq,rs) + A(qp,rs)] * (2 -\delta_{rs})
440  void symm_col_packed4(const SharedTensor2d &a);
441  // (-)At(p>=q, r>=s) = 1/2 [A(pq,rs) - A(qp,rs)] * (2 -\delta_{pq})
442  void antisymm_row_packed4(const SharedTensor2d &a);
443  // (-)At(p>=q, r>=s) = 1/2 [A(pq,rs) - A(qp,rs)] * (2 -\delta_{rs})
444  void antisymm_col_packed4(const SharedTensor2d &a);
445 
446  // A2d_(pq,rs) = 2 <pq|rs> - <pq|sr>
447  void tei_cs1_anti_symm(const SharedTensor2d &J, const SharedTensor2d &K);
448  // A2d_(pq,rs) = 2 <pq|rs> - <qp|rs>
449  void tei_cs2_anti_symm(const SharedTensor2d &J, const SharedTensor2d &K);
450  // A2d_(pq,rs) = 2 (pq|rs) - (ps|rq)
451  void tei_cs3_anti_symm(const SharedTensor2d &J, const SharedTensor2d &K);
452  // A2d_(pq,rs) = 2 (pq|rs) - (rq|ps)
453  void tei_cs4_anti_symm(const SharedTensor2d &J, const SharedTensor2d &K);
454 
455  // A2d(ij,ab) += P_(ij) * P_(ab) A(ia,jb)
456  void P_ijab(const SharedTensor2d &A);
457 
458 
459  // General tensor contractions over
460  // C(pq,rs) = \sum_{tu} A(pq,tu) B(tu,rs)
461  // t_a1: t; t_a2: u; f_a1: p; f_a2: q
462  void cont444(int t_a1, int t_a2, int f_a1, int f_a2, const SharedTensor2d& A, int t_b1, int t_b2, int f_b1, int f_b2, const SharedTensor2d& B, double alpha, double beta);
463  void cont444(bool delete_a, int t_a1, int t_a2, int f_a1, int f_a2, SharedTensor2d& A,
464  bool delete_b, int t_b1, int t_b2, int f_b1, int f_b2, SharedTensor2d& B,
465  double alpha, double beta);
466  void cont444(string idx_c, string idx_a, string idx_b, bool delete_a, bool delete_b, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
467  // C(pq) = \sum_{rst} A(pr,st) B(rs,tq)
468  void cont442(string idx_c, string idx_a, string idx_b, bool delete_a, bool delete_b, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
469  // C(pq,rs) = \sum_{t} A(pq,rt) B(t,s)
470  void cont424(string idx_c, string idx_a, string idx_b, bool delete_a, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
471  // C(pq,rs) = \sum_{t} A(p,t) B(tq,rs)
472  void cont244(string idx_c, string idx_a, string idx_b, bool delete_b, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
473  // C(Q,pq) = \sum_{rs} A(Q,rs) B(rs,pq)
474  // where dim(idx_c) & dim(idx_a)=2 but dim(idx_b)=4
475  void cont343(string idx_c, string idx_a, string idx_b, bool delete_b, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
476  // C(Q,pq) = \sum_{r} A(p,r) B(Q,rq)
477  void cont233(string idx_c, string idx_a, string idx_b, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
478  // C(Q,pq) = \sum_{r} A(Q,pr) B(r,q)
479  void cont323(string idx_c, string idx_a, string idx_b, bool delete_a, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
480  // C(pq) = \sum_{Qr} A(Q,rp) B(Q,rq)
481  void cont332(string idx_c, string idx_a, string idx_b, bool delete_a, bool delete_b, SharedTensor2d& A, SharedTensor2d& B, double alpha, double beta);
482 
483  friend class Tensor1d;
484  friend class Tensor3d;
485  friend class Tensor1i;
486  friend class Tensor2i;
487 };
488 
489 class Tensor3d
490 {
491 
492  private:
493  double ***A3d_;
494  int dim1_,dim2_,dim3_;
495  string name_; // Name of the array
496 
497  public:
498  Tensor3d(int d1,int d2, int d3);
499  Tensor3d(string name, int d1,int d2, int d3);
500  Tensor3d(); //default constructer
501  ~Tensor3d(); //destructer
502 
503  void init(string name, int d1,int d2, int d3);
504  void init(int d1,int d2, int d3);
505  void memalloc();
506  void zero();
507  void print();
508  void release();
509  void set(int h, int i, int j, double value);
510  double get(int h, int i, int j);
511 
512  friend class Tensor1d;
513  friend class Tensor2d;
514 };
515 
516 class Tensor1i
517 {
518 
519  private:
520  int *A1i_;
521  int dim1_;
522  string name_; // Name of the array
523 
524  public:
525  Tensor1i(int d1);
526  Tensor1i(string name, int d1);
527  Tensor1i(); //default constructer
528  ~Tensor1i(); //destructer
529 
530  void init(string name, int d1);
531  void init(int d1);
532  void memalloc();
533  void zero();
534  void print();
535  void release();
536  void set(int i, int value);
537  int get(int i);
538  void add(const SharedTensor1i& Adum);
539  void add(int i, int value);
540  void subtract(const SharedTensor1i& Adum);
541  void subtract(int i, int value);
542 
543 };
544 
545 class Tensor2i
546 {
547 
548  private:
549  int **A2i_;
550  int dim1_,dim2_;
551  string name_; // Name of the array
552 
553  public:
554  Tensor2i(int d1,int d2);
555  Tensor2i(string name, int d1,int d2);
556  Tensor2i(); //default constructer
557  ~Tensor2i(); //destructer
558 
559  void init(string name, int d1,int d2);
560  void init(int d1,int d2);
561  void memalloc();
562  void zero();
563  void zero_diagonal();
564  void print();
565  void print(std::string OutFileRMR);
566  void release();
567  void set(int i, int j, int value);
568  void set(int **A);
569  double get(int i, int j);
570  void add(const SharedTensor2i& Adum);
571  void add(int i, int j, int value);
572  void subtract(const SharedTensor2i& Adum);
573  void subtract(int i, int j, int value);
574  SharedTensor2i transpose();
575  void copy(const SharedTensor2i& Adum);
576  void copy(int **a);
577  void identity();
578  int trace();
579  int **to_int_matrix();
580  int dim1() const { return dim1_; }
581  int dim2() const { return dim2_; }
582 
583  friend class Tensor1i;
584  friend class Tensor3i;
585  friend class Tensor1d;
586  friend class Tensor2d;
587 };
588 
589 
590 class Tensor3i
591 {
592 
593  private:
594  int ***A3i_;
595  int dim1_,dim2_,dim3_;
596  string name_; // Name of the array
597 
598  public:
599  Tensor3i(int d1,int d2, int d3);
600  Tensor3i(string name, int d1,int d2, int d3);
601  Tensor3i(); //default constructer
602  ~Tensor3i(); //destructer
603 
604  void init(string name, int d1,int d2, int d3);
605  void init(int d1,int d2, int d3);
606  void memalloc();
607  void zero();
608  void print();
609  void release();
610  void set(int h, int i, int j, int value);
611  int get(int h, int i, int j);
612 };
613 }} // End Namespaces
614 #endif // _dfocc_tensors_h_
double * A1d_
Definition: tensors.h:63
double dot(BlockMatrix *A, BlockMatrix *B)
Definition: mcscf/block_matrix.cc:223
std::shared_ptr< Tensor2i > SharedTensor2i
Definition: tensors.h:56
double *** A3d_
Definition: tensors.h:493
string name_
Definition: tensors.h:65
double ** A2d_
Definition: tensors.h:133
Definition: tensors.h:129
string name_
Definition: tensors.h:495
Definition: tensors.h:590
std::shared_ptr< Tensor1i > SharedTensor1i
Definition: tensors.h:55
Definition: tensors.h:59
Definition: matrixtmp.h:41
Definition: pointgrp.h:105
Definition: libpsio/config.h:64
int dim3_
Definition: tensors.h:595
int *** A3i_
Definition: tensors.h:594
string name_
Definition: tensors.h:596
unsigned long int ULI
Definition: libpsio/config.h:62
std::shared_ptr< Tensor2d > SharedTensor2d
Definition: tensors.h:53
int dim2() const
Definition: tensors.h:581
void sort(double *A, double **B, int n)
Definition: libqt/sort.cc:49
int * row2d2_
Definition: tensors.h:136
std::shared_ptr< Tensor1d > SharedTensor1d
Definition: tensors.h:50
std::shared_ptr< Tensor3i > SharedTensor3i
Definition: tensors.h:57
int dim1() const
Definition: tensors.h:116
int dim1_
Definition: tensors.h:521
int dim2_
Definition: tensors.h:550
int * A1i_
Definition: tensors.h:520
int dim2_
Definition: tensors.h:134
int ** row_idx_
Definition: tensors.h:135
Definition: tensors.h:545
int dim1() const
Definition: tensors.h:580
Definition: tensors.h:489
Definition: tensors.h:516
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:50
CCTransform * trans
Definition: psimrcc/main.cc:79
int dim2() const
Definition: tensors.h:282
SharedWavefunction dfoccwave(SharedWavefunction, Options &)
Definition: dfocc/main.cc:35
int dim1() const
Definition: tensors.h:281
string name_
Definition: tensors.h:551
std::shared_ptr< Tensor3d > SharedTensor3d
Definition: tensors.h:54
string name_
Definition: tensors.h:522
Definition: PsiFileImpl.h:38
int dim1_
Definition: tensors.h:64
std::shared_ptr< Vector > SharedVector
Definition: adc.h:52
int dim3_
Definition: tensors.h:494
Definition: psio.hpp:196
string name_
Definition: tensors.h:137
int ** A2i_
Definition: tensors.h:549