packages/Alea/mc_solvers/GmresPolynomial.cc
//---------------------------------------------------------------------------//
// Convert roots of polynomial to coefficients
//---------------------------------------------------------------------------//
Convert roots of polynomial to coefficients
Convert roots of polynomial to coefficients
Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SCALAR> >
GmresPolynomial::computeArnoldiMatrix(int order) const
{
    REQUIRE( b_A  != Teuchos::null );
    REQUIRE( order > 0 );
    if( b_verbosity >= HIGH )
        std::cout << "Computing GMRES roots" << std::endl;
    MV V(b_A->getDomainMap(),order+1);
    bool zero_out = true;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SCALAR> > H(
        new Teuchos::SerialDenseMatrix<LO,SCALAR>(order+1,order,zero_out) );
    std::vector<SCALAR> dots(1), vec_norm(1);
    std::vector<int> v_ind(1), q_ind(1);
    v_ind[0] = 0;
    Teuchos::RCP<MV> v = MVT::CloneViewNonConst(V,v_ind);
    Teuchos::RCP<const MV> q;
    Teuchos::RCP<const MV> b = computeInitialKrylovVector();
    MVT::SetBlock(*b,v_ind,*v);
    for( int n=0; n<order; ++n )
    {
        v_ind[0] = n+1;
        v = MVT::CloneViewNonConst(V,v_ind);
        q_ind[0] = n;
        q = MVT::CloneView(V,q_ind);
        OPT::Apply(*b_A,*q,*v);
        for( int j=0; j<n+1; ++j )
        {
            q_ind[0] = j;
            q = MVT::CloneView(V,q_ind);
            MVT::MvDot(*v,*q,dots);
            (*H)(j,n) = dots[0];
            MVT::MvAddMv(SCALAR_TRAITS::one(),*v,-(*H)(j,n),*q,*v);
        }
        MVT::MvNorm(*v,vec_norm);
        (*H)(n+1,n) = vec_norm[0];
        MVT::MvScale(*v,1.0/vec_norm[0]);
    }
    return H;
}
Teuchos :: RCP < Teuchos :: SerialDenseMatrix < LO , SCALAR > > GmresPolynomial :: computeArnoldiMatrix ( int order ) const { REQUIRE ( b_A != Teuchos :: null ) ; REQUIRE ( order > 0 ) ; if ( b_verbosity >= HIGH ) std :: cout << " " << std :: endl ; MV V ( b_A -> getDomainMap ( ) , order + 1 ) ; bool zero_out = true ; Teuchos :: RCP < Teuchos :: SerialDenseMatrix < LO , SCALAR > > H ( new Teuchos :: SerialDenseMatrix < LO , SCALAR > ( order + 1 , order , zero_out ) ) ; std :: vector < SCALAR > dots ( 1 ) , vec_norm ( 1 ) ; std :: vector < int > v_ind ( 1 ) , q_ind ( 1 ) ; v_ind [ 0 ] = 0 ; Teuchos :: RCP < MV > v = MVT :: CloneViewNonConst ( V , v_ind ) ; Teuchos :: RCP < const MV > q ; Teuchos :: RCP < const MV > b = computeInitialKrylovVector ( ) ; MVT :: SetBlock ( * b , v_ind , * v ) ; for ( int n = 0 ; n < order ; ++ n ) { v_ind [ 0 ] = n + 1 ; v = MVT :: CloneViewNonConst ( V , v_ind ) ; q_ind [ 0 ] = n ; q = MVT :: CloneView ( V , q_ind ) ; OPT :: Apply ( * b_A , * q , * v ) ; for ( int j = 0 ; j < n + 1 ; ++ j ) { q_ind [ 0 ] = j ; q = MVT :: CloneView ( V , q_ind ) ; MVT :: MvDot ( * v , * q , dots ) ; ( * H ) ( j , n ) = dots [ 0 ] ; MVT :: MvAddMv ( SCALAR_TRAITS :: one ( ) , * v , - ( * H ) ( j , n ) , * q , * v ) ; } MVT :: MvNorm ( * v , vec_norm ) ; ( * H ) ( n + 1 , n ) = vec_norm [ 0 ] ; MVT :: MvScale ( * v , 1.0 / vec_norm [ 0 ] ) ; } return H ; }
Convert roots of polynomial to coefficients
Convert roots of polynomial to coefficients
// Do a quick Arnoldi iteration to compute Hessenberg matrix H // Multivector to hold basis vectors // Dense upper Hessenberg matrix // Set initial vector // Get initial vector // Start Arnoldi iteration: Trefethen & Bau Algorithm 33.1 // v = A*q_n // H_jn = q_j^* v // v = v - H_jn * q_j // H_n+1,n = ||v|| // q_n+1 = v / ||v||