<a name="lac"></a>
<h3>lac</h3>
+<ol>
+ <li>
+ <p>
+ New: Based on work by Francisco Alvaro, the existing SLEPcWrappers now
+ have a handle on the generalized eigenvalue problem where B=I.
+ <br>
+ (Toby D. Young 2009/06/25)
+ </p>
+ </li>
+</ol>
+
+
<ol>
<li>
<p>
<ol>
<li>
<p>
- New: Based on work with Rickard Armiento, Francisco Alvaro, and Jose E. Roman,
- SLEPcWrappers that give a handle on some of the features of SLEPc (Scalable Library
- for Eigenvalue Problem Computations): (1) The SLEPcWrappers::SolverBase class can be
- used for specifying an eigenvalue problem, either in standard or generalized form,
- on serial or parallel architectures with support for a few solver types; and
- (2) The SLEPcWrappers::TransformationBase class encapsulates a variety of spectral
- transformations providing some functionality required for acceleration techniques
- based on the transformation of the spectrum.
+ New: Based on work with Rickard Armiento, Francisco Alvaro, and Jose
+ E. Roman, SLEPcWrappers that give a handle on some of the features
+ of SLEPc (Scalable Library for Eigenvalue Problem Computations): (1)
+ The SLEPcWrappers::SolverBase class can be used for specifying an
+ eigenvalue problem, either in standard or generalized form, on
+ serial or parallel architectures with support for a few solver
+ types; and (2) The SLEPcWrappers::TransformationBase class
+ encapsulates a variety of spectral transformations providing some
+ functionality required for acceleration techniques based on the
+ transformation of the spectrum.
<br>
(Toby D. Young 2009/06/25)
</p>
* Base class for solver classes using the SLEPc solvers which are
* selected based on flags passed to the eigenvalue problem solver
* context. Derived classes set the right flags to set the right
- * solver -- in particular, note that the AdditionalData structure is
- * a dummy structure and is there for backward compatibility.
+ * solver. On the other hand, note that: the AdditionalData structure
+ * is a dummy structure and is there for backward/forward
+ * compatibility.
*
* SLEPcWrappers can be implemented in application codes in the
* following way:
mpi_communicator);
system.solve (A, B, lambda, x, n_eigenvectors);
@endverbatim
- * for the generalized eigenvalue problem $Ax=B\lambda x$.
+ * for the generalized eigenvalue problem $Ax=B\lambda x$. See also
+ * @ref step 36 "step-36" for a hands-on example.
*
* An alternative implementation to the one above is to use the API
* internals directly within the application code. In this way the
* calling sequence requires calling several of SolverBase functions
- * rather than just one. This freedom is intended for use by users of
- * the SLEPcWrappers that require a greater handle on the eigenvalue
+ * rather than just one. This freedom is intended for use of the
+ * SLEPcWrappers that require a greater handle on the eigenvalue
* problem solver context. See also:
@verbatim
template <typename OutputVector>
const PETScWrappers::MatrixBase &B,
std::vector<double> &kr,
std::vector<OutputVector> &vr,
- const unsigned int n_eigenvectors
- )
+ const unsigned int n_eigenvectors)
+ {<code>code...</code>}
@endverbatim
* as an example on how to do this.
*
/**
* Composite method that solves the
- * linear system $Ax=\lambda
- * Bx$. The eigenvector sent in has
- * to have at least one element
- * that we can use as a template
- * when resizing, since we do not
- * know the parameters of the
- * specific vector class used
+ * eigensystem $Ax=\lambda x$. The
+ * eigenvector sent in has to have
+ * at least one element that we can
+ * use as a template when resizing,
+ * since we do not know the
+ * parameters of the specific
+ * vector class used
* (i.e. local_dofs for MPI
* vectors). However, while copying
* eigenvectors, at least twice the
* reset.
*
* Note that the number of
- * converged eigenstates can be
+ * converged eigenvectors can be
* larger than the number of
- * eigenstates requested; this is
+ * eigenvectors requested; this is
* due to a round off error
* (success) of the eigenvalue
* solver context. If this is found
* requested but handle that it may
* be more by ignoring any extras.
*/
+ template <typename OutputVector>
+ void
+ solve (const PETScWrappers::MatrixBase &A,
+ std::vector<double> &kr,
+ std::vector<OutputVector> &vr,
+ const unsigned int n_eigenvectors);
+
+ /**
+ * Same as above, but here is a
+ * composite method for solving the
+ * system $A x=\lambda B x$.
+ */
template <typename OutputVector>
void
solve (const PETScWrappers::MatrixBase &A,
/**
* Initialize solver for the linear
- * system $Ax=\lambda
- * Bx$. (required before calling
- * solve)
+ * system $Ax=\lambda x$. (Note:
+ * this is required before calling
+ * solve ())
+ */
+ void
+ set_matrices (const PETScWrappers::MatrixBase &A);
+
+ /**
+ * Same as above, but here is a
+ * composite method for solving the
+ * system $A x=\lambda B x$.
*/
void
set_matrices (const PETScWrappers::MatrixBase &A,
const PETScWrappers::MatrixBase &B);
/**
- * Set the initial vector for the solver.
+ * Set the initial vector for the
+ * solver.
*/
void
set_initial_vector (const PETScWrappers::VectorBase &initial_vec);
/**
* Set the spectral transformation
- * to be used. By default SLEPc
+ * to be used.
*/
void
set_transformation (SLEPcWrappers::TransformationBase &trans);
* Solve the linear system for
* n_eigenvectors
* eigenstates. Parameter
- * n_converged contains the actual
- * number of eigenstates that have
- * . converged; this can be both
+ * <code>n_converged</code>
+ * contains the actual number of
+ * eigenstates that have .
+ * converged; this can be both
* fewer or more than
* n_eigenvectors, depending on the
* SLEPc eigensolver used.
* Access to object that controls
* convergence.
*/
- SolverControl & control() const;
+ SolverControl &control() const;
/**
* Exceptions.
// --------------------------- inline and template functions -----------
-
- /**
- * This is declared here to make it
- * possible to take a std::vector
- * of different PETScWrappers vector
- * types
+ /**
+ * This is declared here to make it possible to take a std::vector
+ * of different PETScWrappers vector types
*/
+
+ template <typename OutputVector>
+ void
+ SolverBase::solve (const PETScWrappers::MatrixBase &A,
+ std::vector<double> &kr,
+ std::vector<OutputVector> &vr,
+ const unsigned int n_eigenvectors = 1)
+ {
+ unsigned int n_converged;
+
+ set_matrices(A);
+ solve(n_eigenvectors,&n_converged);
+
+ if (n_converged > n_eigenvectors)
+ {
+ n_converged = n_eigenvectors;
+ }
+
+ AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError());
+ vr.resize(n_converged, vr.front());
+ kr.resize(n_converged);
+
+ for (unsigned int index=0; index < n_converged;
+ ++index)
+ {
+ get_eigenpair(index, kr[index], vr[index]);
+ }
+ }
+
+
template <typename OutputVector>
void
SolverBase::solve (const PETScWrappers::MatrixBase &A,
const PETScWrappers::MatrixBase &B,
std::vector<double> &kr,
std::vector<OutputVector> &vr,
- const unsigned int n_eigenvectors = 0)
+ const unsigned int n_eigenvectors = 1)
{
unsigned int n_converged;
-
+
set_matrices(A,B);
-
solve(n_eigenvectors,&n_converged);
-
+
if (n_converged > n_eigenvectors)
{
n_converged = n_eigenvectors;
}
-
+
AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError());
vr.resize(n_converged, vr.front());
kr.resize(n_converged);
-
+
for (unsigned int index=0; index < n_converged;
++index)
{
get_eigenpair(index, kr[index], vr[index]);
}
}
-
-
+
+
}
DEAL_II_NAMESPACE_CLOSE
if( solver_data != 0 )
solver_data.reset ();
}
-
+
+ void
+ SolverBase::set_matrices (const PETScWrappers::MatrixBase &A)
+ {
+ opA = &A;
+ opB = NULL;
+ }
+
void
SolverBase::set_matrices (const PETScWrappers::MatrixBase &A,
const PETScWrappers::MatrixBase &B)
ierr = EPSCreate (mpi_communicator, &solver_data->eps);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
- AssertThrow (opA && opB, ExcSLEPcWrappersUsageError());
- ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
+ AssertThrow (opA, ExcSLEPcWrappersUsageError());
+ if (opB)
+ ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
+ else
+ ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
- if( ini_vec && ini_vec->size() != 0 )
+
+ if (ini_vec && ini_vec->size() != 0)
{
ierr = EPSSetInitialVector(solver_data->eps, *ini_vec);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
- if( transform )
+ if (transform)
transform->set_context(solver_data->eps);
-
+
// set runtime options.
set_solver_type (solver_data->eps);
// get number of converged
// eigenstates
- ierr = EPSGetConverged (solver_data->eps, reinterpret_cast<int *>(n_converged));
+ ierr = EPSGetConverged (solver_data->eps,
+ reinterpret_cast<int *>(n_converged));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
}
-
-
DEAL_II_NAMESPACE_CLOSE
#else