eigen_solver.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Non-inline functions for the eigensolvers
31 
32 #ifdef OOMPH_HAS_MPI
33 #include "mpi.h"
34 #endif
35 
36 //Include cfortran.h and the header for the FORTRAN ARPACK routines
37 #include "cfortran.h"
38 #include "arpack.h"
39 #include "lapack_qz.h"
40 
41 //Oomph-lib headers
42 #include "eigen_solver.h"
43 #include "linear_solver.h"
44 #include "problem.h"
45 
46 
47 namespace oomph
48 {
49 //===============================================================
50 /// Constructor, set default values and set the initial
51 /// linear solver to be superlu
52 //===============================================================
53  ARPACK::ARPACK() : EigenSolver(), Spectrum(1), NArnoldi(30),
54  Small(true), Compute_eigenvectors(true)
56 
57 //===============================================================
58 /// Destructor, delete the default linear solver
59 //===============================================================
61 
62 
63 //==========================================================================
64 /// Use ARPACK to solve an eigen problem that is assembled by elements in
65 /// a mesh in a Problem object.
66 //==========================================================================
67 void ARPACK::solve_eigenproblem(Problem* const &problem_pt,
68  const int &n_eval,
69  Vector<std::complex<double> > &eigenvalue,
70  Vector<DoubleVector > &eigenvector)
71 {
72  bool Verbose = false;
73 
74  //Control parameters
75  int ido=0; //Reverse communication flag
76  char bmat[2]; //Specifies the type of matrix on the RHS
77  //Must be 'I' for standard eigenvalue problem
78  // 'G' for generalised eigenvalue problem
79  int n; //Dimension of the problem
80  char which[3]; //Set which eigenvalues are required.
81  int nev; //Number of eigenvalues desired
82  double tol=0.0; //Stopping criteria
83  int ncv; //Number of columns of the matrix V (Dimension of Arnoldi subspace)
84  int iparam[11]={0,0,0,0,0,0,0,0,0,0,0}; //Control parameters
85  //Pointers to starting locations in the workspace vectors
86  int ipntr[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
87  int info=0; //Setting to zero gives random initial guess for vector to start
88  //Arnoldi iteration
89 
90  //Set up the sizes of the matrix
91  n = problem_pt->ndof(); // Total size of matrix
92  nev = n_eval; //Number of desired eigenvalues
93  ncv = NArnoldi < n ? NArnoldi : n; //Number of Arnoldi vectors allowed
94  //Maximum possible value is max
95  //dimension of matrix
96 
97  //If we don't have enough Arnoldi vectors to compute the desired number
98  //of eigenvalues then complain
99  if(nev > ncv)
100  {
101  std::ostringstream warning_stream;
102  warning_stream << "Number of requested eigenvalues " << nev << "\n"
103  << "is greater than the number of Arnoldi vectors:"
104  << ncv << "\n";
105  //Increase the number of Arnoldi vectors
106  ncv = nev + 10;
107  if(ncv > n) {ncv = n;}
108 
109  warning_stream << "Increasing number of Arnoldi vectors to " << ncv
110  << "\n but you may want to increase further using\n"
111  << "ARPACK::narnoldi()\n"
112  << "which will also get rid of this warning.\n";
113 
114  OomphLibWarning(warning_stream.str(),
115  OOMPH_CURRENT_FUNCTION,
116  OOMPH_EXCEPTION_LOCATION);
117  }
118 
119  //Allocate some workspace
120  int lworkl = 3*ncv*ncv + 6*ncv;
121 
122  //Array that holds the final set of Arnoldi basis vectors
123  double **v = new double*;
124  *v = new double[ncv*n];
125  //Leading dimension of v (n)
126  int ldv = n;
127 
128  //Residual vector
129  double *resid = new double[n];
130  //Workspace vector
131  double *workd = new double[3*n];
132  //Workspace vector
133  double *workl = new double[lworkl];
134 
135  bmat[0] = 'G';
136  //If we require eigenvalues to the left of the shift
137  if(Small) {which[0] = 'S';}
138  //Otherwise we require eigenvalues to the right of the shift
139  else {which[0] = 'L';}
140 
141  //Which part of the eigenvalue interests us
142  switch(Spectrum)
143  {
144  //Imaginary part
145  case -1:
146  which[1] = 'I';
147  break;
148 
149  //Magnitude
150  case 0:
151  which[1] = 'M';
152  break;
153 
154  //Real part
155  case 1:
156  which[1] = 'R';
157  break;
158 
159  default:
160  std::ostringstream error_stream;
161  error_stream
162  << "Spectrum is set to an invalid value. It must be 0, 1 or -1\n"
163  << "Not " << Spectrum << std::endl;
164 
165  throw OomphLibError(error_stream.str(),
166  OOMPH_CURRENT_FUNCTION,
167  OOMPH_EXCEPTION_LOCATION);
168  }
169 
170 // %---------------------------------------------------%
171 // | This program uses exact shifts with respect to |
172 // | the current Hessenberg matrix (IPARAM(1) = 1). |
173 // | IPARAM(3) specifies the maximum number of Arnoldi |
174 // | iterations allowed. Mode 3 of DNAUPD is used |
175 // | (IPARAM(7) = 3). All these options can be changed |
176 // | by the user. For details see the documentation in |
177 // | DNAUPD. |
178 // %---------------------------------------------------%
179 
180  int ishifts = 1;
181  int maxitr = 300;
182  int mode = 3; //M is symetric and semi-definite
183 
184  iparam[0] = ishifts; //Exact shifts wrt Hessenberg matrix H
185  iparam[2] = maxitr; //Maximum number of allowed iterations
186  iparam[3] = 1; //Bloacksize to be used in the recurrence
187  iparam[6] = mode; //Mode is shift and invert in real arithmetic
188 
189  //Real and imaginary parts of the shift
190  double sigmar = Sigma_real, sigmai = 0.0;
191 
192  // TEMPORARY
193  // only use non distributed matrice and vectors
194  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n,false);
195  this->build_distribution(dist);
196 
197  //Before we get into the Arnoldi loop solve the shifted matrix problem
198  //Allocated Row compressed matrices for the mass matrix and shifted main
199  //matrix
200  CRDoubleMatrix M(this->distribution_pt()), AsigmaM(this->distribution_pt());
201 
202  //Assemble the matrices
203  problem_pt->get_eigenproblem_matrices(M,AsigmaM,sigmar);
204 
205  //Allocate storage for the vectors to be used in matrix vector products
206  DoubleVector rhs(this->distribution_pt(),0.0);
207  DoubleVector x(this->distribution_pt(),0.0);
208 
209 
210  bool LOOP_FLAG=true;
211  bool First=true;
212 
213  //Enable resolves for the linear solver
215 
216  //Do not report the time taken
218 
219  do
220  {
221 
222 // %---------------------------------------------%
223 // | Repeatedly call the routine DNAUPD and take |
224 // | actions indicated by parameter IDO until |
225 // | either convergence is indicated or maxitr |
226 // | has been exceeded. |
227 // %---------------------------------------------%
228 //
229 
230  DNAUPD(ido, bmat, n, which, nev, tol, resid,
231  ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
232 
233  if(ido == -1)
234  {
235 
236 // %-------------------------------------------%
237 // | Perform matrix vector multiplication |
238 // | y <--- OP*x |
239 // | The user should supply his/her own |
240 // | matrix vector multiplication routine here |
241 // | that takes workd(ipntr(1)) as the input |
242 // | vector, and return the matrix vector |
243 // | product to workd(ipntr(2)). |
244 // %-------------------------------------------%
245 //
246 
247  //Firstly multiply by mx
248  for(int i=0;i<n;i++) {x[i] = workd[ipntr[0]-1+i];}
249  M.multiply(x,rhs);
250 
251  //Now solve the system
252  DoubleVector temp(rhs);
253  if(First) {Linear_solver_pt->solve(&AsigmaM,temp,rhs); First=false;}
254  else {Linear_solver_pt->resolve(temp,rhs);}
255  temp.clear();
256  //AsigmaM.lubksub(rhs);
257  //Put the solution into the workspace
258  for(int i=0;i<n;i++)
259  {
260  workd[ipntr[1]-1+i] = rhs[i];
261  }
262 
263  continue;
264  }
265  else if(ido == 1)
266  {
267  //Already done multiplication by Mx
268  //Need to load the rhs vector
269  for(int i=0;i<n;i++)
270  {
271  rhs[i] = workd[ipntr[2]-1+i];
272  }
273  //Now solve the system
274  //AsigmaM.lubksub(rhs);
275  DoubleVector temp(rhs);
276  if(First) {Linear_solver_pt->solve(&AsigmaM,temp,rhs); First=false;}
277  else {Linear_solver_pt->resolve(temp,rhs);}
278  //Put the solution into the workspace
279  for(int i=0;i<n;i++)
280  {
281  workd[ipntr[1]-1+i] = rhs[i];
282  }
283  continue;
284  }
285  else if(ido == 2)
286  {
287  //Need to multiply by mass matrix
288  //Vector<double> x(n);
289  for(int i=0;i<n;i++) {x[i] = workd[ipntr[0]-1+i];}
290  M.multiply(x,rhs);
291 
292  for(int i=0;i<n;i++)
293  {
294  workd[ipntr[1]-1+i] = rhs[i];
295  }
296  continue;
297  }
298 
299  if(info < 0)
300  {
301  oomph_info << "ERROR" << std::endl;
302  oomph_info << "Error with _naupd, info = '" << info << std::endl;
303  if(info==-7) {oomph_info << "Not enough memory " << std::endl;}
304  }
305 
306 // %-------------------------------------------%
307 // | No fatal errors occurred. |
308 // | Post-Process using DNEUPD. |
309 // | |
310 // | Computed eigenvalues may be extracted. |
311 // | |
312 // | Eigenvectors may also be computed now if |
313 // | desired. (indicated by rvec = .true.) |
314 // %-------------------------------------------%
315 //
316  LOOP_FLAG = false;
317  } while(LOOP_FLAG);
318 
319  //Report any problems
320  if (info == 1)
321  {
322  oomph_info << "Maximum number of iterations reached." << std::endl;
323  }
324  else if(info == 3)
325  {
326  oomph_info
327  << "No shifts could be applied during implicit Arnoldi "
328  << "update, try increasing NCV. "
329  << std::endl;
330  }
331 
332  //Disable resolves for the linear solver
334 
335  //Compute Ritz or Schur vectors, if desired
336  int rvec = Compute_eigenvectors;
337  //Specify the form of the basis computed
338  char howmany[2];
339  howmany[0] = 'A';
340  //Find out the number of converged eigenvalues
341  int nconv = iparam[4];
342 
343  //Note that there is an error (feature) in ARPACK that
344  //means in certain cases, if we request eigenvectors,
345  //consequent Schur factorization of the matrix spanned
346  //by the Arnoldi vectors leads to more converged eigenvalues
347  //than reported by DNAPUD. This is a pain because it
348  //causes a segmentation fault and in every case that I've found
349  //the eigenvalues/vectors are not those desired.
350  //At the moment, I'm just going to leave it, but I note
351  //the problem here to remind myself about it.
352 
353  //Array to select which Ritz vectors are computed
354  int select[ncv];
355  //Array that holds the real part of the eigenvalues
356  double *eval_real = new double[nconv+1];
357  //Array that holds the imaginary part of the eigenvalues
358  double *eval_imag = new double[nconv+1];
359 
360  //Workspace
361  double *workev = new double[3*ncv];
362  //Array to hold the eigenvectors
363  double **z = new double*;
364  *z = new double[(nconv+1)*n];
365 
366  //Error flag
367  int ierr;
368 
369  DNEUPD ( rvec, howmany, select, eval_real, eval_imag, z, ldv,
370  sigmar, sigmai, workev, bmat, n, which, nev, tol,
371  resid, ncv, v, ldv, iparam, ipntr, workd, workl,
372  lworkl, ierr );
373 //
374 // %-----------------------------------------------%
375 // | The real part of the eigenvalue is returned |
376 // | in the first column of the two dimensional |
377 // | array D, and the imaginary part is returned |
378 // | in the second column of D. The corresponding |
379 // | eigenvectors are returned in the first NEV |
380 // | columns of the two dimensional array V if |
381 // | requested. Otherwise, an orthogonal basis |
382 // | for the invariant subspace corresponding to |
383 // | the eigenvalues in D is returned in V. |
384 // %-----------------------------------------------%
385 
386  if(ierr != 0)
387  {
388  oomph_info << "Error with _neupd, info = " << ierr << std::endl;
389  }
390  else
391  //Print the output anyway
392  {
393  //Resize the eigenvalue output vector
394  eigenvalue.resize(nconv);
395  for(int j=0;j<nconv;j++)
396  {
397  //Turn the output from ARPACK into a complex number
398  std::complex<double> eigen(eval_real[j],eval_imag[j]);
399  //Add the eigenvalues to the output vector
400  eigenvalue[j] = eigen;
401  }
402 
404  {
405  //Load the eigenvectors
406  eigenvector.resize(nconv);
407  for(int j=0;j<nconv;j++)
408  {
409  eigenvector[j].build(this->distribution_pt(),0.0);
410  for(int i=0;i<n;i++)
411  {
412  eigenvector[j][i] = z[0][j*n+i];
413  }
414  }
415  }
416  else
417  {
418  eigenvector.resize(0);
419  }
420 
421  //Report the information
422  if(Verbose)
423  {
424  oomph_info << " ARPACK " << std::endl;
425  oomph_info << " ====== " << std::endl;
426  oomph_info << " Size of the matrix is " << n << std::endl;
427  oomph_info << " The number of Ritz values requested is " << nev
428  << std::endl;
429  oomph_info << " The number of Arnoldi vectors generated (NCV) is " << ncv
430  << std::endl;
431  oomph_info << " What portion of the spectrum: " << which << std::endl;
432  oomph_info << " The number of converged Ritz values is "
433  << nconv << std::endl;
434  oomph_info
435  << " The number of Implicit Arnoldi update iterations taken is "
436  << iparam[2] << std::endl;
437  oomph_info << " The number of OP*x is " << iparam[8] << std::endl;
438  oomph_info << " The convergence criterion is " << tol << std::endl;
439  }
440  }
441 
442  //Clean up the allocated memory
443  delete[] *z; *z=0;
444  delete z; z=0;
445  delete[] workev; workev=0;
446  delete[] eval_imag; eval_imag=0;
447  delete[] eval_real; eval_real=0;
448 
449  //Clean up the memory
450  delete[] workl; workl = 0;
451  delete[] workd; workd = 0;
452  delete[] resid; resid = 0;
453  delete[] *v; *v=0;
454  delete v; v=0;
455 }
456 
457 
458 //==========================================================================
459 /// Use LAPACK to solve an eigen problem that is assembled by elements in
460 /// a mesh in a Problem object.
461 //==========================================================================
462  void LAPACK_QZ::solve_eigenproblem(Problem* const &problem_pt,
463  const int &n_eval,
464  Vector<std::complex<double> > &eigenvalue,
465  Vector<DoubleVector> &eigenvector)
466 {
467  //Some character identifiers for use in the LAPACK routine
468  //Do not calculate the left eigenvectors
469  char no_eigvecs[2] = "N";
470  //Do caculate the eigenvectors
471  char eigvecs[2] = "V";
472 
473  //Get the dimension of the matrix
474  int n = problem_pt->ndof(); // Total size of matrix
475 
476  //Initialise a padding integer
477  int padding=0;
478  //If the dimension of the matrix is even, then pad the arrays to
479  //make the size odd. This somehow sorts out a strange run-time behaviour
480  //identified by Rich Hewitt.
481  if(n%2==0) {padding=1;}
482 
483  //Get the real and imaginary parts of the shift
484  double sigmar = Sigma_real;// sigmai = 0.0;
485 
486  //Actual size of matrix that will be allocated
487  int padded_n = n + padding;
488 
489  //Storage for the matrices in the packed form required by the LAPACK
490  //routine
491  double *M = new double[padded_n*padded_n];
492  double *A = new double[padded_n*padded_n];
493 
494  // TEMPORARY
495  // only use non-distributed matrices and vectors
496  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n,false);
497  this->build_distribution(dist);
498 
499  //Enclose in a separate scope so that memory is cleaned after assembly
500  {
501  //Allocated Row compressed matrices for the mass matrix and shifted main
502  //matrix
503  CRDoubleMatrix temp_M(this->distribution_pt()),
504  temp_AsigmaM(this->distribution_pt());
505 
506  //Assemble the matrices
507  problem_pt->get_eigenproblem_matrices(temp_M,temp_AsigmaM,sigmar);
508 
509  //Now convert these matrices into the appropriate packed form
510  unsigned index=0;
511  for(int i=0;i<n;++i)
512  {
513  for(int j=0;j<n;++j)
514  {
515  M[index] = temp_M(j,i);
516  A[index] = temp_AsigmaM(j,i);
517  ++index;
518  }
519  //If necessary pad the columns with zeros
520  if(padding)
521  {
522  M[index] = 0.0;
523  A[index] = 0.0;
524  ++index;
525  }
526  }
527  //No need to pad the final row because it is never accessed by the
528  //routine.
529  }
530 
531  // Temporary eigenvalue storage
532  double* alpha_r = new double[n];
533  double* alpha_i = new double[n];
534  double* beta = new double[n];
535  // Temporary eigenvector storage
536  double* vec_left = new double[1];
537  double* vec_right = new double[n*n];
538 
539  // Workspace for the LAPACK routine
540  std::vector<double> work(1,0.0);
541  //Info flag for the LAPACK routine
542  int info = 0;
543 
544  // Call FORTRAN LAPACK to get the required workspace
545  // Note the use of the padding flag
546  LAPACK_DGGEV( no_eigvecs, eigvecs, n, &A[0], padded_n, &M[0], padded_n,
547  alpha_r, alpha_i, beta, vec_left, 1, vec_right, n,
548  &work[0], -1, info );
549 
550  //Get the amount of requires workspace
551  int required_workspace = ( int ) work[0];
552  //If we need it
553  work.resize( required_workspace );
554 
555 // call FORTRAN LAPACK again with the optimum workspace
556  LAPACK_DGGEV( no_eigvecs, eigvecs, n, &A[0], padded_n, &M[0], padded_n,
557  alpha_r, alpha_i, beta, vec_left, 1, vec_right, n, &work[0],
558  required_workspace, info );
559 
560  //Now resize storage for the eigenvalues and eigenvectors
561  //We get them all!
562  eigenvalue.resize(n);
563  eigenvector.resize(n);
564 
565  //Now convert the output into our format
566  for (int i=0; i<n; ++i )
567  {
568  //N.B. This is naive and dangerous according to the documentation
569  //beta could be very small giving over or under flow
570  //Remember to fix the shift back again
571  eigenvalue[i] =
572  std::complex<double>(sigmar + alpha_r[i]/beta[i],alpha_i[i]/beta[i]);
573 
574  //Resize the eigenvector storage
575  eigenvector[i].build(this->distribution_pt(),0.0);
576  //Load up the eigenvector (assume that it's real)
577  for(int k = 0; k < n; ++k )
578  {
579  eigenvector[i][k] = vec_right[i*n + k ];
580  }
581  }
582 
583  // Delete all allocated storage
584  delete[] vec_right;
585  delete[] vec_left;
586  delete[] beta;
587  delete[] alpha_r;
588  delete[] alpha_i;
589  delete[] A;
590  delete[] M;
591 }
592 
593 
594 //==========================================================================
595 /// Use LAPACK to solve a complex eigen problem specified by the given
596 /// matrices.
597 //==========================================================================
599  const ComplexMatrixBase &M,
600  Vector<std::complex<double> > &eigenvalue,
601  Vector<Vector<std::complex<double> > >
602  &eigenvector)
603 {
604  //Some character identifiers for use in the LAPACK routine
605  //Do not calculate the left eigenvectors
606  char no_eigvecs[2] = "N";
607  //Do caculate the eigenvectors
608  char eigvecs[2] = "V";
609 
610  //Get the dimension of the matrix
611  int n = A.nrow(); // Total size of matrix
612 
613  //Storage for the matrices in the packed form required by the LAPACK
614  //routine
615  double *M_linear = new double[2*n*n];
616  double *A_linear = new double[2*n*n];
617 
618  //Now convert the matrices into the appropriate packed form
619  unsigned index=0;
620  for(int i=0;i<n;++i)
621  {
622  for(int j=0;j<n;++j)
623  {
624  M_linear[index] = M(j,i).real();
625  A_linear[index] = A(j,i).real();
626  ++index;
627  M_linear[index] = M(j,i).imag();
628  A_linear[index] = A(j,i).imag();
629  ++index;
630  }
631  }
632 
633  // Temporary eigenvalue storage
634  double* alpha = new double[2*n];
635  double* beta = new double[2*n];
636  // Temporary eigenvector storage
637  double* vec_left = new double[2];
638  double* vec_right = new double[2*n*n];
639 
640  // Workspace for the LAPACK routine
641  std::vector<double> work(2,0.0);
642  std::vector<double> rwork(8*n,0.0);
643 
644  //Info flag for the LAPACK routine
645  int info = 0;
646 
647  // Call FORTRAN LAPACK to get the required workspace
648  // Note the use of the padding flag
649  LAPACK_ZGGEV( no_eigvecs, eigvecs, n, &A_linear[0], n, &M_linear[0], n,
650  alpha, beta, vec_left, 1, vec_right, n,
651  &work[0], -1, &rwork[0], info );
652 
653  //Get the amount of requires workspace
654  int required_workspace = ( int ) work[0];
655  //If we need it
656  work.resize( 2*required_workspace );
657 
658 // call FORTRAN LAPACK again with the optimum workspace
659  LAPACK_ZGGEV( no_eigvecs, eigvecs, n, &A_linear[0], n, &M_linear[0], n,
660  alpha, beta, vec_left, 1, vec_right, n, &work[0],
661  required_workspace, &rwork[0], info );
662 
663  //Now resize storage for the eigenvalues and eigenvectors
664  //We get them all!
665  eigenvalue.resize(n);
666  eigenvector.resize(n);
667 
668  //Now convert the output into our format
669  for (int i=0; i<n; ++i )
670  {
671  //We have temporary complex numbers
672  std::complex<double> num(alpha[2*i],alpha[2*i+1]);
673  std::complex<double> den(beta[2*i],beta[2*i+1]);
674 
675  //N.B. This is naive and dangerous according to the documentation
676  //beta could be very small giving over or under flow
677  eigenvalue[i] = num/den;
678 
679  //Resize the eigenvector storage
680  eigenvector[i].resize(n);
681  //Load up the eigenvector (assume that it's real)
682  for(int k = 0; k < n; ++k )
683  {
684  eigenvector[i][k] = std::complex<double>(vec_right[2*i*n + 2*k],
685  vec_right[2*i*n + 2*k+1]);
686  }
687  }
688 
689  // Delete all allocated storage
690  delete[] vec_right;
691  delete[] vec_left;
692  delete[] beta;
693  delete[] alpha;
694  delete[] A_linear;
695  delete[] M_linear;
696 }
697 }
bool Compute_eigenvectors
Boolean to indicate whether or not to compute the eigenvectors.
Definition: eigen_solver.h:127
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
Abstract base class for matrices of complex doubles – adds abstract interfaces for solving...
void clear()
wipes the DoubleVector
virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
Get the matrices required by a eigensolver. If the shift parameter is non-zero the second matrix will...
Definition: problem.cc:8337
cstr elem_len * i
Definition: cfortran.h:607
The Problem class.
Definition: problem.h:152
virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the rhs vector. Solution is returned in...
int NArnoldi
Number of Arnoldi vectors to compute.
Definition: eigen_solver.h:119
OomphInfo oomph_info
. SuperLU Project Solver class. This is a combined wrapper for both SuperLU and SuperLU Dist...
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
ARPACK()
Constructor.
Definition: eigen_solver.cc:53
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
LinearSolver * Linear_solver_pt
Pointer to a linear solver.
Definition: eigen_solver.h:108
bool Small
Boolean to set which part of the spectrum left (default) or right of the shifted value.
Definition: eigen_solver.h:124
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
Definition: eigen_solver.h:71
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
void disable_doc_time()
Disable documentation of solve times.
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
virtual ~ARPACK()
Destructor, delete the linear solver.
Definition: eigen_solver.cc:60
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
Definition: eigen_solver.cc:67
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
int Spectrum
Integer to set whether the real, imaginary or magnitude is required to be small or large...
Definition: eigen_solver.h:115
void find_eigenvalues(const ComplexMatrixBase &A, const ComplexMatrixBase &M, Vector< std::complex< double > > &eigenvalue, Vector< Vector< std::complex< double > > > &eigenvector)
LinearSolver * Default_linear_solver_pt
Pointer to a default linear solver.
Definition: eigen_solver.h:111
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872