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