MATLAB demos and test

SPEX/ExampleMats/494_bus.mat.txt

@@ -1664,4 +1664,4 @@
 493 493 111.638400
 304 494 -66.225170
 488 494 -44.722720
-494 494 110.947900
+494 494 110.947900

SPEX/ExampleMats/README.txt

@@ -1,49 +1,49 @@
 Example matrices for SPEX demos and testing.
 
-All these matrices 1-based.
+All matrices 1-based.
 
 Files:
 
-    README.txt
+README.txt
 
-    matrix                  right-hand-side         description
-    ------                  ---------------         -----------
-    10teams.mat.txt         10teams.rhs.txt         Basislib matrix
-    tomography.mat          tomography.mat.rhs      Earth science matrix from 
-                                                    MathWorks 
-    LF10.mat.txt            LF10.rhs.txt            SPD linear 1D beam matrix 
-                                                    from Oberwolfach model 
-                                                    reduction benchmark collection 
-                                                    (collected from the SuiteSparse 
-                                                    Matrix Collection)
-    LFAT5.mat.txt           LFAT5.rhs.txt           SPD linear 1D beam matrix 
-                                                    from Oberwolfach model 
-                                                    reduction benchmark collection 
-                                                    (collected from the SuiteSparse 
-                                                    Matrix Collection)
-    Trefethen_500.mat.txt   Trefethen_500.rhs.txt   Diagonal matrix with primes 
-                                                    (collected from the SuiteSparse 
-                                                    Matrix Collection)
-    494_bus.mat.txt         494_bus.rhs.txt         SPD matrix from the Harwell
-                                                    -Boeing collection (collected 
-                                                    from the SuiteSparse Matrix 
-                                                    Collection)
-    mesh1e1.mat.txt         mesh1e1.rhs.txt         SPD matrix from NASA 
-                                                    (collected from the SuiteSparse 
-                                                    Matrix Collection)
-    example.mat.txt         example.rhs.txt         Integer matrix
-    NSR8K.mat.txt           NSR8K.rhs.txt           Matrix from a real-world 
-                                                    mixed integer program 
-                                                    (collected from MIPLIB)
-    test1.mat.txt           test1.rhs.txt           Unsymmetric double matrix 
-                                                    created to test SPEX Cholesky
-    test2.mat.txt           test2.rhs.txt           Unsymmetric integer matrix 
-                                                    created to test SPEX Cholesky
-    test3.mat.txt           test3.rhs.txt           Unsymmetric matrix created 
-                                                    to test SPEX Cholesky
-    test4.mat.txt           test4.rhs.txt           Symmetric indefinite matrix 
-                                                    created to test SPEX Cholesky
-    test5.mat.txt           test5.rhs.txt           Small rectangular matrix 
-                                                    created to test SPEX
-    test.mat.txt            test.rhs.txt            Integer matrix to test SPEX LU
+matrix                  right-hand-side         description
+------                  ---------------         -----------
+10teams.mat.txt         10teams.rhs.txt         Basislib matrix
+tomography.mat.txt      tomography.mat.rhs.txt  Earth science matrix from
+                                                MathWorks
+LF10.mat.txt            LF10.rhs.txt            SPD linear 1D beam matrix
+                                                from Oberwolfach model
+                                                reduction benchmark collection
+                                                (collected from the SuiteSparse
+                                                Matrix Collection)
+LFAT5.mat.txt           LFAT5.rhs.txt           SPD linear 1D beam matrix
+                                                from Oberwolfach model
+                                                reduction benchmark collection
+                                                (collected from the SuiteSparse
+                                                Matrix Collection)
+Trefethen_500.mat.txt   Trefethen_500.rhs.txt   Diagonal matrix with primes
+                                                (collected from the SuiteSparse
+                                                Matrix Collection)
+494_bus.mat.txt         494_bus.rhs.txt         SPD matrix from the Harwell
+                                                -Boeing collection (collected
+                                                from the SuiteSparse Matrix
+                                                Collection)
+mesh1e1.mat.txt         mesh1e1.rhs.txt         SPD matrix from NASA
+                                                (collected from the SuiteSparse
+                                                Matrix Collection)
+example.mat.txt         example.rhs.txt         Integer matrix
+NSR8K.mat.txt           NSR8K.rhs.txt           Matrix from a real-world
+                                                mixed integer program
+                                                (collected from MIPLIB)
+test1.mat.txt           test1.rhs.txt           Unsymmetric double matrix
+                                                created to test SPEX Cholesky
+test2.mat.txt           test2.rhs.txt           Unsymmetric integer matrix
+                                                created to test SPEX Cholesky
+test3.mat.txt           test3.rhs.txt           Unsymmetric matrix created
+                                                to test SPEX Cholesky
+test4.mat.txt           test4.rhs.txt           Symmetric indefinite matrix
+                                                created to test SPEX Cholesky
+test5.mat.txt           test5.rhs.txt           Small rectangular matrix
+                                                created to test SPEX
+test.mat.txt            test.rhs.txt            Integer matrix to test SPEX LU
 

SPEX/ExampleMats/mesh1e1.mat.txt

@@ -304,4 +304,4 @@
 2 48 1.280289
 42 48 1.130447
 45 48 1.351190
-48 48 5.968440
+48 48 5.968440

SPEX/ExampleMats/tomography.mat → SPEX/ExampleMats/tomography.mat.txt

@@ -28724,4 +28724,4 @@
 477 500 444.602886
 498 500 5295.602821
 499 500 32053.297605
-500 500 12741928.914287
+500 500 12741928.914287

SPEX/MATLAB/spex_backslash.m

@@ -123,7 +123,7 @@
 end
 
 % SPEX Backslash expects sparse input.
-% So, if A is not sprase it is sprasified.
+% So, if A is not sparse it is sprasified.
 if (~issparse(A))
     A = sparse(A);
 end

SPEX/MATLAB/spex_lu_backslash.m

@@ -1,4 +1,4 @@
-function x = SPEX_lu_backslash (A,b,option)
+function x = spex_lu_backslash (A,b,option)
 % SPEX_LU_BACKSLASH: solve Ax=b via sparse left-looking integer-preserving LU
 % spex_lu_backslash computes the exact solution to the sparse linear system
 % Ax = b where A and b are stored as doubles. A must be stored as a sparse

SPEX/MATLAB/spex_mex_demo.m

@@ -1,24 +1,27 @@
-%% spex_mex_demo a demo of SPEX MATLAB interface
-% SPEX is a package for solving sparse linear systems of equations
-% with a roundoff-free integer-preserving method.  The result is
-% always exact, unless the matrix A is perfectly singular.
+%% spex_mex_demo a demo of the SPEX MATLAB interface
+% SPEX is a package for solving sparse linear systems of
+% equations with a roundoff-free integer-preserving method.
+% The result is always exact, unless the matrix A is perfectly
+% singular.
+%
+% See also vpa, spex_backslash, spex_lu_backslash,
+%   spex_cholesky_backslash, spex_ldl_backslash,
+%   spex_mex_install, spex_mex_test.
 %
-% See also vpa, spex_backslash, spex_lu_backslash, spex_cholesky_backslash,
-%   spex_ldl_backslash, spex_mex_install, spex_mex_test.
-
 % SPEX: (c) 2022-2024, Christopher Lourenco, Jinhao Chen,
-% Lorena Mejia Domenzain, Erick Moreno-Centeno, and Timothy A. Davis.
-% All Rights Reserved.
+% Lorena Mejia Domenzain, Erick Moreno-Centeno, and Timothy A.
+% Davis.  All Rights Reserved.
 % SPDX-License-Identifier: GPL-2.0-or-later or LGPL-3.0-or-later
 
-
-format compact
+%#ok<*NOPTS>
+%#ok<*NASGU>
 
 %% SPEX vs MATLAB backslash: first example
-% In this first example, x = spex_backslash (A, b) ; returns an approximate
-% solution, but not because it was computed incorrectly in spex_backslash.
-% It is computed exactly as a rational result in spex_lu_backslash with
-% arbitrary precision, but then converted to double precision on output.
+% In this first example, x = spex_backslash (A, b) returns an
+% approximate solution, but not because it was computed
+% incorrectly in spex_backslash.  It is computed exactly as a
+% rational result in spex_backslash with arbitrary precision,
+% but then converted to double precision on output.
 
 format long g
 load west0479
@@ -27,18 +30,21 @@
 xtrue = rand (n,1) ;
 b = A*xtrue ;
 x = spex_backslash (A, b) ;
-% error is nonzero: x is computed exactly in rational arbitrary-precision,
-% but then lost precision when returned to MATLAB:
+% error is nonzero: x is computed exactly in rational
+% arbitrary-precision, but then loses precision when
+% returned as a double vector to MATLAB:
 err_spex = norm (x-xtrue)
+
 x = A\b ;
-% error is nonzero: MATLAB x=A\b experiences floating-point error
-% throughout its computations:
+% error is nonzero: MATLAB x=A\b experiences floating-point
+% error throughout its computations:
 err_matlab = norm (x-xtrue)
 
 %% SPEX backslash: exact, vs MATLAB backslash: approximate
-% In this example, x = spex_backslash (A, b) ; is returned exactly in the
-% MATLAB vector x, because x contains only integers representable exactly
-% in double precision.  x = A\b results in floating-point roundoff error.
+% In this example, x = spex_backslash (A, b) is returned
+% exactly in the MATLAB vector x, because x contains only
+% integers representable exactly in double precision.  x = A\b
+% results in floating-point roundoff error.
 
 amax = max (abs (A), [ ], 'all') ;
 A = floor (2^20 * (A / amax)) + n * speye (n) ;
@@ -52,109 +58,101 @@
 err_matlab = norm (x-xtrue)
 
 %% SPEX Backslash on singular problems
-% The matrix 2008 is a rank deficient rectangular matrix
-% If we compute A = A'*A we obtain this 9*9 integer matrix:
-%
-%     4    -1    -1     0    -1     0     0    -1     0
-%    -1     4     0    -1     0    -1     0    -1     0
-%    -1     0     4    -1    -1     0    -1     0     0
-%     0    -1    -1     4     0    -1    -1     0     0
-%    -1     0    -1     0     4    -1     0     0    -1
-%     0    -1     0    -1    -1     4     0     0    -1
-%     0     0    -1    -1     0     0     4    -1    -1
-%    -1    -1     0     0     0     0    -1     4    -1
-%     0     0     0     0    -1    -1    -1    -1     4
-%
-% Since this is an integer matrix, SPEX can obtain it exactly.
-% This matrix is ill-conditioned and singular, SPEX correctly determines
-% that it is. However, MATLAB Cholesky factorization both
-% succeeds and returns a useless answer.
-% Note that to continue the demo, the backslash line is wrapped in a try catch
-% the default behavior of SPEX is to return an error to MATLAB and exit
-
-fprintf("Testing singular matrix");
-prob = ssget(2008); A = prob.A;
-B = full(A'*A);
-b = ones(9,1);
-
-try x = spex_backslash (B, b) ;
-catch fprintf("\nError, spex determines matrix is singular\n");
+% SPEX correctly determines that the following integer matrix is
+% singular.  The MATLAB numerical Cholesky factorization fails to
+% do so, returning a useless solution.
+
+A = sparse ([ % matrix JGD_Homology/ch3-3-b1:
+      4    -1    -1     0    -1     0     0    -1     0
+     -1     4     0    -1     0    -1     0    -1     0
+     -1     0     4    -1    -1     0    -1     0     0
+      0    -1    -1     4     0    -1    -1     0     0
+     -1     0    -1     0     4    -1     0     0    -1
+      0    -1     0    -1    -1     4     0     0    -1
+      0     0    -1    -1     0     0     4    -1    -1
+     -1    -1     0     0     0     0    -1     4    -1
+      0     0     0     0    -1    -1    -1    -1     4 ]) ;
+b = ones (9,1);
+try
+    x = spex_backslash (A, b) ;
+catch me
+    fprintf ("\nSPEX Error: %s\n", me.message) ;
 end
 
-R = chol(B);
-x_chol = R \ (R' \ b)
-
-R = vpa(chol(B));
-x_vpa = vpa ( R \ vpa( (R'\b)))
+R = chol(A);
+x_chol_matlab = R \ (R' \ b)
 
 %% SPEX on an ill-conditioned problem
-% x = spex_backslash (A,b) is able to accurately solve problems that x=A\b cannot.
-% Consider the following matrix in the MATLAB gallery:
+% x = spex_backslash (A,b) is able to accurately solve problems
+% that x=A\b cannot.  Consider the Wilkinson gallery matrix:
 
-[U, b] = gallery ('wilk', 3)
+[A, b] = gallery ('wilk', 3)
 
-%%    vpa can find a good but not perfect solution:
-xvpa = vpa (U) \ b
+%% vpa can find a good but not perfect solution:
+xvpa = vpa (A) \ b
 
-%     but MATLAB's numerical x = U\b computes a poor solution:
-xapprox = U \ b
+% but MATLAB's numerical x = A\b computes a poor solution:
+xapprox = A \ b
 
 %% spex_backslash computes the exact answer
-% It returns it to MATLAB as a double vector, obtaining the exact results
+% SPEX returns the exact solution a double vector:
 
-xspex = spex_backslash (U, b)
+xspex = spex_backslash (A, b)
 err = xvpa - xspex
 relerr = double (err (2:3) ./ xvpa (2:3))
 
-%% spex_lu_backslash with exact results
-% spex_lu_backslash can also return x as a cell array of strings, which
-% preserves the exact rational result.  The printing option is also
-% enabled in this example.  The floating-point matrices U and b are
-% converted into a scaled integer matrix before solving U*x=b with
-% SPEX Left LU.
+%% spex_backslash with exact results
+% spex_backslash can also return x as a cell array of strings,
+% which preserves the exact rational result.  The printing option
+% is also enabled in this example.  The floating-point matrices A
+% and b are converted into a scaled integer matrix before solving
+% A*x=b with SPEX Left LU.
 %
-% Note that, importantly, SPEX obtains an integer matrix by scaling the
-% input. SPEX scales all input by 1e16. This is because consider the number
-% U(1,2). The value U(1,2)=0.9 is a floating point number and cannot be
-% represented exactly in IEEE floating-point. Specifically, the rational
-% represenation of it is fl(0.9) = 45000000000000001 / 50000000000000000.
+% Note that, importantly, SPEX obtains an integer matrix by
+% scaling the input. SPEX scales all input by 1e16. This is
+% because consider the number A(1,2). The value A(1,2)=0.9 is a
+% floating point number and cannot be represented exactly in IEEE
+% floating-point. Specifically, the rational represenation of it
+% is fl(0.9) = 45000000000000001 / 50000000000000000.
 %
-% SPEX assumes the user wants what they typed in. Scaling this matrix exactly
-% gives the above rational. Conversely scaling with 16 digits of precision
-% gives the exact fraction 9/10.
+% SPEX assumes the user wants what they typed in. Scaling this
+% matrix exactly gives the above rational. Conversely scaling
+% with 16 digits of precision gives the exact fraction 9/10.
 %
-% If one wishes to obtain FULL floating-point precision and/or support
-% for floating point values smaller than 1e-16 there are two options:
+% If one wishes to obtain FULL floating-point precision and/or
+% support for floating point values smaller than 1e-16 there are
+% two options:
 %
-%   1) Within MATLAB the user scales the matrix themselves. If SPEX is
-%      given an integer matrix it is preserved exactly
+%   1) Within MATLAB the user scales the matrix themselves. If
+%      SPEX is given an integer matrix it is preserved exactly.
 %
 %   2) Within C convert the matrix to a SPEX_MPFR. MPFR numbers
-%      can exactly store doubles. The conversion from MPFR to integer
-%      is fully exact and all one would obtain the rational representation
-%      of the floating-point number itself. This can be done with 1 call
-%      to the SPEX matrix copy function.
+%      can exactly store doubles. The conversion from MPFR to
+%      integer is fully exact and all one would obtain the
+%      rational representation of the floating-point number
+%      itself. This can be done with 1 call to the SPEX matrix
+%      copy function.
 %
-% In any case, below is the exact solution with U(1,2) = 9/10
+% In any case, below is the exact solution with A(1,2) = 9/10
 
 option.print = 3 ;          % also print the details
 option.solution = 'char' ;  % return x as a cell array of strings
 
 %%
 
-xspex = spex_lu_backslash (U, b, option)
+xspex = spex_backslash (A, b, option)
 
 %% Converting an exact rational result to vpa or double
-% If spex_lu_backslash returns x as a cell array of strings, it cannot
-% be immediately used in computations in MATLAB.  It can be converted
-% into a vpa or double matrix, as illustrated below.
+% If spex_backslash returns x as a cell array of strings, it
+% cannot be immediately used in computations in MATLAB.  It can
+% be converted into a vpa or double matrix, as illustrated below.
 
 xspex_as_vpa = vpa (xspex)
 xspex_as_double = double (vpa (xspex))
 xvpa_as_double = double (xvpa)
 
 %% Comparing the VPA and SPEX_BACKSLASH solutions in double
-% Both vpa(U)\b and spex_backslash(U,b) compute the same result
+% Both vpa(A)\b and spex_backslash(A,b) compute the same result
 % in the end, when their results are converted to double.
 err = xvpa_as_double - xspex_as_double