SPI
SLEPc PETSc Interface is used for rapid development and intuitive matrix operations similar to MatLab or NumPy style of operations
SPIMat.cpp
Go to the documentation of this file.
1 #include "SPIMat.hpp"
2 #include "SPIprint.hpp"
3 #include <petscviewerhdf5.h>
4 #include <petscmath.h>
5 #include <slepcsvd.h>
6 // use 1 to use the GPU and 0 or anything else to only use CPU and MPI. Be sure this matches what is in SPIVec.cpp
7 #define USE_GPU 0
8 
9 namespace SPI{
10 
11  // constructors
13  SPIMat::SPIMat(
14  std::string _name
15  ){name=_name; }
17  SPIMat::SPIMat(
18  const SPIMat &A,
19  std::string _name
20  ){
21  name=_name;
22  (*this) = A;
23  }
25  SPIMat::SPIMat(
26  const std::vector<SPIVec> &A,
27  std::string _name
28  ){
29  PetscInt m=A[0].rows;
30  PetscInt n=A.size();
31  Init(m,n,_name); // initialize the matrix
32  for(PetscInt j=0; j<n; ++j){
33  set_col(j,A[j]);
34  }
35  (*this)();
36  }
38  SPIMat::SPIMat(
39  PetscInt rowscols,
40  std::string _name
41  ){
42  Init(rowscols,rowscols,_name);
43  }
45  SPIMat::SPIMat(
46  PetscInt rowsm,
47  PetscInt colsn,
48  std::string _name
49  ){
50  Init(rowsm,colsn,_name);
51  }
52 
53  // Initialize matrix
55  PetscInt SPIMat::Init(
56  PetscInt m,
57  PetscInt n,
58  std::string _name
59  ){
60  name=_name;
61  rows=m;
62  cols=n;
63  ierr = MatCreate(PETSC_COMM_WORLD,&mat);CHKERRXX(ierr);
64  ierr = MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRXX(ierr);
65  //ierr = MatSetFromOptions(mat);CHKERRXX(ierr);
66 #if UUSE_GPU == 1
67  ierr = MatSetType(mat,MATMPIAIJCUSPARSE);CHKERRXX(ierr);
68 #else
69  ierr = MatSetType(mat,MATMPIAIJ);CHKERRXX(ierr);
70 #endif
71  ierr = MatSetUp(mat);CHKERRXX(ierr);
72  flag_init=PETSC_TRUE;
73  return 0;
74  }
75 
77  SPIMat& SPIMat::set(
78  PetscInt m,
79  PetscInt n,
80  const PetscScalar v
81  ){
82  ierr = MatSetValue(mat,m,n,v,INSERT_VALUES);CHKERRXX(ierr);
83  return (*this);
84  }
86  SPIMat& SPIMat::set_col(
87  const PetscInt col,
88  const SPIVec &v
89  ){
90  for(PetscInt i=0; i<rows; ++i){
91  (*this)(i,col,v(i,PETSC_TRUE));
92  }
93  //(*this)();
94  return (*this);
95  }
97  SPIMat& SPIMat::add(
98  PetscInt m,
99  PetscInt n,
100  const PetscScalar v
101  ){
102  ierr = MatSetValue(mat,m,n,v,ADD_VALUES);CHKERRXX(ierr);
103  return (*this);
104  }
106  PetscScalar SPIMat::operator()(
107  PetscInt m,
108  PetscInt n,
109  PetscBool global
110  )const{
111  PetscErrorCode ierr2;
112  PetscScalar v,v_global=0.;
113  PetscInt low,high;
114  ierr2 = MatGetOwnershipRange(mat,&low, &high);CHKERRXX(ierr2);
115  if ((low<=m) && (m<high)){
116  ierr2 = MatGetValues(mat,1,&m, 1,&n, &v);
117  }
118  if (global){ // if broadcast to all processors
119  MPIU_Allreduce(&v,&v_global,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
120  }
121  else{
122  v_global=v; // return local value
123  }
124  return v_global;
125  }
126 
127  // overloaded operators, get
129  PetscScalar SPIMat::operator()(
130  PetscInt m,
131  PetscInt n,
132  PetscBool global
133  ){
134  PetscScalar v,v_global=0.;
135  PetscInt low,high;
136  ierr = MatGetOwnershipRange(mat,&low, &high);CHKERRXX(ierr);
137  if ((low<=m) && (m<high)){
138  ierr = MatGetValues(mat,1,&m, 1,&n, &v);
139  }
140  if (global){ // if broadcast to all processors
141  MPIU_Allreduce(&v,&v_global,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
142  }
143  else{
144  v_global=v; // return local value
145  }
146  return v_global;
147  }
148  // overloaded operator, set
150  SPIMat& SPIMat::operator()(
151  PetscInt m,
152  PetscInt n,
153  const PetscScalar v
154  ){
155  ierr = MatSetValue(mat,m,n,v,INSERT_VALUES);CHKERRXX(ierr);
156  //(*this)(); // assemble after every insertion
157  return (*this);
158  }
160  SPIMat& SPIMat::operator()(
161  PetscInt m,
162  PetscInt n,
163  const double v
164  ){
165  //ierr = (*this)(m,n,(PetscScalar)v);CHKERRXX(ierr);
166  return (*this)(m,n,(PetscScalar)(v+0.0*PETSC_i));
167  }
169  SPIMat& SPIMat::operator()(
170  PetscInt m,
171  PetscInt n,
172  const int v
173  ){
174  //ierr = (*this)(m,n,(PetscScalar)v);CHKERRXX(ierr);
175  return (*this)(m,n,(PetscScalar)((double)v+0.0*PETSC_i));
176  }
177 
178  // overloaded operator, set
180  SPIMat& SPIMat::operator()(
181  PetscInt m,
182  PetscInt n,
183  const SPIMat &Asub,
184  InsertMode addv
185  ){
186  //InsertMode addv = INSERT_VALUES;
187  PetscInt rowoffset = m;
188  PetscInt coloffset = n;
189  PetscInt nsub = Asub.rows;
190  PetscInt Isubstart,Isubend;
191  PetscInt ncols;
192  const PetscInt *cols;
193  const PetscScalar *vals;
194 
195  // get ranges for each matrix
196  MatGetOwnershipRange(Asub.mat,&Isubstart,&Isubend);
197 
198  PetscErrorCode ierr;
199  for (PetscInt i=Isubstart; i<Isubend && i<nsub; i++){
200  ierr = MatGetRow(Asub.mat,i,&ncols,&cols,&vals);CHKERRXX(ierr);
201  PetscInt offcols[ncols];
202  PetscScalar avals[ncols];
203  for (PetscInt j=0; j<ncols; j++) {
204  offcols[j] = cols[j]+coloffset;
205  avals[j] = vals[j];
206  }
207  //printScalar(vals,ncols);
208  PetscInt rowoffseti = i+rowoffset;
209  ierr = MatSetValues(mat,1,&rowoffseti,ncols,offcols,avals,addv);CHKERRXX(ierr);
210 
211  MatRestoreRow(Asub.mat,i,&ncols,&cols,&vals);
212  }
213  //(*this)();//assemble
214  return (*this);
215  }
216  // overloaded operator, assemble
218  SPIMat& SPIMat::operator()(){
219  ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRXX(ierr);
220  ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRXX(ierr);
221  return (*this);
222  }
223  // overloaded operator, MatAXPY
225  SPIMat& SPIMat::operator+=(
226  const SPIMat &X
227  ){
228  ierr = MatAXPY(this->mat,1.,X.mat,DIFFERENT_NONZERO_PATTERN);CHKERRXX(ierr);
229  return *this;
230  }
232  SPIMat& SPIMat::axpy(
233  const PetscScalar a,
234  const SPIMat &X
235  ){
236  ierr = MatAXPY(this->mat,a,X.mat,DIFFERENT_NONZERO_PATTERN);CHKERRXX(ierr);
237  return (*this);
238  }
239  // overloaded operator, MatAXPY
241  SPIMat SPIMat::operator+(
242  const SPIMat &X
243  ){
244  SPIMat A;
245  A=*this;
246  ierr = MatAXPY(A.mat,1.,X.mat,DIFFERENT_NONZERO_PATTERN);CHKERRXX(ierr);
247 #if USE_GPU == 1
248  ierr = MatSetType(A.mat,MATMPIAIJCUSPARSE);CHKERRXX(ierr);
249 #else
250  ierr = MatSetType(A.mat,MATMPIAIJ);CHKERRXX(ierr);
251 #endif
252  return A;
253  }
254  // overloaded operator, MatAXPY
256  SPIMat& SPIMat::operator-=(
257  const SPIMat &X
258  ){
259  ierr = MatAXPY(this->mat,-1.,X.mat,DIFFERENT_NONZERO_PATTERN);CHKERRXX(ierr);
260  return *this;
261  }
262  // overloaded operator, MatAXPY
264  SPIMat SPIMat::operator-(
265  const SPIMat &X
266  ){
267  SPIMat A;
268  A=*this;
269  ierr = MatAXPY(A.mat,-1.,X.mat,DIFFERENT_NONZERO_PATTERN);CHKERRXX(ierr);
270 #if USE_GPU == 1
271  ierr = MatSetType(A.mat,MATMPIAIJCUSPARSE);CHKERRXX(ierr);
272 #else
273  ierr = MatSetType(A.mat,MATMPIAIJ);CHKERRXX(ierr);
274 #endif
275  return A;
276  }
278  SPIMat SPIMat::operator-() const {
279  SPIMat A;
280  A=-1.*(*this);
281  return A;
282  }
283  // overload operator, scale with scalar
285  SPIMat SPIMat::operator*(
286  const PetscScalar a
287  ){
288  SPIMat A;
289  A=*this;
290  ierr = MatScale(A.mat,a);CHKERRXX(ierr);
291  return A;
292  }
294  SPIMat SPIMat::operator*(
295  const double a
296  ){
297  SPIMat A;
298  A=*this;
299  ierr = MatScale(A.mat,a);CHKERRXX(ierr);
300  return A;
301  }
303  SPIVec SPIMat::operator*(
304  const SPIVec &x
305  ){
306  //SPIVec b(x.rows); // only works for square matrices
307  SPIVec b(rows); // works for non-square matrices
308  ierr = MatMult(mat,x.vec,b.vec);CHKERRXX(ierr);
309  return b;
310  }
311  // overload operator, scale with scalar
313  SPIMat& SPIMat::operator*=(
314  const double a
315  ){
316  ierr = MatScale(this->mat,(PetscScalar)(a+0.*PETSC_i));CHKERRXX(ierr);
317  return *this;
318  }
320  SPIMat& SPIMat::operator*=(
321  const PetscScalar a
322  ){
323  ierr = MatScale(this->mat,a);CHKERRXX(ierr);
324  return *this;
325  }
327  SPIMat& SPIMat::operator/=(
328  const PetscScalar a
329  ){
330  ierr = MatScale(this->mat,1./a);CHKERRXX(ierr);
331  return *this;
332  }
334  SPIMat SPIMat::operator/(
335  const PetscScalar a
336  ){
337  return (1./a)*(*this);
338  }
340  SPIMat SPIMat::operator/(
341  const SPIMat &A
342  ){
343  SPIMat Z(A.rows,A.cols);
344  for(PetscInt i=0; i<Z.rows; ++i){
345  for(PetscInt j=0; j<Z.cols; ++j){
346  Z(i,j,((*this)(i,j,PETSC_TRUE))/A(i,j,PETSC_TRUE));
347  }
348  }
349  Z();
350  return Z;
351  }
352  // overload operator, matrix multiply
354  SPIMat SPIMat::operator*(
355  const SPIMat &A
356  ){
357  SPIMat C;
358  C.rows=rows;
359  C.cols=A.cols;
360  ierr = MatMatMult(mat,A.mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C.mat); CHKERRXX(ierr);
361  C.flag_init=PETSC_TRUE;
362  //ierr = MatSetType(C.mat,MATMPIAIJ);CHKERRXX(ierr);
363  return C;
364  }
365  // overload operator, copy and initialize
367  SPIMat& SPIMat::operator=(
368  const SPIMat &A
369  ){
370  if(flag_init){
371  this->~SPIMat();
372  //ierr = MatCopy(A.mat,mat,DIFFERENT_NONZERO_PATTERN);CHKERRXX(ierr);
373  }
374  //else
375  rows=A.rows;
376  cols=A.cols;
377  ierr = MatConvert(A.mat,MATSAME,MAT_INITIAL_MATRIX,&mat);CHKERRXX(ierr);
378 #if USE_GPU == 1
379  ierr = MatSetType(mat,MATMPIAIJCUSPARSE);CHKERRXX(ierr);
380 #else
381  ierr = MatSetType(mat,MATMPIAIJ);CHKERRXX(ierr);
382 #endif
383  flag_init=PETSC_TRUE;
384  //
385  return *this;
386  }
387  // overload % for element wise multiplication
388  //SPIMat operator%(SPIMat A)
389  //return *this;
390  //]
392  PetscInt SPIMat::T(
393  SPIMat &A
394  ){
395  //ierr = MatTranspose(mat,MAT_INITIAL_MATRIX,&A.mat);CHKERRXX(ierr);
396  //A();
397  //A.Init(cols,rows);
398  ierr = MatTranspose(mat,MAT_INITIAL_MATRIX,&A.mat);CHKERRXX(ierr);
399  A.rows=this->cols;
400  A.cols=this->rows;
401  A.flag_init=PETSC_TRUE;
402  return 0;
403  }
405  SPIMat& SPIMat::T(){
406  ierr = MatTranspose(mat,MAT_INPLACE_MATRIX,&mat);CHKERRXX(ierr);
407  return (*this);
408  //SPIMat T1;
409  //ierr = MatCreateTranspose(mat,&T1.mat);CHKERRXX(ierr);
410  //SPIMat T1;
411  //ierr = MatTranspose(mat,MAT_INITIAL_MATRIX,&T1.mat);CHKERRXX(ierr);
412  //return T1;
413  }
415  PetscInt SPIMat::H(
416  SPIMat &A
417  ){ // A = Hermitian Transpose(*this.mat) operation with initialization of A (tranpose and complex conjugate)
418  ierr = MatHermitianTranspose(mat,MAT_INITIAL_MATRIX,&A.mat);CHKERRXX(ierr);
419  return 0;
420  }
422  SPIMat& SPIMat::H(){ // Hermitian Transpose the current mat
423  ierr = MatHermitianTranspose(mat,MAT_INPLACE_MATRIX,&mat);CHKERRXX(ierr);
424  return (*this);
425  }
427  SPIVec SPIMat::H(
428  const SPIVec &x
429  ){ // Hermitian Transpose the current mat
430  SPIVec y(cols);
431  ierr = MatMultHermitianTranspose(mat,x.vec,y.vec);CHKERRXX(ierr);
432  return y;
433  }
435  SPIMat& SPIMat::conj(){
436  ierr = MatConjugate(mat);CHKERRXX(ierr);
437  return (*this);
438  }
440  SPIMat& SPIMat::real(){
441  ierr = MatRealPart(mat); CHKERRXX(ierr);
442  return (*this);
443  }
445  SPIVec SPIMat::diag(){ // get diagonal of matrix
446  SPIVec d(rows);
447  ierr = MatGetDiagonal(mat,d.vec); CHKERRXX(ierr);
448  return d;
449  }
451  SPIMat& SPIMat::zero_row(
452  const PetscInt row
453  ){
454  ierr = MatZeroRows(mat,1,&row,0.,0,0); CHKERRXX(ierr);
455  return (*this);
456  }
458  SPIMat& SPIMat::eye_row(
459  const PetscInt row
460  ){
461  ierr = MatZeroRows(mat,1,&row,1.,0,0); CHKERRXX(ierr);
462  return (*this);
463  }
465  SPIMat& SPIMat::zero_row_full(
466  const PetscInt row
467  ){
468  for(PetscInt j=0; j<cols; j++){
469  (*this)(row,j,0.0);
470  }
471  //ierr = MatZeroRows(mat,1,&row,0.,0,0); CHKERRXX(ierr);
472  return (*this);
473  }
475  SPIMat& SPIMat::zero_rows(
476  std::vector<PetscInt> rows
477  ){
478  ierr = MatZeroRows(mat,rows.size(),rows.data(),0.,0,0); CHKERRXX(ierr);
479  return (*this);
480  }
482  SPIMat& SPIMat::eye_rows(
483  std::vector<PetscInt> rows
484  ){
485  ierr = MatZeroRows(mat,rows.size(),rows.data(),1.,0,0); CHKERRXX(ierr);
486  return (*this);
487  }
489  SPIVec SPIMat::col(
490  const PetscInt i
491  ){
492  SPIVec yy(rows);
493  yy.ierr = MatGetColumnVector(mat,yy.vec,i); CHKERRXX(yy.ierr);
494  return yy;
495  }
496  // print matrix to screen
498  PetscInt SPIMat::print(){
499  (*this)();
500  PetscPrintf(PETSC_COMM_WORLD,("\n---------------- "+name+"---start------\n").c_str());
501  SPI::printf("shape = "+std::to_string(this->rows)+" x "+std::to_string(this->cols));
502  ierr = MatView(mat,PETSC_VIEWER_STDOUT_WORLD);CHKERRXX(ierr);
503  PetscPrintf(PETSC_COMM_WORLD,("---------------- "+name+"---done-------\n\n").c_str());
504  return 0;
505  }
506 
508  SPIMat::~SPIMat(){
509  if(flag_init){
510  flag_init=PETSC_FALSE;
511  ierr = MatDestroy(&mat);CHKERRXX(ierr);
512  }
513  }
514 
515  // overload operator, scale with scalar
517  SPIMat operator*(
518  const PetscScalar a,
519  const SPIMat A
520  ){
521  SPIMat B;
522  B=A;
523  B.ierr = MatScale(B.mat,a);CHKERRXX(B.ierr);
524  return B;
525  }
527  SPIMat operator*(const SPIMat A, const PetscScalar a){
528  SPIMat B;
529  B=A;
530  B.ierr = MatScale(B.mat,a);CHKERRXX(B.ierr);
531  return B;
532  }
533  // overload operator, Linear System solve Ax=b
535  SPIVec operator/(
536  const SPIVec &b,
537  const SPIMat &A
538  ){
539  SPIVec x;
540  x.rows=b.rows;
541  KSP ksp; // Linear solver context
542  PetscErrorCode ierr;
543  ierr = VecDuplicate(b.vec,&x.vec);CHKERRXX(ierr);
544 
545  // Create the linear solver and set various options
546  // Create linear solver context
547  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRXX(ierr);
548  // Set operators. Here the matrix that defines the linear system
549  // also serves as the preconditioning matrix.
550  ierr = KSPSetOperators(ksp,A.mat,A.mat);CHKERRXX(ierr);
551  // Set runtime options, e.g.,
552  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> -ksp_type <type> -pc_type <type>
553  //ierr = KSPSetFromOptions(ksp);CHKERRXX(ierr);
554  //ierr = KSPSetType(ksp,KSPPREONLY);CHKERRXX(ierr);
555 
556  PC pc;
557  ierr = KSPGetPC(ksp,&pc);CHKERRXX(ierr);
558  ierr = PCSetType(pc,PCLU);CHKERRXX(ierr);
559  ierr = KSPSetType(ksp,KSPPREONLY);CHKERRXX(ierr);
560  //ierr = PCSetOperators(pc,A.mat,A.mat); CHKERRXX(ierr);
561 
562 
563  // Solve the linear system
564  ierr = KSPSolve(ksp,b.vec,x.vec);CHKERRXX(ierr);
565  x.flag_init = PETSC_TRUE;
566 
567  // output iterations
568  //PetscInt its;
569  //ierr = KSPGetIterationNumber(ksp,&its);CHKERRXX(ierr);
570  //PetscPrintf(PETSC_COMM_WORLD,"KSP Solved in %D iterations \n",its);
571  // Free work space. All PETSc objects should be destroyed when they
572  //set_Vec(x);
573  // are no longer needed.
574  ierr = KSPDestroy(&ksp);CHKERRXX(ierr);
575  return x;
576  }
578  SPIMat operator^(
579  const PetscScalar a,
580  const SPIMat &A
581  ){
582  SPIMat Y(A.rows,A.cols);
583  for(PetscInt i=0; i<Y.rows; ++i){
584  for(PetscInt j=0; j<Y.cols; ++j){
585  Y(i,j,std::pow<PetscReal>(a,A(i,j,PETSC_TRUE)));
586  //std::cout<<"A("<<i<<","<<j<<") = "<<A(i,j,PETSC_TRUE)<<std::endl;
587  //std::cout<<"Y("<<i<<","<<j<<") = "<<Y(i,j,PETSC_TRUE)<<std::endl;
588  //Y();
589  }
590  }
591  Y();
592  //Y.real();
593  return Y;
594  }
596  SPIVec solve(
597  const SPIMat &A,
598  const SPIVec &b
599  ){
600  return b/A;
601  }
603  SPIMat inv(
604  const SPIMat &A
605  ){
606  if(0){ // only works on single processor and with Dense matrices
607  SPIMat Acopy;
608  Acopy.rows=A.rows;
609  Acopy.cols=A.cols;
610  Acopy.name=A.name;
611  Mat Amat;
612  Acopy.ierr = MatCreate(PETSC_COMM_WORLD,&Amat);CHKERRXX(Acopy.ierr);
613  Acopy.ierr = MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,A.rows,A.cols);CHKERRXX(Acopy.ierr);
614  Acopy.ierr = MatSetType(Amat,MATDENSE);CHKERRXX(Acopy.ierr);
615  Acopy.ierr = MatSetUp(Amat);CHKERRXX(Acopy.ierr);
616  Acopy.flag_init=PETSC_TRUE;
617  Acopy.mat = Amat;
618  Acopy(0,0,A);
619  //Acopy.ierr = MatSetType(Acopy.mat,MATMPIDENSE);CHKERRXX(Acopy.ierr);
620  Acopy();CHKERRXX(Acopy.ierr);
621  Acopy.ierr = MatLUFactor(Acopy.mat,NULL,NULL,NULL); CHKERRXX(Acopy.ierr);
622  SPIMat B;
623  Mat Bmat;
624  B.ierr = MatCreate(PETSC_COMM_WORLD,&Bmat);CHKERRXX(B.ierr);
625  B.ierr = MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,A.rows,A.cols);CHKERRXX(B.ierr);
626  B.ierr = MatSetType(Bmat,MATDENSE);CHKERRXX(B.ierr);
627  B.ierr = MatSetUp(Bmat);CHKERRXX(B.ierr);
628  B.mat = Bmat;
629  B(0,0,eye(A.rows));
630  SPIMat Ainv;
631  Mat Ainvmat;
632  Ainv.ierr = MatCreate(PETSC_COMM_WORLD,&Ainvmat);CHKERRXX(Ainv.ierr);
633  Ainv.ierr = MatSetSizes(Ainvmat,PETSC_DECIDE,PETSC_DECIDE,A.rows,A.cols);CHKERRXX(Ainv.ierr);
634  Ainv.ierr = MatSetType(Ainvmat,MATDENSE);CHKERRXX(Ainv.ierr);
635  Ainv.ierr = MatSetUp(Ainvmat);CHKERRXX(Ainv.ierr);
636  Ainv.mat = Ainvmat;
637  Ainv.rows=A.rows;
638  Ainv.cols=A.cols;
639  Ainv.name = A.name + "^-1";
640  //Ainv.ierr = MatSetType(Ainv.mat,MATSEQDENSE);CHKERRXX(Ainv.ierr);
641  //Ainv();
642  Ainv.ierr = MatMatSolve(Acopy.mat,B.mat,Ainv.mat); CHKERRXX(Ainv.ierr);
643  return Ainv;
644  }
645  else{ // solve using ksp multiple times
646  SPIMat B(eye(A.rows));
647  SPIVec b(B.col(0));
648  SPIVec x;
649  x.rows=b.rows;
650  KSP ksp; // Linear solver context
651  PetscErrorCode ierr;
652  ierr = VecDuplicate(b.vec,&x.vec);CHKERRXX(ierr);
653  x.flag_init=PETSC_TRUE;
654 
655  // Create the linear solver and set various options
656  // Create linear solver context
657  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRXX(ierr);
658  // Set operators. Here the matrix that defines the linear system
659  // also serves as the preconditioning matrix.
660  ierr = KSPSetOperators(ksp,A.mat,A.mat);CHKERRXX(ierr);
661  // Set runtime options, e.g.,
662  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> -ksp_type <type> -pc_type <type>
663  //ierr = KSPSetFromOptions(ksp);CHKERRXX(ierr);
664  //ierr = KSPSetType(ksp,KSPPREONLY);CHKERRXX(ierr);
665 
666  PC pc;
667  ierr = KSPGetPC(ksp,&pc);CHKERRXX(ierr);
668  ierr = PCSetType(pc,PCLU);CHKERRXX(ierr);
669  ierr = KSPSetType(ksp,KSPPREONLY);CHKERRXX(ierr);
670  //ierr = PCSetOperators(pc,A.mat,A.mat); CHKERRXX(ierr);
671 
672 
673  // Solve the linear system multiple times to get Ainv
674  SPIMat Ainv(A.rows,A.cols,A.name+"^-1");
675  for(PetscInt j=0; j<A.cols; ++j){
676  ierr = KSPSolve(ksp,B.col(j).vec,x.vec);CHKERRXX(ierr);
677  for(PetscInt i=0; i<A.cols; ++i){
678  PetscScalar value=x(i,PETSC_TRUE);
679  if(value != 0.) Ainv(i,j,value); // try to preserve zero structure
680  }
681  }
682  Ainv();
683 
684  // output iterations
685  //PetscInt its;
686  //ierr = KSPGetIterationNumber(ksp,&its);CHKERRXX(ierr);
687  //PetscPrintf(PETSC_COMM_WORLD,"KSP Solved in %D iterations \n",its);
688  // Free work space. All PETSc objects should be destroyed when they
689  //set_Vec(x);
690  // are no longer needed.
691  ierr = KSPDestroy(&ksp);CHKERRXX(ierr);
692  return Ainv;
693  }
694  }
695  // identity matrix formation
697  SPIMat eye(
698  const PetscInt n
699  ){
700  SPIMat I(n,"I");
701  SPIVec one(n);
702  I.ierr = VecSet(one.vec,1.);CHKERRXX(I.ierr);
703  I.ierr = MatDiagonalSet(I.mat,one.vec,INSERT_VALUES);CHKERRXX(I.ierr);
704 
705  return I;
706  }
708  SPIMat zeros(
709  const PetscInt m,
710  const PetscInt n
711  ){
712  SPIMat O(m,n,"zero");
713  O();
714  //O.ierr = MatZeroEntries(O.mat); CHKERRXX(O.ierr);
715  // set main diagonal to zero... but these didn't quite work.
716  //SPIVec zero(m);
717  //O.ierr = VecSet(zero.vec,1.);CHKERRXX(O.ierr);
718  //zero *= 0.;
719  //O.ierr = MatDiagonalSet(O.mat,zero.vec,INSERT_VALUES);CHKERRXX(O.ierr);
720  //O();
721  //for(PetscInt i=0; i<m; i++)
722  //O(i,i,1.0);
723  //
724  return O;
725  }
726  // diagonal matrix
728  SPIMat diag(
729  const SPIVec &d,
730  const PetscInt k
731  ){ // set diagonal of matrix
732  if(k==0){
733  SPIMat A(d.rows);
734  A.ierr = MatDiagonalSet(A.mat,d.vec,INSERT_VALUES);CHKERRXX(A.ierr);
735  return A;
736  }
737  else if(k>0){
738  PetscInt r0=d.rows, r1=k;
739  PetscInt c0=k, c1=d.rows;
740  SPIMat A00(zeros(r0,c0)), A01(r0,c1),
741  A10(zeros(r1,c0)), A11(zeros(r1,c1));
742  A01.ierr = MatDiagonalSet(A01.mat,d.vec,INSERT_VALUES);CHKERRXX(A01.ierr);
743  //SPIMat A(SPI::block(A00,A01,
744  //A10,A11));
745  SPIMat A(r0+r1,c0+c1);
746  A(0,c0,A01);
747  A();
748  return A;
749 
750  }
751  else if(k<0){
752  PetscInt r0=-k, r1=d.rows;
753  PetscInt c0=d.rows, c1=-k;
754  SPIMat A00(zeros(r0,c0)); CHKERRXX(A00.ierr);
755  SPIMat A01(zeros(r0,c1)); CHKERRXX(A01.ierr);
756  SPIMat A10(r1,c0); CHKERRXX(A10.ierr);
757  SPIMat A11(zeros(r1,c1)); CHKERRXX(A11.ierr);
758  //CHKERRXX(A00.ierr);
759  //CHKERRXX(A10.ierr);
760  //CHKERRXX(A01.ierr);
761  //CHKERRXX(A11.ierr);
762  A10.ierr = MatDiagonalSet(A10.mat,d.vec,INSERT_VALUES);CHKERRXX(A10.ierr);
763  //A00();CHKERRXX(A00.ierr);
764  //A01();CHKERRXX(A01.ierr);
765  //A10();CHKERRXX(A10.ierr);
766  //A11();CHKERRXX(A11.ierr);
767  SPIMat A(SPI::block({{A00,A01},
768  {A10,A11}}));
769  A();
770  return A();
771 
772  }
773  else{
774  exit(0);
775  }
776  }
777 
778  // kron inner product
780  SPIMat kron(
781  const SPIMat &A,
782  const SPIMat &B
783  ){
784  PetscErrorCode ierr;
785 
786  // get A,B information
787  PetscInt m,n,p,q;
788  MatGetSize(A.mat,&m,&n);
789  MatGetSize(B.mat,&p,&q);
790 
791  // assume square matrices A and B, so we can use set_Mat for the square submatrices
792  PetscInt na=m, nb=p,nc;
793  nc=m*p;
794 
795  // init C
796  SPIMat C(nc);
797 
798  // kron function C=kron(A,B)
799  PetscInt ncols;
800  const PetscInt *cols;
801  const PetscScalar *vals;
802  PetscInt Isubstart,Isubend;
803  ierr = MatGetOwnershipRange(A.mat,&Isubstart,&Isubend);CHKERRXX(ierr);
804  for (PetscInt rowi=0; rowi<na; rowi++){
805  //PetscPrintf(PETSC_COMM_WORLD,"kron rowi=%i of %i\n",rowi,m);
806  bool onprocessor=(Isubstart<=rowi) and (rowi<Isubend);
807  if(onprocessor){
808  // extract row of one A
809  ierr = MatGetRow(A.mat,rowi,&ncols,&cols,&vals);CHKERRXX(ierr); // extract the one row of A if owned by processor
810  }
811  else{
812  ncols=0;
813  }
814  PetscInt ncols2=0;
815  MPI_Allreduce(&ncols,&ncols2,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
816 
817  // set global vals2 array
818  PetscScalar vals2[ncols2];
819  PetscInt cols2[ncols2];
820  PetscScalar *vals_temp=new PetscScalar[ncols2] ();// new array and set to 0 with ()
821  PetscInt *cols_temp=new PetscInt[ncols2] ();
822  if(onprocessor){
823  for(PetscInt i=0; i<ncols2; i++){
824  vals_temp[i]=vals[i];
825  cols_temp[i]=cols[i];
826  }
827 
828  }
829  MPIU_Allreduce(vals_temp,vals2,ncols2,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
830  MPI_Allreduce(cols_temp,cols2,ncols2,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
831 
832  // every processor calls set_Mat for each col in row
833  for(PetscInt i=0; i<ncols2; i++){
834  // ierr = set_Mat(vals2[i],B,nb,C,nc,rowi*nb,cols2[i]*nb,INSERT_VALUES);CHKERRXX(ierr);
835  C(rowi*nb,cols2[i]*nb,B*vals2[i],INSERT_VALUES);
836  }
837 
838  if(onprocessor){
839  // restore row
840  ierr = MatRestoreRow(A.mat,rowi,&ncols,&cols,&vals); CHKERRXX(ierr);
841  }
842  delete[] vals_temp;
843  delete[] cols_temp;
844 
845  }
846  C();
847  return C;
848  }
849 
851  std::tuple<std::vector<PetscReal>, std::vector<SPIVec>, std::vector<SPIVec>> svd(
852  const SPIMat &A
853  ){
854  // set m to be the number of rows and k to be number of columns
855  PetscInt m=A.rows; // number of rows
856  PetscInt n=A.cols; // number of columns
857  SVD svd; // SVD context
858  //PetscReal sigmai;
859  //SPIVec ui(m,"u");
860  //SPIVec vi(n,"v");
861  std::vector<PetscReal> sigma(n);
862  std::vector<SPIVec> u(n);
863  std::vector<SPIVec> v(n);
864  PetscInt nconv;
865  //PetscReal error;
866 
867  SVDCreate(PETSC_COMM_WORLD, &svd);
868  SVDSetOperator(svd, A.mat);
869  //SVDSetProblemType(svd, SVD_STANDARD);
870  SVDSetFromOptions(svd);
871  SVDSetDimensions(svd,n,PETSC_DEFAULT,PETSC_DEFAULT); // force it to solve all of the SVD values
872  SVDSolve(svd);
873  //MatCreateVecs(A.mat,&vi.vec,&ui.vec);
874  //ui.flag_init=PETSC_TRUE;
875  //vi.flag_init=PETSC_TRUE;
876  SVDGetConverged(svd, &nconv);
877  std::cout<<"n="<<n<<" nconv = "<<nconv<<std::endl;
878  for(PetscInt j=0; j<nconv; j++){
879  //std::cout<<"j = "<<j<<std::endl;
880  MatCreateVecs(A.mat,&v[j].vec,&u[j].vec);
881  SVDGetSingularTriplet(svd,j,&(sigma[j]),u[j].vec,v[j].vec);
882  u[j].flag_init=PETSC_TRUE;
883  //u[j].name = "u";
884  u[j].rows = m;
885 
886  v[j].flag_init=PETSC_TRUE;
887  //v[j].name = "v";
888  v[j].rows = n;
889  //SVDGetSingularTriplet(svd,j,&sigmai,ui.vec,vi.vec);
890  //SVDComputeError(svd,j,SVD_ERROR_RELATIVE,&error);
891  //sigma.push_back(sigmai);
892  //u.push_back(ui);
893  //v.push_back(vi);
894  //sigma[j] = sigmai;
895  //u[j] = ui;
896  //v[j] = vi;
897  //std::cout<<"j = "<<j<<" sigma = "<<sigma[j]<<" error = "<<error<<std::endl;
898  }
899  //std::cout<<"made it here"<<std::endl;
900  SVDDestroy(&svd);
901  //std::cout<<"made it here deallocate"<<std::endl;
902  return std::make_tuple(sigma,u,v);
903  }
905  SPIVec lstsq(
906  const SPIMat &A,
907  SPIVec &y
908  ){
909  PetscInt n=A.cols;
910  std::vector<PetscReal> sigma;
911  std::vector<SPIVec> u,v;
912  SPIVec x(zeros(n));
913  std::tie(sigma,u,v) = svd(A);
914  for(PetscInt j=0; j<n; ++j){
915  x += ((y.dot(u[j]))/sigma[j]) * v[j];
916  }
917  return x;
918  }
920  SPIVec lstsq(
921  const std::vector<SPIVec> &A,
922  SPIVec &y
923  ){
924  SPIMat S(A);
925  return lstsq(S,y);
926  }
928  std::tuple<PetscScalar, SPIVec, SPIVec> eig(
929  const SPIMat &A,
930  const SPIMat &B,
931  const PetscScalar target,
932  const PetscReal tol,
933  const PetscInt max_iter
934  ){
935  //std::cout<<"target = "<<target<<std::endl;
936  PetscInt rows=A.rows;
937  EPS eps; /* eigenproblem solver context slepc */
938  //ST st;
939  //EPSType type;
940  //KSP ksp; /* linear solver context petsc */
941  PetscErrorCode ierr;
942  PetscScalar ki,alpha;
943  SPIVec eigl_vec(rows),eigr_vec(rows);
944 
945  //PetscScalar kr_temp, ki_temp;
946  PetscScalar kr_temp;
947 
948  // Create the eigenvalue solver and set various options
949  // Create solver contexts
950  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRXX(ierr);
951  // Set operators. Here the matrix that defines the eigenvalue system
952  // swap operators such that LHS matrix is singular, and RHS matrix can be inverted (for slepc)
953  // this will make the Ax=lambda Bx into the problem Bx = (1./lambda) Ax, thus our eigenvalues are inverted
954  ierr = EPSSetOperators(eps,A.mat,B.mat);CHKERRXX(ierr);
955  //ierr = EPSSetOperators(eps,B,A);CHKERRXX(ierr);
956  // Set runtime options, e.g.,
957  // -
958  ierr = EPSSetFromOptions(eps);CHKERRXX(ierr);
959  //std::cout<<"After KSPSetFromOptions"<<std::endl;
960  //
961  // set convergence type
962  EPSWhich which=EPS_TARGET_MAGNITUDE;
963  PetscInt nev=1;
964  EPSSetWhichEigenpairs(eps,which);
965  EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
966  if ((max_iter==-1) && (tol==-1.)){
967  EPSSetTolerances(eps,PETSC_DEFAULT,PETSC_DEFAULT);
968  }
969  else if(tol==-1.){
970  EPSSetTolerances(eps,PETSC_DEFAULT,max_iter);
971  }
972  else if(max_iter==-1){
973  EPSSetTolerances(eps,tol,PETSC_DEFAULT);
974  }
975  else{
976  EPSSetTolerances(eps,tol,max_iter);
977  }
978  if (
979  which==EPS_TARGET_REAL ||
980  which==EPS_TARGET_IMAGINARY ||
981  which==EPS_TARGET_MAGNITUDE){
982  // PetscScalar target=0.-88.5*PETSC_i;
983  EPSSetTarget(eps,target);
984  //EPSSetTolerances(eps,1.E-8,100000);
985  }
986 
987 
988  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
989  ST st;
990  EPSGetST(eps,&st);
991  STSetType(st,STSINVERT);
992  // Solve the system with left and right
993  EPSSetTwoSided(eps,PETSC_TRUE);
994  // Solve the system
995  ierr = EPSSolve(eps);CHKERRXX(ierr);
996  //std::cout<<"After KSPSolve"<<std::endl;
997 
998  // output iterations
999  //PetscInt its, maxit, i;
1000  PetscInt nconv;
1001  //PetscReal error, tol, re, im;
1002  //PetscReal tol2;
1003  /* Optional: Get some information from the solver and display it */
1004  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1005  // ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRXX(ierr);
1006  // ierr = EPSGetType(eps,&type);CHKERRXX(ierr);
1007  // ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRXX(ierr);
1008  // ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRXX(ierr);
1009  // ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRXX(ierr);
1010  // ierr = EPSGetTolerances(eps,&tol2,&maxit);CHKERRXX(ierr);
1011  // ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol2,maxit);CHKERRXX(ierr);
1012 
1013  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1014  Display solution and clean up
1015  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1016  /*
1017  Get number of converged approximate eigenpairs
1018  */
1019  ierr = EPSGetConverged(eps,&nconv);CHKERRXX(ierr);
1020  // ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);CHKERRXX(ierr);
1021 
1022  if (nconv>0) {
1023  /*
1024  Display eigenvalues and relative errors
1025  */
1026  // ierr = PetscPrintf(PETSC_COMM_WORLD,
1027  // " k ||Ax-kx||/||kx||\n"
1028  // " ----------------- ------------------\n");CHKERRXX(ierr);
1029 
1030  // for (PetscInt i=0;i<nconv;i++)
1031  // /*
1032  // Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
1033  // ki (imaginary part)
1034  // */
1035  // ierr = EPSGetEigenpair(eps,i,&alpha,&ki,eig_vec.vec,xi.vec);CHKERRXX(ierr);
1036  // /*
1037  // Compute the relative error associated to each eigenpair
1038  // */
1039  // PetscReal error, re, im;
1040  // ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRXX(ierr);
1041  // re = PetscRealPart(alpha);
1042  // im = PetscImaginaryPart(alpha);
1043  // if (im!=0.0)
1044  // ierr = PetscPrintf(PETSC_COMM_WORLD," (%9e+%9ei) %12g\n",(double)re,(double)im,(double)error);CHKERRXX(ierr);
1045  // else
1046  // ierr = PetscPrintf(PETSC_COMM_WORLD," %12e %12g\n",(double)re,(double)error);CHKERRXX(ierr);
1047  //
1048  //
1049  // ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRXX(ierr);
1050 
1051  ierr = EPSGetEigenpair(eps,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1052  ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1053  }
1054 
1055  // PetscInt its;
1056  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1057  // //ierr = PetscPrintf(PETSC_COMM_WORLD,"ksp iterations %D\n",its);CHKERRXX(ierr);
1058  // ierr = PetscPrintf(PETSC_COMM_WORLD,"EPS Solved in %D iterations \n",its); CHKERRXX(ierr);
1059  // Free work space. All PETSc objects should be destroyed when they
1060  // are no longer needed.
1061  //set_Vec(x);
1062  ierr = EPSDestroy(&eps);CHKERRXX(ierr);
1063  //ierr = VecDestroy(&x);CHKERRXX(ierr);
1064  //ierr = VecDestroy(&b);CHKERRXX(ierr);
1065  //ierr = MatDestroy(&A);CHKERRXX(ierr);
1066  //ierr = PetscFinalize();
1067 
1068  return std::make_tuple(alpha,eigl_vec,eigr_vec);
1069  //return std::make_tuple(alpha,alpha);
1070  }
1072  std::tuple<PetscScalar, SPIVec> eig_right(
1073  const SPIMat &A,
1074  const SPIMat &B,
1075  const PetscScalar target,
1076  const PetscReal tol,
1077  const PetscInt max_iter
1078  ){
1079  //std::cout<<"target = "<<target<<std::endl;
1080  PetscInt rows=A.rows;
1081  EPS eps; /* eigenproblem solver context slepc */
1082  //ST st;
1083  //EPSType type;
1084  //KSP ksp; /* linear solver context petsc */
1085  PetscErrorCode ierr;
1086  PetscScalar ki,alpha;
1087  SPIVec eigr_vec(rows);
1088 
1089  //PetscScalar kr_temp, ki_temp;
1090  PetscScalar kr_temp;
1091 
1092  // Create the eigenvalue solver and set various options
1093  // Create solver contexts
1094  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRXX(ierr);
1095  // Set operators. Here the matrix that defines the eigenvalue system
1096  // swap operators such that LHS matrix is singular, and RHS matrix can be inverted (for slepc)
1097  // this will make the Ax=lambda Bx into the problem Bx = (1./lambda) Ax, thus our eigenvalues are inverted
1098  ierr = EPSSetOperators(eps,A.mat,B.mat);CHKERRXX(ierr);
1099  //ierr = EPSSetOperators(eps,B,A);CHKERRXX(ierr);
1100  // Set runtime options, e.g.,
1101  // -
1102  ierr = EPSSetFromOptions(eps);CHKERRXX(ierr);
1103  //std::cout<<"After KSPSetFromOptions"<<std::endl;
1104  //
1105  // set convergence type
1106  EPSWhich which=EPS_TARGET_MAGNITUDE;
1107  PetscInt nev=1;
1108  EPSSetWhichEigenpairs(eps,which);
1109  EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
1110  if ((max_iter==-1) && (tol==-1.)){
1111  EPSSetTolerances(eps,PETSC_DEFAULT,PETSC_DEFAULT);
1112  }
1113  else if(tol==-1.){
1114  EPSSetTolerances(eps,PETSC_DEFAULT,max_iter);
1115  }
1116  else if(max_iter==-1){
1117  EPSSetTolerances(eps,tol,PETSC_DEFAULT);
1118  }
1119  else{
1120  EPSSetTolerances(eps,tol,max_iter);
1121  }
1122  if (
1123  which==EPS_TARGET_REAL ||
1124  which==EPS_TARGET_IMAGINARY ||
1125  which==EPS_TARGET_MAGNITUDE){
1126  // PetscScalar target=0.-88.5*PETSC_i;
1127  EPSSetTarget(eps,target);
1128  //EPSSetTolerances(eps,1.E-8,100000);
1129  }
1130 
1131 
1132  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
1133  ST st;
1134  EPSGetST(eps,&st);
1135  STSetType(st,STSINVERT);
1136  // Solve the system with left and right
1137  //EPSSetTwoSided(eps,PETSC_TRUE);
1138  // Solve the system
1139  ierr = EPSSolve(eps);CHKERRXX(ierr);
1140  //std::cout<<"After KSPSolve"<<std::endl;
1141 
1142  // output iterations
1143  //PetscInt its, maxit, i;
1144  PetscInt nconv;
1145  //PetscReal error, tol, re, im;
1146  //PetscReal tol2;
1147  /* Optional: Get some information from the solver and display it */
1148  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1149  // ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRXX(ierr);
1150  // ierr = EPSGetType(eps,&type);CHKERRXX(ierr);
1151  // ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRXX(ierr);
1152  // ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRXX(ierr);
1153  // ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRXX(ierr);
1154  // ierr = EPSGetTolerances(eps,&tol2,&maxit);CHKERRXX(ierr);
1155  // ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol2,maxit);CHKERRXX(ierr);
1156 
1157  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1158  Display solution and clean up
1159  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1160  /*
1161  Get number of converged approximate eigenpairs
1162  */
1163  ierr = EPSGetConverged(eps,&nconv);CHKERRXX(ierr);
1164  // ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);CHKERRXX(ierr);
1165 
1166  if (nconv>0) {
1167  /*
1168  Display eigenvalues and relative errors
1169  */
1170  // ierr = PetscPrintf(PETSC_COMM_WORLD,
1171  // " k ||Ax-kx||/||kx||\n"
1172  // " ----------------- ------------------\n");CHKERRXX(ierr);
1173 
1174  // for (PetscInt i=0;i<nconv;i++)
1175  // /*
1176  // Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
1177  // ki (imaginary part)
1178  // */
1179  // ierr = EPSGetEigenpair(eps,i,&alpha,&ki,eig_vec.vec,xi.vec);CHKERRXX(ierr);
1180  // /*
1181  // Compute the relative error associated to each eigenpair
1182  // */
1183  // PetscReal error, re, im;
1184  // ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRXX(ierr);
1185  // re = PetscRealPart(alpha);
1186  // im = PetscImaginaryPart(alpha);
1187  // if (im!=0.0)
1188  // ierr = PetscPrintf(PETSC_COMM_WORLD," (%9e+%9ei) %12g\n",(double)re,(double)im,(double)error);CHKERRXX(ierr);
1189  // else
1190  // ierr = PetscPrintf(PETSC_COMM_WORLD," %12e %12g\n",(double)re,(double)error);CHKERRXX(ierr);
1191  //
1192  //
1193  // ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRXX(ierr);
1194 
1195  ierr = EPSGetEigenpair(eps,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1196  eigr_vec.rows=rows;
1197  //ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1198  }
1199 
1200  // PetscInt its;
1201  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1202  // //ierr = PetscPrintf(PETSC_COMM_WORLD,"ksp iterations %D\n",its);CHKERRXX(ierr);
1203  // ierr = PetscPrintf(PETSC_COMM_WORLD,"EPS Solved in %D iterations \n",its); CHKERRXX(ierr);
1204  // Free work space. All PETSc objects should be destroyed when they
1205  // are no longer needed.
1206  //set_Vec(x);
1207  ierr = EPSDestroy(&eps);CHKERRXX(ierr);
1208  //ierr = VecDestroy(&x);CHKERRXX(ierr);
1209  //ierr = VecDestroy(&b);CHKERRXX(ierr);
1210  //ierr = MatDestroy(&A);CHKERRXX(ierr);
1211  //ierr = PetscFinalize();
1212 
1213  return std::make_tuple(alpha,eigr_vec);
1214  //return std::make_tuple(alpha,alpha);
1215  }
1217  std::tuple<PetscScalar, SPIVec, SPIVec> eig_init(
1218  const SPIMat &A,
1219  const SPIMat &B,
1220  const PetscScalar target,
1221  const SPIVec &ql,
1222  const SPIVec &qr,
1223  PetscReal tol,
1224  const PetscInt max_iter
1225  ){
1226  //std::cout<<"target = "<<target<<std::endl;
1227  PetscInt rows=A.rows;
1228  EPS eps; /* eigenproblem solver context slepc */
1229  //ST st;
1230  //EPSType type;
1231  //KSP ksp; /* linear solver context petsc */
1232  PetscErrorCode ierr;
1233  PetscScalar alpha;
1234  SPIVec eigl_vec(rows),eigr_vec(rows);
1235 
1236  PetscScalar kr_temp, ki_temp;
1237 
1238  // Create the eigenvalue solver and set various options
1239  // Create solver contexts
1240  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRXX(ierr);
1241  // Set operators. Here the matrix that defines the eigenvalue system
1242  // swap operators such that LHS matrix is singular, and RHS matrix can be inverted (for slepc)
1243  // this will make the Ax=lambda Bx into the problem Bx = (1./lambda) Ax, thus our eigenvalues are inverted
1244  ierr = EPSSetOperators(eps,A.mat,B.mat);CHKERRXX(ierr);
1245  //ierr = EPSSetOperators(eps,B,A);CHKERRXX(ierr);
1246  // Set runtime options, e.g.,
1247  // -
1248  ierr = EPSSetFromOptions(eps);CHKERRXX(ierr);
1249  //std::cout<<"After KSPSetFromOptions"<<std::endl;
1250  //
1251  // set convergence type
1252  EPSWhich which=EPS_TARGET_MAGNITUDE;
1253  PetscInt nev=1;
1254  EPSSetWhichEigenpairs(eps,which);
1255  EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
1256  if ((max_iter==-1) && (tol==-1.)){
1257  EPSSetTolerances(eps,PETSC_DEFAULT,PETSC_DEFAULT);
1258  }
1259  else if(tol==-1.){
1260  EPSSetTolerances(eps,PETSC_DEFAULT,max_iter);
1261  }
1262  else if(max_iter==-1){
1263  EPSSetTolerances(eps,tol,PETSC_DEFAULT);
1264  }
1265  else{
1266  EPSSetTolerances(eps,tol,max_iter);
1267  }
1268  if (
1269  which==EPS_TARGET_REAL ||
1270  which==EPS_TARGET_IMAGINARY ||
1271  which==EPS_TARGET_MAGNITUDE){
1272  // PetscScalar target=0.-88.5*PETSC_i;
1273  EPSSetTarget(eps,target);
1274  //EPSSetTolerances(eps,1.E-8,100000);
1275  }
1276 
1277 
1278  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
1279  ST st;
1280  EPSGetST(eps,&st);
1281  STSetType(st,STSINVERT);
1282  // set initial guess
1283  std::vector<Vec> qrvec(1);
1284  qrvec[0] = qr.vec;
1285  ierr = EPSSetInitialSpace(eps,1,qrvec.data());CHKERRXX(ierr);
1286  std::vector<Vec> qlvec(1);
1287  qlvec[0] = ql.vec;
1288  ierr = EPSSetLeftInitialSpace(eps,1,qlvec.data());CHKERRXX(ierr); // TODO doesn't work using 3.12.1 version, but should work in 3.14 documentation
1289 
1290  // Solve the system with left and right
1291  EPSSetTwoSided(eps,PETSC_TRUE);
1292  ierr = EPSSolve(eps);CHKERRXX(ierr);
1293  //std::cout<<"After KSPSolve"<<std::endl;
1294 
1295  // output iterations
1296  //PetscInt its, maxit;
1297  PetscInt nconv;
1298  //PetscReal error, tol, re, im;
1299  /* Optional: Get some information from the solver and display it */
1300  //ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1301  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRXX(ierr);
1302  //ierr = EPSGetType(eps,&type);CHKERRXX(ierr);
1303  //ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRXX(ierr);
1304  //ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRXX(ierr);
1305  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRXX(ierr);
1306  //ierr = EPSGetTolerances(eps,&tol,&maxit);CHKERRXX(ierr);
1307  //ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);CHKERRXX(ierr);
1308 
1309  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1310  Display solution and clean up
1311  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1312  /*
1313  Get number of converged approximate eigenpairs
1314  */
1315  ierr = EPSGetConverged(eps,&nconv);CHKERRXX(ierr);
1316  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);CHKERRXX(ierr);
1317 
1318  if (nconv>0) {
1319  /*
1320  Display eigenvalues and relative errors
1321  */
1322  // ierr = PetscPrintf(PETSC_COMM_WORLD,
1323  // " k ||Ax-kx||/||kx||\n"
1324  // " ----------------- ------------------\n");CHKERRXX(ierr);
1325 
1326  // for (PetscInt i=0;i<nconv;i++)
1327  // /*
1328  // Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
1329  // ki (imaginary part)
1330  // */
1331  // ierr = EPSGetEigenpair(eps,i,&alpha,&ki,eigr_vec.vec,eigl_vec.vec);CHKERRXX(ierr);
1332  // /*
1333  // Compute the relative error associated to each eigenpair
1334  // */
1335  // PetscReal error, re, im;
1336  // ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRXX(ierr);
1337  // re = PetscRealPart(alpha);
1338  // im = PetscImaginaryPart(alpha);
1339  // if (im!=0.0)
1340  // ierr = PetscPrintf(PETSC_COMM_WORLD," (%9e+%9ei) %12g\n",(double)re,(double)im,(double)error);CHKERRXX(ierr);
1341  // else
1342  // ierr = PetscPrintf(PETSC_COMM_WORLD," %12e %12g\n",(double)re,(double)error);CHKERRXX(ierr);
1343  //
1344  //
1345  // ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRXX(ierr);
1346 
1347  ierr = EPSGetEigenpair(eps,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1348  ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1349  }
1350 
1351  // PetscInt its;
1352  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1353  // //ierr = PetscPrintf(PETSC_COMM_WORLD,"ksp iterations %D\n",its);CHKERRXX(ierr);
1354  // ierr = PetscPrintf(PETSC_COMM_WORLD,"EPS Solved in %D iterations \n",its); CHKERRXX(ierr);
1355  // Free work space. All PETSc objects should be destroyed when they
1356  // are no longer needed.
1357  //set_Vec(x);
1358  ierr = EPSDestroy(&eps);CHKERRXX(ierr);
1359  //ierr = VecDestroy(&x);CHKERRXX(ierr);
1360  //ierr = VecDestroy(&b);CHKERRXX(ierr);
1361  //ierr = MatDestroy(&A);CHKERRXX(ierr);
1362  //ierr = PetscFinalize();
1363 
1364  return std::make_tuple(alpha,eigl_vec,eigr_vec);
1365  //return std::make_tuple(alpha,alpha);
1366  }
1368  std::tuple<PetscScalar, SPIVec> eig_init_right(
1369  const SPIMat &A,
1370  const SPIMat &B,
1371  const PetscScalar target,
1372  const SPIVec &qr,
1373  PetscReal tol,
1374  const PetscInt max_iter
1375  ){
1376  //std::cout<<"target = "<<target<<std::endl;
1377  PetscInt rows=A.rows;
1378  EPS eps; /* eigenproblem solver context slepc */
1379  //ST st;
1380  //EPSType type;
1381  //KSP ksp; /* linear solver context petsc */
1382  PetscErrorCode ierr;
1383  PetscScalar alpha;
1384  SPIVec eigr_vec(rows);
1385 
1386  PetscScalar kr_temp, ki_temp;
1387 
1388  // Create the eigenvalue solver and set various options
1389  // Create solver contexts
1390  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRXX(ierr);
1391  // Set operators. Here the matrix that defines the eigenvalue system
1392  // swap operators such that LHS matrix is singular, and RHS matrix can be inverted (for slepc)
1393  // this will make the Ax=lambda Bx into the problem Bx = (1./lambda) Ax, thus our eigenvalues are inverted
1394  ierr = EPSSetOperators(eps,A.mat,B.mat);CHKERRXX(ierr);
1395  //ierr = EPSSetOperators(eps,B,A);CHKERRXX(ierr);
1396  // Set runtime options, e.g.,
1397  // -
1398  ierr = EPSSetFromOptions(eps);CHKERRXX(ierr);
1399  //std::cout<<"After KSPSetFromOptions"<<std::endl;
1400  //
1401  // set convergence type
1402  EPSWhich which=EPS_TARGET_MAGNITUDE;
1403  PetscInt nev=1;
1404  EPSSetWhichEigenpairs(eps,which);
1405  EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
1406  if ((max_iter==-1) && (tol==-1.)){
1407  EPSSetTolerances(eps,PETSC_DEFAULT,PETSC_DEFAULT);
1408  }
1409  else if(tol==-1.){
1410  EPSSetTolerances(eps,PETSC_DEFAULT,max_iter);
1411  }
1412  else if(max_iter==-1){
1413  EPSSetTolerances(eps,tol,PETSC_DEFAULT);
1414  }
1415  else{
1416  EPSSetTolerances(eps,tol,max_iter);
1417  }
1418  if (
1419  which==EPS_TARGET_REAL ||
1420  which==EPS_TARGET_IMAGINARY ||
1421  which==EPS_TARGET_MAGNITUDE){
1422  // PetscScalar target=0.-88.5*PETSC_i;
1423  EPSSetTarget(eps,target);
1424  //EPSSetTolerances(eps,1.E-8,100000);
1425  }
1426 
1427 
1428  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
1429  ST st;
1430  EPSGetST(eps,&st);
1431  STSetType(st,STSINVERT);
1432  // set initial guess
1433  std::vector<Vec> qrvec(1);
1434  qrvec[0] = qr.vec;
1435  ierr = EPSSetInitialSpace(eps,1,qrvec.data());CHKERRXX(ierr);
1436  //std::vector<Vec> qlvec(1);
1437  //qlvec[0] = ql.vec;
1438  //ierr = EPSSetLeftInitialSpace(eps,1,qlvec.data());CHKERRXX(ierr); // TODO doesn't work using 3.12.1 version, but should work in 3.14 documentation
1439 
1440  // Solve the system with left and right
1441  //EPSSetTwoSided(eps,PETSC_TRUE);
1442  ierr = EPSSolve(eps);CHKERRXX(ierr);
1443  //std::cout<<"After KSPSolve"<<std::endl;
1444 
1445  // output iterations
1446  //PetscInt its, maxit;
1447  PetscInt nconv;
1448  //PetscReal error, tol, re, im;
1449  /*
1450 Optional: Get some information from the solver and display it
1451 */
1452  //ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1453  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRXX(ierr);
1454  //ierr = EPSGetType(eps,&type);CHKERRXX(ierr);
1455  //ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRXX(ierr);
1456  //ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRXX(ierr);
1457  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRXX(ierr);
1458  //ierr = EPSGetTolerances(eps,&tol,&maxit);CHKERRXX(ierr);
1459  //ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);CHKERRXX(ierr);
1460 
1461  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1462  Display solution and clean up
1463  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1464  /*
1465  Get number of converged approximate eigenpairs
1466  */
1467  ierr = EPSGetConverged(eps,&nconv);CHKERRXX(ierr);
1468  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);CHKERRXX(ierr);
1469 
1470  if (nconv>0) {
1471  /*
1472  Display eigenvalues and relative errors
1473  */
1474  // ierr = PetscPrintf(PETSC_COMM_WORLD,
1475  // " k ||Ax-kx||/||kx||\n"
1476  // " ----------------- ------------------\n");CHKERRXX(ierr);
1477 
1478  // for (PetscInt i=0;i<nconv;i++)
1479  // /*
1480  // Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
1481  // ki (imaginary part)
1482  // */
1483  // ierr = EPSGetEigenpair(eps,i,&alpha,&ki,eigr_vec.vec,eigl_vec.vec);CHKERRXX(ierr);
1484  // /*
1485  // Compute the relative error associated to each eigenpair
1486  // */
1487  // PetscReal error, re, im;
1488  // ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRXX(ierr);
1489  // re = PetscRealPart(alpha);
1490  // im = PetscImaginaryPart(alpha);
1491  // if (im!=0.0)
1492  // ierr = PetscPrintf(PETSC_COMM_WORLD," (%9e+%9ei) %12g\n",(double)re,(double)im,(double)error);CHKERRXX(ierr);
1493  // else
1494  // ierr = PetscPrintf(PETSC_COMM_WORLD," %12e %12g\n",(double)re,(double)error);CHKERRXX(ierr);
1495  //
1496  //
1497  // ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRXX(ierr);
1498 
1499  ierr = EPSGetEigenpair(eps,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1500  eigr_vec.rows=rows;
1501  //ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1502  }
1503 
1504  // PetscInt its;
1505  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1506  // //ierr = PetscPrintf(PETSC_COMM_WORLD,"ksp iterations %D\n",its);CHKERRXX(ierr);
1507  // ierr = PetscPrintf(PETSC_COMM_WORLD,"EPS Solved in %D iterations \n",its); CHKERRXX(ierr);
1508  // Free work space. All PETSc objects should be destroyed when they
1509  // are no longer needed.
1510  //set_Vec(x);
1511  ierr = EPSDestroy(&eps);CHKERRXX(ierr);
1512  //ierr = VecDestroy(&x);CHKERRXX(ierr);
1513  //ierr = VecDestroy(&b);CHKERRXX(ierr);
1514  //ierr = MatDestroy(&A);CHKERRXX(ierr);
1515  //ierr = PetscFinalize();
1516 
1517  return std::make_tuple(alpha,eigr_vec);
1518  //return std::make_tuple(alpha,alpha);
1519  }
1521  std::tuple<std::vector<PetscScalar>, std::vector<SPIVec>> eig_init_rights(
1522  const SPIMat &A,
1523  const SPIMat &B,
1524  const std::vector<PetscScalar> targets,
1525  const std::vector<SPIVec> &qrs,
1526  PetscReal tol,
1527  const PetscInt max_iter
1528  ){
1529  //std::cout<<"target = "<<target<<std::endl;
1530  PetscInt rows=A.rows;
1531  EPS eps; /* eigenproblem solver context slepc */
1532  //ST st;
1533  //EPSType type;
1534  //KSP ksp; /* linear solver context petsc */
1535  PetscErrorCode ierr;
1536  std::vector<PetscScalar> alphas(qrs.size());
1537  PetscScalar alpha;
1538  std::vector<SPIVec> eigr_vecs(qrs.size());
1539  SPIVec eigr_vec(rows);
1540 
1541  PetscScalar kr_temp, ki_temp;
1542 
1543  // Create the eigenvalue solver and set various options
1544  // Create solver contexts
1545  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRXX(ierr);
1546  // Set operators. Here the matrix that defines the eigenvalue system
1547  // swap operators such that LHS matrix is singular, and RHS matrix can be inverted (for slepc)
1548  // this will make the Ax=lambda Bx into the problem Bx = (1./lambda) Ax, thus our eigenvalues are inverted
1549  ierr = EPSSetOperators(eps,A.mat,B.mat);CHKERRXX(ierr);
1550  //ierr = EPSSetOperators(eps,B,A);CHKERRXX(ierr);
1551  // Set runtime options, e.g.,
1552  // -
1553  ierr = EPSSetFromOptions(eps);CHKERRXX(ierr);
1554  //std::cout<<"After KSPSetFromOptions"<<std::endl;
1555  //
1556  // set convergence type
1557  EPSWhich which=EPS_TARGET_MAGNITUDE;
1558  PetscInt nev=1;
1559  EPSSetWhichEigenpairs(eps,which);
1560  EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
1561  if ((max_iter==-1) && (tol==-1.)){
1562  EPSSetTolerances(eps,PETSC_DEFAULT,PETSC_DEFAULT);
1563  }
1564  else if(tol==-1.){
1565  EPSSetTolerances(eps,PETSC_DEFAULT,max_iter);
1566  }
1567  else if(max_iter==-1){
1568  EPSSetTolerances(eps,tol,PETSC_DEFAULT);
1569  }
1570  else{
1571  EPSSetTolerances(eps,tol,max_iter);
1572  }
1573  if (
1574  which==EPS_TARGET_REAL ||
1575  which==EPS_TARGET_IMAGINARY ||
1576  which==EPS_TARGET_MAGNITUDE){
1577  // PetscScalar target=0.-88.5*PETSC_i;
1578  EPSSetTarget(eps,targets[0]);
1579  //EPSSetTolerances(eps,1.E-8,100000);
1580  }
1581 
1582 
1583  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
1584  ST st;
1585  EPSGetST(eps,&st);
1586  STSetType(st,STSINVERT);
1587  // set initial guess
1588  PetscInt k=qrs.size();
1589  std::vector<Vec> qrvec(k);
1590  for(PetscInt i=0; i<k; ++i){
1591  qrvec[i] = qrs[i].vec;
1592  }
1593  ierr = EPSSetInitialSpace(eps,qrvec.size(),qrvec.data());CHKERRXX(ierr);
1594  //std::vector<Vec> qlvec(1);
1595  //qlvec[0] = ql.vec;
1596  //ierr = EPSSetLeftInitialSpace(eps,1,qlvec.data());CHKERRXX(ierr); // TODO doesn't work using 3.12.1 version, but should work in 3.14 documentation
1597 
1598  // Solve the system with left and right
1599  //EPSSetTwoSided(eps,PETSC_TRUE);
1600  PetscInt nconv;
1601  for(PetscInt i=0; i<k; ++i){
1602  EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL); // have to add this to switch the eigenvalue solver to a new target
1603  EPSSetWhichEigenpairs(eps,which);
1604  EPSSetTarget(eps,targets[i]);
1605  ierr = EPSSolve(eps);CHKERRXX(ierr);
1606  //std::cout<<"After KSPSolve"<<std::endl;
1607 
1608  // output iterations
1609  //PetscInt its, maxit;
1610  //PetscReal error, tol, re, im;
1611  /*
1612 Optional: Get some information from the solver and display it
1613 */
1614  //ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1615  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRXX(ierr);
1616  //ierr = EPSGetType(eps,&type);CHKERRXX(ierr);
1617  //ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRXX(ierr);
1618  //ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRXX(ierr);
1619  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRXX(ierr);
1620  //ierr = EPSGetTolerances(eps,&tol,&maxit);CHKERRXX(ierr);
1621  //ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);CHKERRXX(ierr);
1622 
1623  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1624  Display solution and clean up
1625  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1626  /*
1627  Get number of converged approximate eigenpairs
1628  */
1629  ierr = EPSGetConverged(eps,&nconv);CHKERRXX(ierr);
1630  //ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);CHKERRXX(ierr);
1631 
1632  if (nconv>0) {
1633  /*
1634  Display eigenvalues and relative errors
1635  */
1636  // ierr = PetscPrintf(PETSC_COMM_WORLD,
1637  // " k ||Ax-kx||/||kx||\n"
1638  // " ----------------- ------------------\n");CHKERRXX(ierr);
1639 
1640  // for (PetscInt i=0;i<nconv;i++)
1641  // /*
1642  // Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
1643  // ki (imaginary part)
1644  // */
1645  // ierr = EPSGetEigenpair(eps,i,&alpha,&ki,eigr_vec.vec,eigl_vec.vec);CHKERRXX(ierr);
1646  // /*
1647  // Compute the relative error associated to each eigenpair
1648  // */
1649  // PetscReal error, re, im;
1650  // ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRXX(ierr);
1651  // re = PetscRealPart(alpha);
1652  // im = PetscImaginaryPart(alpha);
1653  // if (im!=0.0)
1654  // ierr = PetscPrintf(PETSC_COMM_WORLD," (%9e+%9ei) %12g\n",(double)re,(double)im,(double)error);CHKERRXX(ierr);
1655  // else
1656  // ierr = PetscPrintf(PETSC_COMM_WORLD," %12e %12g\n",(double)re,(double)error);CHKERRXX(ierr);
1657  //
1658  //
1659  // ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRXX(ierr);
1660 
1661  ierr = EPSGetEigenpair(eps,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1662  alphas[i] = alpha;
1663  //EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
1664  eigr_vecs[i] = eigr_vec;
1665  //ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1666  }
1667  }
1668 
1669  // PetscInt its;
1670  // ierr = EPSGetIterationNumber(eps,&its);CHKERRXX(ierr);
1671  // //ierr = PetscPrintf(PETSC_COMM_WORLD,"ksp iterations %D\n",its);CHKERRXX(ierr);
1672  // ierr = PetscPrintf(PETSC_COMM_WORLD,"EPS Solved in %D iterations \n",its); CHKERRXX(ierr);
1673  // Free work space. All PETSc objects should be destroyed when they
1674  // are no longer needed.
1675  //set_Vec(x);
1676  ierr = EPSDestroy(&eps);CHKERRXX(ierr);
1677  //ierr = VecDestroy(&x);CHKERRXX(ierr);
1678  //ierr = VecDestroy(&b);CHKERRXX(ierr);
1679  //ierr = MatDestroy(&A);CHKERRXX(ierr);
1680  //ierr = PetscFinalize();
1681 
1682  return std::make_tuple(alphas,eigr_vecs);
1683  //return std::make_tuple(alpha,alpha);
1684  }
1686  std::tuple<PetscScalar, SPIVec> polyeig(
1687  const std::vector<SPIMat> &As,
1688  const PetscScalar target,
1689  const PetscReal tol,
1690  const PetscInt max_iter
1691  ){
1692  // if linear eigenvalue problem, use EPS
1693  if(0){//As.size()==1)
1694  //return eig(As[0],SPI::eye(As[0].rows),target,tol,max_iter);
1695  }
1696  else if(0){//As.size()==2)
1697  //return eig(As[0],-As[1],target,tol,max_iter);
1698  }
1699  else{// else polynomial eigenvalue problem use PEP
1700  PetscInt rows=As[0].rows;
1701  PEP pep; /* polynomial eigenproblem solver context slepc */
1702  PetscErrorCode ierr;
1703  PetscScalar ki,alpha;
1704  //SPIVec eigl_vec(rows),eigr_vec(rows);
1705  SPIVec eigr_vec(rows);
1706  PetscScalar kr_temp;
1707  // Create the eigenvalue solver and set various options
1708  ierr = PEPCreate(PETSC_COMM_WORLD,&pep);CHKERRXX(ierr);
1709  // Set operators. Here the matrix that defines the eigenvalue system
1710  std::vector<Mat> AsMat(As.size());
1711  for (unsigned i=0; i<As.size(); ++i){ AsMat[i]=As[i].mat; }
1712  ierr = PEPSetOperators(pep,AsMat.size(),AsMat.data());CHKERRXX(ierr);
1713  // Set runtime options, e.g.,
1714  ierr = PEPSetFromOptions(pep);CHKERRXX(ierr);
1715  // set convergence type
1716  PEPWhich which=PEP_TARGET_MAGNITUDE;
1717  PetscInt nev=1;
1718  PEPSetWhichEigenpairs(pep,which);
1719  PEPSetDimensions(pep,nev,PETSC_DEFAULT,PETSC_DEFAULT);
1720  if ((max_iter==-1) && (tol==-1.)){ PEPSetTolerances(pep,PETSC_DEFAULT, PETSC_DEFAULT); }
1721  else if(tol==-1.){ PEPSetTolerances(pep,PETSC_DEFAULT, max_iter); }
1722  else if(max_iter==-1){ PEPSetTolerances(pep,tol, PETSC_DEFAULT); }
1723  else{ PEPSetTolerances(pep,tol, max_iter); }
1724  PEPSetTarget(pep,target);
1725  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
1726  ST st;
1727  PEPGetST(pep,&st);
1728  STSetType(st,STSINVERT);
1729  // get EPS and set two sided
1730  //EPS eps;
1731  //ierr = PEPLinearGetEPS(pep,&eps);// CHKERRXX(ierr);
1732  //EPSSetTwoSided(eps,PETSC_TRUE); // this doesn't work TODO
1733  // Solve the system
1734  ierr = PEPSolve(pep);CHKERRXX(ierr);
1735  //ierr = EPSSolve(eps);CHKERRXX(ierr);
1736  // output iterations
1737  PetscInt nconv;
1738  ierr = PEPGetConverged(pep,&nconv);CHKERRXX(ierr);
1739  //ierr = EPSGetConverged(eps,&nconv);CHKERRXX(ierr);
1740  if (nconv>0) {
1741  ierr = PEPGetEigenpair(pep,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1742  eigr_vec.flag_init=PETSC_TRUE;
1743  //ierr = EPSGetEigenpair(eps,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1744  //ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1745 
1746  }
1747  // destroy pep
1748  ierr = PEPDestroy(&pep);CHKERRXX(ierr);
1749  // return
1750  return std::make_tuple(alpha,eigr_vec);
1751  }
1752  }
1754  std::tuple<PetscScalar, SPIVec> polyeig_init(
1755  const std::vector<SPIMat> &As,
1756  const PetscScalar target,
1757  const SPIVec &qr,
1758  const PetscReal tol,
1759  const PetscInt max_iter
1760  ){
1761  // if linear eigenvalue problem, use EPS
1762  if(0){//As.size()==1)
1763  //return eig(As[0],SPI::eye(As[0].rows),target,tol,max_iter);
1764  }
1765  else if(0){//As.size()==2)
1766  //return eig(As[0],-As[1],target,tol,max_iter);
1767  }
1768  else{// else polynomial eigenvalue problem use PEP
1769  PetscInt rows=As[0].rows;
1770  PEP pep; /* polynomial eigenproblem solver context slepc */
1771  PetscErrorCode ierr;
1772  PetscScalar alpha;
1773  //SPIVec eigl_vec(rows),eigr_vec(rows);
1774  SPIVec eigr_vec(rows);
1775  PetscScalar kr_temp;
1776  // Create the eigenvalue solver and set various options
1777  ierr = PEPCreate(PETSC_COMM_WORLD,&pep);CHKERRXX(ierr);
1778  // Set operators. Here the matrix that defines the eigenvalue system
1779  std::vector<Mat> AsMat(As.size());
1780  for (unsigned i=0; i<As.size(); ++i){ AsMat[i]=As[i].mat; }
1781  ierr = PEPSetOperators(pep,AsMat.size(),AsMat.data());CHKERRXX(ierr);
1782  // Set runtime options, e.g.,
1783  ierr = PEPSetFromOptions(pep);CHKERRXX(ierr);
1784  // set convergence type
1785  PEPWhich which=PEP_TARGET_MAGNITUDE;
1786  PetscInt nev=1;
1787  PEPSetWhichEigenpairs(pep,which);
1788  PEPSetDimensions(pep,nev,PETSC_DEFAULT,PETSC_DEFAULT);
1789  if ((max_iter==-1) && (tol==-1.)){ PEPSetTolerances(pep,PETSC_DEFAULT, PETSC_DEFAULT); }
1790  else if(tol==-1.){ PEPSetTolerances(pep,PETSC_DEFAULT, max_iter); }
1791  else if(max_iter==-1){ PEPSetTolerances(pep,tol, PETSC_DEFAULT); }
1792  else{ PEPSetTolerances(pep,tol, max_iter); }
1793  PEPSetTarget(pep,target);
1794  // set spectral transform to shift-and-invert (seems to work best for LST_spatial)
1795  ST st;
1796  PEPGetST(pep,&st);
1797  STSetType(st,STSINVERT);
1798  // set initial guess for right eigenvalue problem
1799  std::vector<Vec> qrvec(1);
1800  qrvec[0] = qr.vec;
1801  ierr = PEPSetInitialSpace(pep,1,qrvec.data());CHKERRXX(ierr);
1802  // now for left (adjoint) eigenvalue problem initial guess
1803  //EPS eps;
1804  //ierr = PEPLinearGetEPS(pep,&eps); CHKERRXX(ierr);
1805  //std::vector<Vec> qlvec(1);
1806  //qlvec[0] = ql.vec;
1807  //ierr = EPSSetLeftInitialSpace(eps,1,qlvec.data()); CHKERRXX(ierr);
1808 
1809  // Solve the system
1810  ierr = PEPSolve(pep);CHKERRXX(ierr);
1811  // output iterations
1812  PetscInt nconv;
1813  ierr = PEPGetConverged(pep,&nconv);CHKERRXX(ierr);
1814  if (nconv>0) {
1815  ierr = PEPGetEigenpair(pep,0,&alpha,PETSC_NULL,eigr_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1816  //ierr = EPSGetLeftEigenvector(eps,0,eigl_vec.vec,PETSC_NULL);CHKERRXX(ierr);
1817  }
1818  // destroy pep
1819  ierr = PEPDestroy(&pep);CHKERRXX(ierr);
1820  // return
1821  return std::make_tuple(alpha,eigr_vec);
1822  }
1823  }
1824  // /** \brief set block matrices using an input array of size rows*cols. Fills rows first \return new matrix with inserted blocks */
1825  //SPIMat block(
1826  // const SPIMat Blocks[], ///< [in] array of matrices to set into larger matrix e.g. [ A, B, C, D ]
1827  // const PetscInt rows, ///< [in] number of rows of submatrices e.g. 2
1828  // const PetscInt cols ///< [in] number of columns of submatrices e.g. 2
1829  // )
1830  // PetscInt m[rows];
1831  // PetscInt msum=Blocks[0].rows;
1832  // PetscInt n[cols];
1833  // PetscInt nsum=Blocks[0].cols;
1834  // m[0]=n[0]=0;
1835  // for (PetscInt i=1; i<rows; ++i)[
1836  // m[i] = m[i-1]+Blocks[i-1].rows;
1837  // msum += m[i];
1838  //
1839  // for (PetscInt j=1; j<cols; ++j)[
1840  // n[j] = n[j-1]+Blocks[j*rows].cols;
1841  // nsum += n[j];
1842  //
1843 
1844  // // TODO check if all rows and columns match for block matrix....
1845 
1846  // SPIMat A(msum,nsum);
1847 
1848  // for (PetscInt j=0; j<cols; ++j)[
1849  // for(PetscInt i=0; i<rows; ++i)[
1850  // A(m[i],n[j],Blocks[i+j*rows]);
1851  //
1852  //
1853 
1854  // return A;
1855  //
1857  SPIMat block(
1858  //std::vector<std::vector<SPIMat>> Blocks ///< [in] array of matrices to set into larger matrix e.g. [ A, B, C, D ]
1859  const Block2D<SPIMat> Blocks
1860  ){
1861  const PetscInt rows=Blocks.size(); // number of rows of submatrices e.g. 2
1862  const PetscInt cols=Blocks[0].size(); // number of columns of submatrices e.g. 2
1863  PetscInt m[rows];
1864  PetscInt msum=Blocks[0][0].rows;
1865  PetscInt n[cols];
1866  PetscInt nsum=Blocks[0][0].cols;
1867  m[0]=n[0]=0;
1868  for (PetscInt i=1; i<rows; ++i) m[i] = m[i-1]+Blocks[i-1][0].rows;
1869  for (PetscInt i=1; i<rows; ++i) msum += Blocks[i][0].rows;
1870  for (PetscInt j=1; j<cols; ++j) nsum += Blocks[0][j].cols;
1871  for (PetscInt j=1; j<cols; ++j) n[j] = n[j-1]+Blocks[0][j-1].cols;
1872 
1873  // TODO check if all rows and columns match for block matrix.... user error catch
1874 
1875  SPIMat A(msum,nsum);
1876  //printf("msum,nsum = %d,%d",msum,nsum);
1877  //for (PetscInt j=0; j<cols; ++j)[
1878  //for(PetscInt i=0; i<rows; ++i)[
1879  //printf("m[%d],n[%d] = %d,%d",i,j,m[i],n[j]);
1880  //]
1881  //]
1882 
1883  for(PetscInt i=0; i<rows; ++i){
1884  for (PetscInt j=0; j<cols; ++j){
1885  //printf("setting A[%d,%d] with shape=%dx%d",m[i],n[j],Blocks[i][j].rows,Blocks[i][j].cols);
1886  A(m[i],n[j],Blocks[i][j],INSERT_VALUES);
1887  }
1888  }
1889  //A();
1890 
1891  return A;
1892  }
1893  /* brief create meshgrid with ij indexing \brief tuple of X and Y matrices */
1894  std::tuple<SPIMat,SPIMat> meshgrid(
1895  SPIVec &x,
1896  SPIVec &y
1897  ){
1898  PetscInt m=x.rows;
1899  PetscInt n=y.rows;
1900  SPIMat X(m,n);
1901  SPIMat Y(m,n);
1902  for(PetscInt i=0; i<m; ++i){
1903  for(PetscInt j=0; j<n; ++j){
1904  X(i,j,x(i,PETSC_TRUE));
1905  Y(i,j,y(j,PETSC_TRUE));
1906  }
1907  }
1908  X();
1909  Y();
1910  return std::make_tuple(X,Y);
1911  }
1912 
1924  PetscInt save(
1925  const SPIMat &A,
1926  const std::string filename
1927  ){
1928  PetscViewer viewer;
1929  //PetscViewerASCIIOpen(PETSC_COMM_WORLD,name.c_str(),&viewer);
1930  PetscErrorCode ierr;
1931  std::ifstream f(filename.c_str());
1932  if(f.good()){
1933  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename.c_str(),FILE_MODE_APPEND,&viewer);CHKERRXX(ierr);
1934  }
1935  else{
1936  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename.c_str(),FILE_MODE_WRITE,&viewer);CHKERRXX(ierr);
1937  }
1938  //ierr = PetscViewerPushFormat(viewer,format);CHKERRXX(ierr);
1939  ierr = MatView(A.mat,viewer);CHKERRXX(ierr);
1940  ierr = PetscViewerDestroy(&viewer);CHKERRXX(ierr);
1941  return 0;
1942  }
1944  PetscInt save(
1945  const std::vector<SPIMat> &As,
1946  const std::string filename
1947  ){
1948  PetscViewer viewer;
1949  //PetscViewerASCIIOpen(PETSC_COMM_WORLD,name.c_str(),&viewer);
1950  PetscErrorCode ierr;
1951  std::ifstream f(filename.c_str());
1952  if(f.good()){
1953  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename.c_str(),FILE_MODE_APPEND,&viewer);CHKERRXX(ierr);
1954  }
1955  else{
1956  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename.c_str(),FILE_MODE_WRITE,&viewer);CHKERRXX(ierr);
1957  }
1958  //ierr = PetscViewerPushFormat(viewer,format);CHKERRXX(ierr);
1959  for(unsigned i=0; i<As.size(); ++i){
1960  ierr = MatView(As[i].mat,viewer);CHKERRXX(ierr);
1961  }
1962  ierr = PetscViewerDestroy(&viewer);CHKERRXX(ierr);
1963  return 0;
1964  }
1966  PetscInt load(
1967  SPIMat &A,
1968  const std::string filename
1969  ){
1970  PetscViewer viewer;
1971  //std::ifstream f(filename.c_str());
1972  A.ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename.c_str(), FILE_MODE_READ, &viewer); CHKERRXX(A.ierr);
1973  A.ierr = MatLoad(A.mat,viewer); CHKERRXX(A.ierr);
1974  A.ierr = PetscViewerDestroy(&viewer); CHKERRXX(A.ierr);
1975  return 0;
1976  }
1978  PetscInt load(
1979  std::vector<SPIMat> &As,
1980  const std::string filename
1981  ){
1982  PetscViewer viewer;
1983  //std::ifstream f(filename.c_str());
1984  PetscErrorCode ierr;
1985  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename.c_str(), FILE_MODE_READ, &viewer); CHKERRXX(ierr);
1986  for(unsigned i=0; i<As.size(); ++i){
1987  ierr = MatLoad(As[i].mat,viewer); CHKERRXX(ierr);
1988  }
1989  ierr = PetscViewerDestroy(&viewer); CHKERRXX(ierr);
1990  return 0;
1991  }
1992 
1994  PetscInt draw(
1995  const SPIMat &A
1996  ){
1997  PetscViewer viewer;
1998  PetscErrorCode ierr;
1999  ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,NULL,A.name.c_str(),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&viewer);CHKERRXX(ierr);
2000  ierr = MatView(A.mat,viewer);CHKERRXX(ierr);
2001 
2002  // pause until user inputs at command line
2003  SPI::printf(" draw(SPIMat) with title=%s, hit ENTER to continue",A.name.c_str());
2004  std::cin.ignore();
2005 
2006  ierr = PetscViewerDestroy(&viewer);CHKERRXX(ierr);
2007  return 0;
2008  }
2010  template <class T>
2011  SPIMat _Function_on_each_element(
2012  T (*f)(T const&),
2013  const SPIMat &A
2014  ){
2015  SPIMat out(A.rows,A.cols);
2016  for (PetscInt i=0; i<out.rows; ++i){
2017  for (PetscInt j=0; j<out.cols; ++j){
2018  out(i,j,(*f)(A(i,j,PETSC_TRUE))); // TODO speed up by getting all values at once on local processor and looping through those
2019  }
2020  }
2021  out();
2022  return out;
2023  }
2025  SPIMat sin(const SPIMat &A){ return _Function_on_each_element(&std::sin<PetscReal>, A); }
2027  SPIMat cos(const SPIMat &A){ return _Function_on_each_element(&std::cos<PetscReal>, A); }
2029  SPIMat acos(const SPIMat &A){ return _Function_on_each_element(&std::acos<PetscReal>, A); }
2031  SPIMat tan(const SPIMat &A){ return _Function_on_each_element(&std::tan<PetscReal>, A); }
2033  SPIMat operator%(const SPIMat &A, const SPIMat &B){
2034  SPIMat out(A.rows,A.cols);
2035  for (PetscInt i=0; i<out.rows; ++i){
2036  for (PetscInt j=0; j<out.cols; ++j){
2037  out(i,j,A(i,j,PETSC_TRUE)*B(i,j,PETSC_TRUE)); // TODO speed up by getting all values at once on local processor and looping through those
2038  }
2039  }
2040  out();
2041  return out;
2042  //return _Function_on_each_element2(&std::multiplies<PetscReal>, A,B);
2043  }
2045  SPIMat abs(const SPIMat &A){
2046  SPIMat out(A.rows,A.cols);
2047  for (PetscInt i=0; i<out.rows; ++i){
2048  for (PetscInt j=0; j<out.cols; ++j){
2049  out(i,j,std::abs(A(i,j,PETSC_TRUE))); // TODO speed up by getting all values at once on local processor and looping through those
2050  }
2051  }
2052  out();
2053  return out;
2054  //return _Function_on_each_element(&std::abs<PetscScalar>, A);
2055  }
2056  /* \brief orthogonalize a basis dense matrix from an array of vec using SLEPc BV */
2057  SPIMat orthogonalize(
2058  const std::vector<SPIVec> &x
2059  ){
2060  SPI::printf("starting orthogonalize");
2061  PetscInt m=x[0].rows; // number of rows
2062  PetscInt n=x.size(); // number of columns
2063  Vec xvec[n];
2064  for(PetscInt i=0; i<n; ++i) xvec[i]=x[i].vec;
2065  SPI::printf(" orthogonalize set vectors");
2066  SPIMat E("E");
2067  // create and initialize BV
2068  BV bv;
2069  E.ierr = BVCreate(PETSC_COMM_WORLD,&bv); CHKERRXX(E.ierr);
2070  E.ierr = BVSetSizesFromVec(bv,x[0].vec,n); CHKERRXX(E.ierr);
2071  E.ierr = BVSetFromOptions(bv);CHKERRXX(E.ierr);
2072  E.ierr = BVInsertVecs(bv,0,&n,xvec,PETSC_TRUE);
2073  SPI::printf(" orthogonalize inserted vecs");
2074  //SPI::SPIMat AorthH("AorthH");
2075  E.ierr = BVOrthogonalize(bv,PETSC_NULL); CHKERRXX(E.ierr);
2076  SPI::printf(" orthogonalize BVOrthogonalize");
2077  E.ierr = BVCreateMat(bv,&E.mat); CHKERRXX(E.ierr);
2078  SPI::printf(" orthogonalize BVCreateMat");
2079 #if USE_GPU == 1
2080  E.ierr = MatConvert(E.mat,MATMPIAIJCUSPARSE,MAT_INPLACE_MATRIX,&E.mat); CHKERRXX(E.ierr);
2081 #else
2082  E.ierr = MatConvert(E.mat,MATMPIAIJ,MAT_INPLACE_MATRIX,&E.mat); CHKERRXX(E.ierr);
2083 #endif
2084  SPI::printf(" orthogonalize MatConvert");
2085  E.rows=m;
2086  E.cols=n;
2087  E();
2088  SPI::printf("done orthogonalize");
2089  return E;
2090  }
2091 
2092 }
2093 
2094 
SPIprint.hpp
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::printf
PetscInt printf(std::string msg,...)
Definition: SPIprint.cpp:5
SPI::SPIVec::flag_init
PetscBool flag_init
flag if it has been initialized
Definition: SPIVec.hpp:25
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::SPIVec::ierr
PetscErrorCode ierr
ierr for various routines and operators
Definition: SPIVec.hpp:22
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
SPIMat.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::SPIVec::rows
PetscInt rows
number of rows in vec
Definition: SPIVec.hpp:14
SPI::SPIMat::set
SPIMat & set(PetscInt m, PetscInt n, const PetscScalar v)
Definition: SPIMat.cpp:77
SPI::SPIVec::dot
PetscScalar dot(SPIVec y)
take the inner product of two vectors
Definition: SPIVec.cpp:459
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::SPIVec::vec
Vec vec
petsc Vec data
Definition: SPIVec.hpp:21
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