Method factorizeDenseMatrix

src/IB/DirectMobilitySolver.cpp:602–694

} // factorizeBodyMobilityMatrix

void
DirectMobilitySolver::factorizeDenseMatrix(double* mat_data,
                                           const int mat_size,
                                           const MobilityMatrixInverseType& inv_type,
                                           int* ipiv,
                                           const std::string& mat_name,
                                           const std::string& err_msg)
{
    using MatrixType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
    Eigen::Map<MatrixType> mat_view(mat_data, Eigen::Index(mat_size), Eigen::Index(mat_size));

    if (inv_type != LAPACK_LU)
    {
        // For compatibility with the old code, we copy the lower triangle into
        // the upper triangle, even if they aren't actually equal
        for (int i = 0; i < mat_size; ++i)
        {
            for (int j = i + 1; j < mat_size; ++j)
            {
                mat_view(i, j) = mat_view(j, i);
            }
        }
    }

    if (inv_type == LAPACK_CHOLESKY)
    {
        Eigen::LLT<MatrixType> cholesky_factorization(mat_view);
        mat_view = cholesky_factorization.matrixL();
    }
    else if (inv_type == LAPACK_LU)
    {
        Eigen::PartialPivLU<MatrixType> plu_factorization(mat_view);
        mat_view = plu_factorization.matrixLU();
        const auto& indices = plu_factorization.permutationP().indices();
        TBOX_ASSERT(indices.rows() * indices.cols() == mat_size);
        std::copy_n(indices.data(), mat_size, ipiv);
    }
    else if (inv_type == LAPACK_SVD)
    {
        // Use the symmetric eigenvalue decomposition as a stand-in for the SVD.
        // In particular, since A = V D V^T, where V is the matrix of eigenvectors
        // and D is the matrix of eigenvalues, we can factorize A as
        //
        //   A = V sqrt(D) sqrt(D) V^T
        //     = V sqrt(D) (V sqrt(D))^T
        //
        // and instead store A <- V sqrt(D)
        Eigen::SelfAdjointEigenSolver<MatrixType> eigensolver(mat_view);
        Eigen::Matrix<double, Eigen::Dynamic, 1> eigenvalues = eigensolver.eigenvalues();
        const MatrixType eigenvectors = eigensolver.eigenvectors();
        // Make negative eigenvalues to be equal to min eigen value from
        // input option
        int counter = 0, counter_zero = 0;
        for (int i = 0; i < mat_size; ++i)
        {
            if (eigenvalues[i] < d_svd_eps)
            {
                eigenvalues[i] = d_svd_replace_value;