packages/Alea/mc_solvers/GmresPolynomial.cc
//---------------------------------------------------------------------------//
// Compute FOM polynomial roots from Arnoldi matrix
//---------------------------------------------------------------------------//
Compute FOM polynomial roots from Arnoldi matrix
Compute FOM polynomial roots from Arnoldi matrix
Teuchos::ArrayRCP<const Anasazi::Value<SCALAR> >
GmresPolynomial::computeFomRoots(
    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::ArrayRCP<SCALAR> real_parts(order), imag_parts(order);
    LO worksize = 6*order;
    Teuchos::ArrayRCP<SCALAR> work(worksize);
    int info;
    Teuchos::LAPACK<LO,SCALAR> lapack;
    lapack.HSEQR('E','N',order,1,order,H_square.values(),H_square.stride(),
        real_parts.get(),imag_parts.get(),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 = real_parts[i];
        roots[i].imagpart = imag_parts[i];
    }
    if( b_verbosity >= HIGH )
    {
        std::cout << "FOM roots" << std::endl;
        for( int i=0; i<order; ++i )
        {
            std::cout << i << " " << real_parts[i] << " "
                << imag_parts[i] << std::endl;
        }
    }
    return roots;
}
Teuchos :: ArrayRCP < const Anasazi :: Value < SCALAR > > GmresPolynomial :: computeFomRoots ( 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 :: ArrayRCP < SCALAR > real_parts ( order ) , imag_parts ( order ) ; LO worksize = 6 * order ; Teuchos :: ArrayRCP < SCALAR > work ( worksize ) ; int info ; Teuchos :: LAPACK < LO , SCALAR > lapack ; lapack . HSEQR ( ' ' , ' ' , order , 1 , order , H_square . values ( ) , H_square . stride ( ) , real_parts . get ( ) , imag_parts . get ( ) , 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 = real_parts [ i ] ; roots [ i ] . imagpart = imag_parts [ i ] ; } if ( b_verbosity >= HIGH ) { std :: cout << " " << std :: endl ; for ( int i = 0 ; i < order ; ++ i ) { std :: cout << i << " " << real_parts [ i ] << " " << imag_parts [ i ] << std :: endl ; } } return roots ; }
Compute FOM polynomial roots from Arnoldi matrix
Compute FOM polynomial roots from Arnoldi matrix
// Compute eigenvalues of leading NxN block of H