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