From 947b4c96dc551e837468dc0ded39a9ffd4847893 Mon Sep 17 00:00:00 2001 From: Tim Davis Date: Mon, 4 Sep 2023 18:02:49 -0500 Subject: [PATCH] test coverage for int32 version of SPQR --- .gitignore | 2 + CHOLMOD/Check/cholmod_check.c | 2 +- CHOLMOD/Cholesky/cholmod_amd.c | 2 +- CHOLMOD/Cholesky/cholmod_colamd.c | 4 +- CHOLMOD/Doc/ChangeLog | 2 + CHOLMOD/Include/cholmod_internal.h | 88 +- CHOLMOD/Include/cholmod_types.h | 139 ++ CHOLMOD/Partition/cholmod_camd.c | 2 +- CHOLMOD/Partition/cholmod_ccolamd.c | 6 +- CHOLMOD/Partition/cholmod_csymamd.c | 4 +- CHOLMOD/Tcov/cm.h | 43 +- CHOLMOD/Tcov/null2.c | 6 +- SPQR/Include/SuiteSparseQR_C.h | 20 +- SPQR/Include/spqr_cholmod_wrappers.hpp | 414 ++++ SPQR/Source/SuiteSparseQR_C.cpp | 1 + SPQR/Source/spqr_cholmod_wrappers.cpp | 526 +++++ SPQR/Tcov/Makefile | 27 +- SPQR/Tcov/qrtest.cpp | 2653 +---------------------- SPQR/Tcov/qrtest_template.hpp | 2686 ++++++++++++++++++++++++ SPQR/Tcov/qrtestc.c | 222 +- SPQR/Tcov/qrtestc_template.c | 127 ++ 21 files changed, 3987 insertions(+), 2989 deletions(-) create mode 100644 CHOLMOD/Include/cholmod_types.h create mode 100644 SPQR/Tcov/qrtest_template.hpp create mode 100644 SPQR/Tcov/qrtestc_template.c diff --git a/.gitignore b/.gitignore index d0d4223fd7..c6499d0a82 100644 --- a/.gitignore +++ b/.gitignore @@ -168,7 +168,9 @@ SPQR/Tcov/gpu_results.txt SPQR/Tcov/gpuqrengine_demo SPQR/Tcov/qrdemo_gpu SPQR/Tcov/qrtest +SPQR/Tcov/qrtest32 SPQR/Tcov/qrtest_out.txt +SPQR/Tcov/qrtest_out32.txt SPQR/Tcov/troll.m SPQR/Tcov/cov.out diff --git a/CHOLMOD/Check/cholmod_check.c b/CHOLMOD/Check/cholmod_check.c index 45972b6f77..857163e409 100644 --- a/CHOLMOD/Check/cholmod_check.c +++ b/CHOLMOD/Check/cholmod_check.c @@ -69,7 +69,7 @@ /* === printing definitions ================================================= */ /* ========================================================================== */ -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) #define I8 "%8ld" #define I_8 "%-8ld" #else diff --git a/CHOLMOD/Cholesky/cholmod_amd.c b/CHOLMOD/Cholesky/cholmod_amd.c index 089830bcf5..ef4eff0068 100644 --- a/CHOLMOD/Cholesky/cholmod_amd.c +++ b/CHOLMOD/Cholesky/cholmod_amd.c @@ -166,7 +166,7 @@ int CHOLMOD(amd) Control [AMD_AGGRESSIVE] = Common->method [Common->current].aggressive ; } -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) amd_l2 (n, C->p, C->i, Len, C->nzmax, cnz, Nv, Next, Perm, Head, Elen, Degree, Wi, Control, Info) ; #else diff --git a/CHOLMOD/Cholesky/cholmod_colamd.c b/CHOLMOD/Cholesky/cholmod_colamd.c index 6980afecb5..c14b7d0d3f 100644 --- a/CHOLMOD/Cholesky/cholmod_colamd.c +++ b/CHOLMOD/Cholesky/cholmod_colamd.c @@ -92,7 +92,7 @@ int CHOLMOD(colamd) s = CHOLMOD(mult_size_t) (nrow, 4, &ok) ; s = CHOLMOD(add_size_t) (s, ncol, &ok) ; -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) alen = colamd_l_recommended (A->nzmax, ncol, nrow) ; colamd_l_set_defaults (knobs) ; #else @@ -151,7 +151,7 @@ int CHOLMOD(colamd) Int stats [COLAMD_STATS] ; Cp = C->p ; -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) colamd_l (ncol, nrow, alen, C->i, Cp, knobs, stats) ; #else colamd (ncol, nrow, alen, C->i, Cp, knobs, stats) ; diff --git a/CHOLMOD/Doc/ChangeLog b/CHOLMOD/Doc/ChangeLog index c78a06c44a..684d6a83be 100644 --- a/CHOLMOD/Doc/ChangeLog +++ b/CHOLMOD/Doc/ChangeLog @@ -4,6 +4,8 @@ FIXME DATE, 2023: version 4.2.0 * bug fix: typecast was broken in cholmod_postorder.c * minor fixes to debug scaffolding and test code: cholmod_internal.h, cholmod_metis_wrapper.c, Tcov/memory.c. + * support for SPQR: minor changes to support SPQR v4.2.0 with 32-bit + integers; no change to user-visible API. June 16, 2023: version 4.0.4 diff --git a/CHOLMOD/Include/cholmod_internal.h b/CHOLMOD/Include/cholmod_internal.h index 8b85695d17..e14eeb48ce 100644 --- a/CHOLMOD/Include/cholmod_internal.h +++ b/CHOLMOD/Include/cholmod_internal.h @@ -133,45 +133,12 @@ /* === int/long and double/float definitions ================================ */ /* ========================================================================== */ -/* CHOLMOD is designed for 3 types of integer variables: - * - * (1) all integers are int - * (2) most integers are int, some are int64_t - * (3) all integers are int64_t - * - * and two kinds of floating-point values: - * - * (1) double - * (2) float - * - * the complex types (ANSI-compatible complex, and MATLAB-compatable zomplex) - * are based on the double or float type, and are not selected here. They - * are typically selected via template routines. - * - * This gives 6 different modes in which CHOLMOD can be compiled (only the - * first two are currently supported): - * - * DINT double, int prefix: cholmod_ - * DLONG double, int64_t prefix: cholmod_l_ - * DMIX double, mixed int/int64_t prefix: cholmod_m_ - * SINT float, int prefix: cholmod_si_ - * SLONG float, int64_t prefix: cholmod_sl_ - * SMIX float, mixed int/log prefix: cholmod_sm_ - * - * These are selected with compile time flags (-DDLONG, for example). If no - * flag is selected, the default is DINT. - * - * All six versions use the same include files. The user-visible include files - * are completely independent of which int/long/double/float version is being - * used. The integer / real types in all data structures (sparse, triplet, - * dense, common, and triplet) are defined at run-time, not compile-time, so - * there is only one "cholmod_sparse" data type. Void pointers are used inside - * that data structure to point to arrays of the proper type. Each data - * structure has an itype and dtype field which determines the kind of basic - * types used. These are defined in Include/cholmod_core.h. - * - * FUTURE WORK: support all six types (float, and mixed int/long) - */ +#include "cholmod_types.h" + +#ifndef DLONG +// GPU acceleration only available for the DLONG case (double, int64) +#undef SUITESPARSE_CUDA +#endif /* -------------------------------------------------------------------------- */ /* routines for doing arithmetic on size_t, and checking for overflow */ @@ -182,49 +149,6 @@ size_t cholmod_mult_size_t (size_t a, size_t k, int *ok) ; size_t cholmod_l_add_size_t (size_t a, size_t b, int *ok) ; size_t cholmod_l_mult_size_t (size_t a, size_t k, int *ok) ; -/* -------------------------------------------------------------------------- */ -/* double (also complex double), int64_t */ -/* -------------------------------------------------------------------------- */ - -#ifdef DLONG -#define Real double -#define Int int64_t -#define UInt uint64_t -#define Int_max INT64_MAX -#define CHOLMOD(name) cholmod_l_ ## name -#define LONG -#define DOUBLE -#define ITYPE CHOLMOD_LONG -#define DTYPE CHOLMOD_DOUBLE -#define ID "%" PRId64 - -/* -------------------------------------------------------------------------- */ -/* double (also complex double), int: this is the default */ -/* -------------------------------------------------------------------------- */ - -#else - -#ifndef DINT -#define DINT -#endif -#define INT -#define DOUBLE - -#define Real double -#define Int int32_t -#define UInt uint32_t -#define Int_max INT32_MAX -#define CHOLMOD(name) cholmod_ ## name -#define ITYPE CHOLMOD_INT -#define DTYPE CHOLMOD_DOUBLE -#define ID "%d" - -/* GPU acceleration is not available for the int version of CHOLMOD */ -#undef SUITESPARSE_CUDA - -#endif - - /* ========================================================================== */ /* === Include/cholmod_complexity.h ========================================= */ /* ========================================================================== */ diff --git a/CHOLMOD/Include/cholmod_types.h b/CHOLMOD/Include/cholmod_types.h new file mode 100644 index 0000000000..e82f783d73 --- /dev/null +++ b/CHOLMOD/Include/cholmod_types.h @@ -0,0 +1,139 @@ +//------------------------------------------------------------------------------ +// CHOLMOD/Include/cholmod_types.h +//------------------------------------------------------------------------------ + +// CHOLMOD/Include/cholmod_types.h. Copyright (C) 2005-2023, +// Timothy A. Davis. All Rights Reserved. +// SPDX-License-Identifier: Apache-2.0 + +//------------------------------------------------------------------------------ + +// CHOLMOD internal include file: defining integer and floating-point types. +// This file is suitable for inclusion in C and C++ codes. It can be +// #include'd more than once. + +// The #include'ing file defines one of four macros (DINT, DLONG, SINT, or +// SLONG) before #include'ing this file. + +/* ========================================================================== */ +/* === int/long and double/float definitions ================================ */ +/* ========================================================================== */ + +/* CHOLMOD is designed for 3 types of integer variables: + * + * (1) all integers are int + * (2) most integers are int, some are int64_t + * (3) all integers are int64_t + * + * and two kinds of floating-point values: + * + * (1) double + * (2) float + * + * the complex types (ANSI-compatible complex, and MATLAB-compatable zomplex) + * are based on the double or float type, and are not selected here. They + * are typically selected via template routines. + * + * This gives 6 different modes in which CHOLMOD can be compiled (only the + * first two are currently supported): + * + * DINT double, int prefix: cholmod_ + * DLONG double, int64_t prefix: cholmod_l_ + * DMIX double, mixed int/int64_t prefix: cholmod_m_ + * SINT float, int prefix: cholmod_si_ + * SLONG float, int64_t prefix: cholmod_sl_ + * SMIX float, mixed int/log prefix: cholmod_sm_ + * + * These are selected with compile time flags (-DDLONG, for example). If no + * flag is selected, the default is DINT. + * + * All six versions use the same include files. The user-visible include files + * are completely independent of which int/long/double/float version is being + * used. The integer / real types in all data structures (sparse, triplet, + * dense, common, and triplet) are defined at run-time, not compile-time, so + * there is only one "cholmod_sparse" data type. Void pointers are used inside + * that data structure to point to arrays of the proper type. Each data + * structure has an itype and dtype field which determines the kind of basic + * types used. These are defined in Include/cholmod_core.h. + * + * FUTURE WORK: support all six types (float, and mixed int/long) + * SINT and SLONG are in progress. + */ + +// ----------------------------------------------------------------------------- + +#undef Real +#undef Int +#undef UInt +#undef Int_max +#undef CHOLMOD +#undef ITYPE +#undef DTYPE +#undef ID + +#if defined ( SLONG ) + + //-------------------------------------------------------------------------- + // SLONG: float (also complex float), int64_t + //-------------------------------------------------------------------------- + + #define Real float + #define Int int64_t + #define UInt uint64_t + #define Int_max INT64_MAX + #define CHOLMOD(name) cholmod_sl_ ## name + #define ITYPE CHOLMOD_LONG + #define DTYPE CHOLMOD_SINGLE + #define ID "%" PRId64 + +#elif defined ( SINT ) + + //-------------------------------------------------------------------------- + // SINT: float (also complex float), int32_t + //-------------------------------------------------------------------------- + + #define Real float + #define Int int32_t + #define UInt uint32_t + #define Int_max INT32_MAX + #define CHOLMOD(name) cholmod_si_ ## name + #define ITYPE CHOLMOD_INT + #define DTYPE CHOLMOD_SINGLE + #define ID "%d" + +#elif defined ( DLONG ) + + //-------------------------------------------------------------------------- + // DLONG: double (also complex double), int64_t + //-------------------------------------------------------------------------- + + #define Real double + #define Int int64_t + #define UInt uint64_t + #define Int_max INT64_MAX + #define CHOLMOD(name) cholmod_l_ ## name + #define ITYPE CHOLMOD_LONG + #define DTYPE CHOLMOD_DOUBLE + #define ID "%" PRId64 + +#else + + //-------------------------------------------------------------------------- + // DINT: double (also complex double), int32 + //-------------------------------------------------------------------------- + + #ifndef DINT + #define DINT + #endif + + #define Real double + #define Int int32_t + #define UInt uint32_t + #define Int_max INT32_MAX + #define CHOLMOD(name) cholmod_ ## name + #define ITYPE CHOLMOD_INT + #define DTYPE CHOLMOD_DOUBLE + #define ID "%d" + +#endif + diff --git a/CHOLMOD/Partition/cholmod_camd.c b/CHOLMOD/Partition/cholmod_camd.c index 7b0b78965c..40747391b0 100644 --- a/CHOLMOD/Partition/cholmod_camd.c +++ b/CHOLMOD/Partition/cholmod_camd.c @@ -177,7 +177,7 @@ int CHOLMOD(camd) Control [CAMD_AGGRESSIVE] = Common->method [Common->current].aggressive; } -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) /* DEBUG (camd_l_debug_init ("cholmod_l_camd")) ; */ camd_l2 (n, C->p, C->i, Len, C->nzmax, cnz, Nv, Next, Perm, Head, Elen, Degree, Wi, Control, Info, Cmember, BucketSet) ; diff --git a/CHOLMOD/Partition/cholmod_ccolamd.c b/CHOLMOD/Partition/cholmod_ccolamd.c index 231cef6ea2..08b79c9b54 100644 --- a/CHOLMOD/Partition/cholmod_ccolamd.c +++ b/CHOLMOD/Partition/cholmod_ccolamd.c @@ -67,7 +67,7 @@ static int ccolamd_interface /* ---------------------------------------------------------------------- */ /* get parameters */ -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) ccolamd_l_set_defaults (knobs) ; #else ccolamd_set_defaults (knobs) ; @@ -90,7 +90,7 @@ static int ccolamd_interface if (ok) { -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) ccolamd_l (ncol, nrow, alen, C->i, C->p, knobs, stats, Cmember) ; #else ccolamd (ncol, nrow, alen, C->i, C->p, knobs, stats, Cmember) ; @@ -164,7 +164,7 @@ int CHOLMOD(ccolamd) /* allocate workspace */ /* ---------------------------------------------------------------------- */ -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) alen = ccolamd_l_recommended (A->nzmax, ncol, nrow) ; #else alen = ccolamd_recommended (A->nzmax, ncol, nrow) ; diff --git a/CHOLMOD/Partition/cholmod_csymamd.c b/CHOLMOD/Partition/cholmod_csymamd.c index 3bdafaa88f..03d46e1f70 100644 --- a/CHOLMOD/Partition/cholmod_csymamd.c +++ b/CHOLMOD/Partition/cholmod_csymamd.c @@ -86,7 +86,7 @@ int CHOLMOD(csymamd) perm = Common->Head ; /* size nrow+1 (i/l/l) */ /* get parameters */ -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) ccolamd_l_set_defaults (knobs) ; #else ccolamd_set_defaults (knobs) ; @@ -103,7 +103,7 @@ int CHOLMOD(csymamd) calloc_func = SuiteSparse_config_calloc_func_get ( ) ; free_func = SuiteSparse_config_free_func_get ( ) ; -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) csymamd_l (nrow, A->i, A->p, perm, knobs, stats, calloc_func, free_func, diff --git a/CHOLMOD/Tcov/cm.h b/CHOLMOD/Tcov/cm.h index b8655544f3..ce83669294 100644 --- a/CHOLMOD/Tcov/cm.h +++ b/CHOLMOD/Tcov/cm.h @@ -16,43 +16,6 @@ #include #include -/* -------------------------------------------------------------------------- */ -/* double, int64_t */ -/* -------------------------------------------------------------------------- */ - -#ifdef DLONG -#define Real double -#define Int int64_t -#define UInt uint64_t -#define Int_max INT64_MAX -#define CHOLMOD(name) cholmod_l_ ## name -#define LONG -#define DOUBLE -#define ITYPE CHOLMOD_LONG -#define DTYPE CHOLMOD_DOUBLE - -/* -------------------------------------------------------------------------- */ -/* double, int: this is the default */ -/* -------------------------------------------------------------------------- */ - -#else - -#ifndef DINT -#define DINT -#endif -#define INT -#define DOUBLE - -#define Real double -#define Int int32_t -#define UInt uint32_t -#define Int_max INT32_MAX -#define CHOLMOD(name) cholmod_ ## name -#define ITYPE CHOLMOD_INT -#define DTYPE CHOLMOD_DOUBLE - -#endif - /* -------------------------------------------------------------------------- */ #include "cholmod_internal.h" @@ -183,9 +146,7 @@ void camdtest (cholmod_sparse *A) ; /* AMD, COLAMD, and CCOLAMD */ /* -------------------------------------------------------------------------- */ -#ifdef LONG - -#define ID "%" PRId64 +#if ( ITYPE == CHOLMOD_LONG ) #define AMD_order amd_l_order #define AMD_defaults amd_l_defaults @@ -238,8 +199,6 @@ void camdtest (cholmod_sparse *A) ; #else -#define ID "%d" - #define AMD_order amd_order #define AMD_defaults amd_defaults #define AMD_control amd_control diff --git a/CHOLMOD/Tcov/null2.c b/CHOLMOD/Tcov/null2.c index 1087fa7d99..4b17bbf67d 100644 --- a/CHOLMOD/Tcov/null2.c +++ b/CHOLMOD/Tcov/null2.c @@ -878,7 +878,7 @@ void null2 (cholmod_triplet *Tok, int do_nantests) cm->print = 1 ; cm->print = 4 ; -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) C->itype = CHOLMOD_INT ; #else C->itype = CHOLMOD_LONG ; @@ -1183,7 +1183,7 @@ void null2 (cholmod_triplet *Tok, int do_nantests) ok = CHOLMOD(print_factor)(L, "L OK", cm) ; OK (ok) ; cm->print = 4 ; -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) L->itype = CHOLMOD_INT ; #else L->itype = CHOLMOD_LONG ; @@ -2369,7 +2369,7 @@ if (do_nantests) ok = CHOLMOD(print_triplet)(T, "T ok", cm) ; OK (ok) ; cm->print = 4 ; -#ifdef LONG +#if ( ITYPE == CHOLMOD_LONG ) T->itype = CHOLMOD_INT ; #else T->itype = CHOLMOD_LONG ; diff --git a/SPQR/Include/SuiteSparseQR_C.h b/SPQR/Include/SuiteSparseQR_C.h index 389562b748..df8559474f 100644 --- a/SPQR/Include/SuiteSparseQR_C.h +++ b/SPQR/Include/SuiteSparseQR_C.h @@ -26,7 +26,7 @@ int64_t SuiteSparseQR_C /* returns rank(A) estimate, (-1) if failure */ /* inputs: */ int ordering, /* all, except 3:given treated as 0:fixed */ double tol, /* columns with 2-norm <= tol treated as 0 */ - int64_t econ, /* e = max(min(m,econ),rank(A)) */ + int64_t econ, /* e = max(min(m,econ),rank(A)) */ int getCTX, /* 0: Z=C (e-by-k), 1: Z=C', 2: Z=X (e-by-k) */ cholmod_sparse *A, /* m-by-n sparse matrix to factorize */ cholmod_sparse *Bsparse, /* sparse m-by-k B */ @@ -35,9 +35,9 @@ int64_t SuiteSparseQR_C /* returns rank(A) estimate, (-1) if failure */ cholmod_sparse **Zsparse, /* sparse Z */ cholmod_dense **Zdense, /* dense Z */ cholmod_sparse **R, /* e-by-n sparse matrix */ - int64_t **E, /* size n column perm, NULL if identity */ + int64_t **E, /* size n column perm, NULL if identity */ cholmod_sparse **H, /* m-by-nh Householder vectors */ - int64_t **HPinv, /* size m row permutation */ + int64_t **HPinv, /* size m row permutation */ cholmod_dense **HTau, /* 1-by-nh Householder coefficients */ cholmod_common *cc /* workspace and parameters */ ) ; @@ -47,7 +47,7 @@ int32_t SuiteSparseQR_i_C /* returns rank(A) estimate, (-1) if failure */ /* inputs: */ int ordering, /* all, except 3:given treated as 0:fixed */ double tol, /* columns with 2-norm <= tol treated as 0 */ - int32_t econ, /* e = max(min(m,econ),rank(A)) */ + int32_t econ, /* e = max(min(m,econ),rank(A)) */ int getCTX, /* 0: Z=C (e-by-k), 1: Z=C', 2: Z=X (e-by-k) */ cholmod_sparse *A, /* m-by-n sparse matrix to factorize */ cholmod_sparse *Bsparse, /* sparse m-by-k B */ @@ -56,9 +56,9 @@ int32_t SuiteSparseQR_i_C /* returns rank(A) estimate, (-1) if failure */ cholmod_sparse **Zsparse, /* sparse Z */ cholmod_dense **Zdense, /* dense Z */ cholmod_sparse **R, /* e-by-n sparse matrix */ - int32_t **E, /* size n column perm, NULL if identity */ + int32_t **E, /* size n column perm, NULL if identity */ cholmod_sparse **H, /* m-by-nh Householder vectors */ - int32_t **HPinv, /* size m row permutation */ + int32_t **HPinv, /* size m row permutation */ cholmod_dense **HTau, /* 1-by-nh Householder coefficients */ cholmod_common *cc /* workspace and parameters */ ) ; @@ -73,12 +73,12 @@ int64_t SuiteSparseQR_C_QR /* returns rank(A) est., (-1) if failure */ /* inputs: */ int ordering, /* all, except 3:given treated as 0:fixed */ double tol, /* columns with 2-norm <= tol treated as 0 */ - int64_t econ, /* e = max(min(m,econ),rank(A)) */ + int64_t econ, /* e = max(min(m,econ),rank(A)) */ cholmod_sparse *A, /* m-by-n sparse matrix to factorize */ /* outputs: */ cholmod_sparse **Q, /* m-by-e sparse matrix */ cholmod_sparse **R, /* e-by-n sparse matrix */ - int64_t **E, /* size n column perm, NULL if identity */ + int64_t **E, /* size n column perm, NULL if identity */ cholmod_common *cc /* workspace and parameters */ ) ; @@ -87,12 +87,12 @@ int32_t SuiteSparseQR_i_C_QR // returns rank(A) estimate, (-1) if failu // inputs: int ordering, // all, except 3:given treated as 0:fixed double tol, // columns with 2-norm <= tol are treated as 0 - int32_t econ, // e = max(min(m,econ),rank(A)) + int32_t econ, // e = max(min(m,econ),rank(A)) cholmod_sparse *A, // m-by-n sparse matrix to factorize // outputs: cholmod_sparse **Q, // m-by-e sparse matrix cholmod_sparse **R, // e-by-n sparse matrix - int32_t **E, // size n column permutation, NULL if identity + int32_t **E, // size n column permutation, NULL if identity cholmod_common *cc // workspace and parameters ); diff --git a/SPQR/Include/spqr_cholmod_wrappers.hpp b/SPQR/Include/spqr_cholmod_wrappers.hpp index b20572b3e6..4daf71cd48 100644 --- a/SPQR/Include/spqr_cholmod_wrappers.hpp +++ b/SPQR/Include/spqr_cholmod_wrappers.hpp @@ -1,5 +1,35 @@ #include "spqr.hpp" +template int spqr_start +( + cholmod_common *Common +) ; + +template <> int spqr_start +( + cholmod_common *Common +) ; + +template <> int spqr_start +( + cholmod_common *Common +) ; + +template int spqr_finish +( + cholmod_common *Common +) ; + +template <> int spqr_finish +( + cholmod_common *Common +) ; + +template <> int spqr_finish +( + cholmod_common *Common +) ; + template void *spqr_malloc (size_t n, size_t size, cholmod_common *Common) ; template <> void *spqr_malloc (size_t n, size_t size, cholmod_common *Common) ; template <> void *spqr_malloc (size_t n, size_t size, cholmod_common *Common) ; @@ -506,3 +536,387 @@ template <> cholmod_sparse *spqr_speye /* --------------- */ cholmod_common *Common ) ; + + +template int spqr_free_work +( + cholmod_common *Common +) ; + +template <> int spqr_free_work +( + cholmod_common *Common +) ; + +template <> int spqr_free_work +( + cholmod_common *Common +) ; + +template cholmod_dense *spqr_zeros +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_dense *spqr_zeros +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_dense *spqr_zeros +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) ; + +template cholmod_dense *spqr_ones +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_dense *spqr_ones +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_dense *spqr_ones +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) ; + +template cholmod_sparse *spqr_ssmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* left matrix to multiply */ + cholmod_sparse *B, /* right matrix to multiply */ + int stype, /* requested stype of C */ + int values, /* TRUE: do numerical values, FALSE: pattern only */ + int sorted, /* if TRUE then return C with sorted columns */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_ssmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* left matrix to multiply */ + cholmod_sparse *B, /* right matrix to multiply */ + int stype, /* requested stype of C */ + int values, /* TRUE: do numerical values, FALSE: pattern only */ + int sorted, /* if TRUE then return C with sorted columns */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_ssmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* left matrix to multiply */ + cholmod_sparse *B, /* right matrix to multiply */ + int stype, /* requested stype of C */ + int values, /* TRUE: do numerical values, FALSE: pattern only */ + int sorted, /* if TRUE then return C with sorted columns */ + /* --------------- */ + cholmod_common *Common +) ; + + +template cholmod_sparse *spqr_ssadd +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to add */ + cholmod_sparse *B, /* matrix to add */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for B */ + int values, /* if TRUE compute the numerical values of C */ + int sorted, /* if TRUE, sort columns of C */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_ssadd +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to add */ + cholmod_sparse *B, /* matrix to add */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for B */ + int values, /* if TRUE compute the numerical values of C */ + int sorted, /* if TRUE, sort columns of C */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_ssadd +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to add */ + cholmod_sparse *B, /* matrix to add */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for B */ + int values, /* if TRUE compute the numerical values of C */ + int sorted, /* if TRUE, sort columns of C */ + /* --------------- */ + cholmod_common *Common +) ; + + +template int spqr_sdmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* sparse matrix to multiply */ + int transpose, /* use A if 0, or A' otherwise */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for Y */ + cholmod_dense *X, /* dense matrix to multiply */ + /* ---- in/out --- */ + cholmod_dense *Y, /* resulting dense matrix */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> int spqr_sdmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* sparse matrix to multiply */ + int transpose, /* use A if 0, or A' otherwise */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for Y */ + cholmod_dense *X, /* dense matrix to multiply */ + /* ---- in/out --- */ + cholmod_dense *Y, /* resulting dense matrix */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> int spqr_sdmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* sparse matrix to multiply */ + int transpose, /* use A if 0, or A' otherwise */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for Y */ + cholmod_dense *X, /* dense matrix to multiply */ + /* ---- in/out --- */ + cholmod_dense *Y, /* resulting dense matrix */ + /* --------------- */ + cholmod_common *Common +) ; + + + +template double spqr_norm_sparse +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> double spqr_norm_sparse +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> double spqr_norm_sparse +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm */ + /* --------------- */ + cholmod_common *Common +) ; + + +template double spqr_norm_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> double spqr_norm_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> double spqr_norm_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */ + /* --------------- */ + cholmod_common *Common +) ; + + +template cholmod_dense *spqr_copy_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to copy */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_dense *spqr_copy_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to copy */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_dense *spqr_copy_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to copy */ + /* --------------- */ + cholmod_common *Common +) ; + +template cholmod_sparse *spqr_copy +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to copy */ + int stype, /* requested stype of C */ + int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_copy +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to copy */ + int stype, /* requested stype of C */ + int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_copy +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to copy */ + int stype, /* requested stype of C */ + int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */ + /* --------------- */ + cholmod_common *Common +) ; + +template int spqr_sparse_xtype +( + /* ---- input ---- */ + int to_xtype, /* requested xtype (pattern, real, complex, zomplex) */ + /* ---- in/out --- */ + cholmod_sparse *A, /* sparse matrix to change */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> int spqr_sparse_xtype +( + /* ---- input ---- */ + int to_xtype, /* requested xtype (pattern, real, complex, zomplex) */ + /* ---- in/out --- */ + cholmod_sparse *A, /* sparse matrix to change */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> int spqr_sparse_xtype +( + /* ---- input ---- */ + int to_xtype, /* requested xtype (pattern, real, complex, zomplex) */ + /* ---- in/out --- */ + cholmod_sparse *A, /* sparse matrix to change */ + /* --------------- */ + cholmod_common *Common +) ; + +template int spqr_gpu_memorysize +( + size_t *total_mem, + size_t *available_mem, + cholmod_common *Common +) ; + +template <> int spqr_gpu_memorysize +( + size_t *total_mem, + size_t *available_mem, + cholmod_common *Common +) ; + +template <> int spqr_gpu_memorysize +( + size_t *total_mem, + size_t *available_mem, + cholmod_common *Common +) ; + + +template cholmod_sparse *spqr_read_sparse +( + /* ---- input ---- */ + FILE *f, /* file to read from, must already be open */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_read_sparse +( + /* ---- input ---- */ + FILE *f, /* file to read from, must already be open */ + /* --------------- */ + cholmod_common *Common +) ; + +template <> cholmod_sparse *spqr_read_sparse +( + /* ---- input ---- */ + FILE *f, /* file to read from, must already be open */ + /* --------------- */ + cholmod_common *Common +) ; + diff --git a/SPQR/Source/SuiteSparseQR_C.cpp b/SPQR/Source/SuiteSparseQR_C.cpp index c373e0e714..18f5932fe1 100644 --- a/SPQR/Source/SuiteSparseQR_C.cpp +++ b/SPQR/Source/SuiteSparseQR_C.cpp @@ -165,6 +165,7 @@ cholmod_dense *SuiteSparseQR_C_backslash // returns X, NULL if failure RETURN_IF_NULL (A, NULL) ; RETURN_IF_NULL (B, NULL) ; cc->status = CHOLMOD_OK ; + printf ("A->itype in C backslash is %d\n", A->itype) ; if (A->itype == CHOLMOD_INT) { return ((A->xtype == CHOLMOD_REAL) ? diff --git a/SPQR/Source/spqr_cholmod_wrappers.cpp b/SPQR/Source/spqr_cholmod_wrappers.cpp index 25a35842cd..37761c9b8b 100644 --- a/SPQR/Source/spqr_cholmod_wrappers.cpp +++ b/SPQR/Source/spqr_cholmod_wrappers.cpp @@ -1,4 +1,60 @@ +// ============================================================================= +// === spqr_cholmod_wrappers =================================================== +// ============================================================================= + +// SPQR, Copyright (c) 2008-2023, Timothy A Davis. All Rights Reserved. +// SPDX-License-Identifier: GPL-2.0+ + +//------------------------------------------------------------------------------ +// wrappers for calling CHOLMOD methods +//------------------------------------------------------------------------------ + #include "spqr.hpp" + +//------------------------------------------------------------------------------ +// spqr_start +//------------------------------------------------------------------------------ + +template <> int spqr_start +( + cholmod_common *Common +) +{ + return cholmod_start (Common) ; +} + +template <> int spqr_start +( + cholmod_common *Common +) +{ + return cholmod_l_start (Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_finish +//------------------------------------------------------------------------------ + +template <> int spqr_finish +( + cholmod_common *Common +) +{ + return cholmod_finish (Common) ; +} + +template <> int spqr_finish +( + cholmod_common *Common +) +{ + return cholmod_l_finish (Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_malloc +//------------------------------------------------------------------------------ + template <> void *spqr_malloc (size_t n, size_t size, cholmod_common *Common) { return cholmod_l_malloc (n, size, Common) ; @@ -8,6 +64,10 @@ template <> void *spqr_malloc (size_t n, size_t size, cholmod_common * return cholmod_malloc (n, size, Common) ; } +//------------------------------------------------------------------------------ +// spqr_calloc +//------------------------------------------------------------------------------ + template <> void *spqr_calloc (size_t n, size_t size, cholmod_common *Common) { return cholmod_l_calloc (n, size, Common) ; @@ -17,6 +77,10 @@ template <> void *spqr_calloc (size_t n, size_t size, cholmod_common * return cholmod_calloc (n, size, Common) ; } +//------------------------------------------------------------------------------ +// spqr_free +//------------------------------------------------------------------------------ + template <> void *spqr_free (size_t n, size_t size, void *p, cholmod_common *Common) { return cholmod_l_free (n, size, p, Common) ; @@ -26,6 +90,10 @@ template <> void *spqr_free (size_t n, size_t size, void *p, cholmod_c return cholmod_free (n, size, p, Common) ; } +//------------------------------------------------------------------------------ +// spqr_realloc +//------------------------------------------------------------------------------ + template <> void *spqr_realloc /* returns pointer to reallocated block */ ( /* ---- input ---- */ @@ -55,6 +123,10 @@ template <> void *spqr_realloc /* returns pointer to reallocated block return cholmod_l_realloc (nnew, size, p, n, Common) ; } +//------------------------------------------------------------------------------ +// spqr_allocate_sparse +//------------------------------------------------------------------------------ + template <> cholmod_sparse *spqr_allocate_sparse ( /* ---- input ---- */ @@ -88,6 +160,10 @@ template <> cholmod_sparse *spqr_allocate_sparse return cholmod_l_allocate_sparse (nrow, ncol, nzmax, sorted, packed, stype, xtype, Common) ; } +//------------------------------------------------------------------------------ +// spqr_free_sparse +//------------------------------------------------------------------------------ + template <> int spqr_free_sparse ( /* ---- in/out --- */ @@ -98,6 +174,7 @@ template <> int spqr_free_sparse { return cholmod_free_sparse (A, Common) ; } + template <> int spqr_free_sparse ( /* ---- in/out --- */ @@ -109,6 +186,10 @@ template <> int spqr_free_sparse return cholmod_l_free_sparse (A, Common) ; } +//------------------------------------------------------------------------------ +// spqr_reallocate_sparse +//------------------------------------------------------------------------------ + template <> int spqr_reallocate_sparse ( /* ---- input ---- */ @@ -134,6 +215,10 @@ template <> int spqr_reallocate_sparse return cholmod_l_reallocate_sparse (nznew, A, Common) ; } +//------------------------------------------------------------------------------ +// spqr_allocate_dense +//------------------------------------------------------------------------------ + template <> cholmod_dense *spqr_allocate_dense ( /* ---- input ---- */ @@ -162,6 +247,10 @@ template <> cholmod_dense *spqr_allocate_dense return cholmod_l_allocate_dense (nrow, ncol, d, xtype, Common) ; } +//------------------------------------------------------------------------------ +// spqr_free_dense +//------------------------------------------------------------------------------ + template <> int spqr_free_dense ( /* ---- in/out --- */ @@ -183,6 +272,10 @@ template <> int spqr_free_dense return cholmod_l_free_dense (X, Common) ; } +//------------------------------------------------------------------------------ +// spqr_allocate_factor +//------------------------------------------------------------------------------ + template <> cholmod_factor *spqr_allocate_factor ( /* ---- input ---- */ @@ -204,6 +297,10 @@ template <> cholmod_factor *spqr_allocate_factor return cholmod_l_allocate_factor (n, Common) ; } +//------------------------------------------------------------------------------ +// spqr_free_factor +//------------------------------------------------------------------------------ + template <> int spqr_free_factor ( /* ---- in/out --- */ @@ -225,6 +322,10 @@ template <> int spqr_free_factor return cholmod_l_free_factor (L, Common) ; } +//------------------------------------------------------------------------------ +// spqr_allocate_work +//------------------------------------------------------------------------------ + template <> int spqr_allocate_work ( /* ---- input ---- */ @@ -250,6 +351,10 @@ template <> int spqr_allocate_work return cholmod_l_allocate_work (nrow, iworksize, xworksize, Common) ; } +//------------------------------------------------------------------------------ +// spqr_amd +//------------------------------------------------------------------------------ + template <> int spqr_amd ( cholmod_sparse *A, int64_t *fset, size_t fsize, int64_t *Perm, cholmod_common *Common @@ -265,6 +370,10 @@ template <> int spqr_amd return cholmod_amd (A, fset, fsize, Perm, Common) ; } +//------------------------------------------------------------------------------ +// spqr_metis +//------------------------------------------------------------------------------ + template <> int spqr_metis ( cholmod_sparse *A, int64_t *fset, size_t fsize, int postorder, int64_t *Perm, cholmod_common *Common @@ -281,6 +390,10 @@ template <> int spqr_metis return cholmod_metis (A, fset, fsize, postorder, Perm, Common) ; } +//------------------------------------------------------------------------------ +// spqr_transpose +//------------------------------------------------------------------------------ + template <> cholmod_sparse *spqr_transpose ( cholmod_sparse *A, int values, cholmod_common *Common @@ -296,6 +409,10 @@ template <> cholmod_sparse *spqr_transpose return cholmod_transpose (A, values, Common) ; } +//------------------------------------------------------------------------------ +// spqr_analyze_p2 +//------------------------------------------------------------------------------ + template <> cholmod_factor *spqr_analyze_p2 ( @@ -331,6 +448,10 @@ cholmod_factor *spqr_analyze_p2 return cholmod_l_analyze_p2 (for_whom, A, UserPerm, fset, fsize, Common) ; } +//------------------------------------------------------------------------------ +// spqr_colamd +//------------------------------------------------------------------------------ + template <> int spqr_colamd ( /* ---- input ---- */ @@ -362,6 +483,10 @@ template <> int spqr_colamd return cholmod_l_colamd (A, fset, fsize, postorder, Perm, Common) ; } +//------------------------------------------------------------------------------ +// spqr_postorder +//------------------------------------------------------------------------------ + template <> int32_t spqr_postorder /* return # of nodes postordered */ ( /* ---- input ---- */ @@ -391,6 +516,10 @@ template <> int64_t spqr_postorder /* return # of nodes postordered */ return cholmod_l_postorder (Parent, n, Weight_p, Post, Common) ; } +//------------------------------------------------------------------------------ +// spqr_nnz +//------------------------------------------------------------------------------ + template <> int64_t spqr_nnz ( cholmod_sparse *A, @@ -408,6 +537,10 @@ template <> int64_t spqr_nnz return cholmod_nnz(A, Common) ; } +//------------------------------------------------------------------------------ +// spqr_dense_to_sparse +//------------------------------------------------------------------------------ + template <> cholmod_sparse *spqr_dense_to_sparse ( /* ---- input ---- */ @@ -431,6 +564,10 @@ template <> cholmod_sparse *spqr_dense_to_sparse return cholmod_l_dense_to_sparse(X, values, Common) ; } +//------------------------------------------------------------------------------ +// spqr_sparse_to_dense +//------------------------------------------------------------------------------ + template <> cholmod_dense *spqr_sparse_to_dense ( /* ---- input ---- */ @@ -452,6 +589,10 @@ template <> cholmod_dense *spqr_sparse_to_dense return cholmod_l_sparse_to_dense (A, Common) ; } +//------------------------------------------------------------------------------ +// spqr_speye +//------------------------------------------------------------------------------ + template <> cholmod_sparse *spqr_speye ( /* ---- input ---- */ @@ -476,3 +617,388 @@ template <> cholmod_sparse *spqr_speye { return cholmod_speye(nrow, ncol, xtype, Common) ; } + +//------------------------------------------------------------------------------ +// spqr_free_work +//------------------------------------------------------------------------------ + +template <> int spqr_free_work +( + cholmod_common *Common +) +{ + return cholmod_free_work (Common) ; +} + +template <> int spqr_free_work +( + cholmod_common *Common +) +{ + return cholmod_l_free_work (Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_zeros +//------------------------------------------------------------------------------ + +template <> cholmod_dense *spqr_zeros +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_zeros (nrow, ncol, xtype, Common) ; +} + +template <> cholmod_dense *spqr_zeros +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_zeros (nrow, ncol, xtype, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_ones +//------------------------------------------------------------------------------ + +template <> cholmod_dense *spqr_ones +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_ones (nrow, ncol, xtype, Common) ; +} + +template <> cholmod_dense *spqr_ones +( + /* ---- input ---- */ + size_t nrow, /* # of rows of matrix */ + size_t ncol, /* # of columns of matrix */ + int xtype, /* CHOLMOD_REAL, _COMPLEX, or _ZOMPLEX */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_ones (nrow, ncol, xtype, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_ssmult +//------------------------------------------------------------------------------ + +template <> cholmod_sparse *spqr_ssmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* left matrix to multiply */ + cholmod_sparse *B, /* right matrix to multiply */ + int stype, /* requested stype of C */ + int values, /* TRUE: do numerical values, FALSE: pattern only */ + int sorted, /* if TRUE then return C with sorted columns */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_ssmult (A, B, stype, values, sorted, Common) ; +} + +template <> cholmod_sparse *spqr_ssmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* left matrix to multiply */ + cholmod_sparse *B, /* right matrix to multiply */ + int stype, /* requested stype of C */ + int values, /* TRUE: do numerical values, FALSE: pattern only */ + int sorted, /* if TRUE then return C with sorted columns */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_ssmult (A, B, stype, values, sorted, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_ssadd +//------------------------------------------------------------------------------ + +template <> cholmod_sparse *spqr_ssadd +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to add */ + cholmod_sparse *B, /* matrix to add */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for B */ + int values, /* if TRUE compute the numerical values of C */ + int sorted, /* if TRUE, sort columns of C */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_add (A, B, alpha, beta, values, sorted, Common) ; +} + +template <> cholmod_sparse *spqr_ssadd +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to add */ + cholmod_sparse *B, /* matrix to add */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for B */ + int values, /* if TRUE compute the numerical values of C */ + int sorted, /* if TRUE, sort columns of C */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_add (A, B, alpha, beta, values, sorted, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_sdmult +//------------------------------------------------------------------------------ + +template <> int spqr_sdmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* sparse matrix to multiply */ + int transpose, /* use A if 0, or A' otherwise */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for Y */ + cholmod_dense *X, /* dense matrix to multiply */ + /* ---- in/out --- */ + cholmod_dense *Y, /* resulting dense matrix */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_sdmult (A, transpose, alpha, beta, X, Y, Common) ; +} + +template <> int spqr_sdmult +( + /* ---- input ---- */ + cholmod_sparse *A, /* sparse matrix to multiply */ + int transpose, /* use A if 0, or A' otherwise */ + double alpha [2], /* scale factor for A */ + double beta [2], /* scale factor for Y */ + cholmod_dense *X, /* dense matrix to multiply */ + /* ---- in/out --- */ + cholmod_dense *Y, /* resulting dense matrix */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_sdmult (A, transpose, alpha, beta, X, Y, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_norm_sparse +//------------------------------------------------------------------------------ + +template <> double spqr_norm_sparse +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_norm_sparse (A, norm, Common) ; +} + + +template <> double spqr_norm_sparse +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_norm_sparse (A, norm, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_norm_dense +//------------------------------------------------------------------------------ + +template <> double spqr_norm_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_norm_dense (X, norm, Common) ; +} + +template <> double spqr_norm_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to compute the norm of */ + int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_norm_dense (X, norm, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_copy_dense +//------------------------------------------------------------------------------ + +template <> cholmod_dense *spqr_copy_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to copy */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_copy_dense (X, Common) ; +} + +template <> cholmod_dense *spqr_copy_dense +( + /* ---- input ---- */ + cholmod_dense *X, /* matrix to copy */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_copy_dense (X, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_copy +//------------------------------------------------------------------------------ + +template <> cholmod_sparse *spqr_copy +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to copy */ + int stype, /* requested stype of C */ + int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_copy (A, stype, mode, Common) ; +} + +template <> cholmod_sparse *spqr_copy +( + /* ---- input ---- */ + cholmod_sparse *A, /* matrix to copy */ + int stype, /* requested stype of C */ + int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_copy (A, stype, mode, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_xtype +//------------------------------------------------------------------------------ + +template <> int spqr_sparse_xtype +( + /* ---- input ---- */ + int to_xtype, /* requested xtype (pattern, real, complex, zomplex) */ + /* ---- in/out --- */ + cholmod_sparse *A, /* sparse matrix to change */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_sparse_xtype (to_xtype, A, Common) ; +} + +template <> int spqr_sparse_xtype +( + /* ---- input ---- */ + int to_xtype, /* requested xtype (pattern, real, complex, zomplex) */ + /* ---- in/out --- */ + cholmod_sparse *A, /* sparse matrix to change */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_sparse_xtype (to_xtype, A, Common) ; +} + +//------------------------------------------------------------------------------ +// spqr_gpu_memorysize +//------------------------------------------------------------------------------ + +#ifdef SUITESPARSE_CUDA +template <> int spqr_gpu_memorysize +( + size_t *total_mem, + size_t *available_mem, + cholmod_common *Common +) +{ + return cholmod_gpu_memorysize (total_mem, available_mem, Common) ; +} + +template <> int spqr_gpu_memorysize +( + size_t *total_mem, + size_t *available_mem, + cholmod_common *Common +) +{ + return cholmod_l_gpu_memorysize (total_mem, available_mem, Common) ; +} +#endif + +//------------------------------------------------------------------------------ +// spqr_read_sparse +//------------------------------------------------------------------------------ + +template <> cholmod_sparse *spqr_read_sparse +( + /* ---- input ---- */ + FILE *f, /* file to read from, must already be open */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_read_sparse (f, Common) ; +} + +template <> cholmod_sparse *spqr_read_sparse +( + /* ---- input ---- */ + FILE *f, /* file to read from, must already be open */ + /* --------------- */ + cholmod_common *Common +) +{ + return cholmod_l_read_sparse (f, Common) ; +} + + diff --git a/SPQR/Tcov/Makefile b/SPQR/Tcov/Makefile index 19a01c3ec5..248fc86dd9 100644 --- a/SPQR/Tcov/Makefile +++ b/SPQR/Tcov/Makefile @@ -60,15 +60,16 @@ NVCC = /usr/local/cuda/bin/nvcc -g --profile --generate-line-info $(NV20) \ CLIB = -lsuitesparseconfig -lcholmod -lamd -lcolamd -lccolamd -lcamd -all: qrtest gpudemo qrdemo_gpu +all: qrtest qrtest32 gpudemo qrdemo_gpu -library: qrtest gpudemo +library: qrtest qrtest32 gpudemo purge: distclean distclean: clean - - $(RM) qrtest qrtest_out.txt pfile tfile cov.out qrtest_out1.txt + - $(RM) qrtest qrtest32 qrtest_out.txt pfile tfile cov.out - $(RM) gpuqrengine_demo troll.m qrdemo_gpu gpu_results.txt X.mtx + - $(RM) qrtest_out32.txt qrtest_out1.txt - $(RM) -r $(PURGE) clean: @@ -76,7 +77,7 @@ clean: INC = ../Include/spqr.hpp ../Include/SuiteSparseQR_C.h \ ../Include/SuiteSparseQR_definitions.h \ - ../Include/SuiteSparseQR.hpp + ../Include/SuiteSparseQR.hpp qrtest_template.hpp OBJ = \ spqr_rmap.o \ @@ -180,6 +181,9 @@ qrtestc.o: qrtestc.c $(INC) qrtest: qrtest.cpp $(INC) $(OBJ) $(C) qrtest.cpp -o qrtest $(OBJ) $(LIBS) -lm +qrtest32: qrtest.cpp $(INC) $(OBJ) + $(C) -DINT32=1 qrtest.cpp -o qrtest32 $(OBJ) $(LIBS) -lm + ifneq ($(GPU_CONFIG),) gpu: gpuqrengine_demo qrdemo_gpu - ./gpuqrengine_demo @@ -202,21 +206,24 @@ ifneq ($(GPU_CONFIG),) $(C) ../Demo/qrdemo_gpu.cpp -o qrdemo_gpu $(OBJ) $(LIBS) endif -go: qrtest gpu qrdemo_gpu +go: qrtest qrtest32 gpu qrdemo_gpu + - ./qrtest32 matrixlist.txt > qrtest_out32.txt - ./qrtest matrixlist.txt > qrtest_out.txt - ./cov -go1: qrtest +go1: qrtest qrtest32 - ./qrtest matrix1.txt > qrtest_out1.txt + - ./qrtest32 matrix1.txt > qrtest_out1.txt - ./cov -vgo1: qrtest +vgo1: qrtest qrtest32 - valgrind ./qrtest matrix1.txt > qrtest_out1.txt - # - valgrind --leak-check=full --show-reachable=yes ./qrtest matrix1.txt > qrtest_out1.txt + - valgrind ./qrtest32 matrix1.txt > qrtest_out1.txt - ./cov -vgo: qrtest - - valgrind --leak-check=full --show-reachable=yes ./qrtest matrixlist.txt > qrtest_out.txt +vgo: qrtest qrtest32 + - valgrind --leak-check=full --show-reachable=yes ./qrtest matrixlist.txt > qrtest_out32.txt + - valgrind --leak-check=full --show-reachable=yes ./qrtest32 matrixlist.txt > qrtest_out.txt - ./cov spqr_1colamd.o: ../Source/spqr_1colamd.cpp diff --git a/SPQR/Tcov/qrtest.cpp b/SPQR/Tcov/qrtest.cpp index 94805c7260..84f1afca58 100644 --- a/SPQR/Tcov/qrtest.cpp +++ b/SPQR/Tcov/qrtest.cpp @@ -121,7 +121,19 @@ void normal_memory_handler (cholmod_common *cc, bool free_work) SuiteSparse_config_free_func_set (free) ; cc->error_handler = my_handler ; - if (free_work) cholmod_l_free_work (cc) ; + + if (free_work) + { + if (cc->itype == CHOLMOD_LONG) + { + spqr_free_work (cc) ; + } + else + { + spqr_free_work (cc) ; + } + } + my_tries = -2 ; my_punt = FALSE ; } @@ -135,2643 +147,38 @@ void test_memory_handler (cholmod_common *cc, bool free_work) SuiteSparse_config_free_func_set (my_free) ; cc->error_handler = NULL ; - if (free_work) cholmod_l_free_work (cc) ; - my_tries = -2 ; - my_punt = FALSE ; -} - - -// ============================================================================= -// === SPQR_qmult ============================================================== -// ============================================================================= - -// wrapper for SuiteSparseQR_qmult (dense), optionally testing memory alloc. - -template cholmod_dense *SPQR_qmult -( - // arguments for SuiteSparseQR_qmult: - int64_t method, - cholmod_sparse *H, - cholmod_dense *Tau, - int64_t *HPinv, - cholmod_dense *X, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_dense *Y = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (Y) ; -} - - -// ============================================================================= -// === SPQR_qmult_sparse ======================================================= -// ============================================================================= - -// wrapper for SuiteSparseQR_qmult (sparse), optionally testing memory alloc. - -template cholmod_sparse *SPQR_qmult -( - // arguments for SuiteSparseQR_qmult: - int64_t method, - cholmod_sparse *H, - cholmod_dense *Tau, - int64_t *HPinv, - cholmod_sparse *X, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_sparse *Y = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc); - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (Y) ; -} - -#ifndef NEXPERT - -// ============================================================================= -// === SPQR_qmult (dense case) ================================================= -// ============================================================================= - -// wrapper for SuiteSparseQR_qmult (dense), optionally testing memory alloc. - -template cholmod_dense *SPQR_qmult -( - // arguments for SuiteSparseQR_qmult: - int64_t method, - SuiteSparseQR_factorization *QR, - cholmod_dense *X, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_dense *Y = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - Y = SuiteSparseQR_qmult (method, QR, X, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - Y = SuiteSparseQR_qmult (method, QR, X, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (Y) ; -} - -// ============================================================================= -// === SPQR_qmult (sparse case) ================================================ -// ============================================================================= - -// wrapper for SuiteSparseQR_qmult (sparse), optionally testing memory alloc. - -template cholmod_sparse *SPQR_qmult -( - // arguments for SuiteSparseQR_qmult: - int64_t method, - SuiteSparseQR_factorization *QR, - cholmod_sparse *X, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_sparse *Y = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - Y = SuiteSparseQR_qmult (method, QR, X, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - Y = SuiteSparseQR_qmult (method, QR, X, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (Y) ; -} - - -// ============================================================================= -// === SPQR_solve (dense case) ================================================= -// ============================================================================= - -// Wrapper for testing SuiteSparseQR_solve - -template cholmod_dense *SPQR_solve -( - // arguments for SuiteSparseQR_solve: - int64_t system, - SuiteSparseQR_factorization *QR, - cholmod_dense *B, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_dense *X = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - X = SuiteSparseQR_solve (system, QR, B, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - X = SuiteSparseQR_solve (system, QR, B, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (X) ; -} - - -// ============================================================================= -// === SPQR_solve (sparse case) ================================================ -// ============================================================================= - -// Wrapper for testing SuiteSparseQR_solve - -template cholmod_sparse *SPQR_solve -( - // arguments for SuiteSparseQR_solve: - int64_t system, - SuiteSparseQR_factorization *QR, - cholmod_sparse *B, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_sparse *X = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - X = SuiteSparseQR_solve (system, QR, B, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - X = SuiteSparseQR_solve (system, QR, B, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (X) ; -} - - -// ============================================================================= -// === SPQR_min2norm (dense case) ============================================== -// ============================================================================= - -// Wrapper for testing SuiteSparseQR_min2norm - -template cholmod_dense *SPQR_min2norm -( - // arguments for SuiteSparseQR_min2norm: - int ordering, - double tol, - cholmod_sparse *A, - cholmod_dense *B, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_dense *X = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (X) ; -} - -// ============================================================================= -// === SPQR_min2norm (sparse case) ============================================= -// ============================================================================= - -// Wrapper for testing SuiteSparseQR_min2norm - -template cholmod_sparse *SPQR_min2norm -( - // arguments for SuiteSparseQR_min2norm: - int ordering, - double tol, - cholmod_sparse *A, - cholmod_sparse *B, - cholmod_common *cc, - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - cholmod_sparse *X = NULL ; - if (!memory_test) - { - // just call the method directly; no memory testing - X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; - if (cc->status == CHOLMOD_OK) break ; - } - normal_memory_handler (cc, true) ; - } - return (X) ; -} - -// ============================================================================= -// === SPQR_factorize ========================================================== -// ============================================================================= - -// Wrapper for testing SuiteSparseQR_factorize or -// SuiteSparseQR_symbolic and SuiteSparseQR_numeric -template SuiteSparseQR_factorization *SPQR_factorize -( - // arguments for SuiteSparseQR_factorize: - int ordering, - double tol, - cholmod_sparse *A, - cholmod_common *cc, - - // method to use - int split, // if 1 use SuiteSparseQR_symbolic followed by - // SuiteSparseQR_numeric, if 0 use - // SuiteSparseQR_factorize, if 2, do the - // numeric factorization twice, just for testing. - // if 3 use SuiteSparseQR_C_factorize - // if 3 use SuiteSparseQR_C_symbolic / _C_numeric - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - SuiteSparseQR_factorization *QR ; - SuiteSparseQR_C_factorization *C_QR ; - if (!memory_test) + if (free_work) { - // just call the method directly; no memory testing - - if (split == 4) - { - - C_QR = SuiteSparseQR_C_symbolic (ordering, tol >= 0, A, cc) ; - SuiteSparseQR_C_numeric (tol, A, C_QR, cc) ; - if (C_QR == NULL) - { - cc->status = CHOLMOD_OUT_OF_MEMORY ; - QR = NULL ; - } - else - { - QR = (SuiteSparseQR_factorization *) - (C_QR->factors) ; - if (QR == NULL || QR->QRnum == NULL) - { - cc->status = CHOLMOD_OUT_OF_MEMORY ; - QR = NULL ; - } - else - { - // QR itself will be kept; free the C wrapper - C_QR->factors = NULL ; - cc->status = CHOLMOD_OK ; - if (QR == NULL) printf ("Hey!\n") ; - } - } - SuiteSparseQR_C_free (&C_QR, cc) ; - - } - else if (split == 3) - { - - C_QR = SuiteSparseQR_C_factorize (ordering, tol, A, cc) ; - int save = cc->status ; - if (C_QR == NULL) - { - QR = NULL ; - } - else - { - QR = (SuiteSparseQR_factorization *) (C_QR->factors) ; - C_QR->factors = NULL ; - SuiteSparseQR_C_free (&C_QR, cc) ; - } - cc->status = save ; - - } - else if (split == 2) - { - QR = SuiteSparseQR_symbolic (ordering, - tol >= 0 || tol <= SPQR_DEFAULT_TOL, A, cc) ; - ASSERT (QR != NULL) ; - SuiteSparseQR_numeric (tol, A, QR, cc) ; - // just for testing - SuiteSparseQR_numeric (tol, A, QR, cc) ; - } - else if (split == 1) + if (cc->itype == CHOLMOD_LONG) { - // split symbolic/numeric, no singletons exploited - QR = SuiteSparseQR_symbolic (ordering, - tol >= 0 || tol <= SPQR_DEFAULT_TOL, A, cc) ; - ASSERT (QR != NULL) ; - SuiteSparseQR_numeric (tol, A, QR, cc) ; + spqr_free_work (cc) ; } else { - QR = SuiteSparseQR_factorize (ordering, tol, A, cc) ; - } - - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - - if (split == 4) - { - - C_QR = SuiteSparseQR_C_symbolic (ordering, tol >= 0, A, cc) ; - SuiteSparseQR_C_numeric (tol, A, C_QR, cc) ; - if (C_QR == NULL) - { - cc->status = CHOLMOD_OUT_OF_MEMORY ; - QR = NULL ; - } - else - { - QR = (SuiteSparseQR_factorization *) - (C_QR->factors) ; - if (QR == NULL || QR->QRnum == NULL) - { - cc->status = CHOLMOD_OUT_OF_MEMORY ; - QR = NULL ; - } - else - { - // QR itself will be kept; free the C wrapper - C_QR->factors = NULL ; - cc->status = CHOLMOD_OK ; - if (QR == NULL) printf ("Hey!!\n") ; - } - } - SuiteSparseQR_C_free (&C_QR, cc) ; - - - } - else if (split == 3) - { - - C_QR = SuiteSparseQR_C_factorize (ordering, tol, A, cc) ; - int save = cc->status ; - if (C_QR == NULL) - { - QR = NULL ; - } - else - { - QR = (SuiteSparseQR_factorization *) C_QR->factors ; - C_QR->factors = NULL ; - SuiteSparseQR_C_free (&C_QR, cc) ; - } - cc->status = save ; - - } - else if (split == 1 || split == 2) - { - // split symbolic/numeric, no singletons exploited - QR = SuiteSparseQR_symbolic (ordering, - tol >= 0 || tol <= SPQR_DEFAULT_TOL, A, cc) ; - if (cc->status < CHOLMOD_OK) - { - continue ; - } - ASSERT (QR != NULL) ; - SuiteSparseQR_numeric (tol, A, QR, cc) ; - if (cc->status < CHOLMOD_OK) - { - SuiteSparseQR_free (&QR, cc) ; - } - } - else - { - QR = SuiteSparseQR_factorize (ordering, tol, A, cc) ; - } - - if (cc->status == CHOLMOD_OK) break ; + spqr_free_work (cc) ; } - normal_memory_handler (cc, true) ; } - if (QR == NULL) printf ("huh?? split: %d ordering %d tol %g\n", split, - ordering, tol) ; - return (QR) ; + my_tries = -2 ; + my_punt = FALSE ; } - -#endif - // ============================================================================= -// === SPQR_qr ================================================================= +// === qrest_template ========================================================== // ============================================================================= -// wrapper for SuiteSparseQR, optionally testing memory allocation - -template int64_t SPQR_qr -( - // arguments for SuiteSparseQR: - int ordering, - double tol, - int64_t econ, - int64_t getCTX, - cholmod_sparse *A, - cholmod_sparse *Bsparse, - cholmod_dense *Bdense, - cholmod_sparse **Zsparse, - cholmod_dense **Zdense, - cholmod_sparse **R, - int64_t **E, - cholmod_sparse **H, - int64_t **HPinv, - cholmod_dense **HTau, - cholmod_common *cc, - - // which version to use (C or C++) - int use_c_version, // if TRUE use C version, otherwise use C++ - - // malloc control - int64_t memory_test, // if TRUE, test malloc error handling - int64_t memory_punt // if TRUE, test punt case -) -{ - int64_t rank ; - if (!memory_test) - { - // just call the method directly; no memory testing - if (use_c_version) - { - rank = SuiteSparseQR_C (ordering, tol, econ, getCTX, A, - Bsparse, Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc) ; - } - else - { - rank = SuiteSparseQR (ordering, tol, econ, getCTX, A, - Bsparse, Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc) ; - } - } - else - { - // test malloc error handling - int64_t tries ; - test_memory_handler (cc, true) ; - my_punt = memory_punt ; - for (tries = 0 ; my_tries < 0 ; tries++) - { - my_tries = tries ; - if (use_c_version) - { - rank = SuiteSparseQR_C (ordering, tol, econ, getCTX, A, - Bsparse, Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc); - } - else - { - rank = SuiteSparseQR (ordering, tol, econ, getCTX, A, - Bsparse, Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc); - } - if (cc->status == CHOLMOD_OK) - { - break ; - } - } - normal_memory_handler (cc, true) ; - } - return (rank) ; -} - +#include "qrtest_template.hpp" // ============================================================================= -// === my_rand ================================================================= +// === qrtest main ============================================================= // ============================================================================= -// The POSIX example of rand, duplicated here so that the same sequence will -// be generated on different machines. - -static unsigned long next = 1 ; - -#define MY_RAND_MAX 32767 - -// RAND_MAX assumed to be 32767 -int64_t my_rand (void) -{ - next = next * 1103515245 + 12345 ; - return ((unsigned)(next/65536) % (MY_RAND_MAX + 1)) ; -} - -void my_srand (unsigned seed) -{ - next = seed ; -} - -unsigned long my_seed (void) -{ - return (next) ; -} - -int64_t nrand (int64_t n) // return a random int64_t between 0 and n-1 -{ - return ((n <= 0) ? 0 : (my_rand ( ) % n)) ; -} - -double xrand ( ) // return a random double between -1 and 1 -{ - double x = ((double) my_rand ( )) / MY_RAND_MAX ; - return (2*x-1) ; -} - -double erand (double range) -{ - return (range * xrand ( )) ; -} - -Complex erand (Complex range) -{ - /* - Complex x ; - x.real ( ) = xrand ( ) ; - x.imag ( ) = xrand ( ) ; - return (range * x) ; - */ - Complex i = Complex (0,1) ; - return (range * (xrand ( ) + i * xrand ( ))) ; -} - - -// ============================================================================= -// === getreal ================================================================= -// ============================================================================= - -// return real part of a scalar x - -inline double getreal (double x) -{ - return (x) ; -} - -inline double getreal (Complex x) -{ - return (x.real ( )) ; -} - -// ============================================================================= -// === getimag ================================================================= -// ============================================================================= - -// return imaginary part of a scalar x -inline double getimag (double x) // x is an unused parameter -{ - return (0) ; -} - -inline double getimag (Complex x) -{ - return (x.imag ( )) ; -} - - -// ============================================================================= -// === dense_wrapper =========================================================== -// ============================================================================= - -// place a column-oriented matrix in a cholmod_dense wrapper - -template cholmod_dense *dense_wrapper -( - cholmod_dense *X, - int64_t nrow, - int64_t ncol, - Entry *Xx -) -{ - X->xtype = spqr_type ( ) ; - X->nrow = nrow ; - X->ncol = ncol ; - X->d = nrow ; // leading dimension = nrow - X->nzmax = nrow * ncol ; - X->x = Xx ; - X->z = NULL ; // ZOMPLEX case not supported - X->dtype = CHOLMOD_DOUBLE ; - return (X) ; -} - -// ============================================================================= -// === sparse_split ============================================================ -// ============================================================================= - -// Return the real or imaginary part of a packed complex sparse matrix - -cholmod_sparse *sparse_split -( - cholmod_sparse *A, - int64_t part, - cholmod_common *cc -) -{ - cholmod_sparse *C ; - if (!A || A->xtype != CHOLMOD_COMPLEX || A->nz != NULL) return (NULL) ; - if (! (part == 0 || part == 1)) return (NULL) ; - - int64_t nz = cholmod_l_nnz (A, cc) ; - C = cholmod_l_allocate_sparse (A->nrow, A->ncol, nz, TRUE, TRUE, 0, - CHOLMOD_REAL, cc) ; - - int64_t *Ap = (int64_t *) A->p ; - int64_t *Ai = (int64_t *) A->i ; - double *Ax = (double *) A->x ; - - int64_t *Cp = (int64_t *) C->p ; - int64_t *Ci = (int64_t *) C->i ; - double *Cx = (double *) C->x ; - - int64_t n = A->ncol ; - - for (int64_t k = 0 ; k < n+1 ; k++) - { - Cp [k] = Ap [k] ; - } - - for (int64_t k = 0 ; k < nz ; k++) - { - Ci [k] = Ai [k] ; - } - - for (int64_t k = 0 ; k < nz ; k++) - { - Cx [k] = Ax [2*k + part] ; - } - - return (C) ; -} - -cholmod_sparse *sparse_real (cholmod_sparse *A, cholmod_common *cc) -{ - return (sparse_split (A, 0, cc)) ; -} - -cholmod_sparse *sparse_imag (cholmod_sparse *A, cholmod_common *cc) -{ - return (sparse_split (A, 1, cc)) ; -} - -// ============================================================================= -// === sparse_merge ============================================================ -// ============================================================================= - -// Add the real and imaginary parts of a matrix (both stored in real form) -// into a single matrix. The two parts must have the same nonzero pattern. -// A is CHOLMOD_REAL on input and holds the real part, it is CHOLMOD_COMPLEX -// on output. A_imag is CHOLMOD_REAL on input; it holds the imaginary part -// of A as a real matrix. - -int sparse_merge -( - cholmod_sparse *A, // input/output - cholmod_sparse *A_imag, // input only - cholmod_common *cc -) -{ - if (A == NULL || A_imag == NULL) - { - return (FALSE) ; - } - int64_t nz1 = cholmod_l_nnz (A, cc) ; - int64_t nz2 = cholmod_l_nnz (A_imag, cc) ; - if (A->xtype != CHOLMOD_REAL || A_imag->xtype != CHOLMOD_REAL || nz1 != nz2) - { - return (FALSE) ; - } - - // change A from real to complex - cholmod_l_sparse_xtype (CHOLMOD_COMPLEX, A, cc) ; - - double *Ax = (double *) A->x ; - double *Az = (double *) A_imag->x ; - - // merge in the imaginary part from A_imag into A - for (int64_t k = 0 ; k < nz1 ; k++) - { - Ax [2*k+1] = Az [k] ; - } - return (TRUE) ; -} - - -// ============================================================================= -// === sparse_diff ============================================================= -// ============================================================================= - -// Compute C = A-B where A and B are either both real, or both complex - -template cholmod_sparse *sparse_diff -( - cholmod_sparse *A, - cholmod_sparse *B, - cholmod_common *cc -) -{ - cholmod_sparse *C ; - double one [2] = {1,0}, minusone [2] = {-1,0} ; - - if (spqr_type ( ) == CHOLMOD_REAL) - { - // C = A - B - C = cholmod_l_add (A, B, one, minusone, TRUE, TRUE, cc) ; - } - else - { - cholmod_sparse *A_real, *A_imag, *B_real, *B_imag, *C_imag ; - - A_real = sparse_real (A, cc) ; - A_imag = sparse_imag (A, cc) ; - - B_real = sparse_real (B, cc) ; - B_imag = sparse_imag (B, cc) ; - - // real(C) = real(A) - real (B) - C = cholmod_l_add (A_real, B_real, one, minusone, TRUE, TRUE, cc) ; - - // imag(C) = imag(A) - imag (B) - C_imag = cholmod_l_add (A_imag, B_imag, one, minusone, TRUE, TRUE, cc) ; - - // C = real(C) + 1i*imag(C) - sparse_merge (C, C_imag, cc) ; - cholmod_l_free_sparse (&C_imag, cc) ; - - cholmod_l_free_sparse (&A_real, cc) ; - cholmod_l_free_sparse (&A_imag, cc) ; - cholmod_l_free_sparse (&B_real, cc) ; - cholmod_l_free_sparse (&B_imag, cc) ; - } - - return (C) ; -} - - -// ============================================================================= -// === permute_columns ========================================================= -// ============================================================================= - -// compute A(:,P) - -template cholmod_sparse *permute_columns -( - cholmod_sparse *A, - int64_t *P, - cholmod_common *cc -) -{ - int64_t m = A->nrow ; - int64_t n = A->ncol ; - int64_t nz = cholmod_l_nnz (A, cc) ; - int64_t xtype = spqr_type ( ) ; - int64_t *Ap = (int64_t *) A->p ; - int64_t *Ai = (int64_t *) A->i ; - Entry *Ax = (Entry *) A->x ; - cholmod_sparse *C ; - - // allocate empty matrix C with space for nz entries - C = cholmod_l_allocate_sparse (m, n, nz, TRUE, TRUE, 0, xtype, cc) ; - int64_t *Cp = (int64_t *) C->p ; - int64_t *Ci = (int64_t *) C->i ; - Entry *Cx = (Entry *) C->x ; - - // construct column pointers for C - for (int64_t k = 0 ; k < n ; k++) - { - // column j of A becomes column k of C - int64_t j = P ? P [k] : k ; - Cp [k] = Ap [j+1] - Ap [j] ; - } - spqr_cumsum (n, Cp) ; - - // copy columns from A to C - for (int64_t k = 0 ; k < n ; k++) - { - // copy column k of A into column j of C - int64_t j = P ? P [k] : k ; - int64_t pdest = Cp [k] ; - int64_t psrc = Ap [j] ; - int64_t len = Ap [j+1] - Ap [j] ; - for (int64_t t = 0 ; t < len ; t++) - { - Ci [pdest + t] = Ai [psrc + t] ; - Cx [pdest + t] = Ax [psrc + t] ; - } - } - - return (C) ; -} - - -// ============================================================================= -// === sparse_multiply ========================================================= -// ============================================================================= - -// compute A*B where A and B can be both real or both complex - -template cholmod_sparse *sparse_multiply -( - cholmod_sparse *A, - cholmod_sparse *B, - cholmod_common *cc -) -{ - cholmod_sparse *C ; - - if (spqr_type ( ) == CHOLMOD_REAL) - { - // C = A*B - C = cholmod_l_ssmult (A, B, 0, TRUE, TRUE, cc) ; - } - else - { - // cholmod_ssmult and cholmod_add only work for real matrices - cholmod_sparse *A_real, *A_imag, *B_real, *B_imag, *C_imag, *T1, *T2 ; - double one [2] = {1,0}, minusone [2] = {-1,0} ; - - A_real = sparse_real (A, cc) ; - A_imag = sparse_imag (A, cc) ; - - B_real = sparse_real (B, cc) ; - B_imag = sparse_imag (B, cc) ; - - // real(C) = real(A)*real(B) - imag(A)*imag(B) - T1 = cholmod_l_ssmult (A_real, B_real, 0, TRUE, TRUE, cc) ; - T2 = cholmod_l_ssmult (A_imag, B_imag, 0, TRUE, TRUE, cc) ; - C = cholmod_l_add (T1, T2, one, minusone, TRUE, TRUE, cc) ; - cholmod_l_free_sparse (&T1, cc) ; - cholmod_l_free_sparse (&T2, cc) ; - - // imag(C) = imag(A)*real(B) + real(A)*imag(B) - T1 = cholmod_l_ssmult (A_imag, B_real, 0, TRUE, TRUE, cc) ; - T2 = cholmod_l_ssmult (A_real, B_imag, 0, TRUE, TRUE, cc) ; - C_imag = cholmod_l_add (T1, T2, one, one, TRUE, TRUE, cc) ; - cholmod_l_free_sparse (&T1, cc) ; - cholmod_l_free_sparse (&T2, cc) ; - - // C = real(C) + 1i*imag(C) - sparse_merge (C, C_imag, cc) ; - cholmod_l_free_sparse (&C_imag, cc) ; - - cholmod_l_free_sparse (&A_real, cc) ; - cholmod_l_free_sparse (&A_imag, cc) ; - cholmod_l_free_sparse (&B_real, cc) ; - cholmod_l_free_sparse (&B_imag, cc) ; - } - - return (C) ; -} - - -// ============================================================================= -// === sparse_resid ============================================================ -// ============================================================================= - -// compute norm (A*x-b,1) for A,x, and b all sparse - -template double sparse_resid -( - cholmod_sparse *A, - double anorm, - cholmod_sparse *X, - cholmod_sparse *B, - cholmod_common *cc -) -{ - cholmod_sparse *AX, *Resid ; - // AX = A*X - AX = sparse_multiply (A, X, cc) ; - // Resid = AX - B - Resid = sparse_diff (AX, B, cc) ; - // resid = norm (Resid,1) - double resid = cholmod_l_norm_sparse (Resid, 1, cc) ; - resid = CHECK_NAN (resid) ; - cholmod_l_free_sparse (&AX, cc) ; - cholmod_l_free_sparse (&Resid, cc) ; - return (CHECK_NAN (resid / anorm)) ; -} - - -// ============================================================================= -// === dense_resid ============================================================= -// ============================================================================= - -// compute norm (A*x-b,1) for A sparse, x and b dense - -template double dense_resid -( - cholmod_sparse *A, - double anorm, - cholmod_dense *X, - int64_t nb, - Entry *Bx, - cholmod_common *cc -) -{ - cholmod_dense *B, Bmatrix, *Resid ; - double one [2] = {1,0}, minusone [2] = {-1,0} ; - - B = dense_wrapper (&Bmatrix, A->nrow, nb, Bx) ; - // Resid = B - Resid = cholmod_l_copy_dense (B, cc) ; - - // Resid = A*X - Resid - cholmod_l_sdmult (A, FALSE, one, minusone, X, Resid, cc) ; - - // resid = norm (Resid,1) - double resid = cholmod_l_norm_dense (Resid, 1, cc) ; - resid = CHECK_NAN (resid) ; - cholmod_l_free_dense (&Resid, cc) ; - return (CHECK_NAN (resid / anorm)) ; -} - -// ============================================================================= -// === check_r_factor ========================================================== -// ============================================================================= - -// compute norm (R'*R - (A(:,P))'*(A(:,P)), 1) / norm (A'*A,1) - -template double check_r_factor -( - cholmod_sparse *R, - cholmod_sparse *A, - int64_t *P, - cholmod_common *cc -) -{ - cholmod_sparse *RTR, *RT, *C, *CT, *CTC, *D ; - - // RTR = R'*R - RT = cholmod_l_transpose (R, 2, cc) ; - RTR = sparse_multiply (RT, R, cc) ; - cholmod_l_free_sparse (&RT, cc) ; - - // C = A(:,P) - C = permute_columns (A, P, cc) ; - - // CTC = C'*C - CT = cholmod_l_transpose (C, 2, cc) ; - CTC = sparse_multiply (CT, C, cc) ; - cholmod_l_free_sparse (&CT, cc) ; - cholmod_l_free_sparse (&C, cc) ; - - double ctcnorm = cholmod_l_norm_sparse (CTC, 1, cc) ; - if (ctcnorm == 0) ctcnorm = 1 ; - ctcnorm = CHECK_NAN (ctcnorm) ; - - // D = RTR - CTC - D = sparse_diff (RTR, CTC, cc) ; - cholmod_l_free_sparse (&CTC, cc) ; - cholmod_l_free_sparse (&RTR, cc) ; - - // err = norm (D,1) - double err = cholmod_l_norm_sparse (D, 1, cc) / ctcnorm ; - err = CHECK_NAN (err) ; - - cholmod_l_free_sparse (&D, cc) ; - return (CHECK_NAN (err)) ; -} - - -// ============================================================================= -// === check_qr ================================================================ -// ============================================================================= - -// compute norm (Q*R - A(:,P)) / norm (A) - -template double check_qr -( - cholmod_sparse *Q, - cholmod_sparse *R, - cholmod_sparse *A, - int64_t *P, - double anorm, - cholmod_common *cc -) -{ - cholmod_sparse *QR, *C, *D ; - - // C = A(:,P) - C = permute_columns (A, P, cc) ; - - // QR = Q*R - QR = sparse_multiply (Q, R, cc) ; - - // D = RTR - CTC - D = sparse_diff (QR, C, cc) ; - cholmod_l_free_sparse (&QR, cc) ; - cholmod_l_free_sparse (&C, cc) ; - - // err = norm (D,1) - double err = cholmod_l_norm_sparse (D, 1, cc) / anorm ; - err = CHECK_NAN (err) ; - - cholmod_l_free_sparse (&D, cc) ; - return (CHECK_NAN (err)) ; -} - - -// ============================================================================= -// === Rsolve ================================================================== -// ============================================================================= - -template int Rsolve -( - // R is n-by-n, upper triangular with zero-free diagonal - int64_t n, - cholmod_sparse *R, - Entry *X, // X is n-by-nx, leading dimension n, overwritten with soln - int64_t nx, - cholmod_common *cc -) -{ - // int64_t n = R->n ; - int64_t *Rp = (int64_t *) R->p ; - int64_t *Ri = (int64_t *) R->i ; - Entry *Rx = (Entry *) R->x ; - - // check the diagonal - for (int64_t j = 0 ; j < n ; j++) - { - if (Rp [j] == Rp [j+1] || Ri [Rp [j+1]-1] != j) - { - printf ("Rsolve: R not upper triangular w/ zero-free diagonal\n") ; - return (FALSE) ; - } - } - - // do the backsolve - for (int64_t k = 0 ; k < nx ; k++) - { - for (int64_t j = n-1 ; j >= 0 ; j--) - { - Entry rjj = Rx [Rp [j+1]-1] ; - if (rjj == (Entry) 0) - { - printf ("Rsolve: R has an explicit zero on the diagonal\n") ; - return (FALSE) ; - } - X [j] /= rjj ; - for (int64_t p = Rp [j] ; p < Rp [j+1]-1 ; p++) - { - X [Ri [p]] -= Rx [p] * X [j] ; - } - } - X += n ; - } - - return (TRUE) ; -} - - -// ============================================================================= -// === create_Q ================================================================ -// ============================================================================= - -// create a sparse Q -template cholmod_sparse *create_Q -( - cholmod_sparse *H, - cholmod_dense *HTau, - int64_t *HPinv, - cholmod_common *cc -) -{ - cholmod_sparse *Q, *I ; - int64_t m = H->nrow ; - int64_t xtype = spqr_type ( ) ; - I = cholmod_l_speye (m, m, xtype, cc) ; - Q = SPQR_qmult (1, H, HTau, HPinv, I, cc, m<300, nrand (2)) ; - cholmod_l_free_sparse (&I, cc) ; - return (Q) ; -} - - -// ============================================================================= -// === QRsolve ================================================================= -// ============================================================================= - -// solve Ax=b using H, R, and E - -template double QRsolve -( - cholmod_sparse *A, - double anorm, - int64_t rank, - int64_t method, - cholmod_sparse *H, - cholmod_dense *HTau, - int64_t *HPinv, - cholmod_sparse *R, - int64_t *Qfill, - cholmod_dense *Bdense, - cholmod_common *cc -) -{ - double one [2] = {1,0}, zero [2] = {0,0}, resid = EMPTY ; - int64_t xtype = spqr_type ( ) ; - int64_t n = A->ncol ; - int64_t m = A->nrow ; - int64_t nrhs = Bdense->ncol ; - Entry *X, *Y = NULL, *B ; - cholmod_dense *Ydense = NULL ; - cholmod_dense *Xdense ; - - B = (Entry *) Bdense->x ; - Xdense = cholmod_l_zeros (n, nrhs, xtype, cc) ; - X = (Entry *) Xdense->x ; - - if (method == 0) - { - // solve Ax=b using H, R, and E - // Y = Q'*B - Ydense = SPQR_qmult (0, H, HTau, HPinv, Bdense, cc, - m < 300, nrand (2)) ; - } - else - { - cholmod_sparse *Q, *QT ; - // solve using Q instead of qmult - Q = create_Q (H, HTau, HPinv, cc) ; - QT = cholmod_l_transpose (Q, 2, cc) ; - // Y = Q'*B - Ydense = cholmod_l_zeros (m, nrhs, xtype, cc) ; - cholmod_l_sdmult (QT, FALSE, one, zero, Bdense, Ydense, cc) ; - cholmod_l_free_sparse (&Q, cc) ; - cholmod_l_free_sparse (&QT, cc) ; - } - - // Y (1:rank) = R (1:rank,1:rank) \ Y (1:rank) - Y = (Entry *) Ydense->x ; - int64_t ok = Rsolve (rank, R, Y, nrhs, cc) ; - // X = E*Y - if (ok) - { - for (int64_t kk = 0 ; kk < nrhs ; kk++) - { - for (int64_t k = 0 ; k < rank ; k++) - { - int64_t j = Qfill ? Qfill [k] : k ; - X [j] = Y [k] ; - } - X += n ; - Y += m ; - } - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, Xdense, nrhs, B, cc) ; - } - - cholmod_l_free_dense (&Ydense, cc) ; - cholmod_l_free_dense (&Xdense, cc) ; - return (CHECK_NAN (resid)) ; -} - - -// ============================================================================= -// === check_qmult ============================================================= -// ============================================================================= - -// Test qmult - -template double check_qmult -( - cholmod_sparse *H, - cholmod_dense *HTau, - int64_t *HPinv, - int64_t test_errors, - cholmod_common *cc -) -{ - cholmod_sparse *Q, *QT, *Xsparse, *Ssparse, *Zsparse ; - cholmod_dense *Xdense, *Zdense, *Sdense, *Ydense ; - int64_t xtype = spqr_type ( ) ; - Entry *X, *Y, *Z, *S ; - double err, maxerr = 0 ; - double one [2] = {1,0}, zero [2] = {0,0} ; - Entry range = (Entry) 1.0 ; - int64_t k ; - - int64_t m = H->nrow ; - Q = create_Q (H, HTau, HPinv, cc) ; // construct Q from H - - QT = cholmod_l_transpose (Q, 2, cc) ; // QT = Q' - - // compare Q with qmult for sparse and dense X - for (int64_t nx = 0 ; nx < 5 ; nx++) // # of columns of X - { - int64_t xsize = m * nx ; // size of X - for (int64_t nz = 1 ; nz <= xsize+1 ; nz *= 16) // # of nonzeros in X - { - - // ----------------------------------------------------------------- - // create X as m-by-nx, both sparse and dense - // ----------------------------------------------------------------- - - Xdense = cholmod_l_zeros (m, nx, xtype, cc) ; - X = (Entry *) Xdense->x ; - for (k = 0 ; k < nz ; k++) - { - X [nrand (xsize)] += erand (range) ; - } - Xsparse = cholmod_l_dense_to_sparse (Xdense, TRUE, cc) ; - - // ----------------------------------------------------------------- - // Y = Q'*X for method 0, Y = Q*X for method 1 - // ----------------------------------------------------------------- - - for (int method = 0 ; method <= 1 ; method++) - { - // Y = Q'*X or Q*X - Ydense = SPQR_qmult (method, H, HTau, HPinv, Xdense, cc, - m < 300, nrand (2)) ; - Y = (Entry *) Ydense->x ; - - Zdense = cholmod_l_zeros (m, nx, xtype, cc) ; - cholmod_l_sdmult (Q, (method == 0), one, zero, - Xdense, Zdense, cc) ; - Z = (Entry *) Zdense->x ; - - for (err = 0, k = 0 ; k < xsize ; k++) - { - double e1 = spqr_abs (Y [k] - Z [k], cc) ; - e1 = CHECK_NAN (e1) ; - err = MAX (err, e1) ; - } - maxerr = MAX (maxerr, err) ; - - // S = Q'*Xsparse or Q*Xsparse - Ssparse = SPQR_qmult ( - method, H, HTau, HPinv, Xsparse, cc, m < 300, nrand (2)) ; - Sdense = cholmod_l_sparse_to_dense (Ssparse, cc) ; - S = (Entry *) Sdense->x ; - - for (err = 0, k = 0 ; k < xsize ; k++) - { - double e1 = spqr_abs (S [k] - Y [k], cc) ; - e1 = CHECK_NAN (e1) ; - err = MAX (err, e1) ; - } - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_dense (&Ydense, cc) ; - cholmod_l_free_dense (&Zdense, cc) ; - cholmod_l_free_sparse (&Ssparse, cc) ; - cholmod_l_free_dense (&Sdense, cc) ; - } - - cholmod_l_free_dense (&Xdense, cc) ; - cholmod_l_free_sparse (&Xsparse, cc) ; - - // ----------------------------------------------------------------- - // create X as nx-by-m, both sparse and dense - // ----------------------------------------------------------------- - - Xdense = cholmod_l_zeros (nx, m, xtype, cc) ; - X = (Entry *) Xdense->x ; - for (k = 0 ; k < nz ; k++) - { - X [nrand (xsize)] += erand (range) ; - } - Xsparse = cholmod_l_dense_to_sparse (Xdense, TRUE, cc) ; - - // ----------------------------------------------------------------- - // Y = X*Q' for method 2, Y = X*Q for method 3 - // ----------------------------------------------------------------- - - for (int method = 2 ; method <= 3 ; method++) - { - // Y = X*Q' or X*Q - Ydense = SPQR_qmult (method, H, HTau, HPinv, Xdense, cc, - m < 300, nrand (2)) ; - Y = (Entry *) Ydense->x ; - - if (method == 2) - { - // Zsparse = (X*Q') - Zsparse = sparse_multiply (Xsparse, QT, cc) ; - } - else - { - // Zsparse = (X*Q) - Zsparse = sparse_multiply (Xsparse, Q, cc) ; - } - Zdense = cholmod_l_sparse_to_dense (Zsparse, cc) ; - - Z = (Entry *) Zdense->x ; - - for (err = 0, k = 0 ; k < xsize ; k++) - { - double e1 = spqr_abs (Y [k] - Z [k], cc) ; - e1 = CHECK_NAN (e1) ; - err = MAX (err, e1) ; - } - maxerr = MAX (maxerr, err) ; - - // S = X*Q' or X*Q - Ssparse = SPQR_qmult ( - method, H, HTau, HPinv, Xsparse, cc, m < 300, nrand (2)) ; - Sdense = cholmod_l_sparse_to_dense (Ssparse, cc) ; - S = (Entry *) Sdense->x ; - - for (err = 0, k = 0 ; k < xsize ; k++) - { - double e1 = spqr_abs (S [k] - Y [k], cc) ; - e1 = CHECK_NAN (e1) ; - err = MAX (err, e1) ; - } - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_dense (&Ydense, cc) ; - cholmod_l_free_dense (&Zdense, cc) ; - cholmod_l_free_sparse (&Ssparse, cc) ; - cholmod_l_free_sparse (&Zsparse, cc) ; - cholmod_l_free_dense (&Sdense, cc) ; - } - - cholmod_l_free_dense (&Xdense, cc) ; - cholmod_l_free_sparse (&Xsparse, cc) ; - - } - } - cholmod_l_free_sparse (&Q, cc) ; - cholmod_l_free_sparse (&QT, cc) ; - - // ------------------------------------------------------------------------- - // qmult error conditions - // ------------------------------------------------------------------------- - - // These should fail; expect 6 error messages in the output - - if (test_errors) - { - err = 0 ; - printf ("The following six errors are expected:\n") ; - Xdense = cholmod_l_zeros (m+1, 1, xtype,cc) ; - Ydense = SuiteSparseQR_qmult (0, H, HTau, HPinv, Xdense, cc) ; - if (Ydense != NULL || cc->status != CHOLMOD_INVALID) err++ ; - cholmod_l_free_dense (&Xdense, cc) ; - - Xdense = cholmod_l_zeros (1, m+1, xtype,cc) ; - Ydense = SuiteSparseQR_qmult (2, H, HTau, HPinv, Xdense, cc) ; - if (Ydense != NULL || cc->status != CHOLMOD_INVALID) err++ ; - Ydense = SuiteSparseQR_qmult (42, H, HTau, HPinv, Xdense, cc) ; - if (Ydense != NULL || cc->status != CHOLMOD_INVALID) err++ ; - cholmod_l_free_dense (&Xdense, cc) ; - - Xsparse = cholmod_l_speye (m+1, m+1, xtype, cc) ; - - Q = SuiteSparseQR_qmult (0, H, HTau, HPinv, Xsparse, cc); - if (Q != NULL || cc->status != CHOLMOD_INVALID) err++ ; - Q = SuiteSparseQR_qmult (2, H, HTau, HPinv, Xsparse, cc); - if (Q != NULL || cc->status != CHOLMOD_INVALID) err++ ; - Q = SuiteSparseQR_qmult (9, H, HTau, HPinv, Xsparse, cc); - if (Q != NULL || cc->status != CHOLMOD_INVALID) err++ ; - - printf (" ... error handling done\n\n") ; - cholmod_l_free_sparse (&Xsparse, cc) ; - - maxerr = MAX (maxerr, err) ; - } - - return (CHECK_NAN (maxerr)) ; -} - - -// ============================================================================= -// === check_rc ================================================================ -// ============================================================================= - -// X = Q'*B has been done, continue with C=R\C and - -template double check_rc -( - int64_t rank, - cholmod_sparse *R, - cholmod_sparse *A, - Entry *B, - cholmod_dense *X, - int64_t nrhs, - double anorm, - int64_t *Qfill, - cholmod_common *cc -) -{ - double resid = EMPTY ; - int64_t xtype = spqr_type ( ) ; - cholmod_dense *W ; - int64_t n, ok ; - Entry *W1, *X1 ; - if (!R || !X) - { - ok = 0 ; - } - else - { - n = X->nrow ; - X1 = (Entry *) X->x ; - // solve X = R\X, overwriting X with solution - ok = Rsolve (rank, R, X1, nrhs, cc) ; - } - if (ok) - { - // W = E*X - // W = (Entry *) cholmod_l_calloc (A->ncol * nrhs, sizeof (Entry), cc) ; - W = cholmod_l_zeros (A->ncol, nrhs, xtype, cc) ; - W1 = (Entry *) W->x ; - for (int64_t col = 0 ; col < nrhs ; col++) - { - for (int64_t k = 0 ; k < rank ; k++) - { - int64_t j = Qfill ? Qfill [k] : k ; - if (j < (int64_t) A->ncol) W1 [j] = X1 [k] ; - } - W1 += A->ncol ; - X1 += n ; - } - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, W, nrhs, B, cc) ; - // cholmod_l_free (A->ncol * nrhs, sizeof (Entry), W, cc) ; - cholmod_l_free_dense (&W, cc) ; - } - return (CHECK_NAN (resid)) ; -} - - -// ============================================================================= -// === transpose =============================================================== -// ============================================================================= - -// Transpose a dense matrix. - -template cholmod_dense *transpose -( - cholmod_dense *Xdense, - cholmod_common *cc -) -{ - Entry *X, *Y ; - cholmod_dense *Ydense ; - if (Xdense == NULL) - { - printf ("transpose failed!\n") ; - return (NULL) ; - } - int64_t m = Xdense->nrow ; - int64_t n = Xdense->ncol ; - int64_t ldx = Xdense->d ; - int64_t xtype = spqr_type ( ) ; - Ydense = cholmod_l_allocate_dense (n, m, n, xtype, cc) ; - X = (Entry *) Xdense->x ; - Y = (Entry *) Ydense->x ; - for (int64_t i = 0 ; i < m ; i++) - { - for (int64_t j = 0 ; j < n ; j++) - { - Y [j+i*n] = spqr_conj (X [i+j*ldx]) ; - } - } - return (Ydense) ; -} - -// ============================================================================= -// === qrtest ================================================================== -// ============================================================================= - -template void qrtest -( - cholmod_sparse *A, - double errs [5], - cholmod_common *cc -) -{ - cholmod_sparse *H, *I, *R, *Q, *Csparse, *Xsparse, *AT, *Bsparse ; - cholmod_dense *Cdense, *Xdense, *Bdense, *HTau ; ; - double tol = DBL_EPSILON, err, resid, maxerr, maxresid [2][2] ; - double tols [ ] = { SPQR_DEFAULT_TOL, -1, 0, DBL_EPSILON } ; - int64_t n, m, nz, *HPinv, ntol, *Ai, *Ap, k, *Qfill, rank, nb, *Cp, *Ci, econ, - which ; - Entry *B, *Ax, *Cx ; - int64_t xtype = spqr_type ( ) ; - int ordering ; - Entry range = (Entry) 1.0 ; - - errs [0] = EMPTY ; - errs [1] = EMPTY ; - errs [2] = EMPTY ; - errs [3] = EMPTY ; - errs [4] = EMPTY ; - if (A == NULL) - { - fprintf (stderr, "qrtest: no input matrix\n") ; - return ; - } - - m = A->nrow ; - n = A->ncol ; - Ap = (int64_t *) A->p ; - Ai = (int64_t *) A->i ; - Ax = (Entry *) A->x ; - double anorm = cholmod_l_norm_sparse (A, 1, cc) ; - anorm = CHECK_NAN (anorm) ; - printf ("\n===========================================================\n") ; - printf ("Matrix: %ld by %ld, nnz(A) = %ld, norm(A,1) = %g\n", - m, n, cholmod_l_nnz (A, cc), anorm) ; - printf ( "===========================================================\n") ; - - if (anorm == 0) anorm = 1 ; - - // these should all be zero, no matter what the matrix - maxerr = 0 ; - - // residuals for well-determined or under-determined systems. If - // not rank deficient, these should be zero - maxresid [0][0] = 0 ; // for m <= n, default tol, ordering not fixed - maxresid [0][1] = 0 ; // for m <= n, all other cases - - // residuals for least-squares systems (these will not be zero) - maxresid [1][0] = 0 ; // for m > n, default tol, ordering not fixed - maxresid [1][1] = 0 ; // for m > n, all other cases - - my_srand (m+n) ; - - if (m > 45000) - { - // The only thing returned are the statistics (rank, etc) - printf ("\n=== QR huge (expect int overflow):\n") ; - rank = SPQR_qr ( - 0, 0, 0, 0, A, - NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, - cc, FALSE, FALSE, FALSE) ; - return ; - } - - if (MAX (m,n) >= 300) fprintf (stderr, "please wait ... ") ; - - AT = cholmod_l_transpose (A, 2, cc) ; - - // ------------------------------------------------------------------------- - // test QR routines - // ------------------------------------------------------------------------- - - econ = m ; - - for (ordering = 0 ; ordering <= 9 ; ordering++) - { - - // skip SPQR_ORDERING_GIVEN, unless the matrix is tiny - if (ordering == SPQR_ORDERING_GIVEN && MAX (m,n) > 10) continue ; - - for (ntol = 0 ; ntol < NTOL ; ntol++) - { - - tol = tols [ntol] ; - if // (ntol == 0) // this old test is fragile ... - (tol == SPQR_DEFAULT_TOL) // use this instead - { - // with default tolerance, the fixed ordering can sometimes - // fail if the matrix is rank deficient (R cannot be permuted - // to upper trapezoidal form). - which = (ordering == 0) ; - } - else - { - // with non-default tolerance, the solution can sometimes be - // poor; this is expected. - which = 1 ; - } - printf ("\n=== QR with ordering %d tol %g:\n", ordering, tol) ; - - // ----------------------------------------------------------------- - // create dense and sparse right-hand-sides - // ----------------------------------------------------------------- - - nb = 5 ; - Bdense = cholmod_l_zeros (m, nb, xtype, cc) ; - B = (Entry *) Bdense->x ; - for (k = 0 ; k < m*nb ; k++) - { - B [k] = ((double) (my_rand ( ) % 2)) * erand (range) ; - } - Bsparse = cholmod_l_dense_to_sparse (Bdense, TRUE, cc) ; - - // ----------------------------------------------------------------- - // X = qrsolve(A,B) where X and B are dense - // ----------------------------------------------------------------- - - // X = A\B - if (ordering == SPQR_ORDERING_DEFAULT && tol == SPQR_DEFAULT_TOL) - { - printf ("[ backslach, A and B and X dense: defaults\n") ; - Xdense = SuiteSparseQR (A, Bdense, cc) ; - printf ("] done backslach, A and B and X dense: defaults\n") ; - } - else - { - printf ("[ backslach, A and B and X dense: tol %g order %d\n", - tol, ordering) ; - Xdense = SuiteSparseQR (ordering, tol, A, Bdense, cc) ; - printf ("] done backslach, A B X dense: tol %g order %d\n", - tol, ordering) ; - } - - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid0b %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_dense (&Xdense, cc) ; - - if (cc->useGPU) - { - // error testing for infeasible GPU memory - int64_t save = cc->gpuMemorySize ; - cc->gpuMemorySize = 1 ; - printf ("[ Pretend GPU memory is too small:\n") ; - Xdense = SuiteSparseQR (ordering, tol, A, Bdense, cc) ; - cc->gpuMemorySize = save ; - printf ("] test done infeasible GPU, status %2d, useGPU: %d\n", - cc->status, cc->useGPU) ; - cholmod_l_free_dense (&Xdense, cc) ; - } - - cholmod_l_free_dense (&Bdense, cc) ; - - // ----------------------------------------------------------------- - // X = qrsolve(A,B) where X and B are sparse - // ----------------------------------------------------------------- - - // X = A\B - printf ("[ backslash with sparse B: tol %g\n", tol) ; - Xsparse = SuiteSparseQR (ordering, tol, A, Bsparse, cc) ; - printf ("] did tol %g\n", tol) ; - - // check norm (A*x-b), x and b sparse - resid = sparse_resid (A, anorm, Xsparse, Bsparse, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid0 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_sparse (&Xsparse, cc) ; - cholmod_l_free_sparse (&Bsparse, cc) ; - - // ----------------------------------------------------------------- - // X = qrsolve (A,B) where X and B are sparse, with memory test - // ----------------------------------------------------------------- - - // use B = A and solve AX=B where X is sparse - cc->SPQR_shrink = 0 ; // shrink = 0 ; - rank = SPQR_qr ( - ordering, tol, econ, 2, A, - A, NULL, &Xsparse, NULL, NULL, &Qfill, NULL, NULL, NULL, - cc, TRUE, m < 300, nrand (2)) ; - cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; - - if // (ntol == 0) // old test is fragile ... - (tol == SPQR_DEFAULT_TOL) // use this instead. - { - printf ("using default tol: %g\n", cc->SPQR_tol_used) ; - } - - // check norm (A*x-b), x and b sparse - resid = sparse_resid (A, anorm, Xsparse, A, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid1 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_sparse (&Xsparse, cc) ; - cholmod_l_free (n+n, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // X = qrsolve (A,B) where X and B are dense, with memory test - // ----------------------------------------------------------------- - - // use B = dense m-by-2 matrix with some zeros - nb = 5 ; - Bdense = cholmod_l_zeros (m, nb, xtype, cc) ; - B = (Entry *) Bdense->x ; - for (k = 0 ; k < m*nb ; k++) - { - B [k] = (k+1) % 7 ; - } - - rank = SPQR_qr ( - ordering, tol, econ, 2, A, - NULL, Bdense, NULL, &Xdense, NULL, &Qfill, NULL, NULL, NULL, - cc, FALSE, m < 300, nrand (2)) ; - - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid2 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_dense (&Xdense, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // X = qrsolve (A,B) where X and B are full and H is kept - // ----------------------------------------------------------------- - - cc->SPQR_shrink = 2 ; // shrink = 2 ; - rank = SPQR_qr ( - ordering, tol, econ, 2, A, - NULL, Bdense, NULL, &Xdense, NULL, &Qfill, &H, &HPinv, &HTau, - cc, FALSE, m < 300, nrand (2)) ; - cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; - - cholmod_l_free_dense (&HTau, cc) ; - cholmod_l_free (m, sizeof (int64_t), HPinv, cc) ; - cholmod_l_free_sparse (&H, cc) ; - - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid3 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_dense (&Xdense, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [C,R,E] = qr (A,B) where C is sparse and B is full - // ----------------------------------------------------------------- - - cc->SPQR_shrink = 2 ; // shrink = 2 ; - rank = SPQR_qr ( - ordering, tol, econ, 0, A, - NULL, Bdense, &Csparse, NULL, &R, &Qfill, NULL, NULL, NULL, - cc, FALSE, m < 300, nrand (2)) ; - cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; - - // compute x=R\C and check norm (A*x-b) - Cdense = cholmod_l_sparse_to_dense (Csparse, cc) ; - resid = check_rc (rank, R, A, B, Cdense, nb, anorm, Qfill, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid4 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - cholmod_l_free_dense (&Cdense, cc) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err1: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&Csparse, cc) ; - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [C,R,E] = qr (A,B) where C and B are full - // ----------------------------------------------------------------- - - cc->SPQR_shrink = 0 ; // shrink = 0 ; - rank = SPQR_qr ( - ordering, tol, econ, 0, A, - NULL, Bdense, NULL, &Cdense, &R, &Qfill, NULL, NULL, NULL, - cc, FALSE, m < 300, nrand (2)) ; - cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err2: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_dense (&Cdense, cc) ; - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [C,R,E] = qr (A,B) where C and B are full, simple wrapper - // ----------------------------------------------------------------- - - SuiteSparseQR (ordering, tol, econ, A, Bdense, - &Cdense, &R, &Qfill, cc) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err3: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_dense (&Cdense, cc) ; - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [C,R,E] = qr (A,B) where C and B are sparse, simple wrapper - // ----------------------------------------------------------------- - - Bsparse = cholmod_l_dense_to_sparse (Bdense, TRUE, cc) ; - - SuiteSparseQR (ordering, tol, econ, A, Bsparse, - &Csparse, &R, &Qfill, cc) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err4: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&Csparse, cc) ; - cholmod_l_free_sparse (&Bsparse, cc) ; - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [CT,R,E] = qr (A,B), but do not return R - // ----------------------------------------------------------------- - - rank = SPQR_qr ( - ordering, tol, econ, 1, A, - NULL, Bdense, &Csparse, NULL, NULL, &Qfill, NULL, NULL, NULL, - cc, FALSE, m < 300, nrand (2)) ; - - cholmod_l_free_sparse (&Csparse, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [Q,R,E] = qr (A), Q in Householder form - // ----------------------------------------------------------------- - - rank = SPQR_qr ( - ordering, tol, econ, -1, A, - NULL, NULL, NULL, NULL, &R, &Qfill, &H, &HPinv, &HTau, - cc, FALSE, m < 300, nrand (2)) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err5: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - // solve Ax=b using Householder form - resid = QRsolve (A, anorm, rank, 0, H, HTau, HPinv, R, - Qfill, Bdense, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid5 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - // solve Ax=b using Q matrix form - resid = QRsolve (A, anorm, rank, 1, H, HTau, HPinv, R, - Qfill, Bdense, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid6 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_dense (&HTau, cc) ; - cholmod_l_free_sparse (&H, cc) ; - cholmod_l_free (m, sizeof (int64_t), HPinv, cc) ; - cholmod_l_free (n, sizeof (int64_t), Qfill, cc) ; - cholmod_l_free_sparse (&R, cc) ; - - // ----------------------------------------------------------------- - // [Q,R,E] = qr (A), non-economy - // ----------------------------------------------------------------- - - cc->SPQR_shrink = 0 ; // shrink = 0 ; - I = cholmod_l_speye (m, m, xtype, cc) ; - rank = SPQR_qr ( - ordering, tol, m, 1, A, - I, NULL, &Q, NULL, &R, &Qfill, NULL, NULL, NULL, - cc, FALSE, m<300, nrand (2)) ; - cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; - - // ensure norm (Q*R - A*E) is small - err = check_qr (Q, R, A, Qfill, anorm, cc) ; - printf ("order %d : Q*R-A*E Err6: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&I, cc) ; - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free_sparse (&Q, cc) ; - cholmod_l_free (n+m, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [Q,R,E] = qr (A), non-economy, using simple wrapper - // ----------------------------------------------------------------- - - if (nrand (2)) - { - // use C version - SuiteSparseQR_C_QR (ordering, tol, m, A, &Q, &R, &Qfill, cc) ; - } - else - { - // use C++ version - SuiteSparseQR (ordering, tol, m, A, &Q, &R, &Qfill, cc); - } - - // ensure norm (Q*R - A*E) is small - err = check_qr (Q, R, A, Qfill, anorm, cc) ; - printf ("order %d : Q*R-A*E Err7: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free_sparse (&Q, cc) ; - cholmod_l_free (n+m, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [R,E] = qr (A) - // ----------------------------------------------------------------- - - rank = SPQR_qr ( - ordering, tol, econ, 0, A, - NULL, NULL, NULL, NULL, &R, &Qfill, NULL, NULL, NULL, - cc, FALSE, m < 300, nrand (2)) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err8: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - int64_t rank1 = rank ; - - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free (n, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [R,E] = qr (A) using simple wrapper - // ----------------------------------------------------------------- - - SuiteSparseQR (ordering, tol, econ, A, &R, &Qfill, cc) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err9: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free (n, sizeof (int64_t), Qfill, cc) ; - - // ----------------------------------------------------------------- - // [ ] = qr (A) - // ----------------------------------------------------------------- - - // The only thing returned are the statistics (rank, etc) - cc->SPQR_shrink = 0 ; // shrink = 0 ; - rank = SPQR_qr ( - ordering, tol, econ, 0, A, - NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, - cc, FALSE, m < 300, nrand (2)) ; - cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; - - err = (rank != rank1) ; - printf ("order %d : rank %6ld %6ld Err10: %g\n", ordering, - rank, rank1, err) ; - maxerr = MAX (maxerr, err) ; - - // ----------------------------------------------------------------- - // [C,H,R,E] = qr (A) - // ----------------------------------------------------------------- - - rank = SPQR_qr ( - ordering, tol, econ, 0, A, - NULL, Bdense, &Csparse, NULL, &R, &Qfill, &H, &HPinv, &HTau, - cc, FALSE, m < 300, nrand (2)) ; - - // compute x=R\C and check norm (A*x-b) - Cdense = cholmod_l_sparse_to_dense (Csparse, cc) ; - resid = check_rc (rank, R, A, B, Cdense, nb, anorm, Qfill, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("Resid7 %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - cholmod_l_free_dense (&Cdense, cc) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err11: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&Csparse, cc) ; - cholmod_l_free_sparse (&R, cc) ; - - // compare Q with qmult - err = check_qmult (H, HTau, HPinv, - ordering == 2 && /* ntol == 0 */ tol == SPQR_DEFAULT_TOL, cc) ; - printf ("order %d : check qmult Err12: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_dense (&HTau, cc) ; - cholmod_l_free_sparse (&H, cc) ; - cholmod_l_free (m, sizeof (int64_t), HPinv, cc) ; - cholmod_l_free (n+nb, sizeof (int64_t), Qfill, cc) ; - cholmod_l_free_dense (&Bdense, cc) ; - - // ----------------------------------------------------------------- - // [H,R,E] = qr (A), simple wrapper - // ----------------------------------------------------------------- - - SuiteSparseQR (ordering, tol, econ, A, - &R, &Qfill, &H, &HPinv, &HTau, cc) ; - - // check that R'*R = (A*E)'*(A*E) - err = check_r_factor (R, A, Qfill, cc) ; - printf ("order %d : R'R-(A*E)'*(A*E), Err13: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - // compare Q with qmult - err = check_qmult (H, HTau, HPinv, FALSE, cc) ; - printf ("order %d : check qmult Err14: %g\n", ordering, err) ; - maxerr = MAX (maxerr, err) ; - - cholmod_l_free_sparse (&R, cc) ; - cholmod_l_free_dense (&HTau, cc) ; - cholmod_l_free_sparse (&H, cc) ; - cholmod_l_free (m, sizeof (int64_t), HPinv, cc) ; - cholmod_l_free (n, sizeof (int64_t), Qfill, cc) ; - -#ifndef NEXPERT - - // ================================================================= - // === expert routines ============================================= - // ================================================================= - - SuiteSparseQR_factorization *QR ; - cholmod_dense *XT, *Zdense, *BT, *Ydense ; - cholmod_sparse *Ysparse ; - - // ----------------------------------------------------------------- - // QR = qr (A), then solve - // ----------------------------------------------------------------- - - for (int split = 0 ; split <= 4 ; split++) - { - - QR = SPQR_factorize (ordering, tol, A, cc, - split, m < 300, nrand (2)) ; - - // split == 4 does not use singletons, so it can fail if - // rank < n - int wh = which || (split == 4) ; - - // solve Ax=b - nb = 5 ; - Bdense = cholmod_l_zeros (m, nb, xtype, cc) ; - B = (Entry *) Bdense->x ; - for (k = 0 ; k < m*nb ; k++) - { - B [k] = erand (range) ; - } - - // Y = Q'*B - Ydense = SPQR_qmult (SPQR_QTX, QR, Bdense, cc, - m < 300, nrand (2)) ; - - // X = R\(E*Y) - Xdense = SPQR_solve (SPQR_RETX_EQUALS_B, QR, Ydense, cc, - m < 300, nrand (2)) ; - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; - maxresid [m>n][wh] = MAX (maxresid [m>n][wh], resid) ; - printf ("Resid8_%d %d %ld %d : %g (%d) tol %g\n", - split, m>n, ntol, ordering, resid, wh, tol) ; - cholmod_l_free_dense (&Xdense, cc) ; - cholmod_l_free_dense (&Ydense, cc) ; - - // Y = (B'*Q)' - BT = transpose (Bdense, cc) ; - XT = SPQR_qmult (SPQR_XQ, QR, BT, cc, - m < 300, nrand (2)) ; - - Ydense = transpose (XT, cc) ; - cholmod_l_free_dense (&XT, cc) ; - cholmod_l_free_dense (&BT, cc) ; - - // X = R\(E*Y) - Xdense = SPQR_solve (SPQR_RETX_EQUALS_B, QR, Ydense, cc, - m < 300, nrand (2)) ; - // check norm (A*x-b), x and b dense - resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; - maxresid [m>n][wh] = MAX (maxresid [m>n][wh], resid) ; - printf ("Resid9_%d %d %ld %d : %g (%d) tol %g\n", - split, m>n, ntol, ordering, resid, wh, tol) ; - cholmod_l_free_dense (&Xdense, cc) ; - cholmod_l_free_dense (&Ydense, cc) ; - - // ------------------------------------------------------------- - // error testing - // ------------------------------------------------------------- - - if (ordering == 0 && /* ntol == 0 */ tol == SPQR_DEFAULT_TOL) - { - printf ("Error handling ... expect 3 error messages: \n") ; - err = (SuiteSparseQR_qmult (-1, QR, Bdense, cc) - != NULL) ; - cholmod_l_free_dense (&Bdense, cc) ; - Bdense = cholmod_l_zeros (m+1, 1, xtype, cc) ; - err += (SuiteSparseQR_qmult (SPQR_QX,QR,Bdense,cc) - != NULL); - cholmod_l_free_dense (&Bdense, cc) ; - Bdense = cholmod_l_zeros (1, m+1, xtype, cc) ; - err += (SuiteSparseQR_qmult (SPQR_XQ,QR,Bdense,cc) - != NULL); - if (QR->n1cols > 0) - { - // this will fail; cannot refactorize with singletons - printf ("Error handling ... expect error message:\n") ; - err += (SuiteSparseQR_numeric (tol, A, QR, cc) - != FALSE) ; - } - printf ("order %d : error handling Err15: %g\n", - ordering, err) ; - maxerr = MAX (maxerr, err) ; - printf (" ... error handling done\n\n") ; - } - - cholmod_l_free_dense (&Bdense, cc) ; - - // ------------------------------------------------------------- - - // solve A'x=b - nb = 5 ; - Bdense = cholmod_l_zeros (n, nb, xtype, cc) ; - B = (Entry *) Bdense->x ; - for (k = 0 ; k < n*nb ; k++) - { - B [k] = erand (range) ; - } - // Y = R'\(E'*B) - Ydense = SPQR_solve (SPQR_RTX_EQUALS_ETB, QR, Bdense, cc, - m < 300, nrand (2)) ; - // X = Q*Y - Xdense = SPQR_qmult (SPQR_QX, QR, Ydense, cc, - m < 300, nrand (2)) ; - // check norm (A'*x-b), x and b dense - resid = dense_resid (AT, anorm, Xdense, nb, B, cc) ; - maxresid [m (-1, QR, Bdense, cc) - != NULL) ; - cholmod_l_free_dense (&Bdense, cc) ; - err += (SuiteSparseQR_solve (SPQR_RTX_EQUALS_ETB, - QR, Bdense, cc) != NULL) ; - Bdense = cholmod_l_zeros (n+1, 1, xtype, cc) ; - err += (SuiteSparseQR_solve (SPQR_RTX_EQUALS_ETB, - QR, Bdense, cc) != NULL) ; - printf ("order %d : error handling Err16: %g\n", - ordering, err) ; - maxerr = MAX (maxerr, err) ; - printf (" ... error handling done\n\n") ; - } - - SuiteSparseQR_free (&QR, cc) ; - cholmod_l_free_dense (&Bdense, cc) ; - } - - // ----------------------------------------------------------------- - // QR = qr (A'), then solve - // ----------------------------------------------------------------- - - // use qmult to solve min-2-norm problem - QR = SuiteSparseQR_factorize (ordering, tol, AT, cc) ; - - nb = 5 ; - Bdense = cholmod_l_zeros (m, nb, xtype, cc) ; - B = (Entry *) Bdense->x ; - for (k = 0 ; k < m*nb ; k++) - { - B [k] = erand (range) ; - } - - // solve X = R'\B - Xdense = SPQR_solve (SPQR_RTX_EQUALS_B, QR, Bdense, cc, - m < 300, nrand (2)) ; - cholmod_l_free_dense (&Xdense, cc) ; - - // solve X = R'\(E'*B) - Xdense = SPQR_solve (SPQR_RTX_EQUALS_ETB, QR, Bdense, cc, - m < 300, nrand (2)) ; - - // Y = Q*X - Ydense = SPQR_qmult (SPQR_QX, QR, Xdense, cc, - m < 300, nrand (2)) ; - - // check norm (A*y-b), y and b dense - resid = dense_resid (A, anorm, Ydense, nb, B, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("ResidB %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - cholmod_l_free_dense (&Ydense, cc) ; - - // Y = (X'*Q')' - XT = transpose (Xdense, cc) ; - Zdense = SPQR_qmult (SPQR_XQT, QR, XT, cc, - m < 300, nrand (2)) ; - Ydense = transpose (Zdense, cc) ; - cholmod_l_free_dense (&XT, cc) ; - cholmod_l_free_dense (&Zdense, cc) ; - - // check norm (A*y-b), y and b dense - resid = dense_resid (A, anorm, Ydense, nb, B, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("ResidC %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - cholmod_l_free_dense (&Ydense, cc) ; - cholmod_l_free_dense (&Xdense, cc) ; - - // ----------------------------------------------------------------- - // min 2-norm solution using min2norm - // ----------------------------------------------------------------- - - Xdense = SPQR_min2norm (ordering, tol, A, Bdense, cc, - m < 300, nrand (2)) ; - - // check norm (A*x-b), y and b dense - resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("ResidD %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - cholmod_l_free_dense (&Xdense, cc) ; - - cholmod_l_free_dense (&Bdense, cc) ; - - // ----------------------------------------------------------------- - // sparse case - // ----------------------------------------------------------------- - - nb = 5 ; - Bdense = cholmod_l_zeros (m, nb, xtype, cc) ; - B = (Entry *) Bdense->x ; - for (k = 0 ; k < m*nb ; k++) - { - B [k] = ((double) (my_rand ( ) % 2)) * erand (range) ; - } - Bsparse = cholmod_l_dense_to_sparse (Bdense, TRUE, cc) ; - cholmod_l_free_dense (&Bdense, cc) ; - - // solve X = R'\(E'*B) - Xsparse = NULL ; - Xsparse = SPQR_solve (SPQR_RTX_EQUALS_ETB, QR, Bsparse, cc, - m < 300, nrand (2)) ; - - // Y = Q*X - Ysparse = NULL ; - Ysparse = SPQR_qmult (SPQR_QX, QR, Xsparse, cc, - m < 300, nrand (2)) ; - - // check norm (A*y-b), y and b sparse - resid = sparse_resid (A, anorm, Ysparse, Bsparse, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("ResidE %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_sparse (&Xsparse, cc) ; - cholmod_l_free_sparse (&Ysparse, cc) ; - - Xsparse = SPQR_min2norm (ordering, tol, A, Bsparse, cc, - m < 300, nrand (2)) ; - - // check norm (A*x-b), x and b sparse - resid = sparse_resid (A, anorm, Xsparse, Bsparse, cc) ; - maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; - printf ("ResidF %d %ld %d : %g\n", m>n, ntol, ordering, resid) ; - - cholmod_l_free_sparse (&Xsparse, cc) ; - cholmod_l_free_sparse (&Bsparse, cc) ; - - SuiteSparseQR_free (&QR, cc) ; +#ifdef INT32 +#define Int int32_t +#else +#define Int int64_t #endif - } - } - - // ------------------------------------------------------------------------- - // check error handling - // ------------------------------------------------------------------------- - - printf ("Check error handling, one error message is expected:\n") ; - cholmod_dense *Bgunk = cholmod_l_ones (m+1, 1, xtype, cc) ; - rank = SuiteSparseQR ( - 0, 0, econ, -1, A, - NULL, Bgunk, NULL, NULL, NULL, NULL, NULL, NULL, NULL, - cc) ; - cholmod_l_free_dense (&Bgunk, cc) ; - err = (rank != EMPTY) ; - maxerr = MAX (maxerr, err) ; // rank should be EMPTY - printf (" ... error handling done\n\n") ; - - // ------------------------------------------------------------------------- - // test non-user callable functions - // ------------------------------------------------------------------------- - - // attempt to permute A to upper triangular form - int64_t *Qtrap ; - rank = spqr_trapezoidal (n, Ap, Ai, Ax, 0, NULL, FALSE, &Cp, &Ci, &Cx, - &Qtrap, cc) ; - printf ("Rank of A, if A*P permutable to upper trapezoidal: %ld\n", rank) ; - if (Cp != NULL) - { - nz = Cp [n] ; - cholmod_l_free (n+1, sizeof (int64_t), Cp, cc) ; - cholmod_l_free (nz, sizeof (int64_t), Ci, cc) ; - cholmod_l_free (nz, sizeof (Entry), Cx, cc) ; - cholmod_l_free (n, sizeof (int64_t), Qtrap, cc) ; - } - cholmod_l_free_sparse (&AT, cc) ; - - // ------------------------------------------------------------------------- - // test the C API - // ------------------------------------------------------------------------- - - qrtest_C (A, anorm, errs, maxresid, cc) ; - - // ------------------------------------------------------------------------- - // final results - // ------------------------------------------------------------------------- - - errs [0] = CHECK_NAN (maxerr) ; - errs [1] = CHECK_NAN (maxresid [0][0]) ; - errs [2] = CHECK_NAN (maxresid [0][1]) ; - errs [3] = CHECK_NAN (maxresid [1][0]) ; - errs [4] = CHECK_NAN (maxresid [1][1]) ; -} - - -// ============================================================================= -// === do_matrix =============================================================== -// ============================================================================= - -// Read in a matrix, and use it to test SuiteSparseQR -// If kind == 0, then the first two residuals should be low. - -int do_matrix2 (int kind, cholmod_sparse *A, cholmod_common *cc) ; - -int do_matrix (int kind, FILE *file, cholmod_common *cc) -{ - cholmod_sparse *A ; - - int nfail0 = 0 ; - int nfail1 = 0 ; - int nfail2 = 0 ; - int nfail3 = 0 ; - - // ------------------------------------------------------------------------- - // read in the matrix - // ------------------------------------------------------------------------- - - A = cholmod_l_read_sparse (file, cc) ; - if (A == NULL) - { - fprintf (stderr, "Unable to read matrix\n") ; - return (1) ; - } - int64_t m = A->nrow ; - int64_t n = A->ncol ; - fprintf (stderr, "%5ld by %5ld : ", m, n) ; - if (sizeof (int64_t) > sizeof (int) && (m > 10000 || n > 10000)) - { - fprintf (stderr, "(test skipped on 64-bit systems)\n") ; - cholmod_l_free_sparse (&A, cc) ; - return (0) ; - } - - // defaults - cc->SPQR_grain = 1 ; // no parallel analysis - printf ("\nBeginning CPU tests [\n") ; - fprintf (stderr, " CPU ") ; - nfail0 = do_matrix2 (kind, A, cc) ; - - // non-defaults (will not use the GPU) - cc->SPQR_grain = 4 ; // grain size relative to total work - nfail2 = do_matrix2 (kind, A, cc) ; - cc->SPQR_grain = 1 ; // no parallel analysis - printf ("\nCPU tests done ]\n") ; - - // test the GPU, if installed - #ifdef SUITESPARSE_CUDA - cc->useGPU = TRUE ; - // was 3.5 * ((size_t) 1024 * 1024 * 1024) ; - size_t totmem, availmem ; - double t = SuiteSparse_time ( ) ; - cholmod_l_gpu_memorysize (&totmem, &availmem, cc) ; - t = SuiteSparse_time ( ) - t ; - cc->gpuMemorySize = availmem ; - printf ("\nBeginning GPU tests, GPU memory %g MB warmup time %g[\n", - (double) (cc->gpuMemorySize) / (1024*1024), t) ; - fprintf (stderr, " GPU ") ; - nfail1 = do_matrix2 (kind, A, cc) ; - printf ("\nGPU tests done ]\n") ; - if (m > 200) - { - // try with a tiny GPU memory size, but only for a few matrices - // in the test set. Each front will go in its own stage. - printf ("\nBeginning GPU tests with tiny GPU memory [\n") ; - cc->gpuMemorySize = 0 ; - nfail3 = do_matrix2 (kind, A, cc) ; - // restore defaults - cc->useGPU = FALSE ; - printf ("\nGPU tests done (tiny memory) ]\n") ; - } - #endif - - cholmod_l_free_sparse (&A, cc) ; - - printf ("\n") ; - fprintf (stderr, "\n") ; - return (nfail0 + nfail1 + nfail2 + nfail3) ; -} - - - -int do_matrix2 (int kind, cholmod_sparse *A, cholmod_common *cc) -{ - double errs [5] = {0,0,0,0,0} ; - int64_t m = A->nrow ; - int64_t n = A->ncol ; - - // ------------------------------------------------------------------------- - // use it to test SuiteSparseQR - // ------------------------------------------------------------------------- - - if (A->xtype == CHOLMOD_COMPLEX && A->stype == 0) - { - qrtest (A, errs, cc) ; - } - else if (A->xtype == CHOLMOD_REAL) - { - if (A->stype != 0) - { - cholmod_sparse *A1 ; - A1 = cholmod_l_copy (A, 0, 1, cc) ; - qrtest (A1, errs, cc) ; - cholmod_l_free_sparse (&A1, cc) ; - } - else - { - qrtest (A, errs, cc) ; - } - } - else - { - // cannot handle ZOMPLEX, PATTERN, or symmetric/Hermitian COMPLEX - fprintf (stderr, "invalid matrix\n") ; - errs [0] = 1 ; - } - - // ------------------------------------------------------------------------- - // report the results - // ------------------------------------------------------------------------- - - if (kind == 0) - { - printf ("First Resid and ") ; - } - printf ("Err should be low:\n") ; - printf ("RESULT: Err %8.1e Resid %8.1e %8.1e", errs [0], - errs [1], errs [2]) ; - if (m == n) - { - printf (" ") ; - } - else - { - printf (" %8.1e %8.1e", errs [3], errs [4]) ; - } - - if (errs [0] > 1e-10) - { - printf (" : FAIL\n") ; - fprintf (stderr, "Error: %g FAIL\n", errs [0]) ; - return (1) ; - } - - // if kind == 0, then this full-rank matrix should have low residual - if (kind == 0 && (errs [1] > 1e-10)) - { - printf (" : FAIL\n") ; - fprintf (stderr, "error: %g FAIL\n", errs [1]) ; - return (1) ; - } - - printf (" : OK.") ; - fprintf (stderr, "OK.") ; - return (0) ; -} - - -// ============================================================================= -// === qrtest main ============================================================= -// ============================================================================= #define LEN 200 @@ -2786,7 +193,7 @@ int main (int argc, char **argv) // ------------------------------------------------------------------------- cc = &Common ; - cholmod_l_start (cc) ; + spqr_start (cc) ; normal_memory_handler (cc, true) ; if (argc == 1) @@ -2796,7 +203,7 @@ int main (int argc, char **argv) // Usage: qrtest < input.mtx // --------------------------------------------------------------------- - nfail += do_matrix (1, stdin, cc) ; + nfail += do_matrix (1, stdin, cc) ; } else { @@ -2829,7 +236,7 @@ int main (int argc, char **argv) fprintf (stderr, "Unable to open %s\n", matrix_name) ; nfail++ ; } - nfail += do_matrix (kind, matrix, cc) ; + nfail += do_matrix (kind, matrix, cc) ; fclose (matrix) ; } @@ -2840,7 +247,7 @@ int main (int argc, char **argv) // report the results // ------------------------------------------------------------------------- - cholmod_l_finish (cc) ; + spqr_finish (cc) ; if (cc->malloc_count != 0) { diff --git a/SPQR/Tcov/qrtest_template.hpp b/SPQR/Tcov/qrtest_template.hpp new file mode 100644 index 0000000000..35ff71a83e --- /dev/null +++ b/SPQR/Tcov/qrtest_template.hpp @@ -0,0 +1,2686 @@ +// ============================================================================= +// === qrest_template.hpp ====================================================== +// ============================================================================= + + + +// ============================================================================= +// === SPQR_qmult ============================================================== +// ============================================================================= + +// wrapper for SuiteSparseQR_qmult (dense), optionally testing memory alloc. + +template cholmod_dense *SPQR_qmult +( + // arguments for SuiteSparseQR_qmult: + int method, + cholmod_sparse *H, + cholmod_dense *Tau, + Int *HPinv, + cholmod_dense *X, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_dense *Y = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (Y) ; +} + + +// ============================================================================= +// === SPQR_qmult_sparse ======================================================= +// ============================================================================= + +// wrapper for SuiteSparseQR_qmult (sparse), optionally testing memory alloc. + +template cholmod_sparse *SPQR_qmult +( + // arguments for SuiteSparseQR_qmult: + int method, + cholmod_sparse *H, + cholmod_dense *Tau, + Int *HPinv, + cholmod_sparse *X, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_sparse *Y = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + Y = SuiteSparseQR_qmult (method, H, Tau, HPinv, X, cc); + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (Y) ; +} + +#ifndef NEXPERT + +// ============================================================================= +// === SPQR_qmult (dense case) ================================================= +// ============================================================================= + +// wrapper for SuiteSparseQR_qmult (dense), optionally testing memory alloc. + +template cholmod_dense *SPQR_qmult +( + // arguments for SuiteSparseQR_qmult: + int method, + SuiteSparseQR_factorization *QR, + cholmod_dense *X, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_dense *Y = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + Y = SuiteSparseQR_qmult (method, QR, X, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + Y = SuiteSparseQR_qmult (method, QR, X, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (Y) ; +} + +// ============================================================================= +// === SPQR_qmult (sparse case) ================================================ +// ============================================================================= + +// wrapper for SuiteSparseQR_qmult (sparse), optionally testing memory alloc. + +template cholmod_sparse *SPQR_qmult +( + // arguments for SuiteSparseQR_qmult: + int method, + SuiteSparseQR_factorization *QR, + cholmod_sparse *X, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_sparse *Y = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + Y = SuiteSparseQR_qmult (method, QR, X, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + Y = SuiteSparseQR_qmult (method, QR, X, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (Y) ; +} + + +// ============================================================================= +// === SPQR_solve (dense case) ================================================= +// ============================================================================= + +// Wrapper for testing SuiteSparseQR_solve + +template cholmod_dense *SPQR_solve +( + // arguments for SuiteSparseQR_solve: + int system, + SuiteSparseQR_factorization *QR, + cholmod_dense *B, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_dense *X = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + X = SuiteSparseQR_solve (system, QR, B, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + X = SuiteSparseQR_solve (system, QR, B, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (X) ; +} + + +// ============================================================================= +// === SPQR_solve (sparse case) ================================================ +// ============================================================================= + +// Wrapper for testing SuiteSparseQR_solve + +template cholmod_sparse *SPQR_solve +( + // arguments for SuiteSparseQR_solve: + int system, + SuiteSparseQR_factorization *QR, + cholmod_sparse *B, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_sparse *X = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + X = SuiteSparseQR_solve (system, QR, B, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + X = SuiteSparseQR_solve (system, QR, B, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (X) ; +} + + +// ============================================================================= +// === SPQR_min2norm (dense case) ============================================== +// ============================================================================= + +// Wrapper for testing SuiteSparseQR_min2norm + +template cholmod_dense *SPQR_min2norm +( + // arguments for SuiteSparseQR_min2norm: + int ordering, + double tol, + cholmod_sparse *A, + cholmod_dense *B, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_dense *X = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (X) ; +} + +// ============================================================================= +// === SPQR_min2norm (sparse case) ============================================= +// ============================================================================= + +// Wrapper for testing SuiteSparseQR_min2norm + +template cholmod_sparse *SPQR_min2norm +( + // arguments for SuiteSparseQR_min2norm: + int ordering, + double tol, + cholmod_sparse *A, + cholmod_sparse *B, + cholmod_common *cc, + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + cholmod_sparse *X = NULL ; + if (!memory_test) + { + // just call the method directly; no memory testing + X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + X = SuiteSparseQR_min2norm (ordering, tol, A, B, cc) ; + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + return (X) ; +} + +// ============================================================================= +// === SPQR_factorize ========================================================== +// ============================================================================= + +// Wrapper for testing SuiteSparseQR_factorize or +// SuiteSparseQR_symbolic and SuiteSparseQR_numeric + +template +SuiteSparseQR_factorization *SPQR_factorize +( + // arguments for SuiteSparseQR_factorize: + int ordering, + double tol, + cholmod_sparse *A, + cholmod_common *cc, + + // method to use + int split, // if 1 use SuiteSparseQR_symbolic followed by + // SuiteSparseQR_numeric, if 0 use + // SuiteSparseQR_factorize, if 2, do the + // numeric factorization twice, just for testing. + // if 3 use SuiteSparseQR_C_factorize + // if 3 use SuiteSparseQR_C_symbolic / _C_numeric + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + SuiteSparseQR_factorization *QR ; + SuiteSparseQR_C_factorization *C_QR ; + if (!memory_test) + { + // just call the method directly; no memory testing + + if (split == 4) + { + + C_QR = SuiteSparseQR_C_symbolic (ordering, tol >= 0, A, cc) ; + SuiteSparseQR_C_numeric (tol, A, C_QR, cc) ; + if (C_QR == NULL) + { + cc->status = CHOLMOD_OUT_OF_MEMORY ; + QR = NULL ; + } + else + { + QR = (SuiteSparseQR_factorization *) + (C_QR->factors) ; + if (QR == NULL || QR->QRnum == NULL) + { + cc->status = CHOLMOD_OUT_OF_MEMORY ; + QR = NULL ; + } + else + { + // QR itself will be kept; free the C wrapper + C_QR->factors = NULL ; + cc->status = CHOLMOD_OK ; + if (QR == NULL) printf ("Hey!\n") ; + } + } + SuiteSparseQR_C_free (&C_QR, cc) ; + + } + else if (split == 3) + { + + C_QR = SuiteSparseQR_C_factorize (ordering, tol, A, cc) ; + int save = cc->status ; + if (C_QR == NULL) + { + QR = NULL ; + } + else + { + QR = (SuiteSparseQR_factorization *) + (C_QR->factors) ; + C_QR->factors = NULL ; + SuiteSparseQR_C_free (&C_QR, cc) ; + } + cc->status = save ; + + } + else if (split == 2) + { + QR = SuiteSparseQR_symbolic (ordering, + tol >= 0 || tol <= SPQR_DEFAULT_TOL, A, cc) ; + ASSERT (QR != NULL) ; + SuiteSparseQR_numeric (tol, A, QR, cc) ; + // just for testing + SuiteSparseQR_numeric (tol, A, QR, cc) ; + } + else if (split == 1) + { + // split symbolic/numeric, no singletons exploited + QR = SuiteSparseQR_symbolic (ordering, + tol >= 0 || tol <= SPQR_DEFAULT_TOL, A, cc) ; + ASSERT (QR != NULL) ; + SuiteSparseQR_numeric (tol, A, QR, cc) ; + } + else + { + QR = SuiteSparseQR_factorize (ordering, tol, A, cc) ; + } + + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + + if (split == 4) + { + + C_QR = SuiteSparseQR_C_symbolic (ordering, tol >= 0, A, cc) ; + SuiteSparseQR_C_numeric (tol, A, C_QR, cc) ; + if (C_QR == NULL) + { + cc->status = CHOLMOD_OUT_OF_MEMORY ; + QR = NULL ; + } + else + { + QR = (SuiteSparseQR_factorization *) + (C_QR->factors) ; + if (QR == NULL || QR->QRnum == NULL) + { + cc->status = CHOLMOD_OUT_OF_MEMORY ; + QR = NULL ; + } + else + { + // QR itself will be kept; free the C wrapper + C_QR->factors = NULL ; + cc->status = CHOLMOD_OK ; + if (QR == NULL) printf ("Hey!!\n") ; + } + } + SuiteSparseQR_C_free (&C_QR, cc) ; + + + } + else if (split == 3) + { + + C_QR = SuiteSparseQR_C_factorize (ordering, tol, A, cc) ; + int save = cc->status ; + if (C_QR == NULL) + { + QR = NULL ; + } + else + { + QR = (SuiteSparseQR_factorization *) + (C_QR->factors) ; + C_QR->factors = NULL ; + SuiteSparseQR_C_free (&C_QR, cc) ; + } + cc->status = save ; + + } + else if (split == 1 || split == 2) + { + // split symbolic/numeric, no singletons exploited + QR = SuiteSparseQR_symbolic (ordering, + tol >= 0 || tol <= SPQR_DEFAULT_TOL, A, cc) ; + if (cc->status < CHOLMOD_OK) + { + continue ; + } + ASSERT (QR != NULL) ; + SuiteSparseQR_numeric (tol, A, QR, cc) ; + if (cc->status < CHOLMOD_OK) + { + SuiteSparseQR_free (&QR, cc) ; + } + } + else + { + QR = SuiteSparseQR_factorize (ordering, tol, A, cc) ; + } + + if (cc->status == CHOLMOD_OK) break ; + } + normal_memory_handler (cc, true) ; + } + + if (QR == NULL) printf ("huh?? split: %d ordering %d tol %g\n", split, + ordering, tol) ; + return (QR) ; +} + + +#endif + +// ============================================================================= +// === SPQR_qr ================================================================= +// ============================================================================= + +// wrapper for SuiteSparseQR, optionally testing memory allocation + +template Int SPQR_qr +( + // arguments for SuiteSparseQR: + int ordering, + double tol, + Int econ, + int getCTX, + cholmod_sparse *A, + cholmod_sparse *Bsparse, + cholmod_dense *Bdense, + cholmod_sparse **Zsparse, + cholmod_dense **Zdense, + cholmod_sparse **R, + Int **E, + cholmod_sparse **H, + Int **HPinv, + cholmod_dense **HTau, + cholmod_common *cc, + + // which version to use (C or C++) + int use_c_version, // if TRUE use C version, otherwise use C++ + + // malloc control + int memory_test, // if TRUE, test malloc error handling + int memory_punt // if TRUE, test punt case +) +{ + Int rank ; + if (!memory_test) + { + // just call the method directly; no memory testing + if (use_c_version) + { + if (sizeof (Int) == sizeof (int64_t)) + { + rank = SuiteSparseQR_C (ordering, tol, econ, getCTX, A, + Bsparse, Bdense, Zsparse, Zdense, R, + (int64_t **) E, H, (int64_t **) HPinv, HTau, cc) ; + } + else + { + rank = SuiteSparseQR_i_C (ordering, tol, econ, getCTX, A, + Bsparse, Bdense, Zsparse, Zdense, R, + (int32_t **) E, H, (int32_t **) HPinv, HTau, cc) ; + } + } + else + { + rank = SuiteSparseQR (ordering, tol, econ, getCTX, A, + Bsparse, Bdense, Zsparse, Zdense, R, E, H, HPinv, HTau, cc) ; + } + } + else + { + // test malloc error handling + int64_t tries ; + test_memory_handler (cc, true) ; + my_punt = memory_punt ; + for (tries = 0 ; my_tries < 0 ; tries++) + { + my_tries = tries ; + if (use_c_version) + { + if (sizeof (Int) == sizeof (int64_t)) + { + rank = SuiteSparseQR_C (ordering, tol, econ, getCTX, A, + Bsparse, Bdense, Zsparse, Zdense, R, + (int64_t **) E, H, (int64_t **) HPinv, HTau, cc) ; + } + else + { + rank = SuiteSparseQR_i_C (ordering, tol, econ, getCTX, A, + Bsparse, Bdense, Zsparse, Zdense, R, + (int32_t **) E, H, (int32_t **) HPinv, HTau, cc) ; + } + } + else + { + rank = SuiteSparseQR (ordering, tol, econ, + getCTX, A, Bsparse, Bdense, Zsparse, Zdense, R, + E, H, HPinv, HTau, cc); + } + if (cc->status == CHOLMOD_OK) + { + break ; + } + } + normal_memory_handler (cc, true) ; + } + return (rank) ; +} + + +// ============================================================================= +// === my_rand ================================================================= +// ============================================================================= + +// The POSIX example of rand, duplicated here so that the same sequence will +// be generated on different machines. + +static unsigned long next = 1 ; + +#define MY_RAND_MAX 32767 + +// RAND_MAX assumed to be 32767 +int64_t my_rand (void) +{ + next = next * 1103515245 + 12345 ; + return ((unsigned)(next/65536) % (MY_RAND_MAX + 1)) ; +} + +void my_srand (unsigned seed) +{ + next = seed ; +} + +unsigned long my_seed (void) +{ + return (next) ; +} + +int64_t nrand (int64_t n) // return a random int64_t between 0 and n-1 +{ + return ((n <= 0) ? 0 : (my_rand ( ) % n)) ; +} + +double xrand ( ) // return a random double between -1 and 1 +{ + double x = ((double) my_rand ( )) / MY_RAND_MAX ; + return (2*x-1) ; +} + +double erand (double range) +{ + return (range * xrand ( )) ; +} + +Complex erand (Complex range) +{ + /* + Complex x ; + x.real ( ) = xrand ( ) ; + x.imag ( ) = xrand ( ) ; + return (range * x) ; + */ + Complex i = Complex (0,1) ; + return (range * (xrand ( ) + i * xrand ( ))) ; +} + + +// ============================================================================= +// === getreal ================================================================= +// ============================================================================= + +// return real part of a scalar x + +inline double getreal (double x) +{ + return (x) ; +} + +inline double getreal (Complex x) +{ + return (x.real ( )) ; +} + +// ============================================================================= +// === getimag ================================================================= +// ============================================================================= + +// return imaginary part of a scalar x +inline double getimag (double x) // x is an unused parameter +{ + return (0) ; +} + +inline double getimag (Complex x) +{ + return (x.imag ( )) ; +} + + +// ============================================================================= +// === dense_wrapper =========================================================== +// ============================================================================= + +// place a column-oriented matrix in a cholmod_dense wrapper + +template cholmod_dense *dense_wrapper +( + cholmod_dense *X, + Int nrow, + Int ncol, + Entry *Xx +) +{ + X->xtype = spqr_type ( ) ; + X->nrow = nrow ; + X->ncol = ncol ; + X->d = nrow ; // leading dimension = nrow + X->nzmax = nrow * ncol ; + X->x = Xx ; + X->z = NULL ; // ZOMPLEX case not supported + X->dtype = CHOLMOD_DOUBLE ; + return (X) ; +} + +// ============================================================================= +// === sparse_split ============================================================ +// ============================================================================= + +// Return the real or imaginary part of a packed complex sparse matrix + +template +cholmod_sparse *sparse_split +( + cholmod_sparse *A, + int part, + cholmod_common *cc +) +{ + cholmod_sparse *C ; + if (!A || A->xtype != CHOLMOD_COMPLEX || A->nz != NULL) return (NULL) ; + if (! (part == 0 || part == 1)) return (NULL) ; + + int64_t nz = spqr_nnz (A, cc) ; + C = spqr_allocate_sparse (A->nrow, A->ncol, nz, TRUE, TRUE, 0, + CHOLMOD_REAL, cc) ; + + Int *Ap = (Int *) A->p ; + Int *Ai = (Int *) A->i ; + double *Ax = (double *) A->x ; + + Int *Cp = (Int *) C->p ; + Int *Ci = (Int *) C->i ; + double *Cx = (double *) C->x ; + + Int n = A->ncol ; + + for (Int k = 0 ; k < n+1 ; k++) + { + Cp [k] = Ap [k] ; + } + + for (int64_t k = 0 ; k < nz ; k++) + { + Ci [k] = Ai [k] ; + } + + for (Int k = 0 ; k < nz ; k++) + { + Cx [k] = Ax [2*k + part] ; + } + + return (C) ; +} + +template +cholmod_sparse *sparse_real (cholmod_sparse *A, cholmod_common *cc) +{ + return (sparse_split (A, 0, cc)) ; +} + +template +cholmod_sparse *sparse_imag (cholmod_sparse *A, cholmod_common *cc) +{ + return (sparse_split (A, 1, cc)) ; +} + +// ============================================================================= +// === sparse_merge ============================================================ +// ============================================================================= + +// Add the real and imaginary parts of a matrix (both stored in real form) +// into a single matrix. The two parts must have the same nonzero pattern. +// A is CHOLMOD_REAL on input and holds the real part, it is CHOLMOD_COMPLEX +// on output. A_imag is CHOLMOD_REAL on input; it holds the imaginary part +// of A as a real matrix. + +template int sparse_merge +( + cholmod_sparse *A, // input/output + cholmod_sparse *A_imag, // input only + cholmod_common *cc +) +{ + if (A == NULL || A_imag == NULL) + { + return (FALSE) ; + } + int64_t nz1 = spqr_nnz (A, cc) ; + int64_t nz2 = spqr_nnz (A_imag, cc) ; + if (A->xtype != CHOLMOD_REAL || A_imag->xtype != CHOLMOD_REAL || nz1 != nz2) + { + return (FALSE) ; + } + + // change A from real to complex + spqr_sparse_xtype (CHOLMOD_COMPLEX, A, cc) ; + + double *Ax = (double *) A->x ; + double *Az = (double *) A_imag->x ; + + // merge in the imaginary part from A_imag into A + for (int64_t k = 0 ; k < nz1 ; k++) + { + Ax [2*k+1] = Az [k] ; + } + return (TRUE) ; +} + +// ============================================================================= +// === sparse_diff ============================================================= +// ============================================================================= + +// Compute C = A-B where A and B are either both real, or both complex + +template cholmod_sparse *sparse_diff +( + cholmod_sparse *A, + cholmod_sparse *B, + cholmod_common *cc +) +{ + cholmod_sparse *C ; + double one [2] = {1,0}, minusone [2] = {-1,0} ; + + if (spqr_type ( ) == CHOLMOD_REAL) + { + // C = A - B + C = spqr_ssadd (A, B, one, minusone, TRUE, TRUE, cc) ; + } + else + { + cholmod_sparse *A_real, *A_imag, *B_real, *B_imag, *C_imag ; + + A_real = sparse_real (A, cc) ; + A_imag = sparse_imag (A, cc) ; + + B_real = sparse_real (B, cc) ; + B_imag = sparse_imag (B, cc) ; + + // real(C) = real(A) - real (B) + C = spqr_ssadd (A_real, B_real, one, minusone, TRUE, TRUE, cc) ; + + // imag(C) = imag(A) - imag (B) + C_imag = spqr_ssadd (A_imag, B_imag, one, minusone, TRUE, TRUE, cc); + + // C = real(C) + 1i*imag(C) + sparse_merge (C, C_imag, cc) ; + spqr_free_sparse (&C_imag, cc) ; + + spqr_free_sparse (&A_real, cc) ; + spqr_free_sparse (&A_imag, cc) ; + spqr_free_sparse (&B_real, cc) ; + spqr_free_sparse (&B_imag, cc) ; + } + + return (C) ; +} + + +// ============================================================================= +// === permute_columns ========================================================= +// ============================================================================= + +// compute A(:,P) + +template cholmod_sparse *permute_columns +( + cholmod_sparse *A, + Int *P, + cholmod_common *cc +) +{ + Int m = A->nrow ; + Int n = A->ncol ; + Int nz = spqr_nnz (A, cc) ; + int xtype = spqr_type ( ) ; + Int *Ap = (Int *) A->p ; + Int *Ai = (Int *) A->i ; + Entry *Ax = (Entry *) A->x ; + cholmod_sparse *C ; + + // allocate empty matrix C with space for nz entries + C = spqr_allocate_sparse (m, n, nz, TRUE, TRUE, 0, xtype, cc) ; + Int *Cp = (Int *) C->p ; + Int *Ci = (Int *) C->i ; + Entry *Cx = (Entry *) C->x ; + + // construct column pointers for C + for (Int k = 0 ; k < n ; k++) + { + // column j of A becomes column k of C + Int j = P ? P [k] : k ; + Cp [k] = Ap [j+1] - Ap [j] ; + } + spqr_cumsum (n, Cp) ; + + // copy columns from A to C + for (Int k = 0 ; k < n ; k++) + { + // copy column k of A into column j of C + Int j = P ? P [k] : k ; + Int pdest = Cp [k] ; + Int psrc = Ap [j] ; + Int len = Ap [j+1] - Ap [j] ; + for (Int t = 0 ; t < len ; t++) + { + Ci [pdest + t] = Ai [psrc + t] ; + Cx [pdest + t] = Ax [psrc + t] ; + } + } + + return (C) ; +} + + +// ============================================================================= +// === sparse_multiply ========================================================= +// ============================================================================= + +// compute A*B where A and B can be both real or both complex + +template cholmod_sparse *sparse_multiply +( + cholmod_sparse *A, + cholmod_sparse *B, + cholmod_common *cc +) +{ + cholmod_sparse *C ; + + if (spqr_type ( ) == CHOLMOD_REAL) + { + // C = A*B + C = spqr_ssmult (A, B, 0, TRUE, TRUE, cc) ; + } + else + { + // cholmod_ssmult and cholmod_add only work for real matrices + cholmod_sparse *A_real, *A_imag, *B_real, *B_imag, *C_imag, *T1, *T2 ; + double one [2] = {1,0}, minusone [2] = {-1,0} ; + + A_real = sparse_real (A, cc) ; + A_imag = sparse_imag (A, cc) ; + + B_real = sparse_real (B, cc) ; + B_imag = sparse_imag (B, cc) ; + + // real(C) = real(A)*real(B) - imag(A)*imag(B) + T1 = spqr_ssmult (A_real, B_real, 0, TRUE, TRUE, cc) ; + T2 = spqr_ssmult (A_imag, B_imag, 0, TRUE, TRUE, cc) ; + C = spqr_ssadd (T1, T2, one, minusone, TRUE, TRUE, cc) ; + spqr_free_sparse (&T1, cc) ; + spqr_free_sparse (&T2, cc) ; + + // imag(C) = imag(A)*real(B) + real(A)*imag(B) + T1 = spqr_ssmult (A_imag, B_real, 0, TRUE, TRUE, cc) ; + T2 = spqr_ssmult (A_real, B_imag, 0, TRUE, TRUE, cc) ; + C_imag = spqr_ssadd (T1, T2, one, one, TRUE, TRUE, cc) ; + spqr_free_sparse (&T1, cc) ; + spqr_free_sparse (&T2, cc) ; + + // C = real(C) + 1i*imag(C) + sparse_merge (C, C_imag, cc) ; + spqr_free_sparse (&C_imag, cc) ; + + spqr_free_sparse (&A_real, cc) ; + spqr_free_sparse (&A_imag, cc) ; + spqr_free_sparse (&B_real, cc) ; + spqr_free_sparse (&B_imag, cc) ; + } + + return (C) ; +} + + +// ============================================================================= +// === sparse_resid ============================================================ +// ============================================================================= + +// compute norm (A*x-b,1) for A,x, and b all sparse + +template double sparse_resid +( + cholmod_sparse *A, + double anorm, + cholmod_sparse *X, + cholmod_sparse *B, + cholmod_common *cc +) +{ + cholmod_sparse *AX, *Resid ; + // AX = A*X + AX = sparse_multiply (A, X, cc) ; + // Resid = AX - B + Resid = sparse_diff (AX, B, cc) ; + // resid = norm (Resid,1) + double resid = spqr_norm_sparse (Resid, 1, cc) ; + resid = CHECK_NAN (resid) ; + spqr_free_sparse (&AX, cc) ; + spqr_free_sparse (&Resid, cc) ; + return (CHECK_NAN (resid / anorm)) ; +} + +// ============================================================================= +// === dense_resid ============================================================= +// ============================================================================= + +// compute norm (A*x-b,1) for A sparse, x and b dense + +template double dense_resid +( + cholmod_sparse *A, + double anorm, + cholmod_dense *X, + Int nb, + Entry *Bx, + cholmod_common *cc +) +{ + cholmod_dense *B, Bmatrix, *Resid ; + double one [2] = {1,0}, minusone [2] = {-1,0} ; + + B = dense_wrapper (&Bmatrix, A->nrow, nb, Bx) ; + // Resid = B + Resid = spqr_copy_dense (B, cc) ; + + // Resid = A*X - Resid + spqr_sdmult (A, FALSE, one, minusone, X, Resid, cc) ; + + // resid = norm (Resid,1) + double resid = spqr_norm_dense (Resid, 1, cc) ; + resid = CHECK_NAN (resid) ; + spqr_free_dense (&Resid, cc) ; + return (CHECK_NAN (resid / anorm)) ; +} + +// ============================================================================= +// === check_r_factor ========================================================== +// ============================================================================= + +// compute norm (R'*R - (A(:,P))'*(A(:,P)), 1) / norm (A'*A,1) + +template double check_r_factor +( + cholmod_sparse *R, + cholmod_sparse *A, + Int *P, + cholmod_common *cc +) +{ + cholmod_sparse *RTR, *RT, *C, *CT, *CTC, *D ; + + // RTR = R'*R + RT = spqr_transpose (R, 2, cc) ; + RTR = sparse_multiply (RT, R, cc) ; + spqr_free_sparse (&RT, cc) ; + + // C = A(:,P) + C = permute_columns (A, P, cc) ; + + // CTC = C'*C + CT = spqr_transpose (C, 2, cc) ; + CTC = sparse_multiply (CT, C, cc) ; + spqr_free_sparse (&CT, cc) ; + spqr_free_sparse (&C, cc) ; + + double ctcnorm = spqr_norm_sparse (CTC, 1, cc) ; + if (ctcnorm == 0) ctcnorm = 1 ; + ctcnorm = CHECK_NAN (ctcnorm) ; + + // D = RTR - CTC + D = sparse_diff (RTR, CTC, cc) ; + spqr_free_sparse (&CTC, cc) ; + spqr_free_sparse (&RTR, cc) ; + + // err = norm (D,1) + double err = spqr_norm_sparse (D, 1, cc) / ctcnorm ; + err = CHECK_NAN (err) ; + + spqr_free_sparse (&D, cc) ; + return (CHECK_NAN (err)) ; +} + + +// ============================================================================= +// === check_qr ================================================================ +// ============================================================================= + +// compute norm (Q*R - A(:,P)) / norm (A) + +template double check_qr +( + cholmod_sparse *Q, + cholmod_sparse *R, + cholmod_sparse *A, + Int *P, + double anorm, + cholmod_common *cc +) +{ + cholmod_sparse *QR, *C, *D ; + + // C = A(:,P) + C = permute_columns (A, P, cc) ; + + // QR = Q*R + QR = sparse_multiply (Q, R, cc) ; + + // D = RTR - CTC + D = sparse_diff (QR, C, cc) ; + spqr_free_sparse (&QR, cc) ; + spqr_free_sparse (&C, cc) ; + + // err = norm (D,1) + double err = spqr_norm_sparse (D, 1, cc) / anorm ; + err = CHECK_NAN (err) ; + + spqr_free_sparse (&D, cc) ; + return (CHECK_NAN (err)) ; +} + +// ============================================================================= +// === Rsolve ================================================================== +// ============================================================================= + +template int Rsolve +( + // R is n-by-n, upper triangular with zero-free diagonal + Int n, + cholmod_sparse *R, + Entry *X, // X is n-by-nx, leading dimension n, overwritten with soln + Int nx, + cholmod_common *cc +) +{ + // Int n = R->n ; + Int *Rp = (Int *) R->p ; + Int *Ri = (Int *) R->i ; + Entry *Rx = (Entry *) R->x ; + + // check the diagonal + for (Int j = 0 ; j < n ; j++) + { + if (Rp [j] == Rp [j+1] || Ri [Rp [j+1]-1] != j) + { + printf ("Rsolve: R not upper triangular w/ zero-free diagonal\n") ; + return (FALSE) ; + } + } + + // do the backsolve + for (Int k = 0 ; k < nx ; k++) + { + for (Int j = n-1 ; j >= 0 ; j--) + { + Entry rjj = Rx [Rp [j+1]-1] ; + if (rjj == (Entry) 0) + { + printf ("Rsolve: R has an explicit zero on the diagonal\n") ; + return (FALSE) ; + } + X [j] /= rjj ; + for (Int p = Rp [j] ; p < Rp [j+1]-1 ; p++) + { + X [Ri [p]] -= Rx [p] * X [j] ; + } + } + X += n ; + } + + return (TRUE) ; +} + + +// ============================================================================= +// === create_Q ================================================================ +// ============================================================================= + +// create a sparse Q +template cholmod_sparse *create_Q +( + cholmod_sparse *H, + cholmod_dense *HTau, + Int *HPinv, + cholmod_common *cc +) +{ + cholmod_sparse *Q, *I ; + Int m = H->nrow ; + int xtype = spqr_type ( ) ; + I = spqr_speye (m, m, xtype, cc) ; + Q = SPQR_qmult (1, H, HTau, HPinv, I, cc, m<300, nrand (2)) ; + spqr_free_sparse (&I, cc) ; + return (Q) ; +} + + +// ============================================================================= +// === QRsolve ================================================================= +// ============================================================================= + +// solve Ax=b using H, R, and E + +template double QRsolve +( + cholmod_sparse *A, + double anorm, + Int rank, + int method, + cholmod_sparse *H, + cholmod_dense *HTau, + Int *HPinv, + cholmod_sparse *R, + Int *Qfill, + cholmod_dense *Bdense, + cholmod_common *cc +) +{ + double one [2] = {1,0}, zero [2] = {0,0}, resid = EMPTY ; + int xtype = spqr_type ( ) ; + Int n = A->ncol ; + Int m = A->nrow ; + Int nrhs = Bdense->ncol ; + Entry *X, *Y = NULL, *B ; + cholmod_dense *Ydense = NULL ; + cholmod_dense *Xdense ; + + B = (Entry *) Bdense->x ; + Xdense = spqr_zeros (n, nrhs, xtype, cc) ; + X = (Entry *) Xdense->x ; + + if (method == 0) + { + // solve Ax=b using H, R, and E + // Y = Q'*B + Ydense = SPQR_qmult (0, H, HTau, HPinv, Bdense, cc, + m < 300, nrand (2)) ; + } + else + { + cholmod_sparse *Q, *QT ; + // solve using Q instead of qmult + Q = create_Q (H, HTau, HPinv, cc) ; + QT = spqr_transpose (Q, 2, cc) ; + // Y = Q'*B + Ydense = spqr_zeros (m, nrhs, xtype, cc) ; + spqr_sdmult (QT, FALSE, one, zero, Bdense, Ydense, cc) ; + spqr_free_sparse (&Q, cc) ; + spqr_free_sparse (&QT, cc) ; + } + + // Y (1:rank) = R (1:rank,1:rank) \ Y (1:rank) + Y = (Entry *) Ydense->x ; + Int ok = Rsolve (rank, R, Y, nrhs, cc) ; + // X = E*Y + if (ok) + { + for (Int kk = 0 ; kk < nrhs ; kk++) + { + for (Int k = 0 ; k < rank ; k++) + { + Int j = Qfill ? Qfill [k] : k ; + X [j] = Y [k] ; + } + X += n ; + Y += m ; + } + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, Xdense, nrhs, B, cc) ; + } + + spqr_free_dense (&Ydense, cc) ; + spqr_free_dense (&Xdense, cc) ; + return (CHECK_NAN (resid)) ; +} + +// ============================================================================= +// === check_qmult ============================================================= +// ============================================================================= + +// Test qmult + +template double check_qmult +( + cholmod_sparse *H, + cholmod_dense *HTau, + Int *HPinv, + int test_errors, + cholmod_common *cc +) +{ + cholmod_sparse *Q, *QT, *Xsparse, *Ssparse, *Zsparse ; + cholmod_dense *Xdense, *Zdense, *Sdense, *Ydense ; + int xtype = spqr_type ( ) ; + Entry *X, *Y, *Z, *S ; + double err, maxerr = 0 ; + double one [2] = {1,0}, zero [2] = {0,0} ; + Entry range = (Entry) 1.0 ; + Int k ; + + Int m = H->nrow ; + Q = create_Q (H, HTau, HPinv, cc) ; // construct Q from H + + QT = spqr_transpose (Q, 2, cc) ; // QT = Q' + + // compare Q with qmult for sparse and dense X + for (Int nx = 0 ; nx < 5 ; nx++) // # of columns of X + { + Int xsize = m * nx ; // size of X + for (Int nz = 1 ; nz <= xsize+1 ; nz *= 16) // # of nonzeros in X + { + + // ----------------------------------------------------------------- + // create X as m-by-nx, both sparse and dense + // ----------------------------------------------------------------- + + Xdense = spqr_zeros (m, nx, xtype, cc) ; + X = (Entry *) Xdense->x ; + for (k = 0 ; k < nz ; k++) + { + X [nrand (xsize)] += erand (range) ; + } + Xsparse = spqr_dense_to_sparse (Xdense, TRUE, cc) ; + + // ----------------------------------------------------------------- + // Y = Q'*X for method 0, Y = Q*X for method 1 + // ----------------------------------------------------------------- + + for (int method = 0 ; method <= 1 ; method++) + { + // Y = Q'*X or Q*X + Ydense = SPQR_qmult (method, H, HTau, HPinv, + Xdense, cc, m < 300, nrand (2)) ; + Y = (Entry *) Ydense->x ; + + Zdense = spqr_zeros (m, nx, xtype, cc) ; + spqr_sdmult (Q, (method == 0), one, zero, + Xdense, Zdense, cc) ; + Z = (Entry *) Zdense->x ; + + for (err = 0, k = 0 ; k < xsize ; k++) + { + double e1 = spqr_abs (Y [k] - Z [k], cc) ; + e1 = CHECK_NAN (e1) ; + err = MAX (err, e1) ; + } + maxerr = MAX (maxerr, err) ; + + // S = Q'*Xsparse or Q*Xsparse + Ssparse = SPQR_qmult ( + method, H, HTau, HPinv, Xsparse, cc, m < 300, nrand (2)) ; + Sdense = spqr_sparse_to_dense (Ssparse, cc) ; + S = (Entry *) Sdense->x ; + + for (err = 0, k = 0 ; k < xsize ; k++) + { + double e1 = spqr_abs (S [k] - Y [k], cc) ; + e1 = CHECK_NAN (e1) ; + err = MAX (err, e1) ; + } + maxerr = MAX (maxerr, err) ; + + spqr_free_dense (&Ydense, cc) ; + spqr_free_dense (&Zdense, cc) ; + spqr_free_sparse (&Ssparse, cc) ; + spqr_free_dense (&Sdense, cc) ; + } + + spqr_free_dense (&Xdense, cc) ; + spqr_free_sparse (&Xsparse, cc) ; + + // ----------------------------------------------------------------- + // create X as nx-by-m, both sparse and dense + // ----------------------------------------------------------------- + + Xdense = spqr_zeros (nx, m, xtype, cc) ; + X = (Entry *) Xdense->x ; + for (k = 0 ; k < nz ; k++) + { + X [nrand (xsize)] += erand (range) ; + } + Xsparse = spqr_dense_to_sparse (Xdense, TRUE, cc) ; + + // ----------------------------------------------------------------- + // Y = X*Q' for method 2, Y = X*Q for method 3 + // ----------------------------------------------------------------- + + for (int method = 2 ; method <= 3 ; method++) + { + // Y = X*Q' or X*Q + Ydense = SPQR_qmult (method, H, HTau, HPinv, + Xdense, cc, m < 300, nrand (2)) ; + Y = (Entry *) Ydense->x ; + + if (method == 2) + { + // Zsparse = (X*Q') + Zsparse = sparse_multiply (Xsparse, QT, cc) ; + } + else + { + // Zsparse = (X*Q) + Zsparse = sparse_multiply (Xsparse, Q, cc) ; + } + Zdense = spqr_sparse_to_dense (Zsparse, cc) ; + + Z = (Entry *) Zdense->x ; + + for (err = 0, k = 0 ; k < xsize ; k++) + { + double e1 = spqr_abs (Y [k] - Z [k], cc) ; + e1 = CHECK_NAN (e1) ; + err = MAX (err, e1) ; + } + maxerr = MAX (maxerr, err) ; + + // S = X*Q' or X*Q + Ssparse = SPQR_qmult ( + method, H, HTau, HPinv, Xsparse, cc, m < 300, nrand (2)) ; + Sdense = spqr_sparse_to_dense (Ssparse, cc) ; + S = (Entry *) Sdense->x ; + + for (err = 0, k = 0 ; k < xsize ; k++) + { + double e1 = spqr_abs (S [k] - Y [k], cc) ; + e1 = CHECK_NAN (e1) ; + err = MAX (err, e1) ; + } + maxerr = MAX (maxerr, err) ; + + spqr_free_dense (&Ydense, cc) ; + spqr_free_dense (&Zdense, cc) ; + spqr_free_sparse (&Ssparse, cc) ; + spqr_free_sparse (&Zsparse, cc) ; + spqr_free_dense (&Sdense, cc) ; + } + + spqr_free_dense (&Xdense, cc) ; + spqr_free_sparse (&Xsparse, cc) ; + + } + } + spqr_free_sparse (&Q, cc) ; + spqr_free_sparse (&QT, cc) ; + + // ------------------------------------------------------------------------- + // qmult error conditions + // ------------------------------------------------------------------------- + + // These should fail; expect 6 error messages in the output + + if (test_errors) + { + err = 0 ; + printf ("The following six errors are expected:\n") ; + Xdense = spqr_zeros (m+1, 1, xtype,cc) ; + Ydense = SuiteSparseQR_qmult (0, H, HTau, HPinv, Xdense, cc) ; + if (Ydense != NULL || cc->status != CHOLMOD_INVALID) err++ ; + spqr_free_dense (&Xdense, cc) ; + + Xdense = spqr_zeros (1, m+1, xtype,cc) ; + Ydense = SuiteSparseQR_qmult (2, H, HTau, HPinv, Xdense, cc) ; + if (Ydense != NULL || cc->status != CHOLMOD_INVALID) err++ ; + Ydense = SuiteSparseQR_qmult (42, H, HTau, HPinv, Xdense, cc) ; + if (Ydense != NULL || cc->status != CHOLMOD_INVALID) err++ ; + spqr_free_dense (&Xdense, cc) ; + + Xsparse = spqr_speye (m+1, m+1, xtype, cc) ; + + Q = SuiteSparseQR_qmult (0, H, HTau, HPinv, Xsparse, cc); + if (Q != NULL || cc->status != CHOLMOD_INVALID) err++ ; + Q = SuiteSparseQR_qmult (2, H, HTau, HPinv, Xsparse, cc); + if (Q != NULL || cc->status != CHOLMOD_INVALID) err++ ; + Q = SuiteSparseQR_qmult (9, H, HTau, HPinv, Xsparse, cc); + if (Q != NULL || cc->status != CHOLMOD_INVALID) err++ ; + + printf (" ... error handling done\n\n") ; + spqr_free_sparse (&Xsparse, cc) ; + + maxerr = MAX (maxerr, err) ; + } + + return (CHECK_NAN (maxerr)) ; +} + + +// ============================================================================= +// === check_rc ================================================================ +// ============================================================================= + +// X = Q'*B has been done, continue with C=R\C and + +template double check_rc +( + Int rank, + cholmod_sparse *R, + cholmod_sparse *A, + Entry *B, + cholmod_dense *X, + Int nrhs, + double anorm, + Int *Qfill, + cholmod_common *cc +) +{ + double resid = EMPTY ; + int xtype = spqr_type ( ) ; + cholmod_dense *W ; + Int n ; + int ok ; + Entry *W1, *X1 ; + if (!R || !X) + { + ok = 0 ; + } + else + { + n = X->nrow ; + X1 = (Entry *) X->x ; + // solve X = R\X, overwriting X with solution + ok = Rsolve (rank, R, X1, nrhs, cc) ; + } + if (ok) + { + // W = E*X + // W = (Entry *) spqr_calloc (A->ncol * nrhs, sizeof (Entry), cc); + W = spqr_zeros (A->ncol, nrhs, xtype, cc) ; + W1 = (Entry *) W->x ; + for (Int col = 0 ; col < nrhs ; col++) + { + for (Int k = 0 ; k < rank ; k++) + { + Int j = Qfill ? Qfill [k] : k ; + if (j < (Int) A->ncol) W1 [j] = X1 [k] ; + } + W1 += A->ncol ; + X1 += n ; + } + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, W, nrhs, B, cc) ; + // spqr_free (A->ncol * nrhs, sizeof (Entry), W, cc) ; + spqr_free_dense (&W, cc) ; + } + return (CHECK_NAN (resid)) ; +} + + +// ============================================================================= +// === transpose =============================================================== +// ============================================================================= + +// Transpose a dense matrix. + +template cholmod_dense *transpose +( + cholmod_dense *Xdense, + cholmod_common *cc +) +{ + Entry *X, *Y ; + cholmod_dense *Ydense ; + if (Xdense == NULL) + { + printf ("transpose failed!\n") ; + return (NULL) ; + } + Int m = Xdense->nrow ; + Int n = Xdense->ncol ; + Int ldx = Xdense->d ; + int xtype = spqr_type ( ) ; + Ydense = spqr_allocate_dense (n, m, n, xtype, cc) ; + X = (Entry *) Xdense->x ; + Y = (Entry *) Ydense->x ; + for (Int i = 0 ; i < m ; i++) + { + for (Int j = 0 ; j < n ; j++) + { + Y [j+i*n] = spqr_conj (X [i+j*ldx]) ; + } + } + return (Ydense) ; +} + +// ============================================================================= +// === qrtest ================================================================== +// ============================================================================= + +template void qrtest +( + cholmod_sparse *A, + double errs [5], + cholmod_common *cc +) +{ + cholmod_sparse *H, *I, *R, *Q, *Csparse, *Xsparse, *AT, *Bsparse ; + cholmod_dense *Cdense, *Xdense, *Bdense, *HTau ; ; + double tol = DBL_EPSILON, err, resid, maxerr, maxresid [2][2] ; + double tols [ ] = { SPQR_DEFAULT_TOL, -1, 0, DBL_EPSILON } ; + Int n, m, nz, *HPinv, ntol, *Ai, *Ap, k, *Qfill, rank, nb, *Cp, *Ci, econ ; + int which ; + Entry *B, *Ax, *Cx ; + int xtype = spqr_type ( ) ; + int ordering ; + Entry range = (Entry) 1.0 ; + + errs [0] = EMPTY ; + errs [1] = EMPTY ; + errs [2] = EMPTY ; + errs [3] = EMPTY ; + errs [4] = EMPTY ; + if (A == NULL) + { + fprintf (stderr, "qrtest: no input matrix\n") ; + return ; + } + + m = A->nrow ; + n = A->ncol ; + Ap = (Int *) A->p ; + Ai = (Int *) A->i ; + Ax = (Entry *) A->x ; + double anorm = spqr_norm_sparse (A, 1, cc) ; + anorm = CHECK_NAN (anorm) ; + printf ("\n===========================================================\n") ; + printf ("Matrix: %d by %d, nnz(A) = %ld, norm(A,1) = %g\n", + (int) m, (int) n, spqr_nnz (A, cc), anorm) ; + printf ( "===========================================================\n") ; + + if (anorm == 0) anorm = 1 ; + + // these should all be zero, no matter what the matrix + maxerr = 0 ; + + // residuals for well-determined or under-determined systems. If + // not rank deficient, these should be zero + maxresid [0][0] = 0 ; // for m <= n, default tol, ordering not fixed + maxresid [0][1] = 0 ; // for m <= n, all other cases + + // residuals for least-squares systems (these will not be zero) + maxresid [1][0] = 0 ; // for m > n, default tol, ordering not fixed + maxresid [1][1] = 0 ; // for m > n, all other cases + + my_srand (m+n) ; + + if (m > 45000) + { + // The only thing returned are the statistics (rank, etc) + printf ("\n=== QR huge (expect int overflow):\n") ; + rank = SPQR_qr ( + 0, 0, 0, 0, A, + NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, + cc, FALSE, FALSE, FALSE) ; + return ; + } + + if (MAX (m,n) >= 300) fprintf (stderr, "please wait ... ") ; + + AT = spqr_transpose (A, 2, cc) ; + + // ------------------------------------------------------------------------- + // test QR routines + // ------------------------------------------------------------------------- + + econ = m ; + + for (ordering = 0 ; ordering <= 9 ; ordering++) + { + + // skip SPQR_ORDERING_GIVEN, unless the matrix is tiny + if (ordering == SPQR_ORDERING_GIVEN && MAX (m,n) > 10) continue ; + + for (ntol = 0 ; ntol < NTOL ; ntol++) + { + + tol = tols [ntol] ; + if // (ntol == 0) // this old test is fragile ... + (tol == SPQR_DEFAULT_TOL) // use this instead + { + // with default tolerance, the fixed ordering can sometimes + // fail if the matrix is rank deficient (R cannot be permuted + // to upper trapezoidal form). + which = (ordering == 0) ; + } + else + { + // with non-default tolerance, the solution can sometimes be + // poor; this is expected. + which = 1 ; + } + printf ("\n=== QR with ordering %d tol %g:\n", ordering, tol) ; + + // ----------------------------------------------------------------- + // create dense and sparse right-hand-sides + // ----------------------------------------------------------------- + + nb = 5 ; + Bdense = spqr_zeros (m, nb, xtype, cc) ; + B = (Entry *) Bdense->x ; + for (k = 0 ; k < m*nb ; k++) + { + B [k] = ((double) (my_rand ( ) % 2)) * erand (range) ; + } + Bsparse = spqr_dense_to_sparse (Bdense, TRUE, cc) ; + + // ----------------------------------------------------------------- + // X = qrsolve(A,B) where X and B are dense + // ----------------------------------------------------------------- + + // X = A\B + if (ordering == SPQR_ORDERING_DEFAULT && tol == SPQR_DEFAULT_TOL) + { + printf ("[ backslach, A and B and X dense: defaults\n") ; + Xdense = SuiteSparseQR (A, Bdense, cc) ; + printf ("] done backslach, A and B and X dense: defaults\n") ; + } + else + { + printf ("[ backslach, A and B and X dense: tol %g order %d\n", + tol, ordering) ; + Xdense = SuiteSparseQR (ordering, tol, A, Bdense, cc) ; + printf ("] done backslach, A B X dense: tol %g order %d\n", + tol, ordering) ; + } + + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid0b %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_dense (&Xdense, cc) ; + + if (cc->useGPU) + { + // error testing for infeasible GPU memory + size_t save = cc->gpuMemorySize ; + cc->gpuMemorySize = 1 ; + printf ("[ Pretend GPU memory is too small:\n") ; + Xdense = SuiteSparseQR (ordering, tol, A, Bdense, cc) ; + cc->gpuMemorySize = save ; + printf ("] test done infeasible GPU, status %2d, useGPU: %d\n", + cc->status, cc->useGPU) ; + spqr_free_dense (&Xdense, cc) ; + } + + spqr_free_dense (&Bdense, cc) ; + + // ----------------------------------------------------------------- + // X = qrsolve(A,B) where X and B are sparse + // ----------------------------------------------------------------- + + // X = A\B + printf ("[ backslash with sparse B: tol %g\n", tol) ; + Xsparse = SuiteSparseQR (ordering, tol, A, Bsparse, cc) ; + printf ("] did tol %g\n", tol) ; + + // check norm (A*x-b), x and b sparse + resid = sparse_resid (A, anorm, Xsparse, Bsparse, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid0 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_sparse (&Xsparse, cc) ; + spqr_free_sparse (&Bsparse, cc) ; + + // ----------------------------------------------------------------- + // X = qrsolve (A,B) where X and B are sparse, with memory test + // ----------------------------------------------------------------- + + // use B = A and solve AX=B where X is sparse + cc->SPQR_shrink = 0 ; // shrink = 0 ; + rank = SPQR_qr ( + ordering, tol, econ, 2, A, + A, NULL, &Xsparse, NULL, NULL, &Qfill, NULL, NULL, NULL, + cc, TRUE, m < 300, nrand (2)) ; + cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; + + if // (ntol == 0) // old test is fragile ... + (tol == SPQR_DEFAULT_TOL) // use this instead. + { + printf ("using default tol: %g\n", cc->SPQR_tol_used) ; + } + + // check norm (A*x-b), x and b sparse + resid = sparse_resid (A, anorm, Xsparse, A, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid1 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_sparse (&Xsparse, cc) ; + spqr_free (n+n, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // X = qrsolve (A,B) where X and B are dense, with memory test + // ----------------------------------------------------------------- + + // use B = dense m-by-2 matrix with some zeros + nb = 5 ; + Bdense = spqr_zeros (m, nb, xtype, cc) ; + B = (Entry *) Bdense->x ; + for (k = 0 ; k < m*nb ; k++) + { + B [k] = (k+1) % 7 ; + } + + rank = SPQR_qr ( + ordering, tol, econ, 2, A, + NULL, Bdense, NULL, &Xdense, NULL, &Qfill, NULL, NULL, NULL, + cc, FALSE, m < 300, nrand (2)) ; + + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid2 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_dense (&Xdense, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // X = qrsolve (A,B) where X and B are full and H is kept + // ----------------------------------------------------------------- + + cc->SPQR_shrink = 2 ; // shrink = 2 ; + rank = SPQR_qr ( + ordering, tol, econ, 2, A, + NULL, Bdense, NULL, &Xdense, NULL, &Qfill, &H, &HPinv, &HTau, + cc, FALSE, m < 300, nrand (2)) ; + cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; + + spqr_free_dense (&HTau, cc) ; + spqr_free (m, sizeof (Int), HPinv, cc) ; + spqr_free_sparse (&H, cc) ; + + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid3 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_dense (&Xdense, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [C,R,E] = qr (A,B) where C is sparse and B is full + // ----------------------------------------------------------------- + + cc->SPQR_shrink = 2 ; // shrink = 2 ; + rank = SPQR_qr ( + ordering, tol, econ, 0, A, + NULL, Bdense, &Csparse, NULL, &R, &Qfill, NULL, NULL, NULL, + cc, FALSE, m < 300, nrand (2)) ; + cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; + + // compute x=R\C and check norm (A*x-b) + Cdense = spqr_sparse_to_dense (Csparse, cc) ; + resid = check_rc (rank, R, A, B, Cdense, nb, anorm, Qfill, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid4 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + spqr_free_dense (&Cdense, cc) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err1: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&Csparse, cc) ; + spqr_free_sparse (&R, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [C,R,E] = qr (A,B) where C and B are full + // ----------------------------------------------------------------- + + cc->SPQR_shrink = 0 ; // shrink = 0 ; + rank = SPQR_qr ( + ordering, tol, econ, 0, A, + NULL, Bdense, NULL, &Cdense, &R, &Qfill, NULL, NULL, NULL, + cc, FALSE, m < 300, nrand (2)) ; + cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err2: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_dense (&Cdense, cc) ; + spqr_free_sparse (&R, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [C,R,E] = qr (A,B) where C and B are full, simple wrapper + // ----------------------------------------------------------------- + + SuiteSparseQR (ordering, tol, econ, A, Bdense, + &Cdense, &R, &Qfill, cc) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err3: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_dense (&Cdense, cc) ; + spqr_free_sparse (&R, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [C,R,E] = qr (A,B) where C and B are sparse, simple wrapper + // ----------------------------------------------------------------- + + Bsparse = spqr_dense_to_sparse (Bdense, TRUE, cc) ; + + SuiteSparseQR (ordering, tol, econ, A, Bsparse, + &Csparse, &R, &Qfill, cc) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err4: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&Csparse, cc) ; + spqr_free_sparse (&Bsparse, cc) ; + spqr_free_sparse (&R, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [CT,R,E] = qr (A,B), but do not return R + // ----------------------------------------------------------------- + + rank = SPQR_qr ( + ordering, tol, econ, 1, A, + NULL, Bdense, &Csparse, NULL, NULL, &Qfill, NULL, NULL, NULL, + cc, FALSE, m < 300, nrand (2)) ; + + spqr_free_sparse (&Csparse, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [Q,R,E] = qr (A), Q in Householder form + // ----------------------------------------------------------------- + + rank = SPQR_qr ( + ordering, tol, econ, -1, A, + NULL, NULL, NULL, NULL, &R, &Qfill, &H, &HPinv, &HTau, + cc, FALSE, m < 300, nrand (2)) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err5: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + // solve Ax=b using Householder form + resid = QRsolve (A, anorm, rank, 0, H, HTau, HPinv, R, + Qfill, Bdense, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid5 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + // solve Ax=b using Q matrix form + resid = QRsolve (A, anorm, rank, 1, H, HTau, HPinv, R, + Qfill, Bdense, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid6 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_dense (&HTau, cc) ; + spqr_free_sparse (&H, cc) ; + spqr_free (m, sizeof (Int), HPinv, cc) ; + spqr_free (n, sizeof (Int), Qfill, cc) ; + spqr_free_sparse (&R, cc) ; + + // ----------------------------------------------------------------- + // [Q,R,E] = qr (A), non-economy + // ----------------------------------------------------------------- + + cc->SPQR_shrink = 0 ; // shrink = 0 ; + I = spqr_speye (m, m, xtype, cc) ; + rank = SPQR_qr ( + ordering, tol, m, 1, A, + I, NULL, &Q, NULL, &R, &Qfill, NULL, NULL, NULL, + cc, FALSE, m<300, nrand (2)) ; + cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; + + // ensure norm (Q*R - A*E) is small + err = check_qr (Q, R, A, Qfill, anorm, cc) ; + printf ("order %d : Q*R-A*E Err6: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&I, cc) ; + spqr_free_sparse (&R, cc) ; + spqr_free_sparse (&Q, cc) ; + spqr_free (n+m, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [Q,R,E] = qr (A), non-economy, using simple wrapper + // ----------------------------------------------------------------- + + if (nrand (2)) + { + // use C version + if (sizeof (Int) == sizeof (int64_t)) + { + SuiteSparseQR_C_QR (ordering, tol, m, A, &Q, &R, + (int64_t **) &Qfill, cc) ; + } + else + { + SuiteSparseQR_i_C_QR (ordering, tol, m, A, &Q, &R, + (int32_t **) &Qfill, cc) ; + } + } + else + { + // use C++ version + SuiteSparseQR (ordering, tol, m, A, &Q, &R, + &Qfill, cc) ; + } + + // ensure norm (Q*R - A*E) is small + err = check_qr (Q, R, A, Qfill, anorm, cc) ; + printf ("order %d : Q*R-A*E Err7: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&R, cc) ; + spqr_free_sparse (&Q, cc) ; + spqr_free (n+m, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [R,E] = qr (A) + // ----------------------------------------------------------------- + + rank = SPQR_qr ( + ordering, tol, econ, 0, A, + NULL, NULL, NULL, NULL, &R, &Qfill, NULL, NULL, NULL, + cc, FALSE, m < 300, nrand (2)) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err8: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + Int rank1 = rank ; + + spqr_free_sparse (&R, cc) ; + spqr_free (n, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [R,E] = qr (A) using simple wrapper + // ----------------------------------------------------------------- + + SuiteSparseQR (ordering, tol, econ, A, &R, &Qfill, cc) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err9: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&R, cc) ; + spqr_free (n, sizeof (Int), Qfill, cc) ; + + // ----------------------------------------------------------------- + // [ ] = qr (A) + // ----------------------------------------------------------------- + + // The only thing returned are the statistics (rank, etc) + cc->SPQR_shrink = 0 ; // shrink = 0 ; + rank = SPQR_qr ( + ordering, tol, econ, 0, A, + NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, + cc, FALSE, m < 300, nrand (2)) ; + cc->SPQR_shrink = 1 ; // restore default shrink = 1 ; + + err = (rank != rank1) ; + printf ("order %d : rank %6d %6d Err10: %g\n", ordering, + (int) rank, (int) rank1, err) ; + maxerr = MAX (maxerr, err) ; + + // ----------------------------------------------------------------- + // [C,H,R,E] = qr (A) + // ----------------------------------------------------------------- + + rank = SPQR_qr ( + ordering, tol, econ, 0, A, + NULL, Bdense, &Csparse, NULL, &R, &Qfill, &H, &HPinv, &HTau, + cc, FALSE, m < 300, nrand (2)) ; + + // compute x=R\C and check norm (A*x-b) + Cdense = spqr_sparse_to_dense (Csparse, cc) ; + resid = check_rc (rank, R, A, B, Cdense, nb, anorm, Qfill, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("Resid7 %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + spqr_free_dense (&Cdense, cc) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err11: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&Csparse, cc) ; + spqr_free_sparse (&R, cc) ; + + // compare Q with qmult + err = check_qmult (H, HTau, HPinv, + ordering == 2 && /* ntol == 0 */ tol == SPQR_DEFAULT_TOL, cc) ; + printf ("order %d : check qmult Err12: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_dense (&HTau, cc) ; + spqr_free_sparse (&H, cc) ; + spqr_free (m, sizeof (Int), HPinv, cc) ; + spqr_free (n+nb, sizeof (Int), Qfill, cc) ; + spqr_free_dense (&Bdense, cc) ; + + // ----------------------------------------------------------------- + // [H,R,E] = qr (A), simple wrapper + // ----------------------------------------------------------------- + + SuiteSparseQR (ordering, tol, econ, A, + &R, &Qfill, &H, &HPinv, &HTau, cc) ; + + // check that R'*R = (A*E)'*(A*E) + err = check_r_factor (R, A, Qfill, cc) ; + printf ("order %d : R'R-(A*E)'*(A*E), Err13: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + // compare Q with qmult + err = check_qmult (H, HTau, HPinv, FALSE, cc) ; + printf ("order %d : check qmult Err14: %g\n", ordering, err) ; + maxerr = MAX (maxerr, err) ; + + spqr_free_sparse (&R, cc) ; + spqr_free_dense (&HTau, cc) ; + spqr_free_sparse (&H, cc) ; + spqr_free (m, sizeof (Int), HPinv, cc) ; + spqr_free (n, sizeof (Int), Qfill, cc) ; + +#ifndef NEXPERT + + // ================================================================= + // === expert routines ============================================= + // ================================================================= + + SuiteSparseQR_factorization *QR ; + cholmod_dense *XT, *Zdense, *BT, *Ydense ; + cholmod_sparse *Ysparse ; + + // ----------------------------------------------------------------- + // QR = qr (A), then solve + // ----------------------------------------------------------------- + + for (int split = 0 ; split <= 4 ; split++) + { + + QR = SPQR_factorize (ordering, tol, A, cc, + split, m < 300, nrand (2)) ; + + // split == 4 does not use singletons, so it can fail if + // rank < n + int wh = which || (split == 4) ; + + // solve Ax=b + nb = 5 ; + Bdense = spqr_zeros (m, nb, xtype, cc) ; + B = (Entry *) Bdense->x ; + for (k = 0 ; k < m*nb ; k++) + { + B [k] = erand (range) ; + } + + // Y = Q'*B + Ydense = SPQR_qmult (SPQR_QTX, QR, Bdense, cc, + m < 300, nrand (2)) ; + + // X = R\(E*Y) + Xdense = SPQR_solve (SPQR_RETX_EQUALS_B, QR, + Ydense, cc, m < 300, nrand (2)) ; + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; + maxresid [m>n][wh] = MAX (maxresid [m>n][wh], resid) ; + printf ("Resid8_%d %d %d %d : %g (%d) tol %g\n", + split, m>n, (int) ntol, ordering, resid, wh, tol) ; + spqr_free_dense (&Xdense, cc) ; + spqr_free_dense (&Ydense, cc) ; + + // Y = (B'*Q)' + BT = transpose (Bdense, cc) ; + XT = SPQR_qmult (SPQR_XQ, QR, BT, cc, + m < 300, nrand (2)) ; + + Ydense = transpose (XT, cc) ; + spqr_free_dense (&XT, cc) ; + spqr_free_dense (&BT, cc) ; + + // X = R\(E*Y) + Xdense = SPQR_solve (SPQR_RETX_EQUALS_B, QR, + Ydense, cc, m < 300, nrand (2)) ; + // check norm (A*x-b), x and b dense + resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; + maxresid [m>n][wh] = MAX (maxresid [m>n][wh], resid) ; + printf ("Resid9_%d %d %d %d : %g (%d) tol %g\n", + split, m>n, (int) ntol, ordering, resid, wh, tol) ; + spqr_free_dense (&Xdense, cc) ; + spqr_free_dense (&Ydense, cc) ; + + // ------------------------------------------------------------- + // error testing + // ------------------------------------------------------------- + + if (ordering == 0 && /* ntol == 0 */ tol == SPQR_DEFAULT_TOL) + { + printf ("Error handling ... expect 3 error messages: \n") ; + err = (SuiteSparseQR_qmult (-1, QR, Bdense, cc) + != NULL) ; + spqr_free_dense (&Bdense, cc) ; + Bdense = spqr_zeros (m+1, 1, xtype, cc) ; + err += (SuiteSparseQR_qmult (SPQR_QX,QR,Bdense,cc) + != NULL); + spqr_free_dense (&Bdense, cc) ; + Bdense = spqr_zeros (1, m+1, xtype, cc) ; + err += (SuiteSparseQR_qmult (SPQR_XQ,QR,Bdense,cc) + != NULL); + if (QR->n1cols > 0) + { + // this will fail; cannot refactorize with singletons + printf ("Error handling ... expect error message:\n") ; + err += (SuiteSparseQR_numeric + (tol, A, QR, cc) != FALSE) ; + } + printf ("order %d : error handling Err15: %g\n", + ordering, err) ; + maxerr = MAX (maxerr, err) ; + printf (" ... error handling done\n\n") ; + } + + spqr_free_dense (&Bdense, cc) ; + + // ------------------------------------------------------------- + + // solve A'x=b + nb = 5 ; + Bdense = spqr_zeros (n, nb, xtype, cc) ; + B = (Entry *) Bdense->x ; + for (k = 0 ; k < n*nb ; k++) + { + B [k] = erand (range) ; + } + // Y = R'\(E'*B) + Ydense = SPQR_solve (SPQR_RTX_EQUALS_ETB, QR, + Bdense, cc, m < 300, nrand (2)) ; + // X = Q*Y + Xdense = SPQR_qmult (SPQR_QX, QR, Ydense, cc, + m < 300, nrand (2)) ; + // check norm (A'*x-b), x and b dense + resid = dense_resid (AT, anorm, Xdense, nb, B, cc) ; + maxresid [m (&Xdense, cc) ; + spqr_free_dense (&Ydense, cc) ; + + // ------------------------------------------------------------- + // error testing + // ------------------------------------------------------------- + + if (!split && ordering == 0 && /* ntol == 0 */ + tol == SPQR_DEFAULT_TOL) + { + printf ("Error testing ... expect 3 error messages:\n") ; + err = (SuiteSparseQR_solve (-1, QR, Bdense, cc) + != NULL) ; + spqr_free_dense (&Bdense, cc) ; + err += (SuiteSparseQR_solve + (SPQR_RTX_EQUALS_ETB, QR, Bdense, cc) != NULL) ; + Bdense = spqr_zeros (n+1, 1, xtype, cc) ; + err += (SuiteSparseQR_solve (SPQR_RTX_EQUALS_ETB, + QR, Bdense, cc) != NULL) ; + printf ("order %d : error handling Err16: %g\n", + ordering, err) ; + maxerr = MAX (maxerr, err) ; + printf (" ... error handling done\n\n") ; + } + + SuiteSparseQR_free (&QR, cc) ; + spqr_free_dense (&Bdense, cc) ; + } + + // ----------------------------------------------------------------- + // QR = qr (A'), then solve + // ----------------------------------------------------------------- + + // use qmult to solve min-2-norm problem + QR = SuiteSparseQR_factorize (ordering, tol, AT, cc) ; + + nb = 5 ; + Bdense = spqr_zeros (m, nb, xtype, cc) ; + B = (Entry *) Bdense->x ; + for (k = 0 ; k < m*nb ; k++) + { + B [k] = erand (range) ; + } + + // solve X = R'\B + Xdense = SPQR_solve (SPQR_RTX_EQUALS_B, QR, Bdense, cc, + m < 300, nrand (2)) ; + spqr_free_dense (&Xdense, cc) ; + + // solve X = R'\(E'*B) + Xdense = SPQR_solve (SPQR_RTX_EQUALS_ETB, QR, Bdense, + cc, m < 300, nrand (2)) ; + + // Y = Q*X + Ydense = SPQR_qmult (SPQR_QX, QR, Xdense, cc, + m < 300, nrand (2)) ; + + // check norm (A*y-b), y and b dense + resid = dense_resid (A, anorm, Ydense, nb, B, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("ResidB %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + spqr_free_dense (&Ydense, cc) ; + + // Y = (X'*Q')' + XT = transpose (Xdense, cc) ; + Zdense = SPQR_qmult (SPQR_XQT, QR, XT, cc, + m < 300, nrand (2)) ; + Ydense = transpose (Zdense, cc) ; + spqr_free_dense (&XT, cc) ; + spqr_free_dense (&Zdense, cc) ; + + // check norm (A*y-b), y and b dense + resid = dense_resid (A, anorm, Ydense, nb, B, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("ResidC %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + spqr_free_dense (&Ydense, cc) ; + spqr_free_dense (&Xdense, cc) ; + + // ----------------------------------------------------------------- + // min 2-norm solution using min2norm + // ----------------------------------------------------------------- + + Xdense = SPQR_min2norm (ordering, tol, A, Bdense, cc, + m < 300, nrand (2)) ; + + // check norm (A*x-b), y and b dense + resid = dense_resid (A, anorm, Xdense, nb, B, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("ResidD %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + spqr_free_dense (&Xdense, cc) ; + + spqr_free_dense (&Bdense, cc) ; + + // ----------------------------------------------------------------- + // sparse case + // ----------------------------------------------------------------- + + nb = 5 ; + Bdense = spqr_zeros (m, nb, xtype, cc) ; + B = (Entry *) Bdense->x ; + for (k = 0 ; k < m*nb ; k++) + { + B [k] = ((double) (my_rand ( ) % 2)) * erand (range) ; + } + Bsparse = spqr_dense_to_sparse (Bdense, TRUE, cc) ; + spqr_free_dense (&Bdense, cc) ; + + // solve X = R'\(E'*B) + Xsparse = NULL ; + Xsparse = SPQR_solve (SPQR_RTX_EQUALS_ETB, QR, + Bsparse, cc, m < 300, nrand (2)) ; + + // Y = Q*X + Ysparse = NULL ; + Ysparse = SPQR_qmult (SPQR_QX, QR, Xsparse, cc, + m < 300, nrand (2)) ; + + // check norm (A*y-b), y and b sparse + resid = sparse_resid (A, anorm, Ysparse, Bsparse, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("ResidE %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_sparse (&Xsparse, cc) ; + spqr_free_sparse (&Ysparse, cc) ; + + Xsparse = SPQR_min2norm (ordering, tol, A, Bsparse, cc, + m < 300, nrand (2)) ; + + // check norm (A*x-b), x and b sparse + resid = sparse_resid (A, anorm, Xsparse, Bsparse, cc) ; + maxresid [m>n][which] = MAX (maxresid [m>n][which], resid) ; + printf ("ResidF %d %d %d : %g\n", m>n, (int) ntol, + ordering, resid) ; + + spqr_free_sparse (&Xsparse, cc) ; + spqr_free_sparse (&Bsparse, cc) ; + + SuiteSparseQR_free (&QR, cc) ; +#endif + } + } + + // ------------------------------------------------------------------------- + // check error handling + // ------------------------------------------------------------------------- + + printf ("Check error handling, one error message is expected:\n") ; + cholmod_dense *Bgunk = spqr_ones (m+1, 1, xtype, cc) ; + rank = SuiteSparseQR ( + 0, 0, econ, -1, A, + NULL, Bgunk, NULL, NULL, NULL, NULL, NULL, NULL, NULL, + cc) ; + spqr_free_dense (&Bgunk, cc) ; + err = (rank != EMPTY) ; + maxerr = MAX (maxerr, err) ; // rank should be EMPTY + printf (" ... error handling done\n\n") ; + + // ------------------------------------------------------------------------- + // test non-user callable functions + // ------------------------------------------------------------------------- + + // attempt to permute A to upper triangular form + Int *Qtrap ; + rank = spqr_trapezoidal (n, Ap, Ai, Ax, 0, NULL, FALSE, &Cp, &Ci, &Cx, + &Qtrap, cc) ; + printf ("Rank of A, if A*P permutable to upper trapezoidal: %d\n", + (int) rank) ; + if (Cp != NULL) + { + nz = Cp [n] ; + spqr_free (n+1, sizeof (Int), Cp, cc) ; + spqr_free (nz, sizeof (Int), Ci, cc) ; + spqr_free (nz, sizeof (Entry), Cx, cc) ; + spqr_free (n, sizeof (Int), Qtrap, cc) ; + } + spqr_free_sparse (&AT, cc) ; + + // ------------------------------------------------------------------------- + // test the C API + // ------------------------------------------------------------------------- + + qrtest_C (A, anorm, errs, maxresid, cc) ; + + // ------------------------------------------------------------------------- + // final results + // ------------------------------------------------------------------------- + + errs [0] = CHECK_NAN (maxerr) ; + errs [1] = CHECK_NAN (maxresid [0][0]) ; + errs [2] = CHECK_NAN (maxresid [0][1]) ; + errs [3] = CHECK_NAN (maxresid [1][0]) ; + errs [4] = CHECK_NAN (maxresid [1][1]) ; +} + +// ============================================================================= +// === do_matrix2 ============================================================== +// ============================================================================= + +// Read in a matrix, and use it to test SuiteSparseQR +// If kind == 0, then the first two residuals should be low. + +template +int do_matrix2 (int kind, cholmod_sparse *A, cholmod_common *cc) +{ + double errs [5] = {0,0,0,0,0} ; + Int m = A->nrow ; + Int n = A->ncol ; + + // ------------------------------------------------------------------------- + // use it to test SuiteSparseQR + // ------------------------------------------------------------------------- + + if (A->xtype == CHOLMOD_COMPLEX && A->stype == 0) + { + qrtest (A, errs, cc) ; + } + else if (A->xtype == CHOLMOD_REAL) + { + if (A->stype != 0) + { + cholmod_sparse *A1 ; + A1 = spqr_copy (A, 0, 1, cc) ; + qrtest (A1, errs, cc) ; + spqr_free_sparse (&A1, cc) ; + } + else + { + qrtest (A, errs, cc) ; + } + } + else + { + // cannot handle ZOMPLEX, PATTERN, or symmetric/Hermitian COMPLEX + fprintf (stderr, "invalid matrix\n") ; + errs [0] = 1 ; + } + + // ------------------------------------------------------------------------- + // report the results + // ------------------------------------------------------------------------- + + if (kind == 0) + { + printf ("First Resid and ") ; + } + printf ("Err should be low:\n") ; + printf ("RESULT: Err %8.1e Resid %8.1e %8.1e", errs [0], + errs [1], errs [2]) ; + if (m == n) + { + printf (" ") ; + } + else + { + printf (" %8.1e %8.1e", errs [3], errs [4]) ; + } + + if (errs [0] > 1e-10) + { + printf (" : FAIL\n") ; + fprintf (stderr, "Error: %g FAIL\n", errs [0]) ; + return (1) ; + } + + // if kind == 0, then this full-rank matrix should have low residual + if (kind == 0 && (errs [1] > 1e-10)) + { + printf (" : FAIL\n") ; + fprintf (stderr, "error: %g FAIL\n", errs [1]) ; + return (1) ; + } + + printf (" : OK.") ; + fprintf (stderr, "OK.") ; + return (0) ; +} + +// ============================================================================= +// === do_matrix =============================================================== +// ============================================================================= + +// Read in a matrix, and use it to test SuiteSparseQR +// If kind == 0, then the first two residuals should be low. + +template +int do_matrix (int kind, FILE *file, cholmod_common *cc) +{ + cholmod_sparse *A ; + + int nfail0 = 0 ; + int nfail1 = 0 ; + int nfail2 = 0 ; + int nfail3 = 0 ; + + // ------------------------------------------------------------------------- + // read in the matrix + // ------------------------------------------------------------------------- + + A = spqr_read_sparse (file, cc) ; + if (A == NULL) + { + fprintf (stderr, "Unable to read matrix\n") ; + return (1) ; + } + Int m = A->nrow ; + Int n = A->ncol ; + fprintf (stderr, "%5d by %5d : ", (int) m, (int) n) ; + if (sizeof (int64_t) > sizeof (int) && (m > 10000 || n > 10000)) + { + fprintf (stderr, "(test skipped on 64-bit systems)\n") ; + spqr_free_sparse (&A, cc) ; + return (0) ; + } + + // defaults + cc->SPQR_grain = 1 ; // no parallel analysis + printf ("\nBeginning CPU tests [\n") ; + fprintf (stderr, " CPU ") ; + nfail0 = do_matrix2 (kind, A, cc) ; + + // non-defaults (will not use the GPU) + cc->SPQR_grain = 4 ; // grain size relative to total work + nfail2 = do_matrix2 (kind, A, cc) ; + cc->SPQR_grain = 1 ; // no parallel analysis + printf ("\nCPU tests done ]\n") ; + + // test the GPU, if installed + #ifdef SUITESPARSE_CUDA + cc->useGPU = TRUE ; + // was 3.5 * ((size_t) 1024 * 1024 * 1024) ; + size_t totmem, availmem ; + double t = SuiteSparse_time ( ) ; + spqr_gpu_memorysize (&totmem, &availmem, cc) ; + t = SuiteSparse_time ( ) - t ; + cc->gpuMemorySize = availmem ; + printf ("\nBeginning GPU tests, GPU memory %g MB warmup time %g[\n", + (double) (cc->gpuMemorySize) / (1024*1024), t) ; + fprintf (stderr, " GPU ") ; + nfail1 = do_matrix2 (kind, A, cc) ; + printf ("\nGPU tests done ]\n") ; + if (m > 200) + { + // try with a tiny GPU memory size, but only for a few matrices + // in the test set. Each front will go in its own stage. + printf ("\nBeginning GPU tests with tiny GPU memory [\n") ; + cc->gpuMemorySize = 0 ; + nfail3 = do_matrix2 (kind, A, cc) ; + // restore defaults + cc->useGPU = FALSE ; + printf ("\nGPU tests done (tiny memory) ]\n") ; + } + #endif + + spqr_free_sparse (&A, cc) ; + + printf ("\n") ; + fprintf (stderr, "\n") ; + return (nfail0 + nfail1 + nfail2 + nfail3) ; +} + diff --git a/SPQR/Tcov/qrtestc.c b/SPQR/Tcov/qrtestc.c index e9e152baff..7e1a95bfd1 100644 --- a/SPQR/Tcov/qrtestc.c +++ b/SPQR/Tcov/qrtestc.c @@ -14,221 +14,25 @@ #define MAX(a,b) (((a) > (b)) ? (a) : (b)) //------------------------------------------------------------------------------ -// qrtest_C64: int64_t version +// qrtest_c64: int64_t version //------------------------------------------------------------------------------ -static void qrtest_C64 -( - 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} ; - int64_t m, n ; -#ifndef NEXPERT - cholmod_dense *Y ; - int split ; -#endif - - m = A->nrow ; - n = A->ncol ; - - B = cholmod_l_ones (m, 1, A->xtype, cc) ; - - /* X = A\B */ - X = SuiteSparseQR_C_backslash_default (A, B, cc) ; - - /* Resid = A*X - B */ - Resid = cholmod_l_copy_dense (B, cc) ; - cholmod_l_sdmult (A, 0, one, minusone, X, Resid, cc) ; - - /* resid = norm (Resid,1) */ - resid = cholmod_l_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ; - resid = (resid < 0 || resid != resid) ? 9e99 : resid ; - cholmod_l_free_dense (&Resid, cc) ; - cholmod_l_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_l_dense_to_sparse (B, 1, cc) ; - Xsparse = SuiteSparseQR_C_backslash_sparse (2, -2, A, Bsparse, cc) ; - X = cholmod_l_sparse_to_dense (Xsparse, cc) ; - cholmod_l_free_sparse (&Bsparse, cc) ; - cholmod_l_free_sparse (&Xsparse, cc) ; - - /* Resid = A*X - B */ - Resid = cholmod_l_copy_dense (B, cc) ; - cholmod_l_sdmult (A, 0, one, minusone, X, Resid, cc) ; - - /* resid = norm (Resid,1) */ - resid = cholmod_l_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ; - resid = (resid < 0 || resid != resid) ? 9e99 : resid ; - cholmod_l_free_dense (&Resid, cc) ; - cholmod_l_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_l_copy_dense (B, cc) ; - cholmod_l_sdmult (A, 0, one, minusone, X, Resid, cc) ; - - /* resid = norm (Resid,1) */ - resid = cholmod_l_norm_dense (Resid, 1, cc) / MAX (anorm, 1) ; - resid = (resid < 0 || resid != resid) ? 9e99 : resid ; - cholmod_l_free_dense (&Resid, cc) ; - cholmod_l_free_dense (&X, cc) ; - - maxresid [m>n][0] = MAX (maxresid [m>n][0], resid) ; - printf ("Resid_C3 %d : %g\n", m>n, resid) ; - - cholmod_l_free_dense (&Y, cc) ; - SuiteSparseQR_C_free (&QR, cc) ; - } -#endif - - cholmod_l_free_dense (&B, cc) ; -} +#define DLONG +#define QRTESTC qrtest_c64 +#include "qrtestc_template.c" +#undef DLONG +#undef QRTESTC //------------------------------------------------------------------------------ -// qrtest_C32: int32_t version +// 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) ; -} +#define DINT +#define QRTESTC qrtest_c32 +#include "qrtestc_template.c" //------------------------------------------------------------------------------ -// qrtest_C: handles both int32_t and int64_t versions +// qrtest_C: handles both int32_t and int64_t versions (both double, not float) //------------------------------------------------------------------------------ void qrtest_C @@ -242,11 +46,11 @@ void qrtest_C { if (A->itype == CHOLMOD_INT) { - qrtest_C32 (A, anorm, errs, maxresid, cc) ; + qrtest_c32 (A, anorm, errs, maxresid, cc) ; } else { - qrtest_C64 (A, anorm, errs, maxresid, cc) ; + qrtest_c64 (A, anorm, errs, maxresid, cc) ; } } diff --git a/SPQR/Tcov/qrtestc_template.c b/SPQR/Tcov/qrtestc_template.c new file mode 100644 index 0000000000..b88d12bce3 --- /dev/null +++ b/SPQR/Tcov/qrtestc_template.c @@ -0,0 +1,127 @@ +/* ========================================================================== */ +/* === qrtestc_template ===================================================== */ +/* ========================================================================== */ + +// SPQR, Copyright (c) 2008-2023, Timothy A Davis. All Rights Reserved. +// SPDX-License-Identifier: GPL-2.0+ + +//------------------------------------------------------------------------------ +// QRTESTC: test the SPQR C API +//------------------------------------------------------------------------------ + +// The #include'ing file must define the following macros: +// QRTEST: the name of the function + +// and one of: +// SINT: float (and float complex), int32 +// SLONG: float (and float complex), int64 +// DINT: double (and double complex), int32 +// DLONG: double (and double complex), int64 + +// the SINT and SLONG cases are not yet used. + +#include "cholmod_types.h" + +static void QRTESTC +( + cholmod_sparse *A, + Real anorm, + Real errs [5], + Real maxresid [2][2], + cholmod_common *cc +) +{ + cholmod_dense *B, *X, *Resid ; + cholmod_sparse *Bsparse, *Xsparse ; + SuiteSparseQR_C_factorization *QR ; + Real resid, one [2] = {1,0}, minusone [2] = {-1,0} ; + Int 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)) ; +} +