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