SPI
SLEPc PETSc Interface is used for rapid development and intuitive matrix operations similar to MatLab or NumPy style of operations
SPIMat.hpp
Go to the documentation of this file.
1 #ifndef SPIMAT_H
2 #define SPIMAT_H
3 #include <iostream>
4 #include <vector>
5 #include <petscksp.h>
6 #include <petscmath.h>
7 #include <slepceps.h>
8 #include <slepcpep.h>
9 #include <string>
10 #include <tuple>
11 #include "SPIVec.hpp"
12 
13 namespace SPI{
14  template <class T>
15  using Block2D = std::vector<std::vector<T>>;
16 
17  struct SPIMat{
18  PetscInt rows;
19  PetscInt cols;
20 
21  // Constructors
22  SPIMat(std::string _name="SPIMat"); // constructor with no arguments (no initialization)
23  SPIMat(const std::vector<SPIVec> &A, std::string _name="SPIMat"); // constructor using a vector of column vectors
24  SPIMat(const SPIMat &A, std::string _name="SPIMat"); // constructor using another SPIMat
25  SPIMat(PetscInt rowscols, std::string _name="SPIMat"); // constructor with one arguement to make square matrix
26  SPIMat(PetscInt rowsm, PetscInt colsn, std::string _name="SPIMat"); // constructor of rectangular matrix
27 
28  Mat mat;
29  PetscErrorCode ierr;
30 
31  // flags
32  PetscBool flag_init=PETSC_FALSE;
33  std::string name;
34 
35  PetscInt Init(PetscInt m,PetscInt n, std::string name="SPIMat"); // initialize the matrix of size m by n
36  SPIMat& set(PetscInt m, PetscInt n,const PetscScalar v); // set a scalar value at position row m and column n
37  SPIMat& set_col(const PetscInt col,const SPIVec &v); // set a column into a matrix
38  SPIMat& add(PetscInt m, PetscInt n,const PetscScalar v); // add a scalar value at position row m and column n
39  // () operators
40  PetscScalar operator()(PetscInt m, PetscInt n, PetscBool global=PETSC_FALSE); // get local value at row m, column n
41  PetscScalar operator()(PetscInt m, PetscInt n, PetscBool global=PETSC_FALSE) const; // get local value at row m, column n
42  SPIMat& operator()(PetscInt m, PetscInt n,const PetscScalar v); // set operator the same as set function
43  SPIMat& operator()(PetscInt m, PetscInt n,const double v); // set operator the same as set function
44  SPIMat& operator()(PetscInt m, PetscInt n,const int v); // set operator the same as set function
45  SPIMat& operator()(PetscInt m, PetscInt n,const SPIMat &Asub, InsertMode addv=INSERT_VALUES); // set submatrix into matrix at row m, col n
46  SPIMat& operator()(); // assemble the matrix
47  // +- operators
48  SPIMat& operator+=(const SPIMat &X); // MatAXPY, Y = 1.*X + Y operation
49  SPIMat& axpy(const PetscScalar a, const SPIMat &X); // MatAXPY function call to add a*X to the current mat
50  SPIMat operator+(const SPIMat &X); // Y + X operation
51  SPIMat& operator-=(const SPIMat &X); // Y = -1.*X + Y operation
52  SPIMat operator-(const SPIMat &X); // Y - X operation
53  SPIMat operator-() const; // -X operation
54  // * operators
55  SPIMat operator*(const PetscScalar a); // Y*a operation
56  SPIMat operator*(const double a); // Y*a operation
57  SPIVec operator*(const SPIVec &x); // A*x operation to return a vector
58  SPIMat& operator*=(const PetscScalar a); // Y = Y*a operation
59  SPIMat& operator*=(const double a); // Y = Y*a operation
60  SPIMat& operator/=(const PetscScalar a); // Y = Y/a operation
61  SPIMat operator/(const PetscScalar a); // Z = Y/a operation
62  SPIMat operator/(const SPIMat &A); // Z = Y/A elementwise operation
63  SPIMat operator*(const SPIMat &A); // Y*A operation
64  // = operator
65  SPIMat& operator=(const SPIMat &A); // Y=X with initialization of Y using MatConvert
66  // overload % for element wise multiplication
67  //SPIMat operator%(SPIMat A);
68  // Transpose functions
69  PetscInt T(SPIMat &A); // A = Transpose(*this.mat) operation with initialization of A
70  SPIMat& T(); // Transpose the current mat
71  // conjugate
72  PetscInt H(SPIMat &A); // A = Hermitian Transpose(*this.mat) operation with initialization of A (tranpose and complex conjugate)
73  SPIMat& H(); // Hermitian Transpose the current mat
74  SPIVec H(const SPIVec &q); // Hermitian Transpose and multiply by vector
75  SPIMat& conj(); // elemenwise conjugate current matrix
76  SPIMat& real(); // take the real part of the matrix (alters current matrix)
77  SPIVec diag(); // get diagonal of matrix
78  SPIMat& zero_row(const PetscInt row); // zero a row
79  SPIMat& eye_row(const PetscInt row); // zero a row and put 1 in the diagonal entry
80  SPIMat& zero_row_full(const PetscInt row); // zero a row
81  SPIMat& zero_rows(std::vector<PetscInt> rows); // zero every row
82  SPIMat& eye_rows(std::vector<PetscInt> rows); // zero every row
83  SPIVec col(const PetscInt i); // get column vector using MatGetColumnVector(mat,vec,i);
84  PetscInt print(); // print mat to screen using PETSC_VIEWER_STDOUT_WORLD
85 
86  ~SPIMat(); // destructor to delete memory
87 
88  };
89  SPIMat operator*(const PetscScalar a, const SPIMat A); // a*A operation to be equivalent to A*a
90  SPIMat operator*(const SPIMat A, const PetscScalar a); // A*a operation to be equivalent to A*a
91  SPIVec operator/(const SPIVec &b, const SPIMat &A); // Solve linear system, Ax=b using b/A notation
92  SPIMat operator^(const PetscScalar a, const SPIMat &A); // Y = a^A operation
93  //SPIVec solve(const SPIVec &b, const SPIMat &A); // Solve linear system, Ax=b using solve(A,b) notation
94  SPIVec solve(const SPIMat &A, const SPIVec &b); // Solve linear system, Ax=b using solve(A,b) notation
95  SPIMat eye(const PetscInt n); // create, form, and return identity matrix of size n
96  SPIMat inv(const SPIMat &A); // get inverse of matrix by solving A*Ainv = B using MatMatSolve
97  SPIMat zeros(const PetscInt m,const PetscInt n); // create, form, and return zero matrix of size mxn
98  SPIMat diag(const SPIVec &diag,const PetscInt k=0); // set diagonal of matrix
99  SPIMat kron(const SPIMat &A, const SPIMat &B); // set kronecker inner product of two matrices
100  std::tuple<std::vector<PetscReal>,std::vector<SPIVec>,std::vector<SPIVec>> svd(const SPIMat &A); // solve general SVD problem of A = U*E*V^H
101  SPIVec lstsq(const SPIMat &A, SPIVec &y); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector)
102  SPIVec lstsq(const std::vector<SPIVec> &A, SPIVec &y); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector)
103  std::tuple<PetscScalar,SPIVec,SPIVec> eig(const SPIMat &A, const SPIMat &B, const PetscScalar target,const PetscReal tol=-1,const PetscInt max_iter=-1); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector)
104  std::tuple<PetscScalar,SPIVec> eig_right(const SPIMat &A, const SPIMat &B, const PetscScalar target,const PetscReal tol=-1,const PetscInt max_iter=-1); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector)
105  std::tuple<PetscScalar,SPIVec, SPIVec> eig_init(const SPIMat &A, const SPIMat &B, const PetscScalar target,const SPIVec &ql, const SPIVec &qr, PetscReal tol=-1,const PetscInt max_iter=-1); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector) using initial subspace from q
106  std::tuple<PetscScalar, SPIVec> eig_init_right(const SPIMat &A, const SPIMat &B, const PetscScalar target, const SPIVec &qr, PetscReal tol=-1,const PetscInt max_iter=-1); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector) using initial subspace from q
107  std::tuple<std::vector<PetscScalar>, std::vector<SPIVec>> eig_init_rights(const SPIMat &A, const SPIMat &B, const std::vector<PetscScalar> targets, const std::vector<SPIVec> &qrs, PetscReal tol=-1,const PetscInt max_iter=-1); // solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector) using initial subspace from q
108  std::tuple<PetscScalar,SPIVec> polyeig(const std::vector<SPIMat> &As, const PetscScalar target,const PetscReal tol=-1.,const PetscInt max_iter=-1); // solve general polynomial eigenvalue problem of (A0 + A1x + A2x^2 ...) = 0 and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector)
109  std::tuple<PetscScalar,SPIVec> polyeig_init(const std::vector<SPIMat> &As, const PetscScalar target, const SPIVec &qr, const PetscReal tol=-1.,const PetscInt max_iter=-1); // solve general polynomial eigenvalue problem of (A0 + A1x + A2x^2 ...) = 0 and return a tuple of tie(PetscScalar alpha, SPIVec eig_vector) using initial subspace from q
110  //SPIMat block(const SPIMat Blocks[], const PetscInt rows,const PetscInt cols); // set block matrices using an input array of size rows*cols. Fills rows first
111  //SPIMat block(const std::vector<std::vector<SPIMat>> Blocks); // set block matrices using an input array of size rows*cols. Fills rows first
112  SPIMat block(const Block2D<SPIMat> Blocks); // set block matrices using an input array of size rows*cols. Fills rows first
113  std::tuple<SPIMat,SPIMat> meshgrid(SPIVec &x, SPIVec &y); // create meshgrid from two grids using ij indexing. i.e. X(i,j) = x(i) and Y(i,j) = y(j)
114  PetscInt save(const SPIMat &A, const std::string filename); // save matrix to filename to binary format
115  PetscInt save(const std::vector<SPIMat> &A, const std::string filename); // save matrix to filename to binary format
116  PetscInt load(SPIMat &A, const std::string filename); // load matrix to filename from binary format
117  PetscInt load(std::vector<SPIMat> &A, const std::string filename); // load matrix to filename from binary format
118  PetscInt draw(const SPIMat &A); // draw nonzero structure and wait at command line input
119  template <class T> SPIMat _Function_on_each_element( T (*f)(T const&), const SPIMat &A); // function handle template for operations
120  SPIMat sin(const SPIMat &A); // sin of matrix
121  SPIMat cos(const SPIMat &A); // cos of matrix
122  SPIMat acos(const SPIMat &A); // cos of matrix
123  SPIMat tan(const SPIMat &A); // tan of matrix
124  SPIMat abs(const SPIMat &A); // abs of matrix
125  SPIMat operator%(const SPIMat &A,const SPIMat &B); // A*B pointwise operation
126  SPIMat orthogonalize(const std::vector<SPIVec> &x); // create orthonormal basis from array of vectors
127 }
128 
129 
130 
131 #endif // SPIMAT_H
SPI::SPIMat::flag_init
PetscBool flag_init
flag if it has been initialized
Definition: SPIMat.hpp:32
SPI::operator^
SPIMat operator^(const PetscScalar a, const SPIMat &A)
Y=a^A operation.
Definition: SPIMat.cpp:578
SPI::sin
SPIMat sin(const SPIMat &A)
take the sin of each element in a matrix
Definition: SPIMat.cpp:2025
SPI::polyeig
std::tuple< PetscScalar, SPIVec > polyeig(const std::vector< SPIMat > &As, const PetscScalar target, const PetscReal tol=-1., const PetscInt max_iter=-1)
solve general polynomial eigenvalue problem of (A0 + A1*alpha + A2*alpha^2 + ...)*x = 0 and return a ...
Definition: SPIMat.cpp:1686
SPI::polyeig_init
std::tuple< PetscScalar, SPIVec > polyeig_init(const std::vector< SPIMat > &As, const PetscScalar target, const SPIVec &qr, const PetscReal tol=-1., const PetscInt max_iter=-1)
solve general polynomial eigenvalue problem of (A0 + A1*alpha + A2*alpha^2 + ...)*x = 0 and return a ...
Definition: SPIMat.cpp:1754
SPI::SPIMat::print
PetscInt print()
print mat to screen using PETSC_VIEWER_STDOUT_WORLD
Definition: SPIMat.cpp:498
SPI::SPIMat::rows
PetscInt rows
number of rows in mat
Definition: SPIMat.hpp:18
SPI::SPIMat::zero_rows
SPIMat & zero_rows(std::vector< PetscInt > rows)
set rows to zero
Definition: SPIMat.cpp:475
SPI::SPIMat::Init
PetscInt Init(PetscInt m, PetscInt n, std::string name="SPIMat")
Definition: SPIMat.cpp:55
SPI::SPIMat::operator+=
SPIMat & operator+=(const SPIMat &X)
MatAXPY, Y = 1.*X + Y operation.
Definition: SPIMat.cpp:225
SPI::SPIMat::mat
Mat mat
petsc Mat data
Definition: SPIMat.hpp:28
SPI::eig_right
std::tuple< PetscScalar, SPIVec > eig_right(const SPIMat &A, const SPIMat &B, const PetscScalar target, const PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1072
SPI::lstsq
SPIVec lstsq(const SPIMat &A, SPIVec &y)
solve least squares problem of A*x=y for skinny A matrix using SVD
Definition: SPIMat.cpp:905
SPI::Block2D
std::vector< std::vector< T > > Block2D
Definition: SPIMat.hpp:15
SPI::SPIMat::name
std::string name
Matrix name.
Definition: SPIMat.hpp:33
SPI::SPIMat::operator/=
SPIMat & operator/=(const PetscScalar a)
Y = Y/a operation.
Definition: SPIMat.cpp:327
SPI::SPIVec
Definition: SPIVec.hpp:13
SPI::SPIMat::operator-=
SPIMat & operator-=(const SPIMat &X)
Y = -1.*X + Y operation.
Definition: SPIMat.cpp:256
SPI::SPIMat::operator()
SPIMat & operator()()
assmelbe the matrix
Definition: SPIMat.cpp:218
SPI
Definition: SPIbaseflow.hpp:16
SPI::save
PetscInt save(const SPIMat &A, const std::string filename)
save matrix to filename in binary format (see Petsc documentation for format ) Format is (from Petsc ...
Definition: SPIMat.cpp:1924
SPI::SPIMat::ierr
PetscErrorCode ierr
ierr for various routines and operators
Definition: SPIMat.hpp:29
SPI::SPIMat::add
SPIMat & add(PetscInt m, PetscInt n, const PetscScalar v)
add a scalar value at position row m and column n
Definition: SPIMat.cpp:97
SPI::kron
SPIMat kron(const SPIMat &A, const SPIMat &B)
set kronecker inner product of two matrices
Definition: SPIMat.cpp:780
SPI::SPIMat::SPIMat
SPIMat(std::string _name="SPIMat")
constructor with no arguments (no initialization)
Definition: SPIMat.cpp:13
SPI::SPIMat::axpy
SPIMat & axpy(const PetscScalar a, const SPIMat &X)
MatAXPY function call to add a*X to the current mat.
Definition: SPIMat.cpp:232
SPI::SPIMat::operator-
SPIMat operator-() const
-X operation
Definition: SPIMat.cpp:278
SPI::SPIMat::zero_row_full
SPIMat & zero_row_full(const PetscInt row)
set a row to zero using dense format
Definition: SPIMat.cpp:465
SPI::SPIMat::diag
SPIVec diag()
get diagonal of matrix
Definition: SPIMat.cpp:445
SPI::_Function_on_each_element
SPIMat _Function_on_each_element(T(*f)(T const &), const SPIMat &A)
take the function of each element in a matrix, e.g. (*f)(A(i,j)) for each i,j
Definition: SPIMat.cpp:2011
SPI::SPIMat::operator*=
SPIMat & operator*=(const PetscScalar a)
Y = Y*a operation.
Definition: SPIMat.cpp:320
SPI::eye
SPIMat eye(const PetscInt n)
create, form, and return identity matrix of size n
Definition: SPIMat.cpp:697
SPI::inv
SPIMat inv(const SPIMat &A)
Solve linear system A*Ainv=B using MatMatSolve.
Definition: SPIMat.cpp:603
SPI::orthogonalize
std::vector< SPIVec > orthogonalize(std::vector< SPIVec > &x, SPIgrid1D &grid)
Definition: SPIgrid.cpp:1323
SPI::draw
PetscInt draw(const SPIMat &A)
draw nonzero structure of matrix
Definition: SPIMat.cpp:1994
SPI::SPIMat::~SPIMat
~SPIMat()
destructor to delete memory
Definition: SPIMat.cpp:508
SPI::tan
SPIMat tan(const SPIMat &A)
take the tan of each element in a matrix
Definition: SPIMat.cpp:2031
SPI::diag
SPIMat diag(const SPIVec &diag, const PetscInt k=0)
set diagonal of matrix
Definition: SPIMat.cpp:728
SPI::SPIMat::T
SPIMat & T()
Transpose the current mat.
Definition: SPIMat.cpp:405
SPI::SPIMat::operator+
SPIMat operator+(const SPIMat &X)
Y + X operation.
Definition: SPIMat.cpp:241
SPI::abs
SPIMat abs(const SPIMat &A)
take the abs of each element in a matrix
Definition: SPIMat.cpp:2045
SPI::SPIMat::eye_row
SPIMat & eye_row(const PetscInt row)
set a row to zero and set 1 in diagonal entry
Definition: SPIMat.cpp:458
SPIVec.hpp
SPI::eig_init_rights
std::tuple< std::vector< PetscScalar >, std::vector< SPIVec > > eig_init_rights(const SPIMat &A, const SPIMat &B, const std::vector< PetscScalar > targets, const std::vector< SPIVec > &qrs, PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1521
SPI::solve
SPIVec solve(const SPIMat &A, const SPIVec &b)
Solve linear system, Ax=b using solve(A,b) notation.
Definition: SPIMat.cpp:596
SPI::SPIMat::operator*
SPIMat operator*(const PetscScalar a)
Y*a operation.
Definition: SPIMat.cpp:285
SPI::eig_init
std::tuple< PetscScalar, SPIVec, SPIVec > eig_init(const SPIMat &A, const SPIMat &B, const PetscScalar target, const SPIVec &ql, const SPIVec &qr, PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1217
SPI::SPIMat::real
SPIMat & real()
take the real part of the vector
Definition: SPIMat.cpp:440
SPI::SPIMat::H
SPIMat & H()
Hermitian Transpose the current mat.
Definition: SPIMat.cpp:422
SPI::load
PetscInt load(SPIMat &A, const std::string filename)
load matrix from filename from binary format (works with save(SPIMat,std::string) function
Definition: SPIMat.cpp:1966
SPI::SPIMat
Definition: SPIMat.hpp:17
SPI::operator%
SPIMat operator%(const SPIMat &A, const SPIMat &B)
take the elementwise multiplication of each element in a matrix
Definition: SPIMat.cpp:2033
SPI::block
SPIMat block(const Block2D< SPIMat > Blocks)
set block matrices using an input array of size rows*cols. Fills rows first
Definition: SPIMat.cpp:1857
SPI::SPIMat::conj
SPIMat & conj()
elemenwise conjugate current matrix
Definition: SPIMat.cpp:435
SPI::SPIMat::operator/
SPIMat operator/(const PetscScalar a)
Z = Y/a operation.
Definition: SPIMat.cpp:334
SPI::acos
SPIMat acos(const SPIMat &A)
take the arccos of each element in a matrix
Definition: SPIMat.cpp:2029
SPI::SPIMat::set
SPIMat & set(PetscInt m, PetscInt n, const PetscScalar v)
Definition: SPIMat.cpp:77
SPI::operator*
SPIMat operator*(const PetscScalar a, const SPIMat A)
a*A operation to be equivalent to A*a
Definition: SPIMat.cpp:517
SPI::eig_init_right
std::tuple< PetscScalar, SPIVec > eig_init_right(const SPIMat &A, const SPIMat &B, const PetscScalar target, const SPIVec &qr, PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:1368
SPI::SPIMat::operator=
SPIMat & operator=(const SPIMat &A)
Y=X with initialization of Y using MatConvert.
Definition: SPIMat.cpp:367
SPI::meshgrid
std::tuple< SPIMat, SPIMat > meshgrid(SPIVec &x, SPIVec &y)
Definition: SPIMat.cpp:1894
SPI::SPIMat::col
SPIVec col(const PetscInt i)
get column vector
Definition: SPIMat.cpp:489
SPI::eig
std::tuple< PetscScalar, SPIVec, SPIVec > eig(const SPIMat &A, const SPIMat &B, const PetscScalar target, const PetscReal tol=-1, const PetscInt max_iter=-1)
solve general eigenvalue problem of Ax = kBx and return a tuple of tie(PetscScalar alpha,...
Definition: SPIMat.cpp:928
SPI::cos
SPIMat cos(const SPIMat &A)
take the cos of each element in a matrix
Definition: SPIMat.cpp:2027
SPI::SPIMat::zero_row
SPIMat & zero_row(const PetscInt row)
set a row to zero
Definition: SPIMat.cpp:451
SPI::operator/
SPIVec operator/(const SPIVec &b, const SPIMat &A)
Solve linear system, Ax=b using b/A notation.
Definition: SPIMat.cpp:535
SPI::SPIMat::set_col
SPIMat & set_col(const PetscInt col, const SPIVec &v)
Definition: SPIMat.cpp:86
SPI::svd
std::tuple< std::vector< PetscReal >, std::vector< SPIVec >, std::vector< SPIVec > > svd(const SPIMat &A)
solve SVD problem of A=U*E*V^H for skinny A matrix
Definition: SPIMat.cpp:851
SPI::zeros
SPIMat zeros(const PetscInt m, const PetscInt n)
create, form, and return zeros matrix of size mxn
Definition: SPIMat.cpp:708
SPI::SPIMat::cols
PetscInt cols
number of columns in mat
Definition: SPIMat.hpp:19
SPI::SPIMat::eye_rows
SPIMat & eye_rows(std::vector< PetscInt > rows)
set rows to zero and set main diagonal to 1
Definition: SPIMat.cpp:482
Generated by   doxygen 1.8.17