54 Compute_eigenvectors(true)
74 Vector<std::complex<double>>& eigenvalue,
90 int iparam[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
92 int ipntr[14] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
98 n = problem_pt->
ndof();
108 std::ostringstream warning_stream;
109 warning_stream <<
"Number of requested eigenvalues " << nev <<
"\n"
110 <<
"is greater than the number of Arnoldi vectors:" << ncv
119 warning_stream <<
"Increasing number of Arnoldi vectors to " << ncv
120 <<
"\n but you may want to increase further using\n"
121 <<
"ARPACK::narnoldi()\n"
122 <<
"which will also get rid of this warning.\n";
125 warning_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
129 int lworkl = 3 * ncv * ncv + 6 * ncv;
132 double** v =
new double*;
133 *v =
new double[ncv * n];
138 double* resid =
new double[n];
140 double* workd =
new double[3 * n];
142 double* workl =
new double[lworkl];
175 std::ostringstream error_stream;
177 <<
"Spectrum is set to an invalid value. It must be 0, 1 or -1\n"
181 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
224 bool LOOP_FLAG =
true;
274 for (
int i = 0;
i < n;
i++)
276 x[
i] = workd[ipntr[0] - 1 +
i];
294 for (
int i = 0;
i < n;
i++)
296 workd[ipntr[1] - 1 +
i] = rhs[
i];
305 for (
int i = 0;
i < n;
i++)
307 rhs[
i] = workd[ipntr[2] - 1 +
i];
322 for (
int i = 0;
i < n;
i++)
324 workd[ipntr[1] - 1 +
i] = rhs[
i];
332 for (
int i = 0;
i < n;
i++)
334 x[
i] = workd[ipntr[0] - 1 +
i];
338 for (
int i = 0;
i < n;
i++)
340 workd[ipntr[1] - 1 +
i] = rhs[
i];
348 oomph_info <<
"Error with _naupd, info = '" << info << std::endl;
351 oomph_info <<
"Not enough memory " << std::endl;
371 oomph_info <<
"Maximum number of iterations reached." << std::endl;
375 oomph_info <<
"No shifts could be applied during implicit Arnoldi "
376 <<
"update, try increasing NCV. " << std::endl;
388 int nconv = iparam[4];
403 double* eval_real =
new double[nconv + 1];
405 double* eval_imag =
new double[nconv + 1];
408 double* workev =
new double[3 * ncv];
410 double** z =
new double*;
411 *z =
new double[(nconv + 1) * n];
456 oomph_info <<
"Error with _neupd, info = " << ierr << std::endl;
462 eigenvalue.resize(nconv);
463 for (
int j = 0; j < nconv; j++)
466 std::complex<double> eigen(eval_real[j], eval_imag[j]);
468 eigenvalue[j] = eigen;
474 eigenvector.resize(nconv);
475 for (
int j = 0; j < nconv; j++)
478 for (
int i = 0;
i < n;
i++)
480 eigenvector[j][
i] = z[0][j * n +
i];
486 eigenvector.resize(0);
494 oomph_info <<
" Size of the matrix is " << n << std::endl;
495 oomph_info <<
" The number of Ritz values requested is " << nev
497 oomph_info <<
" The number of Arnoldi vectors generated (NCV) is "
499 oomph_info <<
" What portion of the spectrum: " << which << std::endl;
500 oomph_info <<
" The number of converged Ritz values is " << nconv
503 <<
" The number of Implicit Arnoldi update iterations taken is "
504 << iparam[2] << std::endl;
505 oomph_info <<
" The number of OP*x is " << iparam[8] << std::endl;
506 oomph_info <<
" The convergence criterion is " << tol << std::endl;
542 Vector<std::complex<double>>& eigenvalue,
547 char no_eigvecs[2] =
"N";
549 char eigvecs[2] =
"V";
552 int n = problem_pt->
ndof();
568 int padded_n = n + padding;
572 double* M =
new double[padded_n * padded_n];
573 double* A =
new double[padded_n * padded_n];
592 for (
int i = 0;
i < n; ++
i)
594 for (
int j = 0; j < n; ++j)
596 M[index] = temp_M(j,
i);
597 A[index] = temp_AsigmaM(j,
i);
613 double* alpha_r =
new double[n];
614 double* alpha_i =
new double[n];
615 double* beta =
new double[n];
617 double* vec_left =
new double[1];
618 double* vec_right =
new double[n * n];
621 std::vector<double> work(1, 0.0);
627 LAPACK_DGGEV(no_eigvecs,
646 int required_workspace = (int)work[0];
648 work.resize(required_workspace);
651 LAPACK_DGGEV(no_eigvecs,
671 eigenvalue.resize(n);
672 eigenvector.resize(n);
675 for (
int i = 0;
i < n; ++
i)
680 eigenvalue[
i] = std::complex<double>(sigmar + alpha_r[
i] / beta[
i],
681 alpha_i[
i] / beta[
i]);
686 for (
int k = 0; k < n; ++k)
688 eigenvector[
i][k] = vec_right[
i * n + k];
710 Vector<std::complex<double>>& eigenvalue,
715 char no_eigvecs[2] =
"N";
717 char eigvecs[2] =
"V";
724 double* M_linear =
new double[2 * n * n];
725 double* A_linear =
new double[2 * n * n];
729 for (
int i = 0;
i < n; ++
i)
731 for (
int j = 0; j < n; ++j)
733 M_linear[index] = M(j,
i).real();
734 A_linear[index] = A(j,
i).real();
736 M_linear[index] = M(j,
i).imag();
737 A_linear[index] = A(j,
i).imag();
743 double* alpha =
new double[2 * n];
744 double* beta =
new double[2 * n];
746 double* vec_left =
new double[2];
747 double* vec_right =
new double[2 * n * n];
750 std::vector<double> work(2, 0.0);
751 std::vector<double> rwork(8 * n, 0.0);
758 LAPACK_ZGGEV(no_eigvecs,
777 int required_workspace = (int)work[0];
779 work.resize(2 * required_workspace);
782 LAPACK_ZGGEV(no_eigvecs,
802 eigenvalue.resize(n);
803 eigenvector.resize(n);
806 for (
int i = 0;
i < n; ++
i)
809 std::complex<double> num(alpha[2 *
i], alpha[2 *
i + 1]);
810 std::complex<double> den(beta[2 *
i], beta[2 *
i + 1]);
814 eigenvalue[
i] = num / den;
817 eigenvector[
i].resize(n);
819 for (
int k = 0; k < n; ++k)
821 eigenvector[
i][k] = std::complex<double>(
822 vec_right[2 *
i * n + 2 * k], vec_right[2 *
i * n + 2 * k + 1]);
int Spectrum
Integer to set whether the real, imaginary or magnitude is required to be small or large.
int NArnoldi
Number of Arnoldi vectors to compute.
LinearSolver * Linear_solver_pt
Pointer to a linear solver.
virtual ~ARPACK()
Destructor, delete the linear solver.
bool Small
Boolean to set which part of the spectrum left (default) or right of the shifted value.
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
LinearSolver * Default_linear_solver_pt
Pointer to a default linear solver.
bool Compute_eigenvectors
Boolean to indicate whether or not to compute the eigenvectors.
A class for compressed row matrices. This is a distributable object.
Abstract base class for matrices of complex doubles – adds abstract interfaces for solving,...
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
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....
void clear()
wipes the DoubleVector
Base class for all EigenProblem solves. This simply defines standard interfaces so that different sol...
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
void find_eigenvalues(const ComplexMatrixBase &A, const ComplexMatrixBase &M, Vector< std::complex< double > > &eigenvalue, Vector< Vector< std::complex< double > > > &eigenvector)
Find the eigenvalues of a generalised eigenvalue problem specified by .
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
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...
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
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...
void disable_doc_time()
Disable documentation of solve times.
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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...
unsigned long ndof() const
Return the number of dofs.
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
A slight extension to the standard template vector class so that we can include "graceful" array rang...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...