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::ArrayRCP<SCALAR> GmresPolynomial::computeCoeffsFromRoots(
    Teuchos::ArrayRCP<const Anasazi::Value<SCALAR> > roots) const
{
    Teuchos::ArrayRCP<SCALAR> coeffs(roots.size());
    Teuchos::ArrayRCP<SCALAR> c(1,1.0);
    int i=0;
    while( i<roots.size() )
    {
        SCALAR a = roots[i].realpart;
        SCALAR b = roots[i].imagpart;
        if( b == 0.0 )
        {
            Teuchos::ArrayRCP<SCALAR> c1(c.size()+1,0.0);
            Teuchos::ArrayRCP<SCALAR> c2(c.size()+1,0.0);
            for( int j=0; j<i+1; ++j )
            {
                c1[j+1] = -c[j];
                c2[j] = a * c[j];
            }
            c.resize(c.size()+1,0.0);
            for( int j=0; j<c.size(); ++j )
            {
                c[j] = c1[j] + c2[j];
            }
            i++;
        }
        else
        {
            Teuchos::ArrayRCP<SCALAR> c1(c.size()+2,0.0);
            Teuchos::ArrayRCP<SCALAR> c2(c.size()+2,0.0);
            Teuchos::ArrayRCP<SCALAR> c3(c.size()+2,0.0);
            for( int j=0; j<i+1; ++j )
            {
                c1[j+2] = c[j];
                c2[j+1] = -2.0*a*c[j];
                c3[j]   = (a*a + b*b) * c[j];
            }
            c.resize(c.size()+2,0.0);
            for( int j=0; j<c.size(); ++j )
            {
                c[j] = c1[j] + c2[j] + c3[j];
            }
            i+=2;
        }
    }
    CHECK( c.size() == coeffs.size()+1 );
    for( int i=0; i<coeffs.size(); ++i )
    {
        coeffs[i] = -c[i+1] / c[0];
    }
    return coeffs;
}
Teuchos :: ArrayRCP < SCALAR > GmresPolynomial :: computeCoeffsFromRoots ( Teuchos :: ArrayRCP < const Anasazi :: Value < SCALAR > > roots ) const { Teuchos :: ArrayRCP < SCALAR > coeffs ( roots . size ( ) ) ; Teuchos :: ArrayRCP < SCALAR > c ( 1 , 1.0 ) ; int i = 0 ; while ( i < roots . size ( ) ) { SCALAR a = roots [ i ] . realpart ; SCALAR b = roots [ i ] . imagpart ; if ( b == 0.0 ) { Teuchos :: ArrayRCP < SCALAR > c1 ( c . size ( ) + 1 , 0.0 ) ; Teuchos :: ArrayRCP < SCALAR > c2 ( c . size ( ) + 1 , 0.0 ) ; for ( int j = 0 ; j < i + 1 ; ++ j ) { c1 [ j + 1 ] = - c [ j ] ; c2 [ j ] = a * c [ j ] ; } c . resize ( c . size ( ) + 1 , 0.0 ) ; for ( int j = 0 ; j < c . size ( ) ; ++ j ) { c [ j ] = c1 [ j ] + c2 [ j ] ; } i ++ ; } else { Teuchos :: ArrayRCP < SCALAR > c1 ( c . size ( ) + 2 , 0.0 ) ; Teuchos :: ArrayRCP < SCALAR > c2 ( c . size ( ) + 2 , 0.0 ) ; Teuchos :: ArrayRCP < SCALAR > c3 ( c . size ( ) + 2 , 0.0 ) ; for ( int j = 0 ; j < i + 1 ; ++ j ) { c1 [ j + 2 ] = c [ j ] ; c2 [ j + 1 ] = -2.0 * a * c [ j ] ; c3 [ j ] = ( a * a + b * b ) * c [ j ] ; } c . resize ( c . size ( ) + 2 , 0.0 ) ; for ( int j = 0 ; j < c . size ( ) ; ++ j ) { c [ j ] = c1 [ j ] + c2 [ j ] + c3 [ j ] ; } i += 2 ; } } CHECK ( c . size ( ) == coeffs . size ( ) + 1 ) ; for ( int i = 0 ; i < coeffs . size ( ) ; ++ i ) { coeffs [ i ] = - c [ i + 1 ] / c [ 0 ] ; } return coeffs ; }
Convert roots of polynomial to coefficients
Convert roots of polynomial to coefficients
// Expand polynomial from factored form p(x) = (x-l1)(x-l2)... //  into monomial basis p(x) = c0 + c1*x + c2*x^2 + ... // Real roots: p_{n+1}(x) = (a-x)*p_{n}(x) // Conjugate pair: p_{n+2}(x) = (x^2 - 2*a*x + (a^2+b^2))*p_{n}(x)