packages/Alea/mc_solvers/GmresPolynomial.cc
//---------------------------------------------------------------------------//
// Compute GMRES polynomial roots from Arnoldi matrix
//---------------------------------------------------------------------------//
Compute GMRES polynomial roots from Arnoldi matrix
Compute GMRES polynomial roots from Arnoldi matrix
Teuchos::ArrayRCP<const Anasazi::Value<SCALAR> >
GmresPolynomial::computeGmresRoots(
    Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SCALAR> > H ) const
{
    LO order = H->numCols();
    REQUIRE( H->numRows() == order+1 );
    Teuchos::SerialDenseMatrix<LO,SCALAR> H_square(
        Teuchos::View,*H,order,order);
    Teuchos::SerialDenseMatrix<LO,SCALAR> HtransH(order,order);
    int err = HtransH.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,1.0,*H,*H,0.0);
    CHECK( err == 0 );
    Teuchos::ArrayRCP<SCALAR> alpha_r(order), alpha_i(order), beta_r(order);
    LO worksize = 8*order;
    Teuchos::ArrayRCP<SCALAR> work(worksize);
    int info;
    Teuchos::LAPACK<LO,SCALAR> lapack;
    lapack.GGEV('N','N',order,HtransH.values(),HtransH.stride(),
        H_square.values(),H_square.stride(), alpha_r.get(),alpha_i.get(),
        beta_r.get(),0,order,0,order,work.get(),worksize,&info);
    CHECK( info == 0 );
    Teuchos::ArrayRCP<Anasazi::Value<SCALAR> > roots(order);
    for( int i=0; i<order; ++i )
    {
        roots[i].realpart = alpha_r[i] / beta_r[i];
        roots[i].imagpart = alpha_i[i] / beta_r[i];
    }
    if( b_verbosity >= HIGH )
    {
        std::cout << "GMRES roots" << std::endl;
        for( int i=0; i<order; ++i )
        {
            std::cout << i << " " << roots[i].realpart << " "
                << roots[i].imagpart << std::endl;
        }
    }
    return roots;
}
Teuchos :: ArrayRCP < const Anasazi :: Value < SCALAR > > GmresPolynomial :: computeGmresRoots ( Teuchos :: RCP < Teuchos :: SerialDenseMatrix < LO , SCALAR > > H ) const { LO order = H -> numCols ( ) ; REQUIRE ( H -> numRows ( ) == order + 1 ) ; Teuchos :: SerialDenseMatrix < LO , SCALAR > H_square ( Teuchos :: View , * H , order , order ) ; Teuchos :: SerialDenseMatrix < LO , SCALAR > HtransH ( order , order ) ; int err = HtransH . multiply ( Teuchos :: TRANS , Teuchos :: NO_TRANS , 1.0 , * H , * H , 0.0 ) ; CHECK ( err == 0 ) ; Teuchos :: ArrayRCP < SCALAR > alpha_r ( order ) , alpha_i ( order ) , beta_r ( order ) ; LO worksize = 8 * order ; Teuchos :: ArrayRCP < SCALAR > work ( worksize ) ; int info ; Teuchos :: LAPACK < LO , SCALAR > lapack ; lapack . GGEV ( ' ' , ' ' , order , HtransH . values ( ) , HtransH . stride ( ) , H_square . values ( ) , H_square . stride ( ) , alpha_r . get ( ) , alpha_i . get ( ) , beta_r . get ( ) , 0 , order , 0 , order , work . get ( ) , worksize , & info ) ; CHECK ( info == 0 ) ; Teuchos :: ArrayRCP < Anasazi :: Value < SCALAR > > roots ( order ) ; for ( int i = 0 ; i < order ; ++ i ) { roots [ i ] . realpart = alpha_r [ i ] / beta_r [ i ] ; roots [ i ] . imagpart = alpha_i [ i ] / beta_r [ i ] ; } if ( b_verbosity >= HIGH ) { std :: cout << " " << std :: endl ; for ( int i = 0 ; i < order ; ++ i ) { std :: cout << i << " " << roots [ i ] . realpart << " " << roots [ i ] . imagpart << std :: endl ; } } return roots ; }
Compute GMRES polynomial roots from Arnoldi matrix
Compute GMRES polynomial roots from Arnoldi matrix
// Compute eigenvalues of leading NxN block of H // Matrix to hold (H^T * H)