diff --git a/doc/shared/sunmatrix/SUNMatrix_Sparse.rst b/doc/shared/sunmatrix/SUNMatrix_Sparse.rst index c2458fa33b..d47c43e3bf 100644 --- a/doc/shared/sunmatrix/SUNMatrix_Sparse.rst +++ b/doc/shared/sunmatrix/SUNMatrix_Sparse.rst @@ -375,15 +375,18 @@ following additional user-callable routines: -.. c:function:: int SUNSparseMatrix_Realloc(SUNMatrix A) +.. c:function:: SUNErrCode SUNSparseMatrix_Realloc(SUNMatrix A) This function reallocates internal storage arrays in a sparse matrix so that the resulting sparse matrix has no wasted space (i.e. the space allocated for nonzero entries equals the actual number of - nonzeros, ``indexptrs[NP]``). Returns 0 on success and - 1 on failure (e.g. if the input matrix is not sparse). + nonzeros, ``indexptrs[NP]``). Returns a :c:type:`SUNErrCode`. +.. c:function:: SUNErrCode SUNSparseMatrix_Reallocate(SUNMatrix A) + Function to reallocate internal sparse matrix storage arrays so that the + resulting sparse matrix has storage for a specified number of nonzeros. + Returns a :c:type:`SUNErrCode`. .. c:function:: void SUNSparseMatrix_Print(SUNMatrix A, FILE* outfile) diff --git a/examples/arkode/F2003_custom/test_fsunlinsol_fortran_mod.f90 b/examples/arkode/F2003_custom/test_fsunlinsol_fortran_mod.f90 index 1496141b70..35a4b8344b 100644 --- a/examples/arkode/F2003_custom/test_fsunlinsol_fortran_mod.f90 +++ b/examples/arkode/F2003_custom/test_fsunlinsol_fortran_mod.f90 @@ -143,7 +143,7 @@ program main ! compute B = A*X retval = FSUNMatMatvec(sA, sX, sB) - if (retval /= SUNMAT_SUCCESS) then + if (retval /= SUN_SUCCESS) then print *, 'ERROR: FSUNMatMatvec fail' stop 1 end if diff --git a/examples/arkode/F2003_custom/test_fsunmatrix_fortran_mod.f90 b/examples/arkode/F2003_custom/test_fsunmatrix_fortran_mod.f90 index 424ce2ab36..d4cb1f7aff 100644 --- a/examples/arkode/F2003_custom/test_fsunmatrix_fortran_mod.f90 +++ b/examples/arkode/F2003_custom/test_fsunmatrix_fortran_mod.f90 @@ -231,7 +231,7 @@ program main ! test SUNMatZero retval = FSUNMatZero(sB) if ( (check_matrix_entry(sB, 0.d0, 1.d-14, Nvar, N) /= 0) & - .or. (retval /= SUNMAT_SUCCESS) ) then + .or. (retval /= SUN_SUCCESS) ) then fails = fails + 1 print *, '>>> FAILED test -- FSUNMatZero' else @@ -241,7 +241,7 @@ program main ! test SUNMatCopy retval = FSUNMatCopy(sA, sB) if ( (check_matrix(sA, sB, 1.d-14, Nvar, N) /= 0) & - .or. (retval /= SUNMAT_SUCCESS) ) then + .or. (retval /= SUN_SUCCESS) ) then fails = fails + 1 print *, '>>> FAILED test -- FSUNMatCopy' else @@ -252,7 +252,7 @@ program main retval = FSUNMatCopy(sA, sB) retval = FSUNMatScaleAdd(-1.d0, sB, sB) if ( (check_matrix_entry(sB, 0.d0, 1.d-14, Nvar, N) /= 0) & - .or. (retval /= SUNMAT_SUCCESS) ) then + .or. (retval /= SUN_SUCCESS) ) then fails = fails + 1 print *, '>>> FAILED test -- FSUNMatScaleAdd case 1' else @@ -262,9 +262,9 @@ program main retval = FSUNMatCopy(sA, sD) retval = FSUNMatCopy(sI, sC) retval = FSUNMatScaleAdd(1.d0, sD, sI) - if (retval == SUNMAT_SUCCESS) retval = FSUNMatScaleAdd(1.d0, sC, sA) + if (retval == SUN_SUCCESS) retval = FSUNMatScaleAdd(1.d0, sC, sA) if ( (check_matrix(sD, sC, 1.d-14, Nvar, N) /= 0) & - .or. (retval /= SUNMAT_SUCCESS) ) then + .or. (retval /= SUN_SUCCESS) ) then fails = fails + 1 print *, '>>> FAILED test -- FSUNMatScaleAdd case 2' else @@ -275,7 +275,7 @@ program main retval = FSUNMatCopy(sI, sB) retval = FSUNMatScaleAddI(-1.d0, sB) if ( (check_matrix_entry(sB, 0.d0, 1.d-14, Nvar, N) /= 0) & - .or. (retval /= SUNMAT_SUCCESS) ) then + .or. (retval /= SUN_SUCCESS) ) then fails = fails + 1 print *, '>>> FAILED test -- FSUNMatScaleAddI' else @@ -288,7 +288,7 @@ program main retval = FSUNMatMatvec(sB, sX, sZ) call FN_VLinearSum(3.d0, sY, 1.d0, sX, sW) if ( (check_vector(sW, sZ, 1.d-15*Nvar*Nvar, Nvar, N) /= 0) & - .or. (retval /= SUNMAT_SUCCESS) ) then + .or. (retval /= SUN_SUCCESS) ) then fails = fails + 1 print *, '>>> FAILED test -- FSUNMatMatvec' else diff --git a/src/sunmatrix/band/sunmatrix_band.c b/src/sunmatrix/band/sunmatrix_band.c index 6eda3b854d..e1e7c843e8 100644 --- a/src/sunmatrix/band/sunmatrix_band.c +++ b/src/sunmatrix/band/sunmatrix_band.c @@ -26,6 +26,8 @@ #include #include +#include "sundials/sundials_errors.h" + #define ZERO SUN_RCONST(0.0) #define ONE SUN_RCONST(1.0) @@ -123,6 +125,8 @@ void SUNBandMatrix_Print(SUNMatrix A, FILE* outfile) SUNFunctionBegin(A->sunctx); sunindextype i, j, start, finish; + SUNAssertVoid(SUNMatGetID(A) == SUNMATRIX_BAND, SUN_ERR_ARG_WRONGTYPE); + /* perform operation */ fprintf(outfile, "\n"); for (i = 0; i < SM_ROWS_B(A); i++) @@ -361,7 +365,8 @@ SUNErrCode SUNMatScaleAdd_Band(sunrealtype c, SUNMatrix A, SUNMatrix B) /* Call separate routine in B has larger bandwidth(s) than A */ if ((SM_UBAND_B(B) > SM_UBAND_B(A)) || (SM_LBAND_B(B) > SM_LBAND_B(A))) { - return SMScaleAddNew_Band(c, A, B); + SUNCheckCall(SMScaleAddNew_Band(c, A, B)); + return SUN_SUCCESS; } /* Otherwise, perform operation in-place */ @@ -391,6 +396,9 @@ SUNErrCode SUNMatMatvec_Band(SUNMatrix A, N_Vector x, N_Vector y) yd = N_VGetArrayPointer(y); SUNCheckLastErr(); + SUNAssert(xd, SUN_ERR_MEM_FAIL); + SUNAssert(yd, SUN_ERR_MEM_FAIL); + /* Perform operation */ for (i = 0; i < SM_ROWS_B(A); i++) { yd[i] = ZERO; } for (j = 0; j < SM_COLUMNS_B(A); j++) @@ -407,6 +415,8 @@ SUNErrCode SUNMatSpace_Band(SUNMatrix A, long int* lenrw, long int* leniw) { SUNFunctionBegin(A->sunctx); SUNAssert(SUNMatGetID(A) == SUNMATRIX_BAND, SUN_ERR_ARG_WRONGTYPE); + SUNAssert(lenrw, SUN_ERR_ARG_CORRUPT); + SUNAssert(leniw, SUN_ERR_ARG_CORRUPT); *lenrw = SM_COLUMNS_B(A) * (SM_SUBAND_B(A) + SM_LBAND_B(A) + 1); *leniw = 7 + SM_COLUMNS_B(A); return SUN_SUCCESS; diff --git a/src/sunmatrix/dense/sunmatrix_dense.c b/src/sunmatrix/dense/sunmatrix_dense.c index 3dd43080de..19abe71af2 100644 --- a/src/sunmatrix/dense/sunmatrix_dense.c +++ b/src/sunmatrix/dense/sunmatrix_dense.c @@ -23,6 +23,8 @@ #include #include +#include "sundials/sundials_errors.h" + #define ZERO SUN_RCONST(0.0) #define ONE SUN_RCONST(1.0) @@ -300,10 +302,11 @@ SUNErrCode SUNMatScaleAdd_Dense(sunrealtype c, SUNMatrix A, SUNMatrix B) SUNErrCode SUNMatMatvec_Dense(SUNMatrix A, N_Vector x, N_Vector y) { + SUNFunctionBegin(A->sunctx); sunindextype i, j; sunrealtype *col_j, *xd, *yd; - SUNFunctionBegin(A->sunctx); + SUNAssert(SUNMatGetID(A) == SUNMATRIX_DENSE, SUN_ERR_ARG_WRONGTYPE); SUNCheck(compatibleMatrixAndVectors(A, x, y), SUN_ERR_ARG_DIMSMISMATCH); /* access vector data (return if NULL data pointers) */ @@ -312,6 +315,9 @@ SUNErrCode SUNMatMatvec_Dense(SUNMatrix A, N_Vector x, N_Vector y) yd = N_VGetArrayPointer(y); SUNCheckLastErr(); + SUNAssert(xd, SUN_ERR_MEM_FAIL); + SUNAssert(yd, SUN_ERR_MEM_FAIL); + /* Perform operation y = Ax */ for (i = 0; i < SM_ROWS_D(A); i++) { yd[i] = ZERO; } for (j = 0; j < SM_COLUMNS_D(A); j++) @@ -326,6 +332,8 @@ SUNErrCode SUNMatSpace_Dense(SUNMatrix A, long int* lenrw, long int* leniw) { SUNFunctionBegin(A->sunctx); SUNAssert(SUNMatGetID(A) == SUNMATRIX_DENSE, SUN_ERR_ARG_WRONGTYPE); + SUNAssert(lenrw, SUN_ERR_ARG_CORRUPT); + SUNAssert(leniw, SUN_ERR_ARG_CORRUPT); *lenrw = SM_LDATA_D(A); *leniw = 3 + SM_COLUMNS_D(A); return SUN_SUCCESS; diff --git a/src/sunmatrix/sparse/sunmatrix_sparse.c b/src/sunmatrix/sparse/sunmatrix_sparse.c index 619024fdd8..957b9e5b7b 100644 --- a/src/sunmatrix/sparse/sunmatrix_sparse.c +++ b/src/sunmatrix/sparse/sunmatrix_sparse.c @@ -24,9 +24,12 @@ #include #include #include +#include #include #include +#include "sundials/sundials_errors.h" + #define ZERO SUN_RCONST(0.0) #define ONE SUN_RCONST(1.0) @@ -340,7 +343,11 @@ SUNErrCode SUNSparseMatrix_Realloc(SUNMatrix A) /* perform reallocation */ SM_INDEXVALS_S(A) = (sunindextype*)realloc(SM_INDEXVALS_S(A), nzmax * sizeof(sunindextype)); + SUNAssert(SM_INDEXVALS_S(A), SUN_ERR_MALLOC_FAIL); + SM_DATA_S(A) = (sunrealtype*)realloc(SM_DATA_S(A), nzmax * sizeof(sunrealtype)); + SUNAssert(SM_DATA_S(A), SUN_ERR_MALLOC_FAIL); + SM_NNZ_S(A) = nzmax; return SUN_SUCCESS; @@ -362,8 +369,12 @@ SUNErrCode SUNSparseMatrix_Reallocate(SUNMatrix A, sunindextype NNZ) /* perform reallocation */ SM_INDEXVALS_S(A) = (sunindextype*)realloc(SM_INDEXVALS_S(A), NNZ * sizeof(sunindextype)); + SUNAssert(SM_INDEXVALS_S(A), SUN_ERR_MALLOC_FAIL); + SM_DATA_S(A) = (sunrealtype*)realloc(SM_DATA_S(A), NNZ * sizeof(sunrealtype)); - SM_NNZ_S(A) = NNZ; + SUNAssert(SM_DATA_S(A), SUN_ERR_MALLOC_FAIL); + + SM_NNZ_S(A) = NNZ; return SUN_SUCCESS; } @@ -576,9 +587,13 @@ SUNErrCode SUNMatCopy_Sparse(SUNMatrix A, SUNMatrix B) { SM_INDEXVALS_S(B) = (sunindextype*)realloc(SM_INDEXVALS_S(B), A_nz * sizeof(sunindextype)); - SM_DATA_S(B) = (sunrealtype*)realloc(SM_DATA_S(B), - A_nz * sizeof(sunrealtype)); - SM_NNZ_S(B) = A_nz; + SUNAssert(SM_INDEXVALS_S(B), SUN_ERR_MALLOC_FAIL); + + SM_DATA_S(B) = (sunrealtype*)realloc(SM_DATA_S(B), + A_nz * sizeof(sunrealtype)); + SUNAssert(SM_DATA_S(B), SUN_ERR_MALLOC_FAIL); + + SM_NNZ_S(B) = A_nz; } /* zero out B so that copy works correctly */ @@ -623,9 +638,10 @@ SUNErrCode SUNMatScaleAddI_Sparse(sunrealtype c, SUNMatrix A) } /* access data arrays from A */ - Ap = Ai = NULL; - Ax = NULL; - Ap = SM_INDEXPTRS_S(A); + Ap = NULL; + Ai = NULL; + Ax = NULL; + Ap = SM_INDEXPTRS_S(A); SUNCheck(Ap, SUN_ERR_ARG_CORRUPT); Ai = SM_INDEXVALS_S(A); SUNCheck(Ai, SUN_ERR_ARG_CORRUPT); @@ -753,9 +769,10 @@ SUNErrCode SUNMatScaleAddI_Sparse(sunrealtype c, SUNMatrix A) SUNCheckLastErr(); /* access data from CSR structures (return if failure) */ - Cp = Ci = NULL; - Cx = NULL; - Cp = SM_INDEXPTRS_S(C); + Cp = NULL; + Ci = NULL; + Cx = NULL; + Cp = SM_INDEXPTRS_S(C); SUNCheck(Cp, SUN_ERR_ARG_CORRUPT); Ci = SM_INDEXVALS_S(C); SUNCheck(Ci, SUN_ERR_ARG_CORRUPT); @@ -983,9 +1000,10 @@ SUNErrCode SUNMatScaleAdd_Sparse(sunrealtype c, SUNMatrix A, SUNMatrix B) SUNCheckLastErr(); /* access data from CSR structures (return if failure) */ - Cp = Ci = NULL; - Cx = NULL; - Cp = SM_INDEXPTRS_S(C); + Cp = NULL; + Ci = NULL; + Cx = NULL; + Cp = SM_INDEXPTRS_S(C); SUNCheck(Cp, SUN_ERR_ARG_CORRUPT); Ci = SM_INDEXVALS_S(C); SUNCheck(Ci, SUN_ERR_ARG_CORRUPT); @@ -1083,6 +1101,8 @@ SUNErrCode SUNMatSpace_Sparse(SUNMatrix A, long int* lenrw, long int* leniw) { SUNFunctionBegin(A->sunctx); SUNAssert(SUNMatGetID(A) == SUNMATRIX_SPARSE, SUN_ERR_ARG_WRONGTYPE); + SUNAssert(lenrw, SUN_ERR_ARG_CORRUPT); + SUNAssert(leniw, SUN_ERR_ARG_CORRUPT); *lenrw = SM_NNZ_S(A); *leniw = 10 + SM_NP_S(A) + SM_NNZ_S(A); return SUN_SUCCESS; @@ -1224,12 +1244,13 @@ SUNErrCode Matvec_SparseCSR(SUNMatrix A, N_Vector x, N_Vector y) */ SUNErrCode format_convert(const SUNMatrix A, SUNMatrix B) { + SUNFunctionBegin(A->sunctx); + sunrealtype *Ax, *Bx; sunindextype *Ap, *Aj; sunindextype *Bp, *Bi; sunindextype n_row, n_col, nnz; sunindextype n, col, csum, row, last; - SUNFunctionBegin(A->sunctx); if (SM_SPARSETYPE_S(A) == SM_SPARSETYPE_S(B)) { @@ -1238,8 +1259,11 @@ SUNErrCode format_convert(const SUNMatrix A, SUNMatrix B) } Ap = SM_INDEXPTRS_S(A); + SUNCheck(Ap, SUN_ERR_ARG_CORRUPT); Aj = SM_INDEXVALS_S(A); + SUNCheck(Aj, SUN_ERR_ARG_CORRUPT); Ax = SM_DATA_S(A); + SUNCheck(Ax, SUN_ERR_ARG_CORRUPT); n_row = (SM_SPARSETYPE_S(A) == CSR_MAT) ? SM_ROWS_S(A) : SM_COLUMNS_S(A); n_col = (SM_SPARSETYPE_S(A) == CSR_MAT) ? SM_COLUMNS_S(A) : SM_ROWS_S(A);