EigenMatrix.C

template<class cmptType>
Foam::EigenMatrix<cmptType>::EigenMatrix(const SquareMatrix<cmptType>& A)
:
    n_(A.n()),
    EValsRe_(n_, Zero),
    EValsIm_(n_, Zero),
    EVecs_(n_, Zero),
    H_()
{
    if (n_ <= 0)
    {
        FatalErrorInFunction
            << "Input matrix has zero size."
            << abort(FatalError);
    }

    if (A.symmetric())
    {
        EVecs_ = A;
        tridiagonaliseSymmMatrix();
        symmTridiagonalQL();
    }
    else
    {
        H_ = A;
        Hessenberg();
        realSchur();
    }
}


template<class cmptType>
Foam::EigenMatrix<cmptType>::EigenMatrix
(
    const SquareMatrix<cmptType>& A,
    bool symmetric
)
:
    n_(A.n()),
    EValsRe_(n_, Zero),
    EValsIm_(n_, Zero),
    EVecs_(n_, Zero),
    H_()
{
    if (n_ <= 0)
    {
        FatalErrorInFunction
            << "Input matrix has zero size."
            << abort(FatalError);
    }

    if (symmetric)
    {
        EVecs_ = A;
        tridiagonaliseSymmMatrix();
        symmTridiagonalQL();
    }
    else
    {
        H_ = A;
        Hessenberg();
        realSchur();
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class cmptType>
const Foam::SquareMatrix<Foam::complex>
Foam::EigenMatrix<cmptType>::complexEVecs() const
{
    SquareMatrix<complex> EVecs(EVecs_.n());

    std::transform
    (
        EVecs_.cbegin(),
        EVecs_.cend(),
        EVecs.begin(),
        [&](const cmptType& x){ return complex(x); }
    );

    for (label i = 0; i < EValsIm_.size(); ++i)
    {
        if (mag(EValsIm_[i]) > VSMALL)
        {
            for (label j = 0; j < EVecs.m(); ++j)
            {
                EVecs(j, i) = complex(EVecs_(j, i), EVecs_(j, i+1));
                EVecs(j, i+1) = EVecs(j, i).conjugate();
            }
            ++i;
        }
    }

    return EVecs;
}


// ************************************************************************* //