1function [ stats ] = spqr_rank_assign_stats(... 2 call_from, est_sval_upper_bounds, est_sval_lower_bounds, tol, ... 3 numerical_rank, nsvals_small, nsvals_large, stats, stats_ssi, ... 4 opts, nargout_call, stats_ssp_N, stats_ssp_NT, start_tic) 5%SPQR_RANK_ASSIGN_STATS set flag and other statistics. 6% 7% Called by spqr_basic and spqr_cod as well as, indirectly, by spqr_pinv and 8% spqr_null after these methods have determined estimated bounds on the 9% singular values of A. Set the flag which indicates the success of 10% and other statistics for spqr_basic, spqr_null, spqr_pinv, spqr_cod. 11% Not user callable. 12% 13% Output: 14% stats - statistics for spqr_basic, spqr_null, spqr_pinv and spqr_cod. 15% including the following fields: 16% 17% stats.flag which indicates whether the method has determined the 18% correct numerical rank: 19% stats.flag is 0 if it is likely, although not 20% guaranteed, that stats.rank is the correct numerical rank for 21% tolerance tol (i.e. agrees with the numerical rank 22% determined by the singular values of R). 23% stats.flag is 1 if the calculated numerical rank stats.rank may 24% be correct for the tolerance stats.tol but the estimated error 25% bounds are too large to confirm this. However stats.rank appears 26% to be correct for an alternate tolerance tol_alt. More 27% generally stats.rank appears to be correct for any tolerance 28% between stats.est_sval_lower_bounds(nsvals_large) and 29% stats.est_sval_upper_bounds(nsvals_large+1). 30% stats.flag is 2 if the calculated numerical rank stats.numerical 31% may be correct but estimated error bounds are too large to confirm 32% this. The conditions for stats.flag to be 0 or 1 are not 33% satisfied. 34% stats.flag is 3 if is likely that the numerical rank returned, 35% stats.rank, is too large. 36% 37% stats.tol_alt is an alternate tolerance that corresponds to the 38% calculated numerical rank when stats.flag is 1. 39% 40% other fields of statistics 41% 42% Input: 43% call_from == 44% call_from = 1 indicates a call from spqr_basic 45% call_from = 2 indicates a call from spqr_cod 46% est_sval_upper_bounds -- est_sval_upper_bounds(i) is an 47% estimate of an upper bound on singular value number 48% stats.sval_numbers_for_bounds(i) of A. 49% est_sval_lower_bounds -- stats.est_sval_lower_bounds(i) is an 50% estimate of an lower bound on singular value number 51% sval_numbers_for_bounds(i) of A. 52% tol - the tolerance defining the numerical rank. The true 53% numerical rank is the number of singular values larger than tol. 54% numerical_rank -- the estimated numerical rank of A 55% nsvals_small -- the number of estimated singular values, from 56% the set that have been estimated, that appear to be smaller 57% than or equal to tol. 58% nsvals_large -- the number of estimated singular values, from 59% the set that have been estimated, that appear to be larger 60% than tol. 61% stats -- the statistics returned by spqr_basic, spqr_null, spqr_pinv 62% and spqr_cod. 63% stats_ssi -- the statistics returned by spqr_ssi. 64% opts -- options for the calls to spqr_function. 65% nargout_call -- the value of nargout in the calling function 66% stats_ssp_N -- stats from spqr_ssp applied to A * N 67% stats_ssp_NT -- stats from spqr_ssp applied to A' * NT 68% start_tic -- tic time for start of calling routine 69 70% Copyright 2012, Leslie Foster and Timothy A Davis. 71 72get_details = opts.get_details ; 73 74%------------------------------------------------------------------------------- 75% determine flag which indicates accuracy of the estimated numerical rank 76%------------------------------------------------------------------------------- 77 78if ( ( nsvals_small == 0 || ... 79 est_sval_upper_bounds(nsvals_large+1) <= tol ) ) && ... 80 ( nsvals_large > 0 && est_sval_lower_bounds(nsvals_large) > tol ) 81 % numerical rank is correct, assuming estimated error bounds are correct 82 flag = 0; 83elseif ( nsvals_small > 0 && nsvals_large > 0 ) && ... 84 ( est_sval_lower_bounds(nsvals_large) > ... 85 est_sval_upper_bounds(nsvals_large + 1) && ... 86 est_sval_upper_bounds(nsvals_large + 1) > tol ) 87 % in this case, assuming that the estimated error bounds are correct, 88 % then the numerical rank is correct with a modified tolerance 89 flag = 1; 90 stats.tol_alt = est_sval_upper_bounds(nsvals_large + 1); 91 % Note: satisfactory values of tol_alt are in the range 92 % est_sval_lower_bounds(nsvals_large) > tol_alt 93 % >= est_sval_upper_bounds(nsvals_large + 1) 94elseif stats_ssi.flag == 3 95 % in this case ssi failed and it is often the case that then 96 % calculated numerical rank is too high 97 flag = 3; 98else 99 % in this case, assuming that the estimated error bounds are correct, 100 % the errors in the error bounds are too large to determine the 101 % numerical rank 102 flag = 2; 103end 104 105 106%------------------------------------------------------------------------------- 107% return statistics 108%------------------------------------------------------------------------------- 109 110stats.flag = flag ; 111stats.rank = numerical_rank ; 112stats.est_sval_upper_bounds = est_sval_upper_bounds; 113stats.est_sval_lower_bounds = est_sval_lower_bounds; 114stats.sval_numbers_for_bounds = ... 115 numerical_rank - nsvals_large + 1 : ... 116 numerical_rank + nsvals_small ; 117 118if call_from == 2 119 if ( get_details == 1 || nargout_call >= 3 ) 120 stats.est_norm_A_times_N = stats_ssp_N.est_svals(1); 121 end 122 if (get_details == 1) 123 stats.est_err_bound_norm_A_times_N = stats_ssp_N.est_error_bounds(1); 124 end 125end 126 127 128if ( ( call_from == 1 && nargout_call == 3 ) || ... 129 ( call_from == 2 && nargout_call == 4 ) ) 130 % include estimated norm of A transpose time NT from call to ssp 131 stats.est_norm_A_transpose_times_NT = stats_ssp_NT.est_svals(1); 132 if get_details == 1 133 stats.est_err_bound_norm_A_transpose_times_NT = ... 134 stats_ssp_NT.est_error_bounds(1); 135 end 136end 137 138if stats.tol_alt == -1 139 stats = rmfield(stats, 'tol_alt') ; 140end 141 142if get_details == 1 143 % order the fields of stats in a convenient order (the fields when 144 % get_details is 0 or 2 are already in a good order) 145 stats = spqr_rank_order_fields(stats); 146end 147 148if (get_details == 1) 149 stats.time = toc (start_tic) ; 150end 151 152