1function err = test_function (A, strategy, burble)
2%TEST_FUNCTION test various functions applied to a factorize object
3%
4% Example
5%   test_functions (A) ;    % where A is square or rectangular, sparse or dense
6%
7% See also test_all, factorize, inverse, mldivide
8
9% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com
10
11reset_rand ;
12
13state = warning ('off', 'MATLAB:illConditionedMatrix') ;
14
15if (nargin < 1)
16    A = rand (10) ;
17end
18if (nargin < 2)
19    strategy = 'default' ;
20end
21if (nargin < 3)
22    burble = 0 ;
23end
24
25err = 0 ;
26[m, n] = size (A) ;
27F = factorize (A, strategy) ;
28
29%   C = rand (m,n) ;
30%   if (~isempty (A))
31%       C (1,1) = A (1,1) ;
32%   end
33%   if (issparse (A))
34%       C = sparse (A) ;
35%   end
36%   H = factorize (C, strategy) ; %#ok
37
38%-------------------------------------------------------------------------------
39% implicitly-defined methods that return the correct result:
40%-------------------------------------------------------------------------------
41
42if (m == 10)
43    burble = 1 ;
44    % only do these pedantic tests on a few test matrices
45    assert (~ismethod (F, 'ismethod')) ;
46    assert (~ismethod (F, 'iscell')) ;      assert (~iscell (F)) ;
47    assert (~ismethod (F, 'iscellstr')) ;   assert (~iscellstr (F)) ;
48    assert (~ismethod (F, 'ischar')) ;      assert (~ischar (F)) ;
49    assert (~ismethod (F, 'iscom')) ;       assert (~iscom (F)) ;
50    assert (~ismethod (F, 'isinteger')) ;   assert (~isinteger (F)) ;
51    assert (~ismethod (F, 'ishandle')) ;    assert (~ishandle (F)) ;
52    assert (~ismethod (F, 'isobject')) ;    assert ( isobject (F)) ;
53    assert (~ismethod (F, 'isinterface')) ; assert (~isinterface (F)) ;
54    assert (~ismethod (F, 'isjava')) ;      assert (~isjava (F)) ;
55    assert (~ismethod (F, 'iskeyword')) ;   assert (~iskeyword (F)) ;
56    assert (~ismethod (F, 'isletter')) ;    assert (~isletter (F)) ;
57    assert (~ismethod (F, 'isspace')) ;     assert (~isspace (F)) ;
58    assert (~ismethod (F, 'islogical')) ;   assert (~islogical (F)) ;
59    assert (~ismethod (F, 'isstruct')) ;    assert (~isstruct (F)) ;
60    assert (~ismethod (F, 'isvarname')) ;   assert (~isvarname (F)) ;
61    assert (~ismethod (F, 'isglobal')) ;    % not tested since it's deprecated
62    assert (~ismethod (F, 'ndims')) ;       assert (ndims (F) == 2) ;
63    assert (~ismethod (F, 'isequal')) ;
64    assert (~ismethod (F, 'isequalwithequalnans')) ;
65    % One could imagine that these return length(F.A) and numel(F.A), but if
66    % they are defined that way, F.A fails because the size of F is wrong.
67    assert (~ismethod (F, 'length')) ;      assert (length (F) == 1) ;
68    assert (~ismethod (F, 'numel')) ;       assert (numel (F) == 1) ;
69end
70
71G = factorize (A, strategy) ;
72any_nans = any (any (isnan (F.A))) || any (any (isnan (double (inverse (F))))) ;
73if (~any_nans)
74    assert ( isequal (F, G)) ;
75    assert (~isequal (F, inverse (F))) ;
76end
77assert ( isequalwithequalnans (F, G)) ;
78assert (~isequalwithequalnans (F, inverse (F))) ;
79clear G
80
81%-------------------------------------------------------------------------------
82% explicit methods
83%-------------------------------------------------------------------------------
84
85if (min (m,n) < 2 || m == 10)
86
87    % explicitly-defined methods, tested here:
88    assert (ismethod (F, 'isfloat')) ;  assert ( isfloat (F)) ;
89    assert (ismethod (F, 'isnumeric')); assert ( isnumeric (F)) ;
90    assert (ismethod (F, 'issingle')) ; assert (~issingle (F)) ;
91    assert (ismethod (F, 'isreal')) ;   assert (isreal (F) == isreal (A)) ;
92    assert (ismethod (F, 'isempty')) ;  assert (isempty (F) == isempty (A)) ;
93    assert (ismethod (F, 'isvector')) ; assert (isvector (F) == isvector (A)) ;
94    assert (ismethod (F, 'isscalar')) ; assert (isscalar (F) == isscalar (A)) ;
95    assert (ismethod (F, 'issparse')) ; assert (issparse (F) == issparse (A)) ;
96    assert (ismethod (F, 'size')) ;     assert (isequal (size (A), size (F))) ;
97
98    if (~any_nans)
99        assert (ismethod (F, 'abs')) ;    assert (isequal (abs (F), abs (F.A)));
100        assert (ismethod (F, 'double')) ; assert (isequal (A, double (F))) ;
101    end
102
103    assert (ismethod (F, 'isfield')) ;
104    assert ( isfield (F, 'A')) ;
105    assert ( isfield (F, 'Factors')) ;
106    assert ( isfield (F, 'is_inverse')) ;
107    assert ( isfield (F, 'is_ctrans')) ;
108    assert ( isfield (F, 'alpha')) ;
109    assert ( isfield (F, 'kind')) ;
110    assert (~isfield (F, 'anthing_else')) ;
111
112    assert (ismethod (F, 'isa')) ;
113    assert ( isa (F, 'double')) ;
114    assert (~isa (F, 'logical')) ;
115    assert (~isa (F, 'char')) ;
116    assert (~isa (F, 'single')) ;
117    assert ( isa (F, 'float')) ;
118    assert (~isa (F, 'int8')) ;
119    assert (~isa (F, 'uint8')) ;
120    assert (~isa (F, 'int16')) ;
121    assert (~isa (F, 'uint16')) ;
122    assert (~isa (F, 'int32')) ;
123    assert (~isa (F, 'uint32')) ;
124    assert (~isa (F, 'int64')) ;
125    assert (~isa (F, 'uint64')) ;
126    assert (~isa (F, 'integer')) ;
127    assert ( isa (F, 'numeric')) ;
128    assert (~isa (F, 'cell')) ;
129    assert (~isa (F, 'struct')) ;
130    assert (~isa (F, 'function_handle')) ;
131    assert ( isa (F, 'factorization')) ;
132    assert (~isa (F, 'anything_else')) ;
133
134    % explicitly-defined methods, but tested elsewhere:
135    assert (ismethod (F, 'mldivide')) ;
136    assert (ismethod (F, 'mldivide_subclass')) ;
137    assert (ismethod (F, 'mrdivide')) ;
138    assert (ismethod (F, 'mrdivide_subclass')) ;
139    assert (ismethod (F, 'inverse')) ;
140    assert (ismethod (F, 'ctranspose')) ;
141    assert (ismethod (F, 'end')) ;
142    assert (ismethod (F, 'mtimes')) ;
143    assert (ismethod (F, 'subsref')) ;
144    assert (ismethod (F, 'disp')) ;
145    assert (ismethod (F, 'condest')) ;
146    assert (ismethod (F, 'struct')) ;
147    if (isa (F, 'factorization_chol_dense'))
148        assert (ismethod (F, 'cholupdate')) ;
149        if (burble)
150            methods (F)
151        end
152    end
153    if (isa (F, 'factorization_svd'))
154        assert (ismethod (F, 'cond')) ;
155        assert (ismethod (F, 'norm')) ;
156        assert (ismethod (F, 'rank')) ;
157        assert (ismethod (F, 'null')) ;
158        assert (ismethod (F, 'orth')) ;
159        assert (ismethod (F, 'pinv')) ;
160        assert (ismethod (F, 'svd')) ;
161        if (burble)
162            methods (F)
163        end
164    else
165        assert (~ismethod (F, 'cond')) ;
166        assert (~ismethod (F, 'norm')) ;
167        assert (~ismethod (F, 'rank')) ;
168        assert (~ismethod (F, 'null')) ;
169        assert (~ismethod (F, 'orth')) ;
170        assert (~ismethod (F, 'pinv')) ;
171        assert (~ismethod (F, 'svd')) ;
172    end
173end
174
175%-------------------------------------------------------------------------------
176% 1-norm and condition estimates
177%-------------------------------------------------------------------------------
178
179if (m == n)
180
181    % norm(A,1)
182    reset_rand ;
183    if (isempty (A))
184        nest1 = 0 ;
185    else
186        nest1 = full (normest1 (A)) ;
187    end
188    if (isa (F, 'factorization_svd'))
189        nexact = norm (F, 1) ;
190    else
191        nexact = norm (A, 1) ;
192    end
193    if (burble)
194        fprintf ('\nnorm(A,1), exact:             %g\n', nexact) ;
195        fprintf ('  MATLAB normest1(A)          %g\n', nest1) ;
196    end
197
198    % norm(inv(A),1)
199    if (isempty (A))
200        imine = 0 ;
201    else
202        imine = full (normest1 (inverse (A))) ;
203    end
204    if (isa (F, 'factorization_svd'))
205        iexact = norm (inverse (F), 1) ;
206    else
207        iexact = norm (inv (A), 1) ;
208    end
209    try
210        reset_rand ;
211        iest = full (normest1 (inv (A))) ;
212    catch %#ok
213        iest = -1 ;
214    end
215    if (burble)
216        fprintf ('\nnorm (inv(A),1) exact:        %g\n', iexact) ;
217        fprintf ('  MATLAB normest1 (inv (A)):  %g\n', iest) ;
218        fprintf ('  normest1 (inverse (F)):     %g\n', imine) ;
219    end
220
221    reset_rand ;
222    kest = full (condest (A)) ;
223    reset_rand ;
224    kF = full (condest (F)) ;
225    kS = full (condest (inverse (F))) ;
226    err = max (err, abs (rankest (F) - F.A_rank)) ;
227
228    kexact = -1 ;
229    kexact2 = -1 ;
230    kexact3 = -1 ;
231    if (isa (F, 'factorization_svd'))
232        try
233            kexact = cond (F, 1) ;
234            kexact1 = F.A_cond ;
235            kexact2 = cond (F) ;
236            err = max (err, abs (kexact1-kexact2)) ;
237            kexact3 = cond (full (A)) ;
238            err = max (err, abs (kexact2-kexact3) / max (1, abs (kexact3))) ;
239        catch %#ok
240        end
241    end
242
243    if (burble)
244        fprintf ('\n  cond (A,1), exact:          %g\n', kexact) ;
245        fprintf ('  MATLAB condest(A):          %g\n', kest) ;
246        fprintf ('  condest(F):                 %g\n', kF) ;
247        fprintf ('  condest(inverse(A)):        %g\n', kS) ;
248        fprintf ('  cond (A,2), exact:          %g\n', kexact2) ;
249        fprintf ('  cond (F,2), exact:          %g\n', kexact3) ;
250        fprintf ('  rankest %d %d\n', rankest (F), F.A_rank) ;
251        if (~isempty (F.A_condest))
252            fprintf ('  cheap condest:              %g\n', F.A_condest) ;
253        end
254    end
255
256    if (err > 1e-9)
257        display (A) ;
258        display (full (A)) ;
259        display (strategy) ;
260        display (err) ;
261        error ('error too high!') ;
262    end
263end
264
265%-------------------------------------------------------------------------------
266% compute explicit entries of the inverse
267
268if (~any_nans && ~isempty (A))
269    P1 = pinv (full (A)) ;
270    P2 = inverse (A) ;
271    s1 = P1 (1) ;
272    s2 = P2 (1) ;
273    err = max (err, abs (s1-s2) / max (1, abs (s1))) ;
274    s1 = P1 (:,1) ;
275    s2 = P2 (:,1) ;
276    err = max (err, norm (s1-s2,1) / max (1, norm (s1,1))) ;
277    s1 = P1 (1,:) ;
278    s2 = P2 (1,:) ;
279    err = max (err, norm (s1-s2,1) / max (1, norm (s1,1))) ;
280end
281
282%-------------------------------------------------------------------------------
283% test struct
284%-------------------------------------------------------------------------------
285
286K = struct (F) ;
287if (burble)
288    disp ('K =') ;
289    disp (K) ;
290end
291
292if (~any_nans)
293    assert (isequal (K.A, F.A)) ;
294    assert (isequal (K.Factors, F.Factors)) ;
295end
296assert (isequal (K.kind, F.kind)) ;
297assert (isequal (K.is_inverse, F.is_inverse)) ;
298assert (isequal (K.is_ctrans, F.is_ctrans)) ;
299assert (isequal (K.alpha, F.alpha)) ;
300
301K = struct (inverse (F)') ;
302if (burble)
303    disp ('K =') ;
304    disp (K) ;
305end
306
307if (~any_nans)
308    assert (isequal (K.A, F.A)) ;
309    assert (isequal (K.Factors, F.Factors)) ;
310end
311assert (isequal (K.is_inverse, ~F.is_inverse)) ;
312assert (isequal (K.is_ctrans, ~F.is_ctrans)) ;
313assert (isequal (K.alpha, F.alpha)) ;
314assert (isequal (K.kind, F.kind)) ;
315
316% restore user's warnings
317warning (state) ;
318