42#ifndef STOKHOS_LANCZOS_HPP
43#define STOKHOS_LANCZOS_HPP
45#include "Teuchos_Array.hpp"
46#include "Teuchos_SerialDenseVector.hpp"
47#include "Teuchos_SerialDenseMatrix.hpp"
48#include "Teuchos_SerialDenseHelpers.hpp"
52 template <
typename ord_type,
typename val_type>
57 typedef Teuchos::SerialDenseVector<ordinal_type, value_type>
vector_type;
81 v[
j] = a*u1[
j] + b*u2[
j];
89 v[
j] = a*u1[
j] + b*u2[
j] +c*u3[
j];
120 template <
typename vectorspace_type,
typename operator_type>
126 typedef Teuchos::SerialDenseVector<ordinal_type,value_type>
vector_type;
127 typedef Teuchos::SerialDenseMatrix<ordinal_type,value_type>
matrix_type;
131 const vectorspace_type& vs,
132 const operator_type& A,
135 Teuchos::Array<value_type>& alpha,
136 Teuchos::Array<value_type>& beta,
137 Teuchos::Array<value_type>& nrm_sqrd) {
144 u0 = Teuchos::getCol(Teuchos::View, u, 0);
153 nrm_sqrd[i] = vs.inner_product(u1, u1);
159 nrm = vs.inner_product(u1, v);
162 alpha[i] = nrm / nrm_sqrd[i];
166 beta[i] = nrm_sqrd[i] / nrm_sqrd[i-1];
170 u2 = Teuchos::getCol(Teuchos::View, u, i+1);
174 vs.add3(
value_type(1), v, -alpha[i], u1, -beta[i], u0, u2);
181 TEUCHOS_TEST_FOR_EXCEPTION(nrm_sqrd[i] < 0, std::logic_error,
182 "Stokhos::LanczosProjPCEBasis::lanczos(): "
183 <<
" Polynomial " << i <<
" out of " << k
184 <<
" has norm " << nrm_sqrd[i] <<
"!");
195 const vectorspace_type& vs,
196 const operator_type& A,
199 Teuchos::Array<value_type>& alpha,
200 Teuchos::Array<value_type>& beta,
201 Teuchos::Array<value_type>& nrm_sqrd) {
207 u0 = Teuchos::getCol(Teuchos::View, u, 0);
215 beta[i] = std::sqrt(vs.inner_product(u1, u1));
216 u1.scale(1.0/beta[i]);
223 alpha[i] = vs.inner_product(u1, v);
227 u2 = Teuchos::getCol(Teuchos::View, u, i+1);
231 vs.add3(
value_type(1), v, -alpha[i], u1, -beta[i], u0, u2);
238 TEUCHOS_TEST_FOR_EXCEPTION(nrm_sqrd[i] < 0, std::logic_error,
239 "Stokhos::LanczosProjPCEBasis::lanczos(): "
240 <<
" Polynomial " << i <<
" out of " << k
241 <<
" has norm " << nrm_sqrd[i] <<
"!");
256 u0 = Teuchos::getCol(Teuchos::View, u, i);
257 nrm = vs.inner_product(u0, u0);
258 dp = vs.inner_product(u2, u0);
Applies Lanczos procedure to a given matrix.
operator_type::value_type value_type
operator_type::ordinal_type ordinal_type
static void computeNormalized(ordinal_type k, const vectorspace_type &vs, const operator_type &A, const vector_type &u_init, matrix_type &u, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &nrm_sqrd)
Compute Lanczos basis.
static void compute(ordinal_type k, const vectorspace_type &vs, const operator_type &A, const vector_type &u_init, matrix_type &u, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &nrm_sqrd)
Compute Lanczos basis.
Teuchos::SerialDenseVector< ordinal_type, value_type > vector_type
Teuchos::SerialDenseMatrix< ordinal_type, value_type > matrix_type
static void gramSchmidt(ordinal_type k, const vectorspace_type &vs, matrix_type &u, vector_type &u2)
Gram-Schmidt orthogonalization routine.
void add2(const value_type &a, const vector_type &u1, const value_type &b, const vector_type &u2, vector_type &v) const
value_type inner_product(const vector_type &u, const vector_type &v) const
void add3(const value_type &a, const vector_type &u1, const value_type &b, const vector_type &u2, const value_type &c, const vector_type &u3, vector_type &v) const
vector_type create_vector() const
WeightedVectorSpace(const vector_type &weights)
Teuchos::SerialDenseVector< ordinal_type, value_type > vector_type
Top-level namespace for Stokhos classes and functions.