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 
52 
59 class PSI_API Matrix : public std::enable_shared_from_this<Matrix> {
60  protected:
62  double*** matrix_;
64  int nirrep_;
70  std::string name_;
72  int symmetry_;
73 
75  void alloc();
77  void release();
78 
80  void copy_from(double***);
81 
83  static double** matrix(int nrow, int ncol);
85  static void free(double** Block);
86 
87  void print_mat(const double* const* const a, int m, int n, std::string out) const;
88 
90  std::vector<int> numpy_shape_;
91 
92  public:
94  Matrix();
100  Matrix(const std::string& name, int symmetry = 0);
102  Matrix(const Matrix& copy);
104  explicit Matrix(const SharedMatrix& copy);
106  explicit Matrix(const Matrix* copy);
114  Matrix(int nirrep, const int* rowspi, const int* colspi, int symmetry = 0);
123  Matrix(const std::string& name, int nirrep, const int* rowspi, const int* colspi, int symmetry = 0);
130  Matrix(int nirrep, int rows, const int* colspi);
131 
138  Matrix(int nirrep, const int* rowspi, int cols);
139 
147  Matrix(int rows, int cols);
156  Matrix(const std::string& name, int rows, int cols);
157 
163  Matrix(dpdfile2* inFile);
164 
173  Matrix(const std::string& name, const Dimension& rows, const Dimension& cols, int symmetry = 0);
174 
182  Matrix(const Dimension& rows, const Dimension& cols, int symmetry = 0);
183 
185  virtual ~Matrix();
186 
196  void init(int nirrep, const int* rowspi, const int* colspi, const std::string& name = "", int symmetry = 0);
197 
198  void init(const Dimension& rowspi, const Dimension& colspi, const std::string& name = "", int symmetry = 0);
199 
201  SharedMatrix clone() const;
202 
206  static SharedMatrix create(const std::string& name, const Dimension& rows, const Dimension& cols);
207 
213  void copy(const SharedMatrix& cp);
214  void copy(const Matrix& cp);
215  void copy(const Matrix* cp);
222  static SharedMatrix horzcat(const std::vector<SharedMatrix>& mats);
223 
228  static SharedMatrix vertcat(const std::vector<SharedMatrix>& mats);
229 
240  SharedMatrix matrix_3d_rotation(Vector3 axis, double phi, bool Sn);
241 
243  void copy_to_row(int h, int row, double const* const data);
244 
245  enum SaveType { Full, SubBlocks, LowerTriangle };
246 
257  bool load(psi::PSIO* psio, size_t fileno, const std::string& tocentry, int nso);
258  bool load(std::shared_ptr<psi::PSIO>& psio, size_t fileno, const std::string& tocentry, int nso);
270  void load(psi::PSIO* const psio, size_t fileno, SaveType savetype = LowerTriangle);
271  void load(std::shared_ptr<psi::PSIO>& psio, size_t fileno, SaveType savetype = LowerTriangle);
279  void load(const std::string& filename);
280 
287  void load_mpqc(const std::string& filename);
288 
299  void save(const std::string& filename, bool append = true, bool saveLowerTriangle = true,
300  bool saveSubBlocks = false);
311  void save(psi::PSIO* const psio, size_t fileno, SaveType savetype = LowerTriangle);
312  void save(std::shared_ptr<psi::PSIO>& psio, size_t fileno, SaveType savetype = LowerTriangle);
320  void set(double val);
321 
327  void set(const double* const tri);
328 
335  void set(const double* const* const sq);
345  void set(const double* const* const sq, int irrep);
356  void set(int h, int m, int n, double val) { matrix_[h][m][n] = val; }
357 
365  void set(int m, int n, double val) { matrix_[0][m][n] = val; }
366 
373  void set_diagonal(const Vector* const vec);
374  void set_diagonal(const Vector& vec);
375  void set_diagonal(const std::shared_ptr<Vector>& vec);
386  double get(const int& h, const int& m, const int& n) const { return matrix_[h][m][n]; }
387 
395  double get(const int& m, const int& n) const { return matrix_[0][m][n]; }
396 
404  SharedVector get_row(int h, int m);
405 
413  SharedVector get_column(int h, int m);
414 
422  void set_row(int h, int m, SharedVector vec);
423 
431  void set_column(int h, int m, SharedVector vec);
432 
440  SharedMatrix get_block(const Slice& rows, const Slice& cols);
441 
449  void set_block(const Slice& rows, const Slice& cols, SharedMatrix block);
450 
463  double** pointer(const int& h = 0) const { return matrix_[h]; }
464  const double** const_pointer(const int& h = 0) const { return const_cast<const double**>(matrix_[h]); }
465 
478  double* get_pointer(const int& h = 0) const {
479  if (rowspi_[h] * (size_t)colspi_[h] > 0)
480  return &(matrix_[h][0][0]);
481  else
482  return nullptr;
483  }
484  const double* get_const_pointer(const int& h = 0) const {
485  if (rowspi_[h] * (size_t)colspi_[h] > 0)
486  return const_cast<const double*>(&(matrix_[h][0][0]));
487  else
488  return nullptr;
489  }
490 
491  size_t size(const int& h = 0) const { return colspi_[h] * (size_t)rowspi_[h]; }
492 
494  void apply_denominator(const Matrix* const);
496  void apply_denominator(const Matrix&);
498  void apply_denominator(const SharedMatrix&);
499 
505  double** to_block_matrix() const;
511  SharedMatrix to_block_sharedmatrix() const;
517  double* to_lower_triangle() const;
518 
524  void set_name(const std::string& name) { name_ = name; }
525 
529  const std::string& name() const { return name_; }
530 
532  void print_out() const { print("outfile"); }
533 
540  void print(std::string outfile = "outfile", const char* extra = nullptr) const;
541 
543  void print_atom_vector(std::string out_fname = "outfile");
544 
548  void print_to_mathematica();
549 
556  void eivprint(const Vector* const values, std::string out = "outfile");
558  void eivprint(const Vector& values, std::string out = "outfile");
560  void eivprint(const std::shared_ptr<Vector>& values, std::string out = "outfile");
561 
563  int rowdim(const int& h = 0) const { return rowspi_[h]; }
565  int coldim(const int& h = 0) const { return colspi_[h]; }
566 
568  const Dimension& rowspi() const { return rowspi_; }
570  int rowspi(const int& h) const { return rowdim(h); }
572  const Dimension& colspi() const { return colspi_; }
574  int colspi(const int& h) const { return coldim(h); }
576  const int& nirrep() const { return nirrep_; }
577 
579  int nrow() const {
580  int rows = 0;
581  for (int h = 0; h < nirrep(); ++h) rows += rowdim(h);
582  return rows;
583  }
584 
586  int ncol() const {
587  int cols = 0;
588  for (int h = 0; h < nirrep(); ++h) cols += coldim(h);
589  return cols;
590  }
591 
593  int max_nrow() const {
594  int row = 0;
595  for (int h = 0; h < nirrep(); ++h)
596  if (row < rowdim(h)) row = rowdim(h);
597  return row;
598  }
599 
601  int max_ncol() const {
602  int col = 0;
603  for (int h = 0; h < nirrep(); ++h)
604  if (col < coldim(h)) col = coldim(h);
605  return col;
606  }
607 
613  const int& symmetry() const { return symmetry_; }
614 
619  void symmetrize_gradient(std::shared_ptr<Molecule> mol);
620 
625  void symmetrize_hessian(std::shared_ptr<Molecule> mol);
626 
628  void identity();
630  void zero();
632  void zero_diagonal();
633 
634  // Math routines
636  double trace();
638  SharedMatrix transpose();
639 
641  void transpose_this();
642 
644  void add(const Matrix* const);
646  void add(const Matrix&);
648  void add(const SharedMatrix&);
649 
651  void subtract(const Matrix* const);
653  void subtract(const SharedMatrix&);
655  void accumulate_product(const Matrix* const, const Matrix* const);
656  void accumulate_product(const SharedMatrix&, const SharedMatrix&);
658  void scale(double);
660  double sum_of_squares();
662  double rms();
664  double absmax();
666  void add(int h, int m, int n, double val) {
667 #ifdef PSIDEBUG
668  if (m > rowspi_[h] || n > colspi_[h ^ symmetry_]) {
669  outfile->Printf("out of bounds: symmetry_ = %d, h = %d, m = %d, n = %d\n", symmetry_, h, m, n);
670 
671  throw PSIEXCEPTION("What are you doing, Rob?");
672  }
673 #endif
674  matrix_[h][m][n] += val;
675  }
677  void add(int m, int n, double val) {
678 #ifdef PSIDEBUG
679  if (m > rowspi_[0] || n > colspi_[0 ^ symmetry_]) {
680  outfile->Printf("out of bounds: symmetry_ = %d, h = %d, m = %d, n = %d\n", symmetry_, 0, m, n);
681 
682  return;
683  }
684 #endif
685  matrix_[0][m][n] += val;
686  }
687 
689  for (int h = 0; h < nirrep_; ++h) {
690  for (int i = 0; i < rowspi_[h]; ++i) {
691  for (int j = 0; j < i; ++j) {
692  matrix_[h][i][j] = matrix_[h][j][i] = (matrix_[h][i][j] + matrix_[h][j][i]);
693  }
694  }
695  }
696  }
697 
699  void scale_row(int h, int m, double a);
701  void scale_column(int h, int n, double a);
702 
708  void apply_symmetry(const SharedMatrix& a, const SharedMatrix& transformer);
709 
715  void remove_symmetry(const SharedMatrix& a, const SharedMatrix& transformer);
722  void transform(const SharedMatrix& L, const SharedMatrix& F, const SharedMatrix& R);
723 
726  void transform(const Matrix* const a, const Matrix* const transformer);
727  void transform(const SharedMatrix& a, const SharedMatrix& transformer);
729 
732  void transform(const Matrix* const transformer);
733  void transform(const SharedMatrix& transformer);
735 
738  void back_transform(const Matrix* const a, const Matrix* const transformer);
739  void back_transform(const SharedMatrix& a, const SharedMatrix& transformer);
741 
744  void back_transform(const Matrix* const transformer);
745  void back_transform(const SharedMatrix& transformer);
747 
749  double vector_dot(const Matrix* const rhs);
750  double vector_dot(const SharedMatrix& rhs);
751  double vector_dot(const Matrix& rhs);
752 
754 
762  void gemm(bool transa, bool transb, double alpha, const Matrix* const a, const Matrix* const b, double beta);
763  void gemm(bool transa, bool transb, double alpha, const SharedMatrix& a, const SharedMatrix& b, double beta);
764  void gemm(bool transa, bool transb, double alpha, const SharedMatrix& a, const Matrix& b, double beta);
765  void gemm(bool transa, bool transb, double alpha, const Matrix& a, const SharedMatrix& b, double beta);
767 
769 
771  void gemm(const char& transa, const char& transb, const std::vector<int>& m, const std::vector<int>& n,
772  const std::vector<int>& k, const double& alpha, const SharedMatrix& a, const std::vector<int>& lda,
773  const SharedMatrix& b, const std::vector<int>& ldb, const double& beta, const std::vector<int>& ldc,
774  const std::vector<unsigned long>& offset_a = std::vector<unsigned long>(),
775  const std::vector<unsigned long>& offset_b = std::vector<unsigned long>(),
776  const std::vector<unsigned long>& offset_c = std::vector<unsigned long>());
777  void gemm(const char& transa, const char& transb, const int& m, const int& n, const int& k, const double& alpha,
778  const SharedMatrix& a, const int& lda, const SharedMatrix& b, const int& ldb, const double& beta,
779  const int& ldc, const unsigned long& offset_a = 0, const unsigned long& offset_b = 0,
780  const unsigned long& offset_c = 0);
782 
789  static SharedMatrix doublet(const SharedMatrix& A, const SharedMatrix& B, bool transA = false, bool transB = false);
790 
799  static SharedMatrix triplet(const SharedMatrix& A, const SharedMatrix& B, const SharedMatrix& C,
800  bool transA = false, bool transB = false, bool transC = false);
801 
807  void axpy(double a, SharedMatrix X);
808 
814  SharedMatrix collapse(int dim = 0);
815 
818  void diagonalize(Matrix* eigvectors, Vector* eigvalues, diagonalize_order nMatz = ascending);
819  void diagonalize(SharedMatrix& eigvectors, std::shared_ptr<Vector>& eigvalues, diagonalize_order nMatz = ascending);
820  void diagonalize(SharedMatrix& eigvectors, Vector& eigvalues, diagonalize_order nMatz = ascending);
822 
826  void diagonalize(SharedMatrix& metric, SharedMatrix& eigvectors, std::shared_ptr<Vector>& eigvalues,
827  diagonalize_order nMatz = ascending);
829 
832  void svd(SharedMatrix& U, SharedVector& S, SharedMatrix& V);
834 
839  void svd_a(SharedMatrix& U, SharedVector& S, SharedMatrix& V);
841 
844  std::tuple<SharedMatrix, SharedVector, SharedMatrix> svd_temps();
846 
849  std::tuple<SharedMatrix, SharedVector, SharedMatrix> svd_a_temps();
851 
854  SharedMatrix pseudoinverse(double condition, int& nremoved);
856 
862  SharedMatrix canonical_orthogonalization(double delta = 0.0, SharedMatrix eigvec = SharedMatrix());
863 
890  SharedMatrix partial_cholesky_factorize(double delta = 0.0, bool throw_if_negative = false);
891 
906  std::pair<SharedMatrix, SharedMatrix> partial_square_root(double delta = 0.0);
907 
913  void cholesky_factorize();
914 
919  void invert();
920 
925  void general_invert();
926 
948  Dimension power(double alpha, double cutoff = 1.0E-12);
949 
957  void expm(int n = 2, bool scale = false);
958 
960  void swap_rows(int h, int i, int j);
962  void swap_columns(int h, int i, int j);
963 
965  void hermitivitize();
967  void copy_lower_to_upper();
969  void copy_upper_to_lower();
971  void zero_lower();
973  void zero_upper();
975  void zero_row(int h, int i);
977  void zero_column(int h, int i);
978 
979  // Reference versions of the above functions
981  void transform(const Matrix& a, const Matrix& transformer);
983  void transform(const Matrix& transformer);
985  void back_transform(const Matrix& a, const Matrix& transformer);
987  void back_transform(const Matrix& transformer);
988 
993  bool add_and_orthogonalize_row(const SharedVector v);
994 
1010  bool schmidt_add_row(int h, int rows, Vector& v) noexcept;
1011  bool schmidt_add_row(int h, int rows, double* v) noexcept;
1013 
1015  void schmidt();
1016 
1019  void schmidt_orthog(SharedMatrix S, int n);
1020 
1027  Dimension schmidt_orthog_columns(SharedMatrix S, double tol, double* res = nullptr);
1028 
1036  void project_out(Matrix& v);
1037 
1039  void gemm(bool transa, bool transb, double alpha, const Matrix& a, const Matrix& b, double beta);
1041  void diagonalize(Matrix& eigvectors, Vector& eigvalues, int nMatz = 1);
1042 
1045  double** operator[](int i) { return matrix_[i]; }
1046  double& operator()(int i, int j) { return matrix_[0][i][j]; }
1047  const double& operator()(int i, int j) const { return matrix_[0][i][j]; }
1048  double& operator()(int h, int i, int j) { return matrix_[h][i][j]; }
1049  const double& operator()(int h, int i, int j) const { return matrix_[h][i][j]; }
1051 
1053  void write_to_dpdfile2(dpdfile2* outFile);
1054 
1059  bool equal(const Matrix& rhs, double TOL = 1.0e-10);
1060  bool equal(const SharedMatrix& rhs, double TOL = 1.0e-10);
1061  bool equal(const Matrix* rhs, double TOL = 1.0e-10);
1063 
1068  bool equal_but_for_row_order(const Matrix& rhs, double TOL = 1.0e-10);
1069  bool equal_but_for_row_order(const SharedMatrix& rhs, double TOL = 1.0e-10);
1070  bool equal_but_for_row_order(const Matrix* rhs, double TOL = 1.0e-10);
1072 
1076  void set_numpy_shape(std::vector<int> shape) { numpy_shape_ = shape; }
1077  std::vector<int> numpy_shape() { return numpy_shape_; }
1078 
1086  void rotate_columns(int h, int i, int j, double theta);
1087  friend class Vector;
1088 };
1089 
1090 } // namespace psi
1091 
1092 #endif // MATRIX_H
std::vector< int > numpy_shape_
Numpy Shape.
Definition: libmints/matrix.h:90
Definition: libmints/matrix.h:245
const double & operator()(int h, int i, int j) const
Definition: libmints/matrix.h:1049
std::vector< int > numpy_shape()
Definition: libmints/matrix.h:1077
int rowdim(const int &h=0) const
Returns the rows in irrep h.
Definition: libmints/matrix.h:563
int colspi(const int &h) const
Returns the columns per irrep array.
Definition: libmints/matrix.h:574
int max_nrow() const
Returns the row size of the largest block.
Definition: libmints/matrix.h:593
int * U
Definition: stringlist.cc:66
#define TOL
Definition: detci/opdm.cc:58
void set(int m, int n, double val)
Definition: libmints/matrix.h:365
void set_block(Slice rows, Slice cols, SharedMatrix block)
size_t size(const int &h=0) const
Definition: libmints/matrix.h:491
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:677
void set_name(const std::string &name)
Definition: libmints/matrix.h:524
Definition: matrixtmp.h:42
Definition: vector3.h:38
Definition: pointgrp.h:104
int ncol() const
Returns the total number of columns.
Definition: libmints/matrix.h:586
const Dimension & colspi() const
Returns the columns per irrep array.
Definition: libmints/matrix.h:572
int nirrep_
Number of irreps.
Definition: libmints/matrix.h:64
double & operator()(int i, int j)
Definition: libmints/matrix.h:1046
int symmetry_
Symmetry of this matrix (in most cases this will be 0 [totally symmetric])
Definition: libmints/matrix.h:72
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:102
double ** pointer(const int &h=0) const
Definition: libmints/matrix.h:463
void add(int h, int m, int n, double val)
Add val to an element of this.
Definition: libmints/matrix.h:666
Slicing for Matrices and Vectors objects.
Definition: dimension.h:125
const double * get_const_pointer(const int &h=0) const
Definition: libmints/matrix.h:484
Definition: libmints/matrix.h:51
void element_add_mirror()
Definition: libmints/matrix.h:688
double * get_pointer(const int &h=0) const
Definition: libmints/matrix.h:478
diagonalize_order
Definition: libmints/matrix.h:51
double *** matrix_
Matrix data.
Definition: libmints/matrix.h:62
int rowspi(const int &h) const
Returns the rows per irrep array.
Definition: libmints/matrix.h:570
long int size_t
Definition: sortintegrals.cc:41
void print_out() const
Python compatible printer.
Definition: libmints/matrix.h:532
const double ** const_pointer(const int &h=0) const
Definition: libmints/matrix.h:464
Definition: dimension.h:40
Definition: libmints/matrix.h:51
Makes using matrices just a little easier.
Definition: libmints/matrix.h:59
SaveType
Definition: libmints/matrix.h:245
const double & operator()(int i, int j) const
Definition: libmints/matrix.h:1047
const int & symmetry() const
Definition: libmints/matrix.h:613
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:565
const Dimension & rowspi() const
Returns the rows per irrep array.
Definition: libmints/matrix.h:568
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:70
#define PSI_API
Definition: pragma.h:155
void set(int h, int m, int n, double val)
Definition: libmints/matrix.h:356
void set_numpy_shape(std::vector< int > shape)
Definition: libmints/matrix.h:1076
Dimension colspi_
Columns per irrep array.
Definition: libmints/matrix.h:68
Definition: pointgrp.h:104
Dimension rowspi_
Rows per irrep array.
Definition: libmints/matrix.h:66
Definition: libmints/matrix.h:51
double & operator()(int h, int i, int j)
Definition: libmints/matrix.h:1048
const int & nirrep() const
Returns the number of irreps.
Definition: libmints/matrix.h:576
const std::string & name() const
Definition: libmints/matrix.h:529
#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:579
double ** operator[](int i)
Definition: libmints/matrix.h:1045
int max_ncol() const
Returns the column size of the largest block.
Definition: libmints/matrix.h:601
Definition: psio.hpp:195
Definition: libmints/matrix.h:51