1function [N, stats, stats_ssp, est_sval_upper_bounds] = ...
2    spqr_rank_form_basis(call_from, A, U, V ,Q1, rank_spqr, numerical_rank, ...
3    stats, opts, est_sval_upper_bounds, nsvals_small, nsvals_large,  p1, Q2, p2)
4%SPQR_RANK_FORM_BASIS forms the basis for the null space of a matrix.
5%
6% Called from spqr_basic and spqr_cod after these routines call spqr and
7% spqr_rank_ssi.  The input parameters are as used in spqr_basic and spqr_cod.
8% The parameter call_from indicates the type of call:
9%      call_from = 1 -- a call from spqr_basic, form null space basis for A'
10%      call_from = 2 -- a call from spqr_cod, form null space basis for A
11%      call_from = 3 -- a call from spqr_cod, form null space basis for A'
12% Output:
13%   N -- basis for the null space in implicit or explicit form as required
14%      by opts
15%   stats -- an update of the stats vector in spqr_basic or spqr_cod
16%   stats_ssp -- information about the call to spqr_ssp
17%   est_sval_upper_bounds -- an update to the estimated singular value upper
18%      bounds
19% Not user-callable.
20
21% Copyright 2012, Leslie Foster and Timothy A Davis.
22
23get_details = opts.get_details ;
24implicit_null_space_basis = opts.implicit_null_space_basis ;
25start_with_A_transpose = opts.start_with_A_transpose ;
26if (get_details == 1)
27    t = tic ;
28end
29[m , n] = size(A);
30nullity_R11 = rank_spqr - numerical_rank;
31
32%-------------------------------------------------------------------------------
33% form X which contains basis for null space of trapezoidal matrix
34%-------------------------------------------------------------------------------
35
36if call_from == 1    % call from spqr_basic (for null space basis for A')
37    if nullity_R11 == 0
38        X = [sparse(rank_spqr, m-rank_spqr) ; ...
39             speye(m-rank_spqr, m-rank_spqr)];
40    else
41        X = [ [sparse(U(:, end-nullity_R11+1:end)) ; ...
42             sparse(m-rank_spqr, nullity_R11) ], ...
43             [ sparse(rank_spqr, m-rank_spqr) ; ...
44             speye(m-rank_spqr, m-rank_spqr)]];
45    end
46elseif call_from == 2    % call from spqr_cod (for null space basis of A)
47    if nullity_R11 == 0
48        X = [sparse(rank_spqr, n-rank_spqr) ; ...
49             speye(n-rank_spqr, n-rank_spqr)];
50    else
51        if (start_with_A_transpose)
52            X = [ [sparse(V(:, end-nullity_R11+1:end)) ; ...
53                 sparse(n-rank_spqr, nullity_R11) ], ...
54                 [ sparse(rank_spqr, n-rank_spqr) ; ...
55                 speye(n-rank_spqr, n-rank_spqr)]];
56        else
57            X = [ [sparse(U(:, end-nullity_R11+1:end)) ; ...
58                 sparse(n-rank_spqr, nullity_R11) ], ...
59                 [ sparse(rank_spqr, n-rank_spqr) ; ...
60                 speye(n-rank_spqr, n-rank_spqr)]];
61        end
62    end
63elseif call_from == 3  % call from spqr_cod (for null space basis for A')
64    if nullity_R11 == 0
65        X = [sparse(rank_spqr, m-rank_spqr) ; ...
66          speye(m-rank_spqr, m-rank_spqr)];
67    else
68        if (start_with_A_transpose)
69            X = [ [sparse(U(:, end-nullity_R11+1:end)) ; ...
70              sparse(m-rank_spqr, nullity_R11) ], ...
71              [ sparse(rank_spqr, m-rank_spqr) ; ...
72              speye(m-rank_spqr, m-rank_spqr)]];
73        else
74            X = [ [sparse(V(:, end-nullity_R11+1:end)) ; ...
75              sparse(m-rank_spqr, nullity_R11) ], ...
76              [ sparse(rank_spqr, m-rank_spqr) ; ...
77              speye(m-rank_spqr, m-rank_spqr)]];
78        end
79    end
80end
81
82%-------------------------------------------------------------------------------
83% form null space basis for A or A' using X and Q from QR factorization
84%-------------------------------------------------------------------------------
85
86if call_from == 1    % call from spqr_basic (for null space basis for A')
87
88    if implicit_null_space_basis
89        % store null space in implicit form (store Q and X in N = Q*X)
90        N.Q = Q1;
91        N.X = X;
92        N.kind = 'Q*X' ;
93    else
94        % store null space in explicit form
95        N = spqr_qmult(Q1,X,1);
96    end
97
98 elseif call_from == 2    % call from spqr_cod (for null space basis of A)
99
100    if implicit_null_space_basis
101        % store null space basis in implicit form
102        if (start_with_A_transpose)
103            % store P, Q and X in N = Q*P*X
104            p = 1:n;
105            p(p2) = 1:length(p2);
106            N.Q = Q1;
107            N.P = p;    % a permutation vector
108            N.X = X;
109            N.kind = 'Q*P*X' ;
110        else
111            % store P, Q and X in N = P*Q*X
112            p(p1) = 1:length(p1);
113            N.P = p;    % a permutation vector
114            N.Q = Q2;
115            N.X = X;
116            N.kind = 'P*Q*X' ;
117        end
118    else
119        % store null space basis as an explicit matrix
120        if (start_with_A_transpose)
121            p = 1:n;
122            p(p2) = 1:length(p2);
123            N = spqr_qmult(Q1,X(p,:),1);
124        else
125            N = spqr_qmult(Q2,X,1);
126            p(p1) = 1:length(p1);
127            N = N(p,:);
128        end
129    end
130
131elseif call_from == 3  % call from spqr_cod (for null space basis for A')
132
133    p = 1:m;
134    if (start_with_A_transpose)
135        p(p1) = 1:length(p1);
136    else
137        p(p2) = 1:length(p2);
138    end
139
140    if implicit_null_space_basis
141        % store null space basis in implicit form
142        if (start_with_A_transpose)
143            %          (store P, Q and X in N = P*Q*X)
144            N.P = p;    % a permutation vector
145            N.Q = Q2;
146            N.X = X;
147            N.kind = 'P*Q*X' ;
148        else
149            %       (store P, Q and X in N = Q*P*X)
150            N.Q = Q1;
151            N.P = p;    % a permutation vector
152            N.X = X;
153            N.kind = 'Q*P*X' ;
154        end
155    else
156        % store null space explicitly forming Q*P*X;
157        if (start_with_A_transpose)
158            N = spqr_qmult(Q2,X,1);
159            N = N(p,:) ;
160        else
161            N = spqr_qmult(Q1,X(p,:),1);
162        end
163    end
164end
165
166if (get_details == 1)
167    stats.time_basis = toc (t) ;
168end
169
170%-------------------------------------------------------------------------------
171% call spqr_ssp to enhance, potentially, est_sval_upper_bounds, the estimated
172%    upper bounds on the singular values
173% and / or
174%    to estimate ||A * N || or || A' * N ||
175%-------------------------------------------------------------------------------
176
177% Note: nullity = m - numerical_rank;   % the number of columns in X and N
178
179if call_from == 1    % call from spqr_basic (for null space basis for A')
180
181    % Note: opts.k is not the same as k in the algorithm description above.
182
183    % Note that, by the interleave theorem for singular values, for
184    % i=1:nullity, singular value i of A'*N will be an upper bound on singular
185    % value numerical_rank + i of A.  S(i,i) is an estimate of singular value i
186    % of A'*N with an estimated accuracy of stats_ssp.est_error_bounds(i).
187    % Therefore let
188
189    [s_ssp,stats_ssp] = spqr_ssp (A', N, max (nsvals_small, opts.k), opts) ;
190
191 elseif call_from == 2    % call from spqr_cod (for null space basis of A)
192
193    % Note that, by the interleave theorem for singular
194    % values, for i = 1, ..., nullity, singular value i of A*N will be
195    % an upper bound on singular value numerical_rank + i of A.
196    % S(i,i) is an estimate of singular value i of A*N with an estimated
197    % accuracy of stats_ssp.est_error_bounds(i).  Therefore let
198
199    [s_ssp,stats_ssp] = spqr_ssp (A, N, max (nsvals_small, opts.k), opts) ;
200
201end
202
203if call_from == 1 || call_from == 2 % call from spqr_basic (for null space
204               %  basis for A') or from spqr_cod (for null space basis of A)
205
206    % By the comments prior to the call to spqr_ssp we have
207    % s_ssp + stats_ssp.est_error_bounds are estimated upper bounds
208    % for singular values (numerical_rank+1):(numerical_rank+nsvals_small)
209    % of A.  We have two estimates for upper bounds on these singular
210    % values of A.  Choose the smaller of the two:
211    if ( nsvals_small > 0 )
212        est_sval_upper_bounds(nsvals_large+1:end) = ...
213            min( est_sval_upper_bounds(nsvals_large+1:end) , ...
214               s_ssp(1:nsvals_small)' + ...
215               stats_ssp.est_error_bounds(1:nsvals_small) );
216    end
217
218elseif call_from == 3  % call from spqr_cod (for null space basis for A')
219
220    % call ssp to estimate nsvals_small sing. values of A' * N
221    %    (useful to estimate the norm of A' * N)
222    [ignore,stats_ssp] = spqr_ssp (A', N, nsvals_small, opts) ;             %#ok
223    clear ignore
224
225end
226
227if (get_details == 1)
228    if call_from == 1 || call_from == 3
229        stats.stats_ssp_NT = stats_ssp ;
230    elseif call_from == 2
231        stats.stats_ssp_N = stats_ssp ;
232    end
233    stats.time_svd = stats.time_svd + stats_ssp.time_svd ;
234end
235end
236