CMSCOPE changed

1 files changed, 1111 insertions, 0 deletions

diff --git a/usr/include/spqr.hpp b/usr/include/spqr.hpp
new file mode 100755
index 000000000..c88c8c28a
--- /dev/null
+++ b/usr/include/spqr.hpp
@@ -0,0 +1,1111 @@
+// =============================================================================
+// === 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->other1 [0] += (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
+
+// #ifndef NPRINT
+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, cholmod_common *cc) ;
+void spqrDebug_print (Complex x, cholmod_common *cc) ;
+void spqrDebug_printf (double x, cholmod_common *cc) ;
+void spqrDebug_printf (Complex x, cholmod_common *cc) ;
+// #endif
+
+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  // for cc->hypotenuse function
+) ;
+
+template <typename Entry> int spqr_rmap
+(
+    SuiteSparseQR_factorization <Entry> *QR,
+    cholmod_common *cc
+) ;
+
+
+// =============================================================================
+// === 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 (cc->hypotenuse (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 ;
+    cc->complex_divide (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