Is it possible to solve the system using the Cholesky decomposition if the matrix is not positive definite?

Question

I cloned CHOLMOD project (DrTimothyAldenDavis / SuiteSparse) on Ubuntu. Then I find simple program for solving system. In my equation Ax = b, matrix A is symmetric but not positive definite, I hoped in this case using form LDLt (instead LLt) will help me. But program returned strange result X = { -nan, -nan, ...., -nan, inf}, and message " CHOLMOD warning: not positive definite ".

#include <stdio.h>
#include <stdint.h>
#include "../Include/cholmod.h"


int main()
{
    double dot=0;
    int nlab =0;
    int check=1;
    int neq =10;
    double K[10][10]=
            {{0.111, 0.222, 0.333, 0.444, 0.555, 0.666, 0.777, 0.888, 0.999, 0.001},
             {0.222, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.333, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.444, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.555, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.666, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.777, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.888, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001},
             {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001}};


    double inputRHS[10] {4.996000, 0.223000, 0.334000, 0.445000, 0.556000, 0.667000,  0.778000, 0.889000, 1.000000, 0.010000};

    cholmod_dense *Y, *X, *r;

    cholmod_factor *L ;

    double one [2] = {1,0}, m1 [2] = {-1,0};

    cholmod_common c;

    cholmod_start (&c);

    // GET THE SIZE OF THE MATRIX
    int nrow=10;
    int ncol=10;
    // CREATE A NULL CHOLMOD MATRIX and RHS
    cholmod_dense* A = cholmod_eye(nrow,ncol,CHOLMOD_REAL,&c);       
    cholmod_dense* B = cholmod_ones(nrow,1,CHOLMOD_REAL,&c);       

    // CONVERT MATRIX AND RHS INTO CHOLMOD MATRIX
    for(int loopA=0; loopA<nrow; loopA++){
        ((double*)B->x)[loopA] = inputRHS[loopA];            
        for(int loopB=0; loopB<ncol; loopB++){
            ((double*)A->x)[loopA + loopB*A->d] = K[loopA][loopB];
        }
    }


    c.print = 5; // for printing
    // PRINT DENSE MATRIX
    cholmod_print_dense(A,"A",&c);

    // PRINT DENSE VECTOR
    cholmod_print_dense(B,"B",&c);

    // CONVERT CHOLMOD MATRIX TO SPARSE MATRIX
    cholmod_sparse* spA = cholmod_dense_to_sparse(A,1,&c);
    spA->stype = 1; // A matrix is symetric positive definite SPD
    // PRINT SPARSE MATRIX
    cholmod_print_sparse(spA,"spA",&c);

    // X = A\B
    // AX=B
    // c.final_ll = 1; // LL'
    L = cholmod_analyze (spA, &c) ;
    printf("L->is_ll %d, L->is_super %d \n", L->is_ll, L->is_super);
    cholmod_factorize (spA, L, &c) ;
    Y = cholmod_solve (CHOLMOD_L, L, B, &c) ;
    X = cholmod_solve (CHOLMOD_DLt, L, Y, &c) ;
    /* Also I try direct solving like: X = cholmod_solve (CHOLMOD_LDLt, L, B, &c), but get the same result */ 

    r = cholmod_copy_dense (B, &c) ;
    cholmod_sdmult (spA, 0, m1, one, X, r, &c) ;
    // checking norm
    printf ("norm(b-Ax) %8.1e\n",cholmod_norm_dense (r, 0, &c));
    cholmod_print_dense(X,"x",&c);

    cholmod_free_factor (&L, &c) ;
    cholmod_free_sparse (&spA, &c) ;
    cholmod_free_dense (&r, &c) ;
    cholmod_free_dense (&X, &c) ;
    cholmod_free_dense (&A, &c) ;
    cholmod_free_dense (&B, &c) ;
    cholmod_finish (&c) ;


    return 0;
}

I try to solve system directly using parameter CHOLMOD_LDLt: X = cholmod_solve (CHOLMOD_LDLt, L, B, &c); then try to solve in two steps, as in the listing above.

I would like to ask 2 question:

1. It is correct to use LDLt factorization to solve system Ax = b, when A is not positive definite?
2. If answer Yes, then what kind of CHOLMOD routines I should to use?

Thanks for any advice!

Answer 1

Short answer - yes.

In my case I had error because, unfortunately, det(A) = 0.0.
So make sure det(A) != 0.0, then set parameter CHOLMOD_A: `X = cholmod_solve (CHOLMOD_A, L, B, &c) ;`