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