From 82f90810f6bc00b2281f8e54f8e5f1fa3a786e57 Mon Sep 17 00:00:00 2001 From: Tim Davis Date: Mon, 4 Sep 2023 11:37:07 -0500 Subject: [PATCH] starting SPQR Tcov tests for int32 --- SPQR/Tcov/qrtest.cpp | 2 +- SPQR/Tcov/qrtestc.c | 137 ++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 137 insertions(+), 2 deletions(-) diff --git a/SPQR/Tcov/qrtest.cpp b/SPQR/Tcov/qrtest.cpp index 6a2165b633..94805c7260 100644 --- a/SPQR/Tcov/qrtest.cpp +++ b/SPQR/Tcov/qrtest.cpp @@ -40,7 +40,7 @@ double yy = 0 ; // ============================================================================= extern "C" { -void qrtest_C +void qrtest_C // handles both int32_t and int64_t versions ( cholmod_sparse *A, double anorm, diff --git a/SPQR/Tcov/qrtestc.c b/SPQR/Tcov/qrtestc.c index 33377bcdbd..e9e152baff 100644 --- a/SPQR/Tcov/qrtestc.c +++ b/SPQR/Tcov/qrtestc.c @@ -13,7 +13,11 @@ #define MAX(a,b) (((a) > (b)) ? (a) : (b)) -void qrtest_C +//------------------------------------------------------------------------------ +// qrtest_C64: int64_t version +//------------------------------------------------------------------------------ + +static void qrtest_C64 ( cholmod_sparse *A, double anorm, @@ -115,3 +119,134 @@ void qrtest_C cholmod_l_free_dense (&B, cc) ; } + +//------------------------------------------------------------------------------ +// qrtest_C32: int32_t version +//------------------------------------------------------------------------------ + +static void qrtest_C32 +( + cholmod_sparse *A, + double anorm, + double errs [5], + double maxresid [2][2], + cholmod_common *cc +) +{ + cholmod_dense *B, *X, *Resid ; + cholmod_sparse *Bsparse, *Xsparse ; + SuiteSparseQR_C_factorization *QR ; + double resid, one [2] = {1,0}, minusone [2] = {-1,0} ; + int32_t m, n ; +#ifndef NEXPERT + cholmod_dense *Y ; + int split ; +#endif + + m = A->nrow ; + n = A->ncol ; + + B = cholmod_ones (m, 1, A->xtype, cc) ; + + /* X = A\B */ + X = SuiteSparseQR_C_backslash_default (A, B, cc) ; + + /* Resid = A*X - B */ + Resid = cholmod_copy_dense (B, cc) ; + cholmod_sdmult (A, 0, one, minusone, X, Resid, cc) ; + + /* resid = norm (Resid,1) */ + resid = cholmod_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ; + resid = (resid < 0 || resid != resid) ? 9e99 : resid ; + cholmod_free_dense (&Resid, cc) ; + cholmod_free_dense (&X, cc) ; + + maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ; + printf ("Resid_C1 %d : %g\n", m>n, resid) ; + + /* X = A\B */ + Bsparse = cholmod_dense_to_sparse (B, 1, cc) ; + Xsparse = SuiteSparseQR_C_backslash_sparse (2, -2, A, Bsparse, cc) ; + X = cholmod_sparse_to_dense (Xsparse, cc) ; + cholmod_free_sparse (&Bsparse, cc) ; + cholmod_free_sparse (&Xsparse, cc) ; + + /* Resid = A*X - B */ + Resid = cholmod_copy_dense (B, cc) ; + cholmod_sdmult (A, 0, one, minusone, X, Resid, cc) ; + + /* resid = norm (Resid,1) */ + resid = cholmod_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ; + resid = (resid < 0 || resid != resid) ? 9e99 : resid ; + cholmod_free_dense (&Resid, cc) ; + cholmod_free_dense (&X, cc) ; + + maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ; + printf ("Resid_C2 %d : %g\n", m>n, resid) ; + +#ifndef NEXPERT + + for (split = 0 ; split <= 1 ; split++) + { + if (split) + { + /* split symbolic/numeric QR factorization */ + QR = SuiteSparseQR_C_symbolic (2, 1, A, cc) ; + SuiteSparseQR_C_numeric (-2, A, QR, cc) ; + } + else + { + /* QR factorization, single pass */ + QR = SuiteSparseQR_C_factorize (2, -2, A, cc) ; + } + + /* Y = Q'*B */ + Y = SuiteSparseQR_C_qmult (0, QR, B, cc) ; + + /* X = E*(R\Y) */ + X = SuiteSparseQR_C_solve (1, QR, Y, cc) ; + + /* Resid = A*X - B */ + Resid = cholmod_copy_dense (B, cc) ; + cholmod_sdmult (A, 0, one, minusone, X, Resid, cc) ; + + /* resid = norm (Resid,1) */ + resid = cholmod_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ; + resid = (resid < 0 || resid != resid) ? 9e99 : resid ; + cholmod_free_dense (&Resid, cc) ; + cholmod_free_dense (&X, cc) ; + + maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ; + printf ("Resid_C3 %d : %g\n", m>n, resid) ; + + cholmod_free_dense (&Y, cc) ; + SuiteSparseQR_C_free (&QR, cc) ; + } +#endif + + cholmod_free_dense (&B, cc) ; +} + +//------------------------------------------------------------------------------ +// qrtest_C: handles both int32_t and int64_t versions +//------------------------------------------------------------------------------ + +void qrtest_C +( + cholmod_sparse *A, + double anorm, + double errs [5], + double maxresid [2][2], + cholmod_common *cc +) +{ + if (A->itype == CHOLMOD_INT) + { + qrtest_C32 (A, anorm, errs, maxresid, cc) ; + } + else + { + qrtest_C64 (A, anorm, errs, maxresid, cc) ; + } +} +