wmlpolynomialrootsdegn.inl

Wild Math Library数值计算库

            fValue0 = fValue1;
            pkF0 = pkF1;
        }
    }
    else
    {
        pkF0 = kSturm[0];
        fValue0 = (*pkF0)(fT0);
        if ( Math<Real>::FAbs(fValue0) < m_fEpsilon )
            fValue0 = (Real)0.0;

        for (i = 1; i < (int)kSturm.size(); i++)
        {
            pkF1 = kSturm[i];
            fValue1 = (*pkF1)(fT0);
            if ( Math<Real>::FAbs(fValue1) < m_fEpsilon )
                fValue1 = (Real)0.0;

            if ( fValue0*fValue1 < (Real)0.0 || fValue0 == (Real)0.0 )
                iSignChanges0++;

            fValue0 = fValue1;
            pkF0 = pkF1;
        }
    }

    // count the sign changes at t1
    int iSignChanges1 = 0;
    if ( fT1 == Math<Real>::MAX_REAL )
    {
        pkF0 = kSturm[0];
        fValue0 = (*pkF0)[pkF0->GetDegree()];
        if ( Math<Real>::FAbs(fValue0) < m_fEpsilon )
            fValue0 = (Real)0.0;

        for (i = 1; i < (int)kSturm.size(); i++)
        {
            pkF1 = kSturm[i];
            fValue1 = (*pkF1)[pkF1->GetDegree()];
            if ( Math<Real>::FAbs(fValue1) < m_fEpsilon )
                fValue1 = (Real)0.0;

            if ( fValue0*fValue1 < (Real)0.0 || fValue0 == (Real)0.0 )
                iSignChanges1++;

            fValue0 = fValue1;
            pkF0 = pkF1;
        }
    }
    else
    {
        pkF0 = kSturm[0];
        fValue0 = (*pkF0)(fT1);
        if ( Math<Real>::FAbs(fValue0) < m_fEpsilon )
            fValue0 = (Real)0.0;

        for (i = 1; i < (int)kSturm.size(); i++)
        {
            pkF1 = kSturm[i];
            fValue1 = (*pkF1)(fT1);
            if ( Math<Real>::FAbs(fValue1) < m_fEpsilon )
                fValue1 = (Real)0.0;

            if ( fValue0*fValue1 < (Real)0.0 || fValue0 == (Real)0.0 )
                iSignChanges1++;

            fValue0 = fValue1;
            pkF0 = pkF1;
        }
    }

    // clean up
    for (i = 0; i < (int)kSturm.size(); i++)
        delete kSturm[i];

    if ( iSignChanges0 >= iSignChanges1 )
        return iSignChanges0 - iSignChanges1;

    // theoretically we should not get here
    assert( false );
    return 0;
}
//----------------------------------------------------------------------------
template <class Real>
bool PolynomialRoots<Real>::Bisection (const Polynomial1<Real>& rkPoly,
    Real fXMin, Real fXMax, int iDigits, Real& rfRoot)
{
    Real fP0 = rkPoly(fXMin);
    if ( Math<Real>::FAbs(fP0) <= Math<Real>::ZERO_TOLERANCE )
    {
        rfRoot = fXMin;
        return true;
    }
    Real fP1 = rkPoly(fXMax);
    if ( Math<Real>::FAbs(fP1) <= Math<Real>::ZERO_TOLERANCE )
    {
        rfRoot = fXMax;
        return true;
    }
    if ( fP0*fP1 > (Real)0.0 )
        return false;

    // determine number of iterations to get 'digits' accuracy.
    Real fTmp0 = Math<Real>::Log(fXMax-fXMin);
    Real fTmp1 = ((Real)iDigits)*LOG10;
    Real fArg = (fTmp0 + fTmp1)*INVLOG2;
    int iMaxIter = (int)(fArg + 0.5f);

    for (int i = 0; i < iMaxIter; i++)
    {
        rfRoot = ((Real)0.5)*(fXMin + fXMax);
        Real fP = rkPoly(rfRoot);
        Real fProduct = fP*fP0;
        if ( fProduct < -Math<Real>::ZERO_TOLERANCE )
        {
            fXMax = rfRoot;
            fP1 = fP;
        }
        else if ( fProduct > Math<Real>::ZERO_TOLERANCE )
        {
            fXMin = rfRoot;
            fP0 = fP;
        }
        else
        {
            break;
        }
    }

    return true;
}
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// FindE support
//----------------------------------------------------------------------------
template <class Real>
void PolynomialRoots<Real>::GetHouseholderVector (int iSize,
    const Vector3<Real>& rkU, Vector3<Real>& rkV)
{
    // Householder vector V:
    // Given a vector U, compute a vector V such that V[0] = 1 and
    // (I-2*V*V^T/|V|^2)*U is zero in all but the first component.  The
    // matrix P = I-2*V*V^T/|V|^2 is a Householder transformation, a
    // reflection matrix.

    Real fLength = rkU[0]*rkU[0];
    int i;
    for (i = 1; i < iSize; i++)
        fLength += rkU[i]*rkU[i];
    fLength = Math<Real>::Sqrt(fLength);

    if ( fLength > m_fEpsilon )
    {
        Real fInv = ((Real)1.0)/(rkU[0]+Math<Real>::Sign(rkU[0])*fLength);
        rkV[0] = (Real)1.0;
        for (i = 1; i < iSize; i++)
            rkV[i] = fInv*rkU[i];
    }
    else
    {
        // U is the zero vector, any vector will do
        rkV[0] = (Real)1.0;
        for (i = 1; i < iSize; i++)
            rkV[i] = (Real)0.0;
    }
}
//----------------------------------------------------------------------------
template <class Real>
void PolynomialRoots<Real>::PremultiplyHouseholder (GMatrix<Real>& rkMat,
    GVector<Real>& rkW, int iRMin, int iRMax, int iCMin, int iCMax,
    int iVSize, const Vector3<Real>& rkV)
{
    // Householder premultiplication:
    // Given a matrix A and an m-by-1 vector V with V[0] = 1, let S be the
    // submatrix of A with m rows rmin <= r <= m+rmin-1 and columns
    // cmin <= c <= cmax.  Overwrite S with P*S where P = I-2*V*V^T/|V|^2.

    int iSubRows = iRMax - iRMin + 1, iSubCols = iCMax - iCMin + 1;
    int iRow, iCol;

    Real fSqrLen = rkV[0]*rkV[0];
    for (int i = 1; i < iVSize; i++)
        fSqrLen += rkV[i]*rkV[i];

    Real fBeta = -((Real)2.0)/fSqrLen;
    for (iCol = 0; iCol < iSubCols; iCol++)
    {
        rkW[iCol] = (Real)0.0;
        for (iRow = 0; iRow < iSubRows; iRow++)
            rkW[iCol] += rkV[iRow]*rkMat[iRMin+iRow][iCMin+iCol];
        rkW[iCol] *= fBeta;
    }

    for (iRow = 0; iRow < iSubRows; iRow++)
    {
        for (iCol = 0; iCol < iSubCols; iCol++)
            rkMat[iRMin+iRow][iCMin+iCol] += rkV[iRow]*rkW[iCol];
    }
}
//----------------------------------------------------------------------------
template <class Real>
void PolynomialRoots<Real>::PostmultiplyHouseholder (GMatrix<Real>& rkMat,
    GVector<Real>& rkW, int iRMin, int iRMax, int iCMin, int iCMax,
    int iVSize, const Vector3<Real>& rkV)
{
    // Householder postmultiplication:
    // Given a matrix A and an n-by-1 vector V with V[0] = 1, let S be the
    // submatrix of A with n columns cmin <= c <= n+cmin-1 and rows
    // rmin <= r <= rmax.  Overwrite S with S*P where P = I-2*V*V^T/|V|^2.

    int iSubRows = iRMax - iRMin + 1, iSubCols = iCMax - iCMin + 1;
    int iRow, iCol;

    Real fSqrLen = rkV[0]*rkV[0];
    for (int i = 1; i < iVSize; i++)
        fSqrLen += rkV[i]*rkV[i];

    Real fBeta = -((Real)2.0)/fSqrLen;
    for (iRow = 0; iRow < iSubRows; iRow++)
    {
        rkW[iRow] = (Real)0.0;
        for (iCol = 0; iCol < iSubCols; iCol++)
            rkW[iRow] += rkMat[iRMin+iRow][iCMin+iCol]*rkV[iCol];
        rkW[iRow] *= fBeta;
    }

    for (iRow = 0; iRow < iSubRows; iRow++)
    {
        for (iCol = 0; iCol < iSubCols; iCol++)
            rkMat[iRMin+iRow][iCMin+iCol] += rkW[iRow]*rkV[iCol];
    }
}
//----------------------------------------------------------------------------
template <class Real>
void PolynomialRoots<Real>::FrancisQRStep (GMatrix<Real>& rkH,
    GVector<Real>& rkW)
{
    // Given an n-by-n unreduced upper Hessenberg matrix H whose trailing
    // 2-by-2 principal submatrix has eigenvalues a1 and a2, overwrite H
    // with Z^T*H*Z where Z = P(1)*...*P(n-2) is a product of Householder
    // matrices and Z^T*(H-a1*I)*(H-a2*I) is upper triangular.

    // assert:  H is unreduced upper Hessenberg and n >= 3

    // compute first column of (H-a1*I)*(H-a2*I)
    int iN = rkH.GetRows();
    Real fTrace = rkH[iN-2][iN-2] + rkH[iN-1][iN-1];
    Real fDet = rkH[iN-2][iN-2]*rkH[iN-1][iN-1] -
        rkH[iN-2][iN-1]*rkH[iN-1][iN-2];
    Vector3<Real> kU;
    kU[0] = rkH[0][0]*rkH[1][1]+rkH[0][1]*rkH[1][0]-fTrace*rkH[0][0]+fDet;
    kU[1] = rkH[1][0]*(rkH[0][0]+rkH[1][1]-fTrace);
    kU[2] = rkH[1][0]*rkH[2][1];

    // overwrite H with P(0)*H*P(0)^T
    Vector3<Real> kV;
    GetHouseholderVector(3,kU,kV);
    PremultiplyHouseholder(rkH,rkW,0,2,0,iN-1,3,kV);
    PostmultiplyHouseholder(rkH,rkW,0,iN-1,0,2,3,kV);

    for (int i = 1; i <= iN-3; i++)
    {
        // overwrite H with P(i)*H*P(i)^T
        kU[0] = rkH[i  ][i-1];
        kU[1] = rkH[i+1][i-1];
        kU[2] = rkH[i+2][i-1];
        GetHouseholderVector(3,kU,kV);

        // The column range does not need to be 0 to iN-1 because of the
        // pattern of zeros that occur in matrix H.
        PremultiplyHouseholder(rkH,rkW,i,i+2,i-1,iN-1,3,kV);

        // The row range does not need to be 0 to iN-1 because of the pattern
        // of zeros that occur in matrix H.
        int iRMax = i+3;
        if ( iRMax >= iN )
            iRMax = iN-1;
        PostmultiplyHouseholder(rkH,rkW,0,iRMax,i,i+2,3,kV);
    }

    // overwrite H with P(n-2)*H*P(n-2)^T
    kU[0] = rkH[iN-2][iN-3];
    kU[1] = rkH[iN-1][iN-3];
    GetHouseholderVector(2,kU,kV);

    // The column range does not need to be 0 to iN-1 because of the pattern
    // of zeros that occur in matrix H.
    PremultiplyHouseholder(rkH,rkW,iN-2,iN-1,iN-3,iN-1,2,kV);

    PostmultiplyHouseholder(rkH,rkW,0,iN-1,iN-2,iN-1,2,kV);
}
//----------------------------------------------------------------------------