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