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