From 5ca223717b2dcd6163d145338bb95da3d47c8e98 Mon Sep 17 00:00:00 2001 From: Tim Davis Date: Tue, 20 Feb 2024 16:33:39 -0600 Subject: [PATCH] . --- .../spex_cholesky_left_triangular_solve.c | 6 ++- .../spex_cholesky_up_triangular_solve.c | 42 ++++++++++++++++++- .../spex_left_lu_ref_triangular_solve.c | 4 ++ 3 files changed, 49 insertions(+), 3 deletions(-) diff --git a/SPEX/SPEX_Cholesky/Source/spex_cholesky_left_triangular_solve.c b/SPEX/SPEX_Cholesky/Source/spex_cholesky_left_triangular_solve.c index 1caa721e2..2514224b4 100644 --- a/SPEX/SPEX_Cholesky/Source/spex_cholesky_left_triangular_solve.c +++ b/SPEX/SPEX_Cholesky/Source/spex_cholesky_left_triangular_solve.c @@ -57,7 +57,11 @@ // c's default qsort static inline int compare (const void * a, const void * b) { - return ( *(int64_t*)a - *(int64_t*)b ); +// return ( *(int64_t*)a - *(int64_t*)b ); + + int64_t x = (* ((int64_t *) a)) ; + int64_t y = (* ((int64_t *) b)) ; + return (x < y ? -1 : ((x == y) ? 0 : 1)) ; } SPEX_info spex_cholesky_left_triangular_solve diff --git a/SPEX/SPEX_Cholesky/Source/spex_cholesky_up_triangular_solve.c b/SPEX/SPEX_Cholesky/Source/spex_cholesky_up_triangular_solve.c index aad6d68b5..414d74634 100644 --- a/SPEX/SPEX_Cholesky/Source/spex_cholesky_up_triangular_solve.c +++ b/SPEX/SPEX_Cholesky/Source/spex_cholesky_up_triangular_solve.c @@ -55,7 +55,11 @@ // c's default qsort static inline int compare (const void * a, const void * b) { - return ( *(int64_t*)a - *(int64_t*)b ); +// return ( *(int64_t*)a - *(int64_t*)b ); + + int64_t x = (* ((int64_t *) a)) ; + int64_t y = (* ((int64_t *) b)) ; + return (x < y ? -1 : ((x == y) ? 0 : 1)) ; } SPEX_info spex_cholesky_up_triangular_solve @@ -107,13 +111,39 @@ SPEX_info spex_cholesky_up_triangular_solve // xi[top..n-1] SPEX_CHECK(spex_cholesky_ereach(&top, xi, A, k, parent, c)); +if (top < 0 || top > n) +{ + HERE ; + fprintf (stderr, "Hey: top is wierd %" PRId64 " n is %" PRId64 "\n", top, n) ; + abort ( ) ; +} + for (i = top; i < n; i++) + { + int64_t j = xi [i] ; +if (j < 0 || j >= n) +{ + HERE ; + fprintf (stderr, "Hey: j is wierd %" PRId64 " n is %" PRId64 "\n", j, n) ; + abort ( ) ; +} + } + // Sort the nonzero pattern using quicksort (required by IPGE unlike in GE) qsort(&xi[top], n-top, sizeof(int64_t*), compare); // Reset x[i] = 0 for all i in nonzero pattern xi [top..n-1] for (i = top; i < n; i++) { - SPEX_MPZ_SET_UI(x->x.mpz[xi[i]],0); + int64_t j = xi [i] ; + +if (j < 0 || j >= n) +{ + HERE ; + fprintf (stderr, "Hey: j is wierd %" PRId64 " n is %" PRId64 "\n", j, n) ; + abort ( ) ; +} + + SPEX_MPZ_SET_UI(x->x.mpz[j],0); } // Reset value of x[k]. If the matrix is nonsingular, x[k] will @@ -149,6 +179,14 @@ SPEX_info spex_cholesky_up_triangular_solve { // Obtain the index of the current nonzero j = xi[p]; + +if (j < 0 || j >= n) +{ + HERE ; + fprintf (stderr, "Hey: j is wierd %" PRId64 " n is %" PRId64 "\n", j, n) ; + abort ( ) ; +} + SPEX_MPZ_SGN(&sgn, x->x.mpz[j]); if (sgn == 0) continue; // If x[j] == 0 no work must be done diff --git a/SPEX/SPEX_LU/Source/spex_left_lu_ref_triangular_solve.c b/SPEX/SPEX_LU/Source/spex_left_lu_ref_triangular_solve.c index f28161020..be1f9768d 100644 --- a/SPEX/SPEX_LU/Source/spex_left_lu_ref_triangular_solve.c +++ b/SPEX/SPEX_LU/Source/spex_left_lu_ref_triangular_solve.c @@ -74,6 +74,10 @@ static inline int compare (const void * a, const void * b) { return ( *(int64_t*)a - *(int64_t*)b ); + + int64_t x = (* ((int64_t *) a)) ; + int64_t y = (* ((int64_t *) b)) ; + return (x < y ? -1 : ((x == y) ? 0 : 1)) ; } SPEX_info spex_left_lu_ref_triangular_solve // sparse REF triangular solve