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