This is my function for a forwards-backwards substition solver. Unfortunately, my cusparseSpSV_solve function calls are setting the first value in the resulting vector to INF. Curiously, this only happens with matrices of dimensions 641x641 or greater. Anyone know why only the first value is being affected, and why this only happens with a certain matrix size?
void triangular_solve(cusparseHandle_t cusparseHandle, cusparseConstSpMatDescr_t matL, cusparseDnVecDescr_t vecZ, cusparseDnVecDescr_t vecRes, cusparseDnVecDescr_t vecTemp_a) {
size_t bufferSize = 0; // Declare variables for buffer size and buffer allocation
void* buffer = nullptr;
double alpha = 1.0; // Initialize alpha for scaling (usually set to 1.0)
// Forward solve: L * a = r
// Step 1: Create SpSV descriptor
cusparseSpSVDescr_t spSVDescr = NULL;
cusparseSpSV_createDescr(&spSVDescr);
// Step 2: Query buffer size for forward solve (L * a = r)
cusparseSpSV_bufferSize(cusparseHandle,
CUSPARSE_OPERATION_NON_TRANSPOSE, // opA
&alpha,
matL, // matrix A
vecRes, // vector r
vecTemp_a, // vector a (output)
CUDA_R_64F, // Use the appropriate data type (double precision in this case)
CUSPARSE_SPSV_ALG_DEFAULT,
spSVDescr,
&bufferSize);
// Step 3: Allocate memory for buffer
cudaMalloc(&buffer, bufferSize);
// Step 4: Perform analysis for forward solve
cusparseSpSV_analysis(cusparseHandle,
CUSPARSE_OPERATION_NON_TRANSPOSE,
&alpha,
matL,
vecRes,
vecTemp_a,
CUDA_R_64F,
CUSPARSE_SPSV_ALG_DEFAULT,
spSVDescr,
buffer);
// Step 5: Solve L * a = r (forward solve)
cusparseSpSV_solve(cusparseHandle,
CUSPARSE_OPERATION_NON_TRANSPOSE,
&alpha,
matL,
vecRes,
vecTemp_a,
CUDA_R_64F,
CUSPARSE_SPSV_ALG_DEFAULT,
spSVDescr);
// Backward solve: L^T * z = a
// Step 6: Query buffer size for backward solve (L^T * z = a)
cusparseSpSV_bufferSize(cusparseHandle,
CUSPARSE_OPERATION_TRANSPOSE,
&alpha,
matL, // matrix A
vecTemp_a, // vector a (input)
vecZ, // vector z (output)
CUDA_R_64F,
CUSPARSE_SPSV_ALG_DEFAULT,
spSVDescr,
&bufferSize);
// Step 7: Reallocate buffer if necessary (if backward solve needs more memory)
cudaFree(buffer);
cudaMalloc(&buffer, bufferSize);
// Step 8: Perform analysis for backward solve
cusparseSpSV_analysis(cusparseHandle,
CUSPARSE_OPERATION_TRANSPOSE,
&alpha,
matL,
vecTemp_a,
vecZ,
CUDA_R_64F,
CUSPARSE_SPSV_ALG_DEFAULT,
spSVDescr,
buffer);
// Step 9: Solve L^T * z = a (backward solve)
cusparseSpSV_solve(cusparseHandle,
CUSPARSE_OPERATION_TRANSPOSE,
&alpha,
matL,
vecTemp_a,
vecZ,
CUDA_R_64F,
CUSPARSE_SPSV_ALG_DEFAULT,
spSVDescr);
// Clean up: free buffer and destroy the descriptor
cudaFree(buffer);
cusparseSpSV_destroyDescr(spSVDescr);
}