File: //usr/include/suitesparse/spqr.hpp
// =============================================================================
// === spqr.hpp ================================================================
// =============================================================================
// Internal definitions and non-user-callable routines. This should not be
// included in the user's code.
#ifndef SPQR_INTERNAL_H
#define SPQR_INTERNAL_H
// -----------------------------------------------------------------------------
// include files
// -----------------------------------------------------------------------------
#include "SuiteSparseQR.hpp"
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <stdio.h>
#include <cstring>
#include <complex>
typedef std::complex<double> Complex ;
// -----------------------------------------------------------------------------
// debugging and printing control
// -----------------------------------------------------------------------------
// force debugging off
#ifndef NDEBUG
#define NDEBUG
#endif
// force printing off
#ifndef NPRINT
#define NPRINT
#endif
// uncomment the following line to turn on debugging (SPQR will be slow!)
/*
#undef NDEBUG
*/
// uncomment the following line to turn on printing (LOTS of output!)
/*
#undef NPRINT
*/
// uncomment the following line to turn on expensive debugging (very slow!)
/*
#define DEBUG_EXPENSIVE
*/
// -----------------------------------------------------------------------------
// Long is defined at SuiteSparse_long, from SuiteSparse_config.h
// -----------------------------------------------------------------------------
#define Long SuiteSparse_long
// -----------------------------------------------------------------------------
// basic macros
// -----------------------------------------------------------------------------
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define EMPTY (-1)
#define TRUE 1
#define FALSE 0
#define IMPLIES(p,q) (!(p) || (q))
// NULL should already be defined, but ensure it is here.
#ifndef NULL
#define NULL ((void *) 0)
#endif
// column-major indexing; A[i,j] is A (INDEX (i,j,lda))
#define INDEX(i,j,lda) ((i) + ((j)*(lda)))
// FLIP is a "negation about -1", and is used to mark an integer i that is
// normally non-negative. FLIP (EMPTY) is EMPTY. FLIP of a number > EMPTY
// is negative, and FLIP of a number < EMTPY is positive. FLIP (FLIP (i)) = i
// for all integers i. UNFLIP (i) is >= EMPTY.
#define EMPTY (-1)
#define FLIP(i) (-(i)-2)
#define UNFLIP(i) (((i) < EMPTY) ? FLIP (i) : (i))
// -----------------------------------------------------------------------------
// additional include files
// -----------------------------------------------------------------------------
#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif
#define ITYPE CHOLMOD_LONG
#define DTYPE CHOLMOD_DOUBLE
#define ID SuiteSparse_long_id
// -----------------------------------------------------------------------------
#define ERROR(status,msg) \
cholmod_l_error (status, __FILE__, __LINE__, msg, cc)
// Check a pointer and return if null. Set status to invalid, unless the
// status is already "out of memory"
#define RETURN_IF_NULL(A,result) \
{ \
if ((A) == NULL) \
{ \
if (cc->status != CHOLMOD_OUT_OF_MEMORY) \
{ \
ERROR (CHOLMOD_INVALID, NULL) ; \
} \
return (result) ; \
} \
}
// Return if Common is NULL or invalid
#define RETURN_IF_NULL_COMMON(result) \
{ \
if (cc == NULL) \
{ \
return (result) ; \
} \
if (cc->itype != ITYPE || cc->dtype != DTYPE) \
{ \
cc->status = CHOLMOD_INVALID ; \
return (result) ; \
} \
}
#define RETURN_IF_XTYPE_INVALID(A,result) \
{ \
if (A->xtype != xtype) \
{ \
ERROR (CHOLMOD_INVALID, "invalid xtype") ; \
return (result) ; \
} \
}
// -----------------------------------------------------------------------------
// debugging and printing macros
// -----------------------------------------------------------------------------
#ifndef NDEBUG
#ifdef MATLAB_MEX_FILE
// #define ASSERT(e) mxAssert (e, "error: ")
extern char spqr_mx_debug_string [200] ;
char *spqr_mx_id (int line) ;
#define ASSERT(e) \
((e) ? (void) 0 : \
mexErrMsgIdAndTxt (spqr_mx_id (__LINE__), \
"assert: (" #e ") file:" __FILE__ ))
#else
#include <assert.h>
#define ASSERT(e) assert (e)
#endif
#define DEBUG(e) e
#ifdef DEBUG_EXPENSIVE
#define DEBUG2(e) e
#define ASSERT2(e) ASSERT(e)
#else
#define DEBUG2(e)
#define ASSERT2(e)
#endif
#else
#define ASSERT(e)
#define ASSERT2(e)
#define DEBUG(e)
#define DEBUG2(e)
#endif
#ifndef NPRINT
#ifdef MATLAB_MEX_FILE
#define PR(e) mexPrintf e
#else
#define PR(e) printf e
#endif
#define PRVAL(e) spqrDebug_print (e)
#else
#define PR(e)
#define PRVAL(e)
#endif
// -----------------------------------------------------------------------------
// For counting flops; disabled if TBB is used
// -----------------------------------------------------------------------------
#define FLOP_COUNT(f) { if (cc->SPQR_grain <= 1) cc->SPQR_flopcount += (f) ; }
// =============================================================================
// === spqr_work ===============================================================
// =============================================================================
// workspace required for each stack in spqr_factorize and spqr_kernel
template <typename Entry> struct spqr_work
{
Long *Stair1 ; // size maxfn if H not kept
Long *Cmap ; // size maxfn
Long *Fmap ; // size n
Entry *WTwork ; // size (fchunk + (keepH ? 0:1)) * maxfn
Entry *Stack_head ; // head of Stack
Entry *Stack_top ; // top of Stack
Long sumfrank ; // sum of ranks of the fronts in this stack
Long maxfrank ; // largest rank of fronts in this stack
// for computing the 2-norm of w, the vector of the dead column norms
double wscale ; // scale factor for norm (w (of this stack))
double wssq ; // sum-of-squares for norm (w (of this stack))
} ;
// =============================================================================
// === spqr_blob ===============================================================
// =============================================================================
// The spqr_blob is a collection of objects that the spqr_kernel requires.
template <typename Entry> struct spqr_blob
{
double tol ;
spqr_symbolic *QRsym ;
spqr_numeric <Entry> *QRnum ;
spqr_work <Entry> *Work ;
Long *Cm ;
Entry **Cblock ;
Entry *Sx ;
Long ntol ;
Long fchunk ;
cholmod_common *cc ;
} ;
// =============================================================================
// === SuiteSparseQR non-user-callable functions ===============================
// =============================================================================
spqr_symbolic *spqr_analyze
(
// inputs, not modified
cholmod_sparse *A,
int ordering, // all ordering options available
Long *Quser, // user provided ordering, if given (may be NULL)
int do_rank_detection, // if TRUE, then rank deficient matrices may be
// considered during numerical factorization,
// with tol >= 0 (tol < 0 is also allowed). If FALSE, then the tol
// parameter is ignored by the numerical factorization, and no rank
// detection is performed.
int keepH, // if nonzero, H is kept
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> spqr_numeric <Entry> *spqr_factorize
(
// input, optionally freed on output
cholmod_sparse **Ahandle,
// inputs, not modified
Long freeA, // if TRUE, free A on output
double tol, // for rank detection
Long ntol, // apply tol only to first ntol columns
spqr_symbolic *QRsym,
// workspace and parameters
cholmod_common *cc
) ;
// returns tol (-1 if error)
template <typename Entry> double spqr_tol
(
// inputs, not modified
cholmod_sparse *A,
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> double spqr_maxcolnorm
(
// inputs, not modified
cholmod_sparse *A,
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> void spqr_kernel
(
Long task,
spqr_blob <Entry> *Blob
) ;
template <typename Entry> void spqr_parallel
(
Long ntasks,
int nthreads,
spqr_blob <Entry> *Blob
) ;
void spqr_freesym
(
spqr_symbolic **QRsym_handle,
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> void spqr_freenum
(
spqr_numeric <Entry> **QRnum_handle,
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> void spqr_freefac
(
SuiteSparseQR_factorization <Entry> **QR_handle,
// workspace and parameters
cholmod_common *cc
) ;
void spqr_stranspose1
(
// input, not modified
cholmod_sparse *A, // m-by-n
Long *Qfill, // size n, fill-reducing column permutation;
// Qfill [k] = j if the kth column of S is the jth
// column of A. Identity permutation is used if
// Qfill is NULL.
// output, contents not defined on input
Long *Sp, // size m+1, row pointers of S
Long *Sj, // size nz, column indices of S
Long *PLinv, // size m, inverse row permutation, PLinv [i] = k
Long *Sleft, // size n+2, Sleft [j] ... Sleft [j+1]-1 is the list of
// rows of S whose leftmost column index is j. The list
// can be empty (that is, Sleft [j] == Sleft [j+1]).
// Sleft [n] is the number of non-empty rows of S, and
// Sleft [n+1] is always m. That is, Sleft [n] ...
// Sleft [n+1]-1 gives the empty rows of S.
// workspace, not defined on input or output
Long *W // size m
) ;
template <typename Entry> void spqr_stranspose2
(
// input, not modified
cholmod_sparse *A, // m-by-n
Long *Qfill, // size n, fill-reducing column permutation;
// Qfill [k] = j
// if the kth column of S is the jth column of A.
// Identity permutation is used if Qfill is NULL.
Long *Sp, // size m+1, row pointers of S
Long *PLinv, // size m, inverse row permutation, PLinv [i] = k
// output, contents not defined on input
Entry *Sx, // size nz, numerical values of S
// workspace, not defined on input or output
Long *W // size m
) ;
// =============================================================================
#ifndef NDEBUG
template <typename Entry> void spqrDebug_dumpdense
(
Entry *A,
Long m,
Long n,
Long lda,
cholmod_common *cc
) ;
template <typename Entry> void spqrDebug_dumpsparse
(
Long *Ap,
Long *Ai,
Entry *Ax,
Long m,
Long n,
cholmod_common *cc
) ;
void spqrDebug_print (double x) ;
void spqrDebug_print (Complex x) ;
void spqrDebug_dump_Parent (Long n, Long *Parent, const char *filename) ;
Long spqrDebug_rhsize // returns # of entries in R+H block
(
// input, not modified
Long m, // # of rows in F
Long n, // # of columns in F
Long npiv, // number of pivotal columns in F
Long *Stair, // size n; column j is dead if Stair [j] == 0.
// Only the first npiv columns can be dead.
cholmod_common *cc
) ;
#endif
#ifdef DEBUG_EXPENSIVE
Long spqrDebug_listcount
(
Long x, Long *List, Long len, Long what,
cholmod_common *cc
) ;
#endif
// =============================================================================
Long spqr_fsize // returns # of rows of F
(
// inputs, not modified
Long f,
Long *Super, // size nf, from QRsym
Long *Rp, // size nf, from QRsym
Long *Rj, // size rjsize, from QRsym
Long *Sleft, // size n+2, from QRsym
Long *Child, // size nf, from QRsym
Long *Childp, // size nf+1, from QRsym
Long *Cm, // size nf, from QRwork
// outputs, not defined on input
Long *Fmap, // size n, from QRwork
Long *Stair // size fn, from QRwork
) ;
template <typename Entry> void spqr_assemble
(
// inputs, not modified
Long f, // front to assemble F
Long fm, // number of rows of F
int keepH, // if TRUE, then construct row pattern of H
Long *Super,
Long *Rp,
Long *Rj,
Long *Sp,
Long *Sj,
Long *Sleft,
Long *Child,
Long *Childp,
Entry *Sx,
Long *Fmap,
Long *Cm,
Entry **Cblock,
#ifndef NDEBUG
char *Rdead,
#endif
Long *Hr,
// input/output
Long *Stair,
Long *Hii, // if keepH, construct list of row indices for F
// input only
Long *Hip,
// outputs, not defined on input
Entry *F,
// workspace, not defined on input or output
Long *Cmap
) ;
template <typename Entry> Long spqr_cpack // returns # of rows in C
(
// input, not modified
Long m, // # of rows in F
Long n, // # of columns in F
Long npiv, // number of pivotal columns in F
Long g, // the C block starts at F (g,npiv)
// input, not modified unless the pack occurs in-place
Entry *F, // m-by-n frontal matrix in column-major order
// output, contents not defined on input
Entry *C // packed columns of C, of size cm-by-cn in upper
// trapezoidal form.
) ;
Long spqr_fcsize // returns # of entries in C of current front F
(
// input, not modified
Long m, // # of rows in F
Long n, // # of columns in F
Long npiv, // number of pivotal columns in F
Long g // the C block starts at F (g,npiv)
) ;
Long spqr_csize // returns # of entries in C of a child
(
// input, not modified
Long c, // child c
Long *Rp, // size nf+1, pointers for pattern of R
Long *Cm, // size nf, Cm [c] = # of rows in child C
Long *Super // size nf, pivotal columns in each front
) ;
template <typename Entry> void spqr_rcount
(
// inputs, not modified
spqr_symbolic *QRsym,
spqr_numeric <Entry> *QRnum,
Long n1rows, // added to each row index of Ra and Rb
Long econ, // only get entries in rows n1rows to econ-1
Long n2, // Ra = R (:,0:n2-1), Rb = R (:,n2:n-1)
int getT, // if true, count Rb' instead of Rb
// input/output
Long *Ra, // size n2; Ra [j] += nnz (R (:,j)) if j < n2
Long *Rb, // If getT is false: size n-n2 and
// Rb [j-n2] += nnz (R (:,j)) if j >= n2.
// If getT is true: size econ, and
// Rb [i] += nnz (R (i, n2:n-1))
Long *Hp, // size rjsize+1. Column pointers for H.
// Only computed if H was kept during factorization.
// Only Hp [0..nh] is used.
Long *p_nh // number of Householder vectors (nh <= rjsize)
) ;
template <typename Entry> void spqr_rconvert
(
// inputs, not modified
spqr_symbolic *QRsym,
spqr_numeric <Entry> *QRnum,
Long n1rows, // added to each row index of Ra, Rb, and H
Long econ, // only get entries in rows n1rows to econ-1
Long n2, // Ra = R (:,0:n2-1), Rb = R (:,n2:n-1)
int getT, // if true, get Rb' instead of Rb
// input/output
Long *Rap, // size n2+1; on input, Rap [j] is the column pointer
// for Ra. Incremented on output by the number of
// entries added to column j of Ra.
// output, not defined on input
Long *Rai, // size rnz1 = nnz(Ra); row indices of Ra
Entry *Rax, // size rnz; numerical values of Ra
// input/output
Long *Rbp, // if getT is false:
// size (n-n2)+1; on input, Rbp [j] is the column
// pointer for Rb. Incremented on output by the number
// of entries added to column j of Rb.
// if getT is true:
// size econ+1; on input, Rbp [i] is the row
// pointer for Rb. Incremented on output by the number
// of entries added to row i of Rb.
// output, not defined on input
Long *Rbi, // size rnz2 = nnz(Rb); indices of Rb
Entry *Rbx, // size rnz2; numerical values of Rb
// input
Long *H2p, // size nh+1; H2p [j] is the column pointer for H.
// H2p, H2i, and H2x are ignored if H was not kept
// during factorization. nh computed by rcount
// output, not defined on input
Long *H2i, // size hnz = nnz(H); indices of H
Entry *H2x, // size hnz; numerical values of H
Entry *H2Tau // size nh; Householder coefficients
) ;
template <typename Entry> Long spqr_rhpack // returns # of entries in R+H
(
// input, not modified
int keepH, // if true, then H is packed
Long m, // # of rows in F
Long n, // # of columns in F
Long npiv, // number of pivotal columns in F
Long *Stair, // size npiv; column j is dead if Stair [j] == 0.
// Only the first npiv columns can be dead.
// input, not modified (unless the pack occurs in-place)
Entry *F, // m-by-n frontal matrix in column-major order
// output, contents not defined on input
Entry *R, // packed columns of R+H
Long *p_rm // number of rows in R block
) ;
template <typename Entry> void spqr_hpinv
(
// input
spqr_symbolic *QRsym,
// input/output
spqr_numeric <Entry> *QRnum,
// workspace
Long *W // size QRnum->m
) ;
template <typename Entry> int spqr_1colamd
(
// inputs, not modified
int ordering, // all available, except 0:fixed and 3:given
// treated as 1:natural
double tol, // only accept singletons above tol
Long bncols, // number of columns of B
cholmod_sparse *A, // m-by-n sparse matrix
// output arrays, neither allocated nor defined on input.
Long **p_Q1fill, // size n+bncols, fill-reducing
// or natural ordering
Long **p_R1p, // size n1rows+1, R1p [k] = # of nonzeros in kth
// row of R1. NULL if n1cols == 0.
Long **p_P1inv, // size m, singleton row inverse permutation.
// If row i of A is the kth singleton row, then
// P1inv [i] = k. NULL if n1cols is zero.
cholmod_sparse **p_Y, // on output, only the first n-n1cols+1 entries of
// Y->p are defined (if Y is not NULL), where
// Y = [A B] or Y = [A2 B2]. If B is empty and
// there are no column singletons, Y is NULL
Long *p_n1cols, // number of column singletons found
Long *p_n1rows, // number of corresponding rows found
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> int spqr_1fixed
(
// inputs, not modified
double tol, // only accept singletons above tol
Long bncols, // number of columns of B
cholmod_sparse *A, // m-by-n sparse matrix
// output arrays, neither allocated nor defined on input.
Long **p_R1p, // size n1rows+1, R1p [k] = # of nonzeros in kth
// row of R1. NULL if n1cols == 0.
Long **p_P1inv, // size m, singleton row inverse permutation.
// If row i of A is the kth singleton row, then
// P1inv [i] = k. NULL if n1cols is zero.
cholmod_sparse **p_Y, // on output, only the first n-n1cols+1 entries of
// Y->p are defined (if Y is not NULL), where
// Y = [A B] or Y = [A2 B2]. If B is empty and
// there are no column singletons, Y is NULL
Long *p_n1cols, // number of column singletons found
Long *p_n1rows, // number of corresponding rows found
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> SuiteSparseQR_factorization <Entry> *spqr_1factor
(
// inputs, not modified
int ordering, // all ordering options available
double tol, // only accept singletons above tol
Long bncols, // number of columns of B
int keepH, // if TRUE, keep the Householder vectors
cholmod_sparse *A, // m-by-n sparse matrix
Long ldb, // leading dimension of B, if dense
Long *Bp, // size bncols+1, column pointers of B
Long *Bi, // size bnz = Bp [bncols], row indices of B
Entry *Bx, // size bnz, numerical values of B
// workspace and parameters
cholmod_common *cc
) ;
Long spqr_cumsum // returns total sum
(
// input, not modified
Long n,
// input/output
Long *X // size n+1. X = cumsum ([0 X])
) ;
void spqr_shift
(
// input, not modified
Long n,
// input/output
Long *X // size n+1
) ;
template <typename Entry> void spqr_larftb
(
// inputs, not modified (V is modified and then restored on output)
int method, // 0,1,2,3
Long m, // C is m-by-n
Long n,
Long k, // V is v-by-k
// for methods 0 and 1, v = m,
// for methods 2 and 3, v = n
Long ldc, // leading dimension of C
Long ldv, // leading dimension of V
Entry *V, // V is v-by-k, unit lower triangular (diag not stored)
Entry *Tau, // size k, the k Householder coefficients
// input/output
Entry *C, // C is m-by-n, with leading dimension ldc
// workspace, not defined on input or output
Entry *W, // for methods 0,1: size k*k + n*k
// for methods 2,3: size k*k + m*k
cholmod_common *cc
) ;
int spqr_happly_work
(
// input
int method, // 0,1,2,3
Long m, // X is m-by-n
Long n,
// FUTURE : make H cholmod_sparse:
Long nh, // number of Householder vectors
Long *Hp, // size nh+1, column pointers for H
Long hchunk,
// outputs; sizes of workspaces needed
Long *p_vmax,
Long *p_vsize,
Long *p_csize
) ;
template <typename Entry> void spqr_happly
(
// input
int method, // 0,1,2,3
Long m, // X is m-by-n
Long n,
Long nh, // number of Householder vectors
Long *Hp, // size nh+1, column pointers for H
Long *Hi, // size hnz = Hp [nh], row indices of H
Entry *Hx, // size hnz, Householder values. Note that the first
// entry in each column must be equal to 1.0
Entry *Tau, // size nh
// input/output
Entry *X, // size m-by-n with leading dimension m
// workspace
Long vmax,
Long hchunk,
Long *Wi, // size vmax
Long *Wmap, // size MAX(mh,1) where H is mh-by-nh
Entry *C, // size csize
Entry *V, // size vsize
cholmod_common *cc
) ;
template <typename Entry> void spqr_panel
(
// input
int method,
Long m,
Long n,
Long v,
Long h, // number of Householder vectors in the panel
Long *Vi, // Vi [0:v-1] defines the pattern of the panel
Entry *V, // v-by-h, panel of Householder vectors
Entry *Tau, // size h, Householder coefficients for the panel
Long ldx,
// input/output
Entry *X, // m-by-n with leading dimension ldx
// workspace
Entry *C, // method 0,1: v-by-n; method 2,3: m-by-v
Entry *W, // method 0,1: k*k+n*k; method 2,3: k*k+m*k
cholmod_common *cc
) ;
template <typename Entry> int spqr_append // TRUE if OK, FALSE otherwise
(
// inputs, not modified
Entry *X, // size m-by-1
Long *P, // size m, or NULL; permutation to apply to X.
// P [k] = i if row k of A is row i of X
// input/output
cholmod_sparse *A, // size m-by-n2 where n2 > n
Long *p_n, // number of columns of A; increased by one
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> Long spqr_trapezoidal // rank of R, or EMPTY
(
// inputs, not modified
Long n, // R is m-by-n (m is not needed here; can be economy R)
Long *Rp, // size n+1, column pointers of R
Long *Ri, // size rnz = Rp [n], row indices of R
Entry *Rx, // size rnz, numerical values of R
Long bncols, // number of columns of B
Long *Qfill, // size n+bncols, fill-reducing ordering. Qfill [k] = j if
// the jth column of A is the kth column of R. If Qfill is
// NULL, then it is assumed to be the identity
// permutation.
int skip_if_trapezoidal, // if R is already in trapezoidal form,
// and skip_if_trapezoidal is TRUE, then
// the matrix T is not created.
// outputs, not allocated on input
Long **p_Tp, // size n+1, column pointers of T
Long **p_Ti, // size rnz, row indices of T
Entry **p_Tx, // size rnz, numerical values of T
Long **p_Qtrap, // size n+bncols, modified Qfill
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> int spqr_type (void) ;
template <typename Entry> void spqr_rsolve
(
// inputs
SuiteSparseQR_factorization <Entry> *QR,
int use_Q1fill,
Long nrhs, // number of columns of B
Long ldb, // leading dimension of B
Entry *B, // size m-by-nrhs with leading dimesion ldb
// output
Entry *X, // size n-by-nrhs with leading dimension n
// workspace
Entry **Rcolp,
Long *Rlive,
Entry *W,
cholmod_common *cc
) ;
// returns rank of F, or 0 on error
template <typename Entry> Long spqr_front
(
// input, not modified
Long m, // F is m-by-n with leading dimension m
Long n,
Long npiv, // number of pivot columns
double tol, // a column is flagged as dead if its norm is <= tol
Long ntol, // apply tol only to first ntol pivot columns
Long fchunk, // block size for compact WY Householder reflections,
// treated as 1 if fchunk <= 1
// input/output
Entry *F, // frontal matrix F of size m-by-n
Long *Stair, // size n, entries F (Stair[k]:m-1, k) are all zero,
// and remain zero on output.
char *Rdead, // size npiv; all zero on input. If k is dead,
// Rdead [k] is set to 1
// output, not defined on input
Entry *Tau, // size n, Householder coefficients
// workspace, undefined on input and output
Entry *W, // size b*(n+b), where b = min (fchunk,n,m)
// input/output
double *wscale,
double *wssq,
cholmod_common *cc
) ;
template <typename Entry> int spqr_rmap
(
SuiteSparseQR_factorization <Entry> *QR,
cholmod_common *cc
) ;
// =============================================================================
// === spqrgpu features ========================================================
// =============================================================================
#ifdef GPU_BLAS
#include "spqrgpu.hpp"
#endif
// =============================================================================
// === spqr_conj ===============================================================
// =============================================================================
inline double spqr_conj (double x)
{
return (x) ;
}
inline Complex spqr_conj (Complex x)
{
return (std::conj (x)) ;
}
// =============================================================================
// === spqr_abs ================================================================
// =============================================================================
inline double spqr_abs (double x, cholmod_common *cc) // cc is unused
{
return (fabs (x)) ;
}
inline double spqr_abs (Complex x, cholmod_common *cc)
{
return (SuiteSparse_config.hypot_func (x.real ( ), x.imag ( ))) ;
}
// =============================================================================
// === spqr_divide =============================================================
// =============================================================================
inline double spqr_divide (double a, double b, cholmod_common *cc) // cc unused
{
return (a/b) ;
}
inline Complex spqr_divide (Complex a, Complex b, cholmod_common *cc)
{
double creal, cimag ;
SuiteSparse_config.divcomplex_func
(a.real(), a.imag(), b.real(), b.imag(), &creal, &cimag) ;
return (Complex (creal, cimag)) ;
}
// =============================================================================
// === spqr_add ================================================================
// =============================================================================
// Add two non-negative Long's, and return the result. Checks for Long overflow
// and sets ok to FALSE if it occurs.
inline Long spqr_add (Long a, Long b, int *ok)
{
Long c = a + b ;
if (c < 0)
{
(*ok) = FALSE ;
return (EMPTY) ;
}
return (c) ;
}
// =============================================================================
// === spqr_mult ===============================================================
// =============================================================================
// Multiply two positive Long's, and return the result. Checks for Long
// overflow and sets ok to FALSE if it occurs.
inline Long spqr_mult (Long a, Long b, int *ok)
{
Long c = a * b ;
if (((double) c) != ((double) a) * ((double) b))
{
(*ok) = FALSE ;
return (EMPTY) ;
}
return (c) ;
}
// =============================================================================
// === BLAS interface ==========================================================
// =============================================================================
// To compile SuiteSparseQR with 64-bit BLAS, use -DBLAS64. See also
// CHOLMOD/Include/cholmod_blas.h
extern "C" {
#include "cholmod_blas.h"
}
#undef CHECK_BLAS_INT
#undef EQ
#define CHECK_BLAS_INT (sizeof (BLAS_INT) < sizeof (Long))
#define EQ(K,k) (((BLAS_INT) K) == ((Long) k))
#ifdef SUN64
#define BLAS_DNRM2 dnrm2_64_
#define LAPACK_DLARF dlarf_64_
#define LAPACK_DLARFG dlarfg_64_
#define LAPACK_DLARFT dlarft_64_
#define LAPACK_DLARFB dlarfb_64_
#define BLAS_DZNRM2 dznrm2_64_
#define LAPACK_ZLARF zlarf_64_
#define LAPACK_ZLARFG zlarfg_64_
#define LAPACK_ZLARFT zlarft_64_
#define LAPACK_ZLARFB zlarfb_64_
#elif defined (BLAS_NO_UNDERSCORE)
#define BLAS_DNRM2 dnrm2
#define LAPACK_DLARF dlarf
#define LAPACK_DLARFG dlarfg
#define LAPACK_DLARFT dlarft
#define LAPACK_DLARFB dlarfb
#define BLAS_DZNRM2 dznrm2
#define LAPACK_ZLARF zlarf
#define LAPACK_ZLARFG zlarfg
#define LAPACK_ZLARFT zlarft
#define LAPACK_ZLARFB zlarfb
#else
#define BLAS_DNRM2 dnrm2_
#define LAPACK_DLARF dlarf_
#define LAPACK_DLARFG dlarfg_
#define LAPACK_DLARFT dlarft_
#define LAPACK_DLARFB dlarfb_
#define BLAS_DZNRM2 dznrm2_
#define LAPACK_ZLARF zlarf_
#define LAPACK_ZLARFG zlarfg_
#define LAPACK_ZLARFT zlarft_
#define LAPACK_ZLARFB zlarfb_
#endif
// =============================================================================
// === BLAS and LAPACK prototypes ==============================================
// =============================================================================
extern "C"
{
void LAPACK_DLARFT (char *direct, char *storev, BLAS_INT *n, BLAS_INT *k,
double *V, BLAS_INT *ldv, double *Tau, double *T, BLAS_INT *ldt) ;
void LAPACK_ZLARFT (char *direct, char *storev, BLAS_INT *n, BLAS_INT *k,
Complex *V, BLAS_INT *ldv, Complex *Tau, Complex *T, BLAS_INT *ldt) ;
void LAPACK_DLARFB (char *side, char *trans, char *direct, char *storev,
BLAS_INT *m, BLAS_INT *n, BLAS_INT *k, double *V, BLAS_INT *ldv,
double *T, BLAS_INT *ldt, double *C, BLAS_INT *ldc, double *Work,
BLAS_INT *ldwork) ;
void LAPACK_ZLARFB (char *side, char *trans, char *direct, char *storev,
BLAS_INT *m, BLAS_INT *n, BLAS_INT *k, Complex *V, BLAS_INT *ldv,
Complex *T, BLAS_INT *ldt, Complex *C, BLAS_INT *ldc, Complex *Work,
BLAS_INT *ldwork) ;
double BLAS_DNRM2 (BLAS_INT *n, double *X, BLAS_INT *incx) ;
double BLAS_DZNRM2 (BLAS_INT *n, Complex *X, BLAS_INT *incx) ;
void LAPACK_DLARFG (BLAS_INT *n, double *alpha, double *X, BLAS_INT *incx,
double *tau) ;
void LAPACK_ZLARFG (BLAS_INT *n, Complex *alpha, Complex *X, BLAS_INT *incx,
Complex *tau) ;
void LAPACK_DLARF (char *side, BLAS_INT *m, BLAS_INT *n, double *V,
BLAS_INT *incv, double *tau, double *C, BLAS_INT *ldc, double *Work) ;
void LAPACK_ZLARF (char *side, BLAS_INT *m, BLAS_INT *n, Complex *V,
BLAS_INT *incv, Complex *tau, Complex *C, BLAS_INT *ldc, Complex *Work) ;
}
#endif