packages/Alea/mc_solvers/GmresPolynomial.cc
//---------------------------------------------------------------------------//
// Build GMRES coefficients from normal equations solution of (A*K)c=b
//---------------------------------------------------------------------------//
Build GMRES coefficients from normal equations solution of (A*K)c=b
Build GMRES coefficients from normal equations solution of ( A * K ) c = b
void GmresPolynomial::buildGmresPolyFromNormalEqns()
{
    Teuchos::RCP<const MV> b = computeInitialKrylovVector();
    Teuchos::RCP<MV> K = computeKrylovStandardBasis(b_m+1, b);
    std::vector<LO> inds(b_m+1);
    for( LO i=0; i<b_m+1; ++i )
        inds[i]=i+1;
    Teuchos::RCP<const MV> AK = MVT::CloneView(*K,inds);
    Teuchos::RCP<SDM> AK_trans_AK( new SDM(b_m+1, b_m+1) );
    Teuchos::RCP<SDM> coeffs(      new SDM(b_m+1, 1) );
    MVT::MvTransMv(1.0,*AK,*AK,*AK_trans_AK);
    MVT::MvTransMv(1.0,*AK,*b,*coeffs);
    SCALAR reg = b_poly_pl->get<SCALAR>("gmres_regularization",0.0);
    for( LO i=0; i<b_m+1; ++i )
    {
        (*AK_trans_AK)(i,i) += reg;
    }
    Teuchos::LAPACK<LO,SCALAR> lapack;
    int info;
    lapack.POSV('L',b_m+1,1,AK_trans_AK->values(),AK_trans_AK->stride(),
        coeffs->values(),coeffs->stride(),&info);
    CHECK( info == 0 );
    b_coeffs.resize(b_m+1);
    std::copy( coeffs->values(), coeffs->values()+b_m+1, b_coeffs.begin() );
}
void GmresPolynomial :: buildGmresPolyFromNormalEqns ( ) { Teuchos :: RCP < const MV > b = computeInitialKrylovVector ( ) ; Teuchos :: RCP < MV > K = computeKrylovStandardBasis ( b_m + 1 , b ) ; std :: vector < LO > inds ( b_m + 1 ) ; for ( LO i = 0 ; i < b_m + 1 ; ++ i ) inds [ i ] = i + 1 ; Teuchos :: RCP < const MV > AK = MVT :: CloneView ( * K , inds ) ; Teuchos :: RCP < SDM > AK_trans_AK ( new SDM ( b_m + 1 , b_m + 1 ) ) ; Teuchos :: RCP < SDM > coeffs ( new SDM ( b_m + 1 , 1 ) ) ; MVT :: MvTransMv ( 1.0 , * AK , * AK , * AK_trans_AK ) ; MVT :: MvTransMv ( 1.0 , * AK , * b , * coeffs ) ; SCALAR reg = b_poly_pl -> get < SCALAR > ( " " , 0.0 ) ; for ( LO i = 0 ; i < b_m + 1 ; ++ i ) { ( * AK_trans_AK ) ( i , i ) += reg ; } Teuchos :: LAPACK < LO , SCALAR > lapack ; int info ; lapack . POSV ( ' ' , b_m + 1 , 1 , AK_trans_AK -> values ( ) , AK_trans_AK -> stride ( ) , coeffs -> values ( ) , coeffs -> stride ( ) , & info ) ; CHECK ( info == 0 ) ; b_coeffs . resize ( b_m + 1 ) ; std :: copy ( coeffs -> values ( ) , coeffs -> values ( ) + b_m + 1 , b_coeffs . begin ( ) ) ; }
Build GMRES coefficients from normal equations solution of (A*K)c=b
Build GMRES coefficients from normal equations solution of ( A * K ) c = b
// Get standard basis for Krylov subspace // Discard first vector (b) // Add regularization to diagonal // Solve using LAPACK Cholesky decomp // Copy coefficients into array in base class