1function [ x ] = ... 2 spqr_rank_deflation(call_from, R, U, V, C, m, n, rank_spqr, ... 3 numerical_rank, nsvals_large, opts, p1, p2, N, Q1, Q2) 4%SPQR_RANK_DEFLATION constructs pseudoinverse or basic solution using deflation. 5% 6% Called from spqr_basic and spqr_cod after these routines call spqr, 7% spqr_rank_ssi and spqr_rank_form_basis. The input parameters are as used 8% in spqr_basic and spqr_cod. 9% The parameter call_from indicates the type of call: 10% call_from = 1 indicates a call from spqr_basic 11% call_from = 2 indicates a call from spqr_cod 12% Output: 13% x -- for call_from = 1 x is a basic solution to 14% min || b - A * x ||. (1) 15% The basic solution has n - (rank returned by spqr) zero components. 16% For call_from = 2 x is an approximate pseudoinverse solution to (1). 17% Not user-callable. 18 19% Algorithm: R * wh = ch or R' * wh = ch is solved where ch is described 20% in the code and R comes from the QR factorizations in 21% spqr_basic or spqr_cod. R is triangular and potentially 22% numerically singular with left and right singular vectors for 23% small singular values stored in U and V. When R is numerically 24% singular deflation (see SIAM SISC, 11:519-530, 1990) to 25% calculate an approximate truncated singular value solution to 26% the triangular system. Orthogonal transformations 27% are applied to wh to obtain the solutions x to (1). 28 29% Copyright 2012, Leslie Foster and Timothy A Davis. 30 31% disable nearly-singular matrix warnings, and save the current state 32warning_state = warning ('off', 'MATLAB:nearlySingularMatrix') ; 33 34start_with_A_transpose = opts.start_with_A_transpose ; 35implicit_null_space_basis = opts.implicit_null_space_basis ; 36 37if (isempty (C)) 38 39 x = zeros (m,0) ; 40 41elseif (start_with_A_transpose || call_from == 1) 42 43 ch = C(1:rank_spqr,:); 44 if numerical_rank == rank_spqr 45 wh = R \ ch ; 46 else 47 % use deflation (see SIAM SISC, 11:519-530, 1990) to calculate an 48 % approximate truncated singular value solution to R * wh = ch 49 U = U(:,nsvals_large+1:end); 50 wh = ch - U*(U'*ch); 51 wh = R \ wh ; 52 V = V(:,nsvals_large+1:end); 53 wh = wh - V*(V'*wh); 54 end 55 if call_from == 2 56 wh(p2,:)=wh; 57 end 58 wh = [wh ; zeros(n - rank_spqr,size(C,2)) ]; 59 if call_from == 1 60 x(p1,:)=wh; 61 else 62 if implicit_null_space_basis 63 x = spqr_qmult(N.Q,wh,1); 64 else 65 x = spqr_qmult(Q1,wh,1); 66 end 67 end 68 69else 70 71 ch = C(p2,:); 72 if numerical_rank == rank_spqr 73 wh = ( ch' / R )' ; % wh = R' \ ch rewritten to save memory 74 else 75 % use deflation (see SIAM SISC, 11:519-530, 1990) to calculate an 76 % approximate truncated singular value solution to R' * wh = ch 77 V = V(:,nsvals_large+1:end); 78 wh = ch - V*(V'*ch); 79 wh = ( wh' / R )' ; % wh = R' \ ch rewritten to save memory 80 U = U(:,nsvals_large+1:end); 81 wh = wh - U*(U'*wh); 82 end 83 wh = [wh ; zeros(n - rank_spqr,size(C,2)) ]; 84 if implicit_null_space_basis 85 wh = spqr_qmult(N.Q,wh,1); 86 else 87 wh = spqr_qmult(Q2,wh,1); 88 end 89 x(p1,:)=wh; 90 91end 92 93% restore the warning back to what it was 94warning (warning_state) ; 95 96