sparseinv_test

SPARSEINV_TEST tests the sparseinv function.

 Example
   sparseinv_test ;        % basic test
   sparseinv_test (1) ;    % extensive test (requires UFget)

 See also sparseinv, sparseinv_install, UFget.

0001 function sparseinv_test (extensive)
0002 %SPARSEINV_TEST tests the sparseinv function.
0003 %
0004 % Example
0005 %   sparseinv_test ;        % basic test
0006 %   sparseinv_test (1) ;    % extensive test (requires UFget)
0007 %
0008 % See also sparseinv, sparseinv_install, UFget.
0009 
0010 % Copyright 2011, Timothy A. Davis, http://www.suitesparse.com
0011 
0012 if (nargin < 1)
0013     extensive = 0 ;
0014 end
0015 
0016 load west0479 ;
0017 A = west0479 ;
0018 
0019 for k = 1:2
0020 
0021     [Z, Zpattern, L, D, U, P, Q, stats] = sparseinv (A) ;
0022 
0023     n = size (A,1) ;
0024     I = speye (n) ;
0025 
0026     S = inv (A) ;
0027     Snorm = norm (S, 1) ;
0028     e1 = norm (Zpattern.*(Z-S), 1) / Snorm ;
0029     e2 = norm ((L+I)*D*(U+I) - P*A*Q, 1) / norm (A,1) ;
0030     c = condest (A) ;
0031 
0032     fprintf ('west0479: errors %g %g condest %g : ', e1, e2, c) ;
0033     if (e1 > 1e-10 || e2 > 1e-10)
0034         error ('west0479 test failed') ;
0035     end
0036     fprintf ('ok\n') ;
0037     disp (stats) ;
0038 
0039     % create a symmetric positive definite matrix to test with
0040     if (k == 1)
0041         A = A+A' + 1e6*I ;
0042     end
0043 end
0044 
0045 % check error-handling
0046 fprintf ('testing error handling (errors below are expected)\n') ;
0047 ok = 1 ;
0048 A = ones (2) ;
0049 try
0050     Z1 = sparseinv (A) ;                                                    %#ok
0051     ok = 0 ;
0052 catch me
0053     fprintf ('    expected error: %s\n', me.message) ;
0054 end
0055 A = sparse (ones (3,2)) ;
0056 try
0057     Z2 = sparseinv (A) ;                                                    %#ok
0058     ok = 0 ;
0059 catch me
0060     fprintf ('    expected error: %s\n', me.message) ;
0061 end
0062 A = sparse (ones (3)) ;
0063 try
0064     Z3 = sparseinv (A) ;                                                    %#ok
0065     ok = 0 ;
0066 catch me
0067     fprintf ('    expected error: %s\n', me.message) ;
0068 end
0069 A = 1i * sparse (ones (3)) ;
0070 try
0071     Z4 = sparseinv (A) ;                                                    %#ok
0072     ok = 0 ;
0073 catch me
0074     fprintf ('    expected error: %s\n', me.message) ;
0075 end
0076 
0077 if (~ok)
0078     error ('error-handling tests failed') ;
0079 end
0080 
0081 % now try with lots of matrices from the UF Sparse Matrix Collection
0082 if (extensive && exist ('UFget', 'file') == 2)
0083 
0084     fprintf ('Now doing extensive tests with UF Sparse Matrix Collection:\n') ;
0085     dofigures = (exist ('cspy', 'file') == 2) ;
0086     if (dofigures)
0087         clf ;
0088     end
0089 
0090     index = UFget ;
0091     f = find ((index.nrows == index.ncols) & (index.isReal == 1)) ;
0092     [ignore,i] = sort (index.nrows (f)) ;   %#ok
0093     f = f (i) ;
0094     nmat = length (f) ;
0095 
0096     s = warning ('off', 'MATLAB:nearlySingularMatrix') ;
0097 
0098     for k = 1:nmat
0099         id = f (k) ;
0100         Prob = UFget (id, index) ;
0101         A = Prob.A ;
0102         n = size (A,1) ;
0103         I = speye (n) ;
0104         fprintf ('id: %4d  n: %4d : %-30s', id, n, Prob.name) ;
0105 
0106         Z = [ ] ; 
0107         try
0108             [Z, Zpattern, L, D, U, P, Q] = sparseinv (A) ;
0109         catch me
0110             fprintf ('%s', me.message) ;
0111         end
0112 
0113         if (~isempty (Z))
0114             e = norm ((L+I) * D * (U+I) - P*A*Q, 1) / norm (A,1) ;
0115             fprintf ('errs:  %12.3e ', e) ;
0116             if (e > 1e-10)
0117                 error ('error in factorization too high') ;
0118             end
0119             S = inv (A) ;           % normally S has MANY nonzero entries
0120             Snorm = norm (S,1) ;
0121             E = abs (Zpattern .* (Z-S)) / Snorm ;
0122             e = norm (E, 1) ;
0123             c = condest (A) ;
0124             fprintf (' %12.3e  condest: %12.2e', e, c) ;
0125             if (e/c  > 1e-8)
0126                 error ('error in sparseinv too high') ;
0127             end
0128             fprintf (' ok') ;
0129 
0130             if (dofigures)
0131                 subplot (2,2,1) ; cspy (A) ;
0132                 title (Prob.name, 'Interpreter', 'none') ;
0133                 subplot (2,2,2) ; cspy (P*A*Q) ;  title ('P*A*Q') ;
0134                 subplot (2,2,3) ; cspy (Z) ;      title ('sparse inverse') ;
0135                 subplot (2,2,4) ; cspy (S) ;      title ('inverse') ;
0136                 drawnow
0137             end
0138 
0139         end
0140         fprintf ('\n') ;
0141         if (n >= 300)
0142             break ;
0143         end
0144     end
0145 
0146     warning (s) ;       % restore warning status
0147 end
0148 
0149 fprintf ('All sparseinv tests passed.\n') ;