/* Testing CHOLMOD Solve b = Ax + noise, with noise covariance R R is sparse, A is sparse */ #include extern "C" { #include "cholmod.h" } using namespace std; int main (void) { cout<<"Solving a Generalized Least Squares using CHOLMOD \n"<nnz = R_nnz; // put in some values for covariance matrix int *Ti, *Tj; double *Tx; Ti = static_cast(T->i); Tj = static_cast(T->j); Tx = static_cast(T->x); Ti[0] = 0; Tj[0] = 0; Tx[0] = 1.0; Ti[1] = 1; Tj[1] = 1; Tx[1] = 2.0; Ti[2] = 2; Tj[2] = 2; Tx[2] = 3.0; Ti[3] = 3; Tj[3] = 3; Tx[3] = 4.0; Ti[4] = 4; Tj[4] = 4; Tx[4] = 5.0; Ti[5] = 5; Tj[5] = 5; Tx[5] = 6.0; // convert triplet to sparse matrix R = cholmod_triplet_to_sparse(T, R_nnz, &c); // free the triplet as we will reuse T cholmod_free_triplet(&T, &c); // Start Building A size_t A_nrows = 6; size_t A_ncols = 3; size_t A_nnz = 11; T = cholmod_allocate_triplet(A_nrows, A_ncols, A_nnz, 0, CHOLMOD_REAL, &c); T->nnz = A_nnz; Ti = static_cast(T->i); Tj = static_cast(T->j); Tx = static_cast(T->x); Ti[0] = 0; Tj[0] = 0; Tx[0] = 1.0; Ti[1] = 1; Tj[1] = 0; Tx[1] = -1.0; Ti[2] = 1; Tj[2] = 1; Tx[2] = 1.0; Ti[3] = 2; Tj[3] = 1; Tx[3] = -1.0; Ti[4] = 2; Tj[4] = 2; Tx[4] = 1.0; Ti[5] = 3; Tj[5] = 0; Tx[5] = -1.0; Ti[6] = 3; Tj[6] = 2; Tx[6] = 2.0; Ti[7] = 4; Tj[7] = 1; Tx[7] = -3.0; Ti[8] = 4; Tj[8] = 2; Tx[8] = 1.0; Ti[9] = 5; Tj[9] = 0; Tx[9] = 7.0; Ti[10] = 5; Tj[10] = 2; Tx[10] = -1.0; // create sparse matrix from triplet A = cholmod_triplet_to_sparse(T, A_nnz, &c); cholmod_print_sparse (A, "A", &c) ; /* print the matrix */ cholmod_print_sparse (R, "R", &c) ; /* print the matrix */ b = cholmod_ones (A->nrow, 1, A->xtype, &c) ; /* b = ones(n,1) */ cholmod_print_dense(b, "b", &c); // ============================================= // Start solving the generalized least squares // x = ( A' * (R \ A) ) \ ( A' * (R \ b) ) // ============================================= // Step 1, compute R \ A L = cholmod_analyze (R, &c) ; /* analyze R */ cholmod_factorize (R, L, &c) ; /* factorize R */ cholmod_sparse *R_inv_A = cholmod_spsolve (CHOLMOD_A, L, A, &c) ; /* solve Rx=A */ // Step 2 compute R \ b cholmod_dense *R_inv_b = cholmod_solve(CHOLMOD_A, L, b, &c); cholmod_free_factor (&L, &c) ; // Step 3, compute At * R_inv_A cholmod_sparse *At = cholmod_transpose(A, 2, &c); // compute A' cholmod_free_sparse (&A, &c) ; // free up A as no longer needed cholmod_sparse *At_R_inv_A = cholmod_ssmult(At, R_inv_A, 1, 1, 1, &c); // Step 4, compute At * R_inv_b double alpha[2] = {1, 1}; double beta[2] = {1,1}; cholmod_dense *At_R_inv_b = cholmod_zeros(At->nrow, 1, At->xtype, &c); cholmod_sdmult(At, 0, alpha, beta, R_inv_b, At_R_inv_b, &c); // Step 5: At_R_inv_A \ At_R_inv_b L = cholmod_analyze (At_R_inv_A, &c) ; cholmod_factorize(At_R_inv_A, L, &c); x = cholmod_solve(CHOLMOD_A, L, At_R_inv_b, &c); cout<<"Solved linear system, x = " <(x->x); for(unsigned int i =0; i < At->nrow ; i++) { cout<