Psi4
libmints/matrix.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-2018 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 _psi_src_lib_libmints_matrix_h_
30 #define _psi_src_lib_libmints_matrix_h_
31 
33 #include "psi4/libmints/typedefs.h"
35 
36 #include <cstdio>
37 #include <string>
38 #include <vector>
39 #include <memory>
40 
41 namespace psi {
42 
43 struct dpdfile2;
44 
45 class PSIO;
46 class Vector;
47 class Dimension;
48 class Molecule;
49 class Vector3;
50 
51 
54  ascending = 1,
57 };
58 
65 class PSI_API Matrix : public std::enable_shared_from_this<Matrix> {
66 protected:
68  double ***matrix_;
70  int nirrep_;
76  std::string name_;
78  int symmetry_;
79 
81  void alloc();
83  void release();
84 
86  void copy_from(double ***);
87 
89  static double** matrix(int nrow, int ncol);
91  static void free(double** Block);
92 
93  void print_mat(const double *const *const a, int m, int n, std::string out) const;
94 
96  std::vector<int> numpy_shape_;
97 
98 public:
99 
101  Matrix();
107  Matrix(const std::string& name, int symmetry = 0);
109  Matrix(const Matrix& copy);
111  explicit Matrix(const SharedMatrix& copy);
113  explicit Matrix(const Matrix* copy);
121  Matrix(int nirrep, const int *rowspi, const int *colspi, int symmetry = 0);
130  Matrix(const std::string& name, int nirrep, const int *rowspi, const int *colspi, int symmetry = 0);
137  Matrix(int nirrep, int rows, const int *colspi);
138 
145  Matrix(int nirrep, const int* rowspi, int cols);
146 
154  Matrix(int rows, int cols);
163  Matrix(const std::string& name, int rows, int cols);
164 
170  Matrix(dpdfile2 *inFile);
171 
180  Matrix(const std::string& name, const Dimension& rows, const Dimension& cols, int symmetry = 0);
181 
189  Matrix(const Dimension& rows, const Dimension& cols, int symmetry = 0);
190 
192  virtual ~Matrix();
193 
203  void init(int nirrep, const int *rowspi, const int *colspi, const std::string& name = "", int symmetry = 0);
204 
205  void init(const Dimension& rowspi, const Dimension& colspi, const std::string& name = "", int symmetry = 0);
206 
208  SharedMatrix clone() const;
209 
213  static SharedMatrix create(const std::string& name,
214  const Dimension& rows,
215  const Dimension& cols);
216 
222  void copy(const SharedMatrix& cp);
223  void copy(const Matrix& cp);
224  void copy(const Matrix* cp);
231  static SharedMatrix horzcat(const std::vector<SharedMatrix >& mats);
232 
237  static SharedMatrix vertcat(const std::vector<SharedMatrix >& mats);
238 
249  SharedMatrix matrix_3d_rotation(Vector3 axis, double phi, bool Sn);
250 
252  void copy_to_row(int h, int row, double const * const data);
253 
254  enum SaveType {
257  LowerTriangle
258  };
259 
270  bool load(psi::PSIO* psio, size_t fileno, const std::string& tocentry, int nso);
271  bool load(std::shared_ptr<psi::PSIO>& psio, size_t fileno, const std::string& tocentry, int nso);
283  void load(psi::PSIO* const psio, size_t fileno, SaveType savetype=LowerTriangle);
284  void load(std::shared_ptr<psi::PSIO>& psio, size_t fileno, SaveType savetype=LowerTriangle);
292  void load(const std::string& filename);
293 
300  void load_mpqc(const std::string& filename);
301 
311  void save(const std::string& filename, bool append=true, bool saveLowerTriangle = true, bool saveSubBlocks=false);
322  void save(psi::PSIO* const psio, size_t fileno, SaveType savetype=LowerTriangle);
323  void save(std::shared_ptr<psi::PSIO>& psio, size_t fileno, SaveType savetype=LowerTriangle);
331  void set(double val);
332 
338  void set(const double * const tri);
339 
346  void set(const double * const * const sq);
356  void set(const double * const * const sq, int irrep);
367  void set(int h, int m, int n, double val) { matrix_[h][m][n] = val; }
368 
376  void set(int m, int n, double val) { matrix_[0][m][n] = val; }
377 
384  void set_diagonal(const Vector * const vec);
385  void set_diagonal(const Vector& vec);
386  void set_diagonal(const std::shared_ptr<Vector>& vec);
397  double get(const int& h, const int& m, const int& n) const { return matrix_[h][m][n]; }
398 
406  double get(const int& m, const int& n) const { return matrix_[0][m][n]; }
407 
415  SharedVector get_row(int h, int m);
416 
424  SharedVector get_column(int h, int m);
425 
433  void set_row(int h, int m, SharedVector vec);
434 
442  void set_column(int h, int m, SharedVector vec);
443 
451  SharedMatrix get_block(const Slice& rows,const Slice& cols);
452 
460  void set_block(const Slice& rows,const Slice& cols,SharedMatrix block);
461 
474  double** pointer(const int& h = 0) const { return matrix_[h]; }
475  const double** const_pointer(const int& h=0) const { return const_cast<const double**>(matrix_[h]); }
476 
489  double* get_pointer(const int& h = 0) const {
490  if(rowspi_[h]*(size_t)colspi_[h] > 0)
491  return &(matrix_[h][0][0]);
492  else
493  return 0;}
494  const double* get_const_pointer(const int& h=0) const {
495  if(rowspi_[h]*(size_t)colspi_[h] > 0)
496  return const_cast<const double*>(&(matrix_[h][0][0]));
497  else
498  return 0;}
499 
500  size_t size(const int &h=0) const { return colspi_[h] * (size_t)rowspi_[h]; }
501 
503  void apply_denominator(const Matrix* const);
505  void apply_denominator(const Matrix&);
507  void apply_denominator(const SharedMatrix&);
508 
514  double **to_block_matrix() const;
520  SharedMatrix to_block_sharedmatrix() const;
526  double *to_lower_triangle() const;
527 
533  void set_name(const std::string& name) { name_ = name; }
534 
538  const std::string& name() const { return name_; }
539 
541  void print_out() const { print("outfile"); }
542 
549  void print(std::string outfile = "outfile", const char *extra=nullptr) const;
550 
552  void print_atom_vector(std::string out_fname = "outfile");
553 
557  void print_to_mathematica();
558 
565  void eivprint(const Vector * const values, std::string out = "outfile");
567  void eivprint(const Vector& values, std::string out = "outfile");
569  void eivprint(const std::shared_ptr<Vector>& values, std::string out = "outfile");
570 
572  int rowdim(const int& h = 0) const { return rowspi_[h]; }
574  int coldim(const int& h = 0) const { return colspi_[h]; }
575 
577  const Dimension& rowspi() const {
578  return rowspi_;
579  }
581  int rowspi(const int& h) const {
582  return rowdim(h);
583  }
585  const Dimension& colspi() const {
586  return colspi_;
587  }
589  int colspi(const int& h) const {
590  return coldim(h);
591  }
593  const int& nirrep() const {
594  return nirrep_;
595  }
596 
598  int nrow() const {
599  int rows = 0;
600  for (int h=0; h<nirrep(); ++h)
601  rows += rowdim(h);
602  return rows;
603  }
604 
606  int ncol() const {
607  int cols = 0;
608  for (int h=0; h<nirrep(); ++h)
609  cols += coldim(h);
610  return cols;
611  }
612 
614  int max_nrow() const {
615  int row = 0;
616  for (int h=0; h<nirrep(); ++h)
617  if (row < rowdim(h))
618  row = rowdim(h);
619  return row;
620  }
621 
623  int max_ncol() const {
624  int col = 0;
625  for (int h=0; h<nirrep(); ++h)
626  if (col < coldim(h))
627  col = coldim(h);
628  return col;
629  }
630 
636  const int& symmetry() const {
637  return symmetry_;
638  }
639 
644  void symmetrize_gradient(std::shared_ptr<Molecule> mol);
645 
650  void symmetrize_hessian(std::shared_ptr<Molecule> mol);
651 
653  void identity();
655  void zero();
657  void zero_diagonal();
658 
659  // Math routines
661  double trace();
663  SharedMatrix transpose();
664 
666  void transpose_this();
667 
669  void add(const Matrix* const);
671  void add(const Matrix&);
673  void add(const SharedMatrix&);
674 
676  void subtract(const Matrix* const);
678  void subtract(const SharedMatrix&);
680  void accumulate_product(const Matrix* const, const Matrix* const);
681  void accumulate_product(const SharedMatrix&, const SharedMatrix&);
683  void scale(double);
685  double sum_of_squares();
687  double rms();
689  double absmax();
691  void add(int h, int m, int n, double val) {
692  #ifdef PSIDEBUG
693  if (m > rowspi_[h] || n > colspi_[h^symmetry_]) {
694  outfile->Printf( "out of bounds: symmetry_ = %d, h = %d, m = %d, n = %d\n",
695  symmetry_, h, m, n);
696 
697  throw PSIEXCEPTION("What are you doing, Rob?");
698  }
699  #endif
700  matrix_[h][m][n] += val;
701  }
703  void add(int m, int n, double val) {
704  #ifdef PSIDEBUG
705  if (m > rowspi_[0] || n > colspi_[0^symmetry_]) {
706  outfile->Printf( "out of bounds: symmetry_ = %d, h = %d, m = %d, n = %d\n",
707  symmetry_, 0, m, n);
708 
709  return;
710  }
711  #endif
712  matrix_[0][m][n] += val;
713  }
714 
716  for (int h=0; h<nirrep_; ++h) {
717  for (int i=0; i<rowspi_[h]; ++i) {
718  for (int j=0; j<i; ++j) {
719  matrix_[h][i][j] = matrix_[h][j][i] = (matrix_[h][i][j] + matrix_[h][j][i]);
720  }
721  }
722  }
723  }
724 
726  void scale_row(int h, int m, double a);
728  void scale_column(int h, int n, double a);
729 
735  void apply_symmetry(const SharedMatrix& a, const SharedMatrix& transformer);
736 
742  void remove_symmetry(const SharedMatrix& a, const SharedMatrix& transformer);
749  void transform(const SharedMatrix& L,
750  const SharedMatrix& F,
751  const SharedMatrix& R);
752 
755  void transform(const Matrix* const a, const Matrix* const transformer);
756  void transform(const SharedMatrix& a, const SharedMatrix& transformer);
758 
761  void transform(const Matrix* const transformer);
762  void transform(const SharedMatrix& transformer);
764 
767  void back_transform(const Matrix* const a, const Matrix* const transformer);
768  void back_transform(const SharedMatrix& a, const SharedMatrix& transformer);
770 
773  void back_transform(const Matrix* const transformer);
774  void back_transform(const SharedMatrix& transformer);
776 
778  double vector_dot(const Matrix* const rhs);
779  double vector_dot(const SharedMatrix& rhs);
780  double vector_dot(const Matrix& rhs);
781 
783 
791  void gemm(bool transa, bool transb, double alpha, const Matrix* const a, const Matrix* const b, double beta);
792  void gemm(bool transa, bool transb, double alpha, const SharedMatrix& a, const SharedMatrix& b, double beta);
793  void gemm(bool transa, bool transb, double alpha, const SharedMatrix& a, const Matrix& b, double beta);
794  void gemm(bool transa, bool transb, double alpha, const Matrix& a, const SharedMatrix& b, double beta);
796 
798 
800  void gemm(const char& transa, const char& transb,
801  const std::vector<int>& m,
802  const std::vector<int>& n,
803  const std::vector<int>& k,
804  const double& alpha,
805  const SharedMatrix& a, const std::vector<int>& lda,
806  const SharedMatrix& b, const std::vector<int>& ldb,
807  const double& beta,
808  const std::vector<int>& ldc,
809  const std::vector<unsigned long>& offset_a = std::vector<unsigned long>(),
810  const std::vector<unsigned long>& offset_b = std::vector<unsigned long>(),
811  const std::vector<unsigned long>& offset_c = std::vector<unsigned long>());
812  void gemm(const char& transa, const char& transb,
813  const int& m,
814  const int& n,
815  const int& k,
816  const double& alpha,
817  const SharedMatrix& a, const int& lda,
818  const SharedMatrix& b, const int& ldb,
819  const double& beta,
820  const int& ldc,
821  const unsigned long& offset_a = 0,
822  const unsigned long& offset_b = 0,
823  const unsigned long& offset_c = 0);
825 
832  static SharedMatrix doublet(const SharedMatrix& A, const SharedMatrix& B, bool transA = false, bool transB = false);
833 
842  static SharedMatrix triplet(const SharedMatrix& A, const SharedMatrix& B, const SharedMatrix& C, bool transA = false, bool transB = false, bool transC = false);
843 
849  void axpy(double a, SharedMatrix X);
850 
856  SharedMatrix collapse(int dim = 0);
857 
860  void diagonalize(Matrix* eigvectors, Vector* eigvalues, diagonalize_order nMatz = ascending);
861  void diagonalize(SharedMatrix& eigvectors, std::shared_ptr<Vector>& eigvalues, diagonalize_order nMatz = ascending);
862  void diagonalize(SharedMatrix& eigvectors, Vector& eigvalues, diagonalize_order nMatz = ascending);
864 
867  void diagonalize(SharedMatrix& metric, SharedMatrix& eigvectors, std::shared_ptr<Vector>& eigvalues, diagonalize_order nMatz = ascending);
869 
872  void svd(SharedMatrix& U, SharedVector& S, SharedMatrix& V);
874 
879  void svd_a(SharedMatrix& U, SharedVector& S, SharedMatrix& V);
881 
884  std::tuple<SharedMatrix,SharedVector,SharedMatrix> svd_temps();
886 
889  std::tuple<SharedMatrix,SharedVector,SharedMatrix> svd_a_temps();
891 
894  SharedMatrix pseudoinverse(double condition, int &nremoved);
896 
902  SharedMatrix canonical_orthogonalization(double delta = 0.0, SharedMatrix eigvec = SharedMatrix());
903 
930  SharedMatrix partial_cholesky_factorize(double delta = 0.0, bool throw_if_negative = false);
931 
946  std::pair<SharedMatrix, SharedMatrix> partial_square_root(double delta = 0.0);
947 
953  void cholesky_factorize();
954 
959  void invert();
960 
965  void general_invert();
966 
988  Dimension power(double alpha, double cutoff = 1.0E-12);
989 
997  void expm(int n = 2, bool scale = false);
998 
1000  void swap_rows(int h, int i, int j);
1002  void swap_columns(int h, int i, int j);
1003 
1005  void hermitivitize();
1007  void copy_lower_to_upper();
1009  void copy_upper_to_lower();
1011  void zero_lower();
1013  void zero_upper();
1015  void zero_row(int h, int i);
1017  void zero_column(int h, int i);
1018 
1019  // Reference versions of the above functions
1021  void transform(const Matrix& a, const Matrix& transformer);
1023  void transform(const Matrix& transformer);
1025  void back_transform(const Matrix& a, const Matrix& transformer);
1027  void back_transform(const Matrix& transformer);
1028 
1033  bool add_and_orthogonalize_row(const SharedVector v);
1034 
1050  bool schmidt_add_row(int h, int rows, Vector& v) throw();
1051  bool schmidt_add_row(int h, int rows, double* v) throw();
1053 
1055  void schmidt();
1056 
1059  void schmidt_orthog(SharedMatrix S, int n);
1060 
1067  Dimension schmidt_orthog_columns(SharedMatrix S, double tol, double*res=0);
1068 
1076  void project_out(Matrix& v);
1077 
1079  void gemm(bool transa, bool transb, double alpha, const Matrix& a, const Matrix& b, double beta);
1081  void diagonalize(Matrix& eigvectors, Vector& eigvalues, int nMatz = 1);
1082 
1085  double** operator[](int i) { return matrix_[i]; }
1086  double& operator()(int i, int j) { return matrix_[0][i][j]; }
1087  const double& operator()(int i, int j) const { return matrix_[0][i][j]; }
1088  double& operator()(int h, int i, int j) { return matrix_[h][i][j]; }
1089  const double& operator()(int h, int i, int j) const { return matrix_[h][i][j]; }
1091 
1093  void write_to_dpdfile2(dpdfile2 *outFile);
1094 
1099  bool equal(const Matrix& rhs, double TOL=1.0e-10);
1100  bool equal(const SharedMatrix& rhs, double TOL=1.0e-10);
1101  bool equal(const Matrix* rhs, double TOL=1.0e-10);
1103 
1108  bool equal_but_for_row_order(const Matrix& rhs, double TOL=1.0e-10);
1109  bool equal_but_for_row_order(const SharedMatrix& rhs, double TOL=1.0e-10);
1110  bool equal_but_for_row_order(const Matrix* rhs, double TOL=1.0e-10);
1112 
1116  void set_numpy_shape(std::vector<int> shape) { numpy_shape_ = shape; }
1117  std::vector<int> numpy_shape() { return numpy_shape_; }
1118 
1126  void rotate_columns(int h, int i, int j, double theta);
1127  friend class Vector;
1128 };
1129 
1130 }
1131 
1132 
1133 #endif // MATRIX_H
std::vector< int > numpy_shape_
Numpy Shape.
Definition: libmints/matrix.h:96
Definition: libmints/matrix.h:256
const double & operator()(int h, int i, int j) const
Definition: libmints/matrix.h:1089
std::vector< int > numpy_shape()
Definition: libmints/matrix.h:1117
int rowdim(const int &h=0) const
Returns the rows in irrep h.
Definition: libmints/matrix.h:572
int colspi(const int &h) const
Returns the columns per irrep array.
Definition: libmints/matrix.h:589
int max_nrow() const
Returns the row size of the largest block.
Definition: libmints/matrix.h:614
int * U
Definition: stringlist.cc:66
#define TOL
Definition: detci/opdm.cc:57
void set(int m, int n, double val)
Definition: libmints/matrix.h:376
void set_block(Slice rows, Slice cols, SharedMatrix block)
size_t size(const int &h=0) const
Definition: libmints/matrix.h:500
void schmidt(double **A, int rows, int cols, std::string out_fname)
Definition: schmidt.cc:60
void add(int m, int n, double val)
Add val to an element of this.
Definition: libmints/matrix.h:703
void set_name(const std::string &name)
Definition: libmints/matrix.h:533
Definition: matrixtmp.h:42
Definition: vector3.h:37
Definition: pointgrp.h:106
int ncol() const
Returns the total number of columns.
Definition: libmints/matrix.h:606
const Dimension & colspi() const
Returns the columns per irrep array.
Definition: libmints/matrix.h:585
int nirrep_
Number of irreps.
Definition: libmints/matrix.h:70
std::regex symmetry_("\\s*symmetry[\\s=]+(\\w+)\\s*", std::regex_constants::icase)
double & operator()(int i, int j)
Definition: libmints/matrix.h:1086
int symmetry_
Symmetry of this matrix (in most cases this will be 0 [totally symmetric])
Definition: libmints/matrix.h:78
void print_mat(double **a, int rows, int cols, std::string out)
Definition: print_mat.cc:52
std::shared_ptr< PsiOutStream > outfile
Definition: core.cc:103
double ** pointer(const int &h=0) const
Definition: libmints/matrix.h:474
void add(int h, int m, int n, double val)
Add val to an element of this.
Definition: libmints/matrix.h:691
Slicing for Matrices and Vectors objects.
Definition: dimension.h:125
const double * get_const_pointer(const int &h=0) const
Definition: libmints/matrix.h:494
Definition: libmints/matrix.h:56
void element_add_mirror()
Definition: libmints/matrix.h:715
double * get_pointer(const int &h=0) const
Definition: libmints/matrix.h:489
diagonalize_order
Definition: libmints/matrix.h:52
double *** matrix_
Matrix data.
Definition: libmints/matrix.h:68
int rowspi(const int &h) const
Returns the rows per irrep array.
Definition: libmints/matrix.h:581
long int size_t
Definition: sortintegrals.cc:40
void print_out() const
Python compatible printer.
Definition: libmints/matrix.h:541
const double ** const_pointer(const int &h=0) const
Definition: libmints/matrix.h:475
Definition: dimension.h:40
Definition: libmints/matrix.h:53
Makes using matrices just a little easier.
Definition: libmints/matrix.h:65
SaveType
Definition: libmints/matrix.h:254
const double & operator()(int i, int j) const
Definition: libmints/matrix.h:1087
const int & symmetry() const
Definition: libmints/matrix.h:636
Definition: libmints/matrix.h:255
Definition: libdpd/dpd.h:140
int nso
Definition: dx_write.cc:56
int coldim(const int &h=0) const
Returns the cols in irrep h.
Definition: libmints/matrix.h:574
const Dimension & rowspi() const
Returns the rows per irrep array.
Definition: libmints/matrix.h:577
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:49
int delta(const int i, const int j)
Definition: bend.cc:175
std::string name_
Name of the matrix.
Definition: libmints/matrix.h:76
#define PSI_API
Definition: pragma.h:128
void set(int h, int m, int n, double val)
Definition: libmints/matrix.h:367
void set_numpy_shape(std::vector< int > shape)
Definition: libmints/matrix.h:1116
Dimension colspi_
Columns per irrep array.
Definition: libmints/matrix.h:74
Definition: pointgrp.h:106
Dimension rowspi_
Rows per irrep array.
Definition: libmints/matrix.h:72
Definition: libmints/matrix.h:55
double & operator()(int h, int i, int j)
Definition: libmints/matrix.h:1088
const int & nirrep() const
Returns the number of irreps.
Definition: libmints/matrix.h:593
const std::string & name() const
Definition: libmints/matrix.h:538
#define PSIEXCEPTION(message)
Definition: exception.h:50
Definition: vector.h:48
std::shared_ptr< Vector > SharedVector
Definition: adc.h:51
int nrow() const
Returns the total number of rows.
Definition: libmints/matrix.h:598
double ** operator[](int i)
Definition: libmints/matrix.h:1085
int max_ncol() const
Returns the column size of the largest block.
Definition: libmints/matrix.h:623
Definition: psio.hpp:195
Definition: libmints/matrix.h:54