1function err = test_factorize (A, strategy)
2%TEST_FACTORIZE test the accuracy of the factorization object
3%
4% Example
5%   test_factorize (A) ;    % where A is square or rectangular, sparse or dense
6%   test_factorize (A, strategy) ;  % forces a particular strategy;
7%                           % works only if the matrix is compatible.
8%
9% See also test_all, factorize, inverse, mldivide
10
11% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com
12
13reset_rand ;
14if (nargin < 1)
15    A = rand (100) ;
16end
17
18if (nargin < 2)
19    strategy = '' ;
20end
21
22% do not check the sparsity of the result when using the SVD
23spcheck = ~(strcmp (strategy, 'svd')) ;
24
25err = 0 ;
26if (strcmp (strategy, 'ldl') && issparse (A) && ~isreal(A))
27    % do not test ldl on sparse complex matrices
28    return ;
29end
30
31[m, n] = size (A) ;
32if (min (m,n) > 0)
33    anorm = norm (A,1) ;
34else
35    anorm = 1 ;
36end
37
38is_symmetric = ((m == n) && (nnz (A-A') == 0)) ;
39
40for nrhs = 1:3
41
42    for bsparse = 0:1
43
44        b = rand (m,nrhs) ;
45        if (bsparse)
46            b = sparse (b) ;
47        end
48
49        %-----------------------------------------------------------------------
50        % test backslash and related methods
51        %-----------------------------------------------------------------------
52
53        for a = [1 pi (pi+1i)]
54
55            % method 0:
56            x = (a*A)\b ;
57            err = check_resid (err, anorm, a*A, x, b, spcheck) ;
58
59            % method 1:
60            S = inverse (A)/a ;
61            x = S*b ;
62            err = check_resid (err, anorm, a*A, x, b, spcheck) ;
63
64            % method 3:
65            F = testfac (A, strategy) ;
66            S = inverse (F)/a ;
67            x = S*b ;
68            err = check_resid (err, anorm, a*A, x, b, spcheck) ;
69
70            % method 4:
71            F = a*testfac (A, strategy) ;
72            x = F\b ;
73            err = check_resid (err, anorm, a*A, x, b, spcheck) ;
74
75            % method 5:
76            S = inverse (F) ;
77            x = S*b ;
78            err = check_resid (err, anorm, a*A, x, b, spcheck) ;
79
80            % method 6:
81            if (m == n)
82                [L, U, p] = lu (A, 'vector') ;
83                x = a * (U \ (L \ (b (p,:)))) ;
84                err = check_resid (err, anorm, A/a, x, b, spcheck) ;
85
86            % method 7:
87                if (is_symmetric)
88                    F = a*factorize (A, 'symmetric') ;
89                    x = F\b ;
90                    err = check_resid (err, anorm, a*A, x, b, spcheck) ;
91                else
92                    F = a*factorize (A, 'unsymmetric') ;
93                    x = F\b ;
94                    err = check_resid (err, anorm, a*A, x, b, spcheck) ;
95                end
96            end
97
98        end
99
100        clear S F
101
102        %------------------------------------------------------------------
103        % test transposed backslash and related methods
104        %------------------------------------------------------------------
105
106        b = rand (n,nrhs) ;
107        if (bsparse)
108            b = sparse (b) ;
109        end
110
111        for a = [1 pi (pi+1i)]
112
113            % method 0:
114            x = (a*A)'\b ;
115            err = check_resid (err, anorm, (a*A)', x, b, spcheck) ;
116
117            % method 1:
118            S = inverse (A) / a ;
119            x = S'*b ;
120            err = check_resid (err, anorm, (a*A)', x, b, spcheck) ;
121
122            % method 2:
123            F = a*testfac (A, strategy) ;
124            x = F'\b ;
125            err = check_resid (err, anorm, (a*A)', x, b, spcheck) ;
126
127            % method 3:
128            S = inverse (F') ;
129            x = S*b ;
130            err = check_resid (err, anorm, (a*A)', x, b, spcheck) ;
131        end
132
133        clear S F
134
135        %------------------------------------------------------------------
136        % test mtimes, times, plus, minus, rdivide, and ldivide
137        %------------------------------------------------------------------
138
139        for a = [1 pi pi+2i]
140
141            F = a*testfac (A, strategy) ;
142            S = inverse (F) ;
143            B = a*A ;   % F is the factorization of a*A
144
145            % test mtimes
146            d = rand (n,1)  ;
147            x = F*d ;
148            y = B*d ;
149            z = S\d ;
150            err = check_error (err, norm (mtx(x)-mtx(y),1)) ;
151            err = check_error (err, norm (mtx(x)-mtx(z),1)) ;
152
153            % test mtimes transpose
154            d = rand (m,1)  ;
155            x = F'*d ;
156            y = B'*d ;
157            z = S'\d ;
158            err = check_error (err, norm (mtx(x)-mtx(y),1)) ;
159            err = check_error (err, norm (mtx(x)-mtx(z),1)) ;
160
161            % test for scalars
162            for s = [1 42 3-2i]
163                E = s\B  - double (s\F) ; err = check_error (err, norm (E,1)) ;
164                E = B/s  - double (F/s) ; err = check_error (err, norm (E,1)) ;
165%               E = B.*s - F.*s ;       err = check_error (err, norm (E,1)) ;
166%               E = s.*B - s.*F ;       err = check_error (err, norm (E,1)) ;
167%               E = s.\B - s.\F ;       err = check_error (err, norm (E,1)) ;
168%               E = B./s - F./s ;       err = check_error (err, norm (E,1)) ;
169                E = B*s  - double (F*s) ; err = check_error (err, norm (E,1)) ;
170                E = s*B  - double (s*F) ; err = check_error (err, norm (E,1)) ;
171            end
172
173        end
174
175        clear S F C
176
177        %------------------------------------------------------------------
178        % test slash and related methods
179        %------------------------------------------------------------------
180
181        b = rand (nrhs,n) ;
182        if (bsparse)
183            b = sparse (b) ;
184        end
185
186        % method 0:
187        x = b/A ;
188        err = check_resid (err, anorm, A', x', b', spcheck) ;
189
190        % method 1:
191        S = inverse (A) ;
192        x = b*S ;
193        err = check_resid (err, anorm, A', x', b', spcheck) ;
194
195        % method 4:
196        F = testfac (A, strategy) ;
197        x = b/F ;
198        err = check_resid (err, anorm, A', x', b', spcheck) ;
199
200        % method 5:
201        S = inverse (F) ;
202        x = b*S ;
203        err = check_resid (err, anorm, A', x', b', spcheck) ;
204
205        % method 6:
206        if (m == n)
207            [L, U, p] = lu (A, 'vector') ;
208            x = (b / U) / L ; x (:,p) = x ;
209            err = check_resid (err, anorm, A', x', b', spcheck) ;
210        end
211
212        %------------------------------------------------------------------
213        % test transpose slash and related methods
214        %------------------------------------------------------------------
215
216        b = rand (nrhs,m) ;
217        if (bsparse)
218            b = sparse (b) ;
219        end
220
221        % method 0:
222        x = b/A' ;
223        err = check_resid (err, anorm, A, x', b', spcheck) ;
224
225        % method 1:
226        S = inverse (A)' ;
227        x = b*S ;
228        err = check_resid (err, anorm, A, x', b', spcheck) ;
229
230        % method 4:
231        F = testfac (A, strategy)' ;
232        x = b/F ;
233        err = check_resid (err, anorm, A, x', b', spcheck) ;
234
235        % method 5:
236        S = inverse (F) ;
237        x = b*S ;
238        err = check_resid (err, anorm, A, x', b', spcheck) ;
239
240        %------------------------------------------------------------------
241        % test double
242        %------------------------------------------------------------------
243
244        Y = double (inverse (A)) ;
245        if (m == n)
246            Z = inv (A) ;
247        else
248            Z = pinv (full (A)) ;
249        end
250        e = norm (Y-Z,1) ;
251        if (n > 0)
252            e = e / norm (Z,1) ;
253        end
254        err = check_error (e, err) ;
255
256        %------------------------------------------------------------------
257        % test subsref
258        %------------------------------------------------------------------
259
260        F = testfac (A, strategy) ;
261        Y = inverse (A) ;
262        if (numel (A) > 1)
263            if (F (end) ~= A (end))
264                error ('factorization subsref error') ;
265            end
266            if (F.A (end) ~= A (end))
267                error ('factorization subsref error') ;
268            end
269        end
270        if (n > 0)
271            if (F (1,1) ~= A (1,1))
272                error ('factorization subsref error') ;
273            end
274            if (F.A (1,1) ~= A (1,1))
275                error ('factorization subsref error') ;
276            end
277            e = abs (Y (1,1) - Z (1,1)) ;
278            err = check_error (e,err) ;
279            if (m > 1 && n > 1)
280                e = norm (Y (1:2,1:2) - Z (1:2,1:2), 1) ;
281                err = check_error (e,err) ;
282            end
283        end
284        if (m > 3 && n > 1)
285            if (any (F (2:end, 1:2) - A (2:end, 1:2)))
286                error ('factorization subsref error') ;
287            end
288            if (any (F (2:4, :) - A (2:4, :)))
289                error ('factorization subsref error') ;
290            end
291            if (any (F (:, 1:2) - A (:, 1:2)))
292                error ('factorization subsref error') ;
293            end
294        end
295
296        %------------------------------------------------------------------
297        % test transposed subsref
298        %------------------------------------------------------------------
299
300        FT = F' ;
301        YT = Y' ;
302        AT = A' ;
303        ZT = Z' ;
304        if (numel (AT) > 1)
305            if (FT (end) ~= AT (end))
306                error ('factorization subsref error') ;
307            end
308            if (F.A (end) ~= A (end))
309                error ('factorization subsref error') ;
310            end
311        end
312
313        if (n > 0)
314            if (FT (1,1) ~= AT (1,1))
315                error ('factorization subsref error') ;
316            end
317            if (F.A (1,1) ~= A (1,1))
318                error ('factorization subsref error') ;
319            end
320            e = abs (YT (1,1) - ZT (1,1)) ;
321            err = check_error (e,err) ;
322            if (m > 1 && n > 1)
323                e = norm (YT (1:2,1:2) - ZT (1:2,1:2), 1) ;
324                err = check_error (e,err) ;
325            end
326        end
327
328        if (m > 1 && n > 3)
329            if (any (FT (2:end, 1:2) - AT (2:end, 1:2)))
330                error ('factorization subsref error') ;
331            end
332            if (any (FT (2:4, :) - AT (2:4, :)))
333                error ('factorization subsref error') ;
334            end
335            if (any (FT (:, 1:2) - AT (:, 1:2)))
336                error ('factorization subsref error') ;
337            end
338        end
339
340        %------------------------------------------------------------------
341        % test update/downdate
342        %------------------------------------------------------------------
343
344        if (isa (F, 'factorization_chol_dense'))
345            w = rand (n,1) ;
346            b = rand (n,1) ;
347            % update
348            G = cholupdate (F,w) ;
349            x = G\b ;
350            err = check_resid (err, anorm, A+w*w', x, b, spcheck) ;
351            % downdate
352            G = cholupdate (G,w,'-') ;
353            x = G\b ;
354            err = check_resid (err, anorm, A, x, b, spcheck) ;
355            clear G
356        end
357
358        %------------------------------------------------------------------
359        % test size
360        %------------------------------------------------------------------
361
362        [m1, n1] = size (F) ;
363        [m, n] = size (A) ;
364        if (m1 ~= m || n1 ~= n)
365            error ('size error') ;
366        end
367        [m1, n1] = size (Y) ;
368        if (m1 ~= n || n1 ~= m)
369            error ('pinv size error') ;
370        end
371        if (size (Y,1) ~= n || size (Y,2) ~= m)
372            error ('pinv size error') ;
373        end
374        if (size (F,1) ~= m || size (F,2) ~= n)
375            error ('size error') ;
376        end
377
378        %------------------------------------------------------------------
379        % test mtimes
380        %------------------------------------------------------------------
381
382        clear S F Y
383
384        for a = [1 pi pi+2i]
385
386            F = a * testfac (A, strategy) ;
387            S = inverse (F) ;
388            d = rand (1,m) ;
389            x = d*F ;
390            y = d*(A*a) ;
391            z = d/S ;
392
393            err = check_error (err, norm (mtx(x)-mtx(y),1)) ;
394            err = check_error (err, norm (mtx(x)-mtx(z),1)) ;
395
396            d = rand (1,n) ;
397            x = d*F' ;
398            y = d*(A*a)' ;
399            z = d/S' ;
400
401            err = check_error (err, norm (mtx(x)-mtx(y),1)) ;
402            err = check_error (err, norm (mtx(x)-mtx(z),1)) ;
403
404        end
405
406        %------------------------------------------------------------------
407        % test inverse
408        %------------------------------------------------------------------
409
410        Y = double (inverse (inverse (A))) ;
411        e = norm (A-Y,1) ;
412        if (e > 0)
413            error ('inverse error') ;
414        end
415
416        %------------------------------------------------------------------
417        % test mldivide and mrdivide with a matrix b, transpose, etc
418        %------------------------------------------------------------------
419
420        if (max (m,n) < 100)
421            F = testfac (A, strategy) ;
422            B = rand (m,n) ;
423            err = check_error (err, norm (B\A - mtx(B\F), 1) / anorm) ;
424            err = check_error (err, norm (A/B - mtx(F/B), 1) / anorm) ;
425            err = check_error (err, norm (B*pinv(full(A))-mtx(B/F), 1) / anorm);
426        end
427    end
428end
429
430fprintf ('.') ;
431
432%--------------------------------------------------------------------------
433
434function err = check_resid (err, anorm, A, x, b, spcheck)
435[m, n] = size (A) ;
436
437x = mtx (x) ;
438
439if (m <= n)
440    e = norm (A*x-b,1) / (anorm + norm (x,1)) ;
441else
442    e = norm (A'*(A*x)-A'*b,1) / (anorm + norm (x,1)) ;
443end
444
445if (min (m,n) > 1 && spcheck)
446    if (issparse (A) && issparse (b))
447        if (~issparse (x))
448            error ('x must be sparse') ;
449        end
450    else
451        if (issparse (x))
452            error ('x must be full') ;
453        end
454    end
455end
456
457err = check_error (e, err) ;
458
459%--------------------------------------------------------------------------
460
461function x = mtx (x)
462% make sure that x is a matrix.  It might be a factorization.
463if (isobject (x))
464    x = double (x) ;
465end
466
467%--------------------------------------------------------------------------
468
469function err = check_error (err1, err2)
470err = max (err1, err2) ;
471if (err > 1e-8)
472    fprintf ('error: %8.3e\n', full (err)) ;
473    error ('error too high') ;
474end
475err = full (err) ;
476
477%--------------------------------------------------------------------------
478
479function F = testfac (A, strategy)
480if (isempty (strategy))
481    F = factorize (A) ;
482else
483    F = factorize (A, strategy) ;
484end
485