Skip to content

Commit

Permalink
address comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@balos1
balos1 committed Dec 12, 2023
1 parent b7491f2 commit 05fc80f
Show file tree
Hide file tree
Showing 6 changed files with 72 additions and 27 deletions.
9 changes: 6 additions & 3 deletions doc/shared/sunmatrix/SUNMatrix_Sparse.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion examples/arkode/F2003_custom/test_fsunlinsol_fortran_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
14 changes: 7 additions & 7 deletions examples/arkode/F2003_custom/test_fsunmatrix_fortran_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
12 changes: 11 additions & 1 deletion src/sunmatrix/band/sunmatrix_band.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@
#include <sundials/sundials_math.h>
#include <sunmatrix/sunmatrix_band.h>

#include "sundials/sundials_errors.h"

#define ZERO SUN_RCONST(0.0)
#define ONE SUN_RCONST(1.0)

Expand Down Expand Up @@ -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++)
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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++)
Expand All @@ -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;
Expand Down
10 changes: 9 additions & 1 deletion src/sunmatrix/dense/sunmatrix_dense.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@
#include <sundials/priv/sundials_errors_impl.h>
#include <sunmatrix/sunmatrix_dense.h>

#include "sundials/sundials_errors.h"

#define ZERO SUN_RCONST(0.0)
#define ONE SUN_RCONST(1.0)

Expand Down Expand Up @@ -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) */
Expand All @@ -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++)
Expand All @@ -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;
Expand Down
52 changes: 38 additions & 14 deletions src/sunmatrix/sparse/sunmatrix_sparse.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,12 @@
#include <stdlib.h>
#include <sundials/priv/sundials_errors_impl.h>
#include <sundials/sundials_math.h>
#include <sunmatrix/sunmatrix_band.h>
#include <sunmatrix/sunmatrix_dense.h>
#include <sunmatrix/sunmatrix_sparse.h>

#include "sundials/sundials_errors.h"

#define ZERO SUN_RCONST(0.0)
#define ONE SUN_RCONST(1.0)

Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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))
{
Expand All @@ -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);
Expand Down

0 comments on commit 05fc80f

Please sign in to comment.