CHOLMOD稀疏-稠密乘法异常：生成b=A*xe结果为零向量求助

Let's break down why your code is producing a zero vector for b and fix it step by step:

The Core Issue: Mismatched CHOLMOD Function Versions #

You're mixing CHOLMOD's long integer (_l suffix) functions (like cholmod_l_start, cholmod_l_read_sparse) with the regular integer version of cholmod_sdmult. CHOLMOD uses separate function sets for int and long index types—using the wrong one causes silent type mismatches under the hood, leading to incorrect calculations (like your zero vector result).

Secondary Issue: Unnecessary Memory Leaks #

Your code reassigns b three times without freeing the previous allocations, which wastes memory. While this isn't causing the zero vector problem, it's a bad practice to fix.

Corrected Code #

Here's the fixed version with explanations embedded:

#include "cholmod.h"
#include <iostream>
#pragma warning(disable:4996)

int main() {
    const char* matFile = "C:/www/CholeskyMatGPU/ex15.mtx";
    const char* bfile = "C:/www/CholeskyMatGPU/b.mtx";
    FILE* bp = fopen(bfile, "w");
    FILE* fp = fopen(matFile, "r");
    
    cholmod_sparse *A;
    cholmod_dense *x, *xe, *b;
    cholmod_factor *L;
    cholmod_common* c = (cholmod_common*)malloc(sizeof(cholmod_common));
    
    // Initialize CHOLMOD's long integer context
    cholmod_l_start(c);
    c->useGPU = 0;
    c->supernodal = CHOLMOD_SUPERNODAL;
    
    // Read sparse matrix (long version)
    A = cholmod_l_read_sparse(fp, c);
    if (!A) {
        std::cerr << "Failed to read matrix A" << std::endl;
        fclose(fp);
        cholmod_l_finish(c);
        free(c);
        return 1;
    }
    cholmod_l_print_sparse(A, "A", c);
    fclose(fp);
    
    double alpha[2] = {1., 0.};  // Real scalar for multiplication
    double beta[2] = {0., 0.};   // Zero out existing b before multiplication
    
    // Create all-ones vector xe (matches A's type and dimensions)
    xe = cholmod_l_ones(A->nrow, 1, A->xtype, c);
    // Allocate b with correct dimensions and type (no redundant reassignments)
    b = cholmod_l_allocate_dense(A->nrow, 1, A->nrow, A->xtype, c);
    if (!b) {
        std::cerr << "Failed to allocate vector b" << std::endl;
        cholmod_l_free_sparse(&A, c);
        cholmod_l_free_dense(&xe, c);
        cholmod_l_finish(c);
        free(c);
        return 1;
    }
    
    // Use THE LONG VERSION of sdmult to compute b = alpha*A*xe + beta*b
    cholmod_l_sdmult(A, 0, alpha, beta, xe, b, c);
    
    // Write b to file
    if (bp) {
        cholmod_l_write_dense(bp, b, "", c);
        fclose(bp);
    }
    
    // Proceed with factorization and solve
    L = cholmod_l_analyze(A, c);
    if (!L) {
        std::cerr << "Failed to analyze matrix A" << std::endl;
        // Cleanup resources
        cholmod_l_free_sparse(&A, c);
        cholmod_l_free_dense(&xe, c);
        cholmod_l_free_dense(&b, c);
        cholmod_l_finish(c);
        free(c);
        return 1;
    }
    
    cholmod_l_factorize(A, L, c);
    x = cholmod_l_solve(CHOLMOD_A, L, b, c);
    
    // Optional: Verify solution by checking x vs xe (should be nearly identical)
    cholmod_dense* diff = cholmod_l_copy_dense(x, c);
    double neg_alpha[2] = {-1.0, 0.0};
    cholmod_l_dense_add(neg_alpha, xe, diff, c);
    double error_norm = cholmod_l_norm_dense(diff, 2, c);
    std::cout << "Norm of x - xe (solution error): " << error_norm << std::endl;
    
    // Cleanup all CHOLMOD resources
    cholmod_l_free_factor(&L, c);
    cholmod_l_free_dense(&x, c);
    cholmod_l_free_dense(&diff, c);
    cholmod_l_free_dense(&b, c);
    cholmod_l_free_dense(&xe, c);
    cholmod_l_free_sparse(&A, c);
    
    cholmod_l_finish(c);
    free(c);
    
    return 0;
}

Key Fixes Explained #

1. Used cholmod_l_sdmult instead of cholmod_sdmult: This matches the long integer context initialized with cholmod_l_start, ensuring index types align correctly during multiplication.
2. Removed redundant b assignments: We allocate b once with the correct dimensions and type, avoiding memory leaks.
3. Added error checking: We validate return values of CHOLMOD functions to catch allocation/reading failures early.
4. Added solution verification: The code calculates the norm of x - xe to confirm the solved x is close to the original all-ones vector (as expected).

Additional Notes #

• Always pair cholmod_l_* functions with each other—never mix regular and long versions in the same context.
• Don't forget to free all CHOLMOD-allocated memory with cholmod_l_free_* functions and call cholmod_l_finish before exiting to clean up the library's internal state.

内容的提问来源于stack exchange，提问作者David