1function [Q,R,P,info] = spqr (A,arg2,arg3)                                  %#ok
2%SPQR multithreaded multifrontal rank-revealing sparse QR.
3%For a sparse m-by-n matrix A, and sparse or full B of size m-by-k:
4%
5%   R = spqr (A)              Q-less QR
6%   R = spqr (A,0)            economy variant (size(R,1) = min(m,n))
7%   R = spqr (A,opts)         as above, with non-default options
8%
9%   [Q,R] = spqr (A)          A=Q*R
10%   [Q,R] = spqr (A,0)        economy variant (size(Q,2) = size(R,1) = min(m,n))
11%   [Q,R] = spqr (A,opts)     A=Q*R, with non-default options
12%
13%   [Q,R,P] = spqr (A)        A*P=Q*R where P reduces fill-in
14%   [Q,R,P] = spqr (A,0)      economy variant (size(Q,2) = size(R,1) = min(m,n))
15%   [Q,R,P] = spqr (A,opts)   as above, with non-default options
16%
17%   [C,R] = spqr (A,B)        as R=spqr(A), also returns C=Q'*B
18%   [C,R] = spqr (A,B,0)      economy variant (size(C,1) = size(R,1) = min(m,n))
19%   [C,R] = spqr (A,B,opts)   as above, with non-default options
20%
21%   [C,R,P] = spqr (A,B)      as R=spqr(A*P), also returns C=Q'*B
22%   [C,R,P] = spqr (A,B,0)    economy variant (size(C,1) = size(R,1) = min(m,n))
23%   [C,R,P] = spqr (A,B,opts) as above, with non-default options
24%
25% P is chosen to reduce fill-in and to return R in upper trapezoidal form if A
26% is estimated to have less than full rank.  opts provides non-default options.
27% Q can be optionally returned in Householder form, which is far sparser than
28% returning Q as a sparse matrix.
29%
30% With 4 output arguments, [Q,R,P,info]=spqr(...) or [C,R,P,info]=spqr(...),
31% a struct "info" is returned with statistics about the QR factorization.  The
32% contents of info are mostly self-explanatory, except for info.norm_E_fro.
33% This is equal to the Frobenius norm of E where E=A*P-Q*R.  If
34% info.norm_E_fro <= info.tol, then this guarantees that the true numerical
35% rank is no larger than the rank r returned by SPQR (r=info.rank_A_estimate).
36% If in addition the smallest singular value of R(1:r,1:r) is larger than
37% info.tol, then info.rank_A_estimate is the true numerical rank of A.  Note
38% that info.norm_E_fro is returned as zero if A is determined by SPQR to have
39% full rank (r = min(m,n) where [m n]=size(A)) or if opts.tol=0.
40%
41% Example:
42%   The least-squares solution of an overdetermined system A*x=b with
43%   m > n can be found in at least one of seven ways (in increasing order of
44%   efficiency):
45%
46%      x = pinv(full(A)) * b ;
47%      [Q,R] = spqr (A) ; x = R\(Q'*b) ;
48%      [Q,R,P] = spqr (A) ; x = P*(R\(Q'*b)) ;
49%      [Q,R,P] = spqr (A,struct('Q','Householder')) ; x=P*(R\spqr_qmult(Q,b,0));
50%      [c,R,P] = spqr (A,b) ; x = P*(R\c) ;
51%      [c,R,p] = spqr (A,b,0) ; x = (R\c) ; x (p) = x ;
52%      x = spqr_solve (A,b) ;
53%
54%   The minimum-norm solution of an underdetermined system A*x=b with
55%   m < n can be found in one of five ways (in increasing order of efficiency):
56%
57%      x = pinv(full(A)) * b ;
58%      [Q,R] = spqr (A') ; x = Q*(R'\b) ;
59%      [Q,R,P] = spqr (A') ; x = Q*(R'\(P'*b)) ;
60%      [Q,R,P] = spqr(A',struct('Q','Householder'));x=spqr_qmult(Q,R'\(P'*b),1);
61%      x = spqr_solve (A,b,struct('solution','min2norm')) ;
62%
63% Entries not present in opts are set to their defaults:
64%
65%   opts.tol:   columns with norm <= tol are treated as zero. The default is
66%       20 * (m+n) * eps * sqrt(max(diag(A'*A))).  Returned as info->tol.
67%       If opts.tol=0, no nonzero columns are treated as zero.  If opts.tol>=0,
68%       and P is present on output, R is returned in upper trapezoidal form
69%       R(1:r,1:r) is upper triangular with zero-free diagonal and where
70%       r=info.rank_A_estimate, unless this is not possible due to structural
71%       rank deficiency.  If opts.tol<0, R is not permuted into this form.
72%
73%   opts.econ:  number of rows of R and columns of Q to return.  m is the
74%   default.  n gives the standard economy form.  A value less than the
75%   estimated rank r is set to r, so opts.econ = 0 gives the "rank-sized"
76%   factorization, where size(R,1) == nnz(diag(R)) == r.
77%
78%   opts.ordering: a string describing what ordering method to use.  Let [m n]
79%   = size (S) where S is obtained by removing singletons from A.  'default':
80%   the default ordering: COLAMD(S).  'amd': AMD(S'*S). 'colamd': COLAMD(S)
81%   'metis': METIS(S'*S), only if METIS is installed. 'best': try all three
82%   (AMD, COLAMD, METIS) and take the best 'bestamd': try AMD and COLAMD and
83%   take the best. 'fixed': P=I; this is the only option if P is not present in
84%   the output. 'natural': singleton removal only.  The singleton pre-ordering
85%   permutes A prior to factorization into the form [A11 A12 ; 0 A22] where A11
86%   is upper triangular with all(abs(diag(A11)) > opts.tol) (see
87%   spqr_singletons).
88%
89%   opts.Q: a string describing how Q is returned.  The default is 'discard' if
90%   Q is not present in the output, or 'matrix' otherwise.  If Q is present and
91%   opts.Q is 'discard', then Q=[] is returned (thus R=spqr(A*P) is
92%   [Q,R,P]=spqr(A) where spqr finds P and Q is discarded instead). 'matrix'
93%   returns Q as a sparse matrix where A=Q*R or A*P=Q*R.  'Householder' returns
94%   Q as a struct containing the Householder reflections applied to A to obtain
95%   R, resulting in a far sparser Q than the 'matrix' option.
96%
97%   opts.permutation: a string describing how P is to be returned.  The default
98%   is 'matrix', so that A*P=Q*R.  'vector' gives A(:,P)=Q*R instead.
99%
100%   opts.spumoni: acts just like spparms('spumoni',k).
101%
102%   opts.grain, opts.small, opts.nthreads: multitasking control (if compiled
103%   with TBB); the scheduler tries to ensure that all parallel tasks have at
104%   least max (total_flops / opts.grain, opts.small) flops.  No TBB parallelism
105%   is exploited if opts.grain = 1.  opts.nthreads gives the number of threads
106%   to use for TBB (which is different than the number of threads used by the
107%   BLAS).  opts.nthreads <= 0 means to let TBB determine the number of threads
108%   (normally equal to the number of cores); otherwise, use exactly
109%   opts.nthreads threads.  Defaults: 1, 1e6, and 0, respectively.  TBB is
110%   disabled by default since it conflicts with BLAS multithreading.  If you
111%   enable TBB, be sure to disable BLAS multithreading with
112%   maxNumCompThreads(1), or choose opts.nthreads * (number of BLAS threads)
113%   equal to the number of cores.  A good value of opts.grain is twice that of
114%   opts.nthreads.  If TBB parallelism is enabled, the METIS ordering normally
115%   gives the best speedup for large problems.
116%
117%   opts.solution: used by spqr_solve; 'basic' (default), or 'min2norm'.
118%   Determines the kind of solution that spqr_solve computes for
119%   underdetermined systems.  Has no effect for least-squares problems; ignored
120%   by spqr itself.
121%
122% See also SPQR_QMULT, SPQR_SOLVE, LU, NULL, ORTH, QRDELETE, QRINSERT,
123% QRUPDATE, SPQR_SINGLETONS.
124
125% Copyright 2008, Timothy A. Davis, http://www.suitesparse.com
126
127error ('spqr mexFunction not found') ;
128