Psi4
dfocc/arrays.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_arrays_h_
30 #define _dfocc_arrays_h_
31 
32 using namespace psi;
33 using namespace std;
34 
35 namespace psi{ namespace dfoccwave{
36 
37 class Array1d;
38 class Array2d;
39 class Array3d;
40 class Array1i;
41 class Array2i;
42 class Array3i;
43 
44 class Array1d
45 {
46 
47  private:
48  double *A1d_;
49  int dim1_;
50  string name_; // Name of the array
51 
52  public:
53  Array1d(int d1);
54  Array1d(string name, int d1);
55  Array1d(); //default constructer
56  ~Array1d(); //destructer
57 
58  Array1d* generate(int d1);
59  Array1d* generate(string name, int d1);
60  void init(string name, int d1);
61  void init(int d1);
62  void memalloc();
63  void zero();
64  void print();
65  void print(std::string OutFileRMR);
66  void release();
67  void set(int i, double value);
68  void set(double *vec);
69  void set(const Array1d *vec);
70  void add(const Array1d* Adum);
71  void add(int i, double value);// add value to ith element of the vector
72  void subtract(const Array1d* Adum);
73  void subtract(int i, double value);
74  double get(int i);
75  // rms: rms of A1d_
76  double rms();
77  // rms: rms of (A1d_ - Atemp)
78  double rms(const Array1d* Atemp);
79  // dot: return result of A1d_' * y
80  double dot(const Array1d *y);
81  // gemv: A1d_ = alpha * A * b + beta, where A is a general matrix
82  void gemv(bool transa, double alpha, const Array2d* a, const Array1d* b, double beta);
83  // gbmv: This function may NOT working correctly!!!!
84  void gbmv(bool transa, double alpha, const Array2d* a, const Array1d* b, double beta);
85  // xay: return result of A1d_' * A * y
86  double xay(const Array2d *a, const Array1d *y);
87  void scale(double a);
88  void copy(double *x);
89  void copy(const Array1d *x);
90  // row_vector: set A1d to nth row of A, dim1_ = A->dim2
91  void row_vector(Array2d *A, int n);
92  // column_vector: set A1d to nth column of A, dim1_ = A->dim1
93  void column_vector(Array2d *A, int n);
94  int dim1() const { return dim1_; }
95  // dirprd: A1d_[i] = a[i] * b[i]
96  void dirprd(Array1d *a, Array1d *b);
97 
98  friend class Array2d;
99  friend class Array3d;
100 };
101 
102 
103 class Array2d
104 {
105 
106  private:
107  double **A2d_;
108  int dim1_,dim2_;
109  string name_; // Name of the array
110 
111  public:
112  Array2d(int d1,int d2);
113  Array2d(string name, int d1,int d2);
114  Array2d(psi::PSIO* psio, unsigned int fileno, string name, int d1,int d2);
115  Array2d(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, string name, int d1,int d2);
116  Array2d(psi::PSIO& psio, unsigned int fileno, string name, int d1,int d2);
117  Array2d(); //default constructer
118  ~Array2d(); //destructer
119 
120  Array2d* generate(int d1,int d2);
121  Array2d* generate(string name, int d1,int d2);
122  void init(string name, int d1,int d2);
123  void init(int d1,int d2);
124  void memalloc();
125  void zero();
126  void zero_diagonal();
127  void print();
128  void print(std::string OutFileRMR);
129  void release();
130  void set(int i, int j, double value);
131  void set(double **A);
132  void set(Array2d *A);
133  void set(SharedMatrix A);
134  double get(int i, int j);
135  void add(const Array2d* Adum);
136  // A2d = alpha * Adum
137  void add(double alpha, const Array2d* Adum);
138  void add(int i, int j, double value);
139  void subtract(const Array2d* Adum);
140  void subtract(int i, int j, double value);
141  Array2d* transpose();
142  void copy(const Array2d* Adum);
143  void copy(double **a);
144  // diagonalize: diagonalize via rsp
145  void diagonalize(Array2d* eigvectors, Array1d* eigvalues, double cutoff);
146  // cdsyev: diagonalize via lapack
147  void cdsyev(char jobz, char uplo, Array2d* eigvectors, Array1d* eigvalues);
148  // davidson: diagonalize via davidson algorithm
149  void davidson(int n_eigval, Array2d* eigvectors, Array1d* eigvalues, double cutoff, int print);
150  // cdgesv: solve a linear equation via lapack
151  void cdgesv(Array1d* Xvec);
152  void cdgesv(double* Xvec);
153  void cdgesv(Array1d* Xvec, int errcod);
154  void cdgesv(double* Xvec, int errcod);
155  // lineq_flin: solve a linear equation via FLIN
156  void lineq_flin(Array1d* Xvec, double *det);
157  void lineq_flin(double* Xvec, double *det);
158  // lineq_pople: solve a linear equation via Pople's algorithm
159  void lineq_pople(Array1d* Xvec, int num_vecs, double cutoff);
160  void lineq_pople(double* Xvec, int num_vecs, double cutoff);
161  // gemm: matrix multiplication C = A * B
162  void gemm(bool transa, bool transb, const Array2d* a, const Array2d* b, double alpha, double beta);
163  // contract: general contraction C(m,n) = \sum_{k} A(m,k) * B(k,n)
164  void contract(bool transa, bool transb, int m, int n, int k, const Array2d* a, const Array2d* b, double alpha, double beta);
165  // contract323: C[Q](m,n) = \sum_{k} A[Q](m,k) * B(k,n)
166  void contract323(bool transa, bool transb, int m, int n, const Array2d* a, const Array2d* b, double alpha, double beta);
167  // contract233: C[Q](m,n) = \sum_{k} A(m,k) * B[Q](k,n)
168  void contract233(bool transa, bool transb, int m, int n, const Array2d* a, const Array2d* b, double alpha, double beta);
169  // level_shift: A[i][i] = A[i][i] - value
170  void level_shift(double value);
171  // outer_product: A = x * y'
172  void outer_product(const Array1d *x, const Array1d *y);
173  void scale(double a);
174  // scale_row: scales mth row with a
175  void scale_row(int m, double a);
176  // scale_column: scales nth column with a
177  void scale_column(int n, double a);
178  // identity: A = I
179  void identity();
180  double trace();
181  // transform: A = L' * B * L
182  void transform(const Array2d* a, const Array2d* transformer);
183  // back_transform: A = L * B * L'
184  void back_transform(const Array2d* a, const Array2d* transformer);
185  // pseudo_transform: A = L * B * L
186  void pseudo_transform(const Array2d* a, const Array2d* transformer);
187  // triple_gemm: A2d_ = a * b * c
188  void triple_gemm(const Array2d* a, const Array2d* b, const Array2d* c);
189  // vector_dot: value = Tr(A' * B)
190  double vector_dot(Array2d *rhs);
191  double vector_dot(double **rhs);
192  double **to_block_matrix();
193  double *to_lower_triangle();
194  void to_shared_matrix(SharedMatrix A);
195  // mgs: orthogonalize with a Modified Gram-Schmid algorithm
196  void mgs();
197  // gs: orthogonalize with a Classical Gram-Schmid algorithm
198  void gs();
199  // row_vector: return nth row as a vector
200  double *row_vector(int n);
201  // column_vector: return nth column as a vector
202  double *column_vector(int n);
203  int dim1() const { return dim1_; }
204  int dim2() const { return dim2_; }
205 
206  void write(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
207  void write(psi::PSIO* const psio, unsigned int fileno);
208  void write(psi::PSIO& psio, unsigned int fileno);
209  void read(psi::PSIO* psio, unsigned int fileno);
210  void read(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
211  void read(psi::PSIO& psio, unsigned int fileno);
212  bool read(PSIO* psio, int itap, const char *label, int dim);
213  bool read(std::shared_ptr<psi::PSIO> psio, int itap, const char *label, int dim);
214  void save(std::shared_ptr<psi::PSIO> psio, unsigned int fileno);
215  void save(psi::PSIO* const psio, unsigned int fileno);
216  void save(psi::PSIO& psio, unsigned int fileno);
217  void load(std::shared_ptr<psi::PSIO> psio, unsigned int fileno, string name, int d1,int d2);
218  void load(psi::PSIO* const psio, unsigned int fileno, string name, int d1,int d2);
219  void load(psi::PSIO& psio, unsigned int fileno, string name, int d1,int d2);
220 
221  // sort1432: A2d_(ps,rq) = A(pq,rs): d1 = num p, d2 = num q, d3 = num r, d4 = num s
222  // col/row_pair_idx are belong to A, while col/row_pairidx2 are belong to A2d_
223  void sort1432(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
224  // sort2134: A2d_(qp,rs) = A(pq,rs): d1 = num p, d2 = num q, d3 = num r, d4 = num s
225  void sort2134(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2);
226  // sort1243: A2d_(pq,sr) = A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
227  void sort1243(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *col_pair_idx2);
228  // sort2413: A2d_(qs,pr) = A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
229  void sort2413(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
230  // sort2143: A2d_(qp,sr) = A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
231  void sort2143(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
232  // sort4231: A2d_(sq,rp) = A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
233  void sort4231(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
234  // sort3142: A2d_(rp,sq) = A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
235  void sort3142(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
236 
237  // sortp1432: A2d_(qp,rs) += A(pq,rs): d1 = num p, d2 = num q, d3 = num r, d4 = num s
238  void sortp1432(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
239  // sortp2134: A2d_(qp,rs) += A(pq,rs): d1 = num p, d2 = num q, d3 = num r, d4 = num s
240  void sortp2134(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2);
241  // sortp1243: A2d_(pq,sr) += A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
242  void sortp1243(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *col_pair_idx2);
243  // sortp2413: A2d_(qs,pr) += A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
244  void sortp2413(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
245  // sortp2143: A2d_(qp,sr) += A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
246  void sortp2143(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
247  // sortp4231: A2d_(sq,rp) += A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
248  void sortp4231(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
249  // sortp3142: A2d_(rp,sq) += A(pq,sr): d1 = num p, d2 = num q, d3 = num r, d4 = num s
250  void sortp3142(int d1, int d2, int d3, int d4, Array2d *A, Array2i *row_pair_idx, Array2i *col_pair_idx, Array2i *row_pair_idx2, Array2i *col_pair_idx2);
251 
252 
253  friend class Array1d;
254  friend class Array3d;
255  friend class Array1i;
256  friend class Array2i;
257 };
258 
259 class Array3d
260 {
261 
262  private:
263  double ***A3d_;
264  int dim1_,dim2_,dim3_;
265  string name_; // Name of the array
266 
267  public:
268  Array3d(int d1,int d2, int d3);
269  Array3d(string name, int d1,int d2, int d3);
270  Array3d(); //default constructer
271  ~Array3d(); //destructer
272 
273  Array3d* generate(int d1,int d2, int d3);
274  Array3d* generate(string name, int d1,int d2, int d3);
275  void init(string name, int d1,int d2, int d3);
276  void init(int d1,int d2, int d3);
277  void memalloc();
278  void zero();
279  void print();
280  void release();
281  void set(int h, int i, int j, double value);
282  double get(int h, int i, int j);
283 
284  friend class Array1d;
285  friend class Array2d;
286 };
287 
288 class Array1i
289 {
290 
291  private:
292  int *A1i_;
293  int dim1_;
294  string name_; // Name of the array
295 
296  public:
297  Array1i(int d1);
298  Array1i(string name, int d1);
299  Array1i(); //default constructer
300  ~Array1i(); //destructer
301 
302  Array1i* generate(int d1);
303  Array1i* generate(string name, int d1);
304  void init(string name, int d1);
305  void init(int d1);
306  void memalloc();
307  void zero();
308  void print();
309  void release();
310  void set(int i, int value);
311  int get(int i);
312  void add(const Array1i* Adum);
313  void add(int i, int value);
314  void subtract(const Array1i* Adum);
315  void subtract(int i, int value);
316 
317 };
318 
319 class Array2i
320 {
321 
322  private:
323  int **A2i_;
324  int dim1_,dim2_;
325  string name_; // Name of the array
326 
327  public:
328  Array2i(int d1,int d2);
329  Array2i(string name, int d1,int d2);
330  Array2i(); //default constructer
331  ~Array2i(); //destructer
332 
333  Array2i* generate(int d1,int d2);
334  Array2i* generate(string name, int d1,int d2);
335  void init(string name, int d1,int d2);
336  void init(int d1,int d2);
337  void memalloc();
338  void zero();
339  void zero_diagonal();
340  void print();
341  void print(std::string OutFileRMR);
342  void release();
343  void set(int i, int j, int value);
344  void set(int **A);
345  double get(int i, int j);
346  void add(const Array2i* Adum);
347  void add(int i, int j, int value);
348  void subtract(const Array2i* Adum);
349  void subtract(int i, int j, int value);
350  Array2i* transpose();
351  void copy(const Array2i* Adum);
352  void copy(int **a);
353  void identity();
354  int trace();
355  int **to_int_matrix();
356  int dim1() const { return dim1_; }
357  int dim2() const { return dim2_; }
358 
359  friend class Array1i;
360  friend class Array3i;
361  friend class Array1d;
362  friend class Array2d;
363 };
364 
365 
366 class Array3i
367 {
368 
369  private:
370  int ***A3i_;
371  int dim1_,dim2_,dim3_;
372  string name_; // Name of the array
373 
374  public:
375  Array3i(int d1,int d2, int d3);
376  Array3i(string name, int d1,int d2, int d3);
377  Array3i(); //default constructer
378  ~Array3i(); //destructer
379 
380  Array3i* generate(int d1,int d2, int d3);
381  Array3i* generate(string name, int d1,int d2, int d3);
382  void init(string name, int d1,int d2, int d3);
383  void init(int d1,int d2, int d3);
384  void memalloc();
385  void zero();
386  void print();
387  void release();
388  void set(int h, int i, int j, int value);
389  int get(int h, int i, int j);
390 };
391 }} // End Namespaces
392 #endif // _dfocc_arrays_h_
int ** A2i_
Definition: dfocc/arrays.h:323
string name_
Definition: dfocc/arrays.h:294
int dim2_
Definition: dfocc/arrays.h:108
double dot(BlockMatrix *A, BlockMatrix *B)
Definition: mcscf/block_matrix.cc:224
string name_
Definition: dfocc/arrays.h:265
int dim1() const
Definition: dfocc/arrays.h:356
Definition: dfocc/arrays.h:103
Definition: matrixtmp.h:42
Definition: pointgrp.h:106
int dim3_
Definition: dfocc/arrays.h:264
Definition: dfocc/arrays.h:366
int dim1() const
Definition: dfocc/arrays.h:94
Definition: dfocc/arrays.h:44
string name_
Definition: dfocc/arrays.h:109
int dim2() const
Definition: dfocc/arrays.h:204
string name_
Definition: dfocc/arrays.h:325
int dim1_
Definition: dfocc/arrays.h:293
double ** A2d_
Definition: dfocc/arrays.h:107
int dim3_
Definition: dfocc/arrays.h:371
int * A1i_
Definition: dfocc/arrays.h:292
int dim1_
Definition: dfocc/arrays.h:49
double * A1d_
Definition: dfocc/arrays.h:48
int dim2() const
Definition: dfocc/arrays.h:357
int dim1() const
Definition: dfocc/arrays.h:203
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:49
Definition: dfocc/arrays.h:319
string name_
Definition: dfocc/arrays.h:50
double *** A3d_
Definition: dfocc/arrays.h:263
Definition: dfocc/arrays.h:259
SharedWavefunction dfoccwave(SharedWavefunction, Options &)
Definition: dfocc/main.cc:36
string name_
Definition: dfocc/arrays.h:372
Definition: dfocc/arrays.h:288
int *** A3i_
Definition: dfocc/arrays.h:370
Definition: PsiFileImpl.h:39
int dim2_
Definition: dfocc/arrays.h:324
Definition: psio.hpp:197