19 #ifndef GEOS_LINEARALGEBRA_UTILITIES_ARNOLDI_HPP_
20 #define GEOS_LINEARALGEBRA_UTILITIES_ARNOLDI_HPP_
22 #include "denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp"
35 template<
typename VECTOR >
40 localIndex const mInternal = std::min( numGlobalRows, m );
47 V[0].create( numLocalRows, op.
comm() );
48 V[0].set( 1.0 / sqrt(
static_cast< real64 >( numGlobalRows ) ) );
53 V[j+1].create( numLocalRows, op.
comm() );
54 op.
apply( V[j], V[j+1] );
58 H( i, j ) = V[i].dot( V[j+1] );
59 V[j+1].axpy( -H( i, j ), V[i] );
61 H( j+1, j ) = V[j+1].norm2();
62 V[j+1].scale( 1.0 / H( j+1, j ) );
67 H.resize( mInternal, mInternal );
71 BlasLapackLA::matrixEigenvalues( H, lambda );
77 lambdaMax = std::max( std::abs( lambda[i] ), lambdaMax );
Abstract base class for linear operators.
virtual globalIndex numGlobalRows() const =0
Get the number of global rows.
virtual localIndex numLocalRows() const =0
Get the number of local rows.
virtual void apply(Vector const &src, Vector &dst) const =0
Apply operator to a vector, dst = this(src).
virtual MPI_Comm comm() const =0
Get the MPI communicator the matrix was created with.
Array< T, 2, PERMUTATION > array2d
Alias for 2D array.
double real64
64-bit floating point type.
real64 ArnoldiLargestEigenvalue(LinearOperator< VECTOR > const &op, localIndex const m=4)
Function implementing the Arnoldi scheme to compute the largest eigenvalue.
GEOSX_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Array< T, 1 > array1d
Alias for 1D array.