Skip to content
Snippets Groups Projects
EigenMatrix.C 28 KiB
Newer Older
  • Learn to ignore specific revisions
  • 
    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;
    }
    
    
    // ************************************************************************* //