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