Chiao/dev (#3)

* changed to enable -NDEBUG in mmex for -O

* Updated license

* Added FGMRES solver for nullspace

* Set rank for backsolve

* Added PIPIT

* updated reference to hifir

* Removed unused codes

* Updated CPP source directory

* Added version and updated ijv2crs

* Updated C++ kernels in HIFIR4M with new HIFIR

* Updated examples using new APIs.

* Updated using field name for params

* Fixed authors

* Fixed checking before recompilation

* updated .gitignore

* updated pipitHifir

* Added download hifir script

* Updated building script to use downloaded hifir

* Removed submodule

* Added option to use local C++ hifir

* Updated building process

* Updated ignore list for user hifir source

* Changed to urlwrite

* Updated varaible name

* Fixed integer size for Octave

* Added update and refactor

* Fixed bugs in apply

* Fixed bugs in create

* Updated building script

* Updated cite bibtex in readme

* Fixed bug in hifir

* Added version file and build

* Added user HIF

* Updated control parameters

* Updated README

* Updated to addparam

* Updated installation

* Added user-provided HIF

* Added changelog

Co-authored-by: Xiangmin Jiao <[email protected]>

.gitignore

@@ -7,4 +7,7 @@
 .DS_Store
 .*swp
 *.o
+*.bak
+*.asv
 .MATLABDriveTag
+hifir-*

CHANGELOG.md

@@ -0,0 +1,7 @@
+# Changelog #
+
+All notable changes to this project will be documented in this file.
+
+## [v0.1.0](https://github.com/hifirworks/hifir4m/releases/tag/v0.1.0) (2021-09-02) ##
+
+This is the initial official release of the HIFIR4M package based on the C++ HIFIR [v0.1.0](https://github.com/hifirworks/hifir/releases/tag/v0.1.0).

README.md

@@ -4,20 +4,36 @@
 
 ## Installation ##
 
-Clone this project to your preferred location. Then start MATLAB or GNU Octave under the directory that contains `hifir4m`, or run the command
+Clone this project to your preferred location, i.e.,
+
+```console
+git clone -b release https://github.com/hifirworks/hifir4m.git
+```
+
+Use `git pull` to download any new changes that have been added since `git clone` or last `git pull`. Alternatively, use `git checkout v[GLOBAL].[MAJOR].[MINOR]` to download a specific version.
+
+Then start MATLAB or GNU Octave under the directory that contains `hifir4m`, or run the command
 
 ```matlab
->> run('/path/to/hifir4m/load_hifir')
+>> run('/path/to/hifir4m/startup_hifir')
 ```
 
 (and replace `/path/to/hifir4m/` to the directory that contains `hifir4m`) bo build `hifir4m` and load its path. It will build the *mex* kernels if needed by linking to MATLAB and Octave's built-in BLAS/LAPACK libraries for the low-level QRCP.
 
+Note that for the first time, `hifir4m` will download C++ package [HIFIR](https://github.com/hifirworks/hifir) while building *mex* kernels. If you don't have access to network during building *mex* kernels, then you can obtain HIFIR beforehand and put its source in `hifir4m/hifir-[hifir-version]` folder; for instance, you can run the following command to download the C++ HIFIR package
+
+```console
+cd /path/to/hifir4m
+wget -qO- https://github.com/hifirworks/hifir/archive/refs/tags/v`cat VERSION`.tar.gz|tar xzf -
+```
+
 ## Usage ##
 
 The easiest way to use `hifir4m` is to call the `gmresHif` interface. For example,
 ```matlab
 >> [x, flag, relres, iter, reshis, times] = gmresHif(A, b);
 ```
+
 where `A` is a MATLAB's built-in sparse matrix or a `MATLAB` `struct` containing the filds of `row_ptr`, `col_ind`, `vals` of a standard CRS storage format, and `b` is a right-hand-side vector (or RHS vectors with two columns).
 
 To access the intermediate-level interfaces of `hifir4m`, please see `gmresHif.m` for the calling sequence of `hifCreate`, `hifApply`, and `hifDestroy`.
@@ -29,6 +45,7 @@ To access the intermediate-level interfaces of `hifir4m`, please see `gmresHif.m
 The software suite is released under a dual-license model. For academic users, individual users, or open-source software developers, you can use HIFIR under the GPLv3+ license free of charge for research and evaluation purpose. For commercial users, separate commercial licenses are available through the Stony Brook University.  For inqueries regarding commercial licenses, please contact Prof. Xiangmin Jiao <[email protected]>.
 
 ## How to Cite `HIFIR` ##
+
 If you use `HIFIR`,  `hifir4m`, or `hifir4py` in your research for nonsingular systems, please cite the `HILUCSI` paper:
 
 ```bibtex
@@ -38,7 +55,6 @@ If you use `HIFIR`,  `hifir4m`, or `hifir4py` in your research for nonsingular s
              large-scale saddle-point problems from {PDE}s},
   journal = {Numer. Linear Algebra Appl.},
   year    = {2021},
-  note    = {To appear},
   doi     = {10.1002/nla.2400},
 ```
 
@@ -47,20 +63,20 @@ If you use them to solve highly ill-conditioned of singular systems, please cite
 ```bibtex
 @Article{jiao2020approximate,
   author  = {Xiangmin Jiao and Qiao Chen},
-  journal = {arxiv},
+  journal = {SIAM J. Matrix Anal. Appl.},
   title   = {Approximate Generalized Inverses with Iterative Refinement for 
              $\epsilon$-Accurate Preconditioning of Singular Systems},
-  year    = {2020},
-  note    = {arXiv:2009.01673},
+  year    = {2021},
+  note    = {To appear},
 }
 
 @Article{chen2021hifir,
-  author  = {Jiao, Xiangmin and Chen, Qiao},
+  author  = {Chen, Qiao and Jiao, Xiangmin},
   title   = {{HIFIR}: Hybrid Incomplete Factorization with Iterative Refinement 
              for Preconditioning Ill-conditioned and Singular Systems},
   journal = {arxiv},
   year    = {2021},
-  note    = {arXiv:21...},
+  note    = {arXiv:2106.09877},
 }
 ```
 

VERSION

@@ -0,0 +1 @@
+0.1.0

api/HifEnum.m

@@ -8,11 +8,8 @@
         DESTROY = int32(4);
         FACTORIZE = int32(5);
         M_SOLVE = int32(6);
-        KSP_SOLVE = int32(7);
-        KSP_NULL = int32(8);
-        EXPORT_DATA = int32(9);
-        M_SOLVE2 = int32(10);
-        M_MULTIPLY = int32(11);
-        QUERY = int32(12);
+        M_MULTIPLY = int32(7);
+        EXPORT_DATA = int32(8);
+        QUERY = int32(9);
     end
 end

api/HifParams.m

@@ -21,7 +21,7 @@
 %           verbose: 0
 %            rf_par: 1
 %           reorder: 1
-%            saddle: 1
+%               spd: 0
 %             check: 1
 %         pre_scale: 0
 %     symm_pre_lvls: 1
@@ -54,33 +54,35 @@
         'verbose', 0, ...
         'rf_par', 1, ...
         'reorder', 1, ...
-        'saddle', 1, ...
+        'spd', 0, ...
         'check', 1, ...
         'pre_scale', 0, ...
         'symm_pre_lvls', 1, ...
         'threads', 0, ...
         'mumps_blr', 1, ...
         'fat_schur_1st', 0, ...
         'rrqr_cond', 0, ...
-        'pivot', 0, ...
+        'pivot', 2, ...
         'gamma', 1, ...
         'beta', 1e3, ...
         'is_symm', 0, ...
         'no_pre', 0, ...
+        'nzp_thres', 0.65, ...
+        'dense_thres', 2000, ...
         'is_mixed', 0, ...
         'is_complex', 0);
     fields = fieldnames(defaults);
 end
 
 p = inputParser;
 positive_params = {'tau_L', 'tau_U', 'kappa_d', 'kappa', 'alpha_L', ...
-    'alpha_U', 'rho', 'c_d', 'c_h'};
+    'alpha_U', 'rho', 'c_d', 'c_h', 'nzp_thres', 'dense_thres'};
 for i = 1:numel(positive_params)
     addParameter(p, positive_params{i}, defaults.(positive_params{i}), ...
         @(x) isscalar(x) && x > 0);
 end
 
-nonnegative_params = {'verbose', 'rf_par', 'reorder', 'saddle', 'check', ...
+nonnegative_params = {'verbose', 'rf_par', 'reorder', 'spd', 'check', ...
     'pre_scale', 'symm_pre_lvls', 'threads', 'mumps_blr', 'fat_schur_1st', ...
     'pivot', 'gamma', 'beta', 'is_symm', 'no_pre', 'is_mixed', 'is_complex'};
 for i = 1:numel(nonnegative_params)
@@ -91,6 +93,7 @@
 addParameter(p, 'N', defaults.N, @isscalar);
 addParameter(p, 'rrqr_cond', defaults.rrqr_cond, @isscalar);
 
+p.KeepUnmatched = true;  % This seems needed for Octave
 parse(p, varargin{:});
 sorted_opts = p.Results;
 

api/Hifir.m

@@ -34,6 +34,16 @@
             [varargout{1:nargout}] = hifApply(h, x, varargin{:});
         end
 
+        function update(h, A)
+            % Update coefficient matrix used in iterative refinement
+            h = hifUpdate(h, A);
+        end
+
+        function varargout = refactorize(h, S, varargin)
+            % Refactorization
+            [varargout{1:nargout}] = hifRefactorize(h, S, varargin{:});
+        end
+
         function delete(h)
             % Destructor for handle
             if ~isempty(h.hdl)

api/hifApply.m

@@ -1,4 +1,4 @@
-function [y, varargout] = hifApply(hif, x, op, rank, nirs)
+function varargout = hifApply(hif, x, op, rank, nirs)
 % hifApply  Applies a HIFIR preconditioner to solve or multiply.
 %
 %     y = hifApply(hif, x [op, rank, nirs])
@@ -25,19 +25,20 @@
     assert(nargin <= 5, 'Solve accepts up to five arguments.');
     if nargin == 5 && nirs > 1
         % with iterative refinement
-        [y, varargout{:}] = hifir4m_mex(HifEnum.M_SOLVE, hif.hdl, ...
-            A.row_ptr, A.col_ind, A.val, x, nirs, ...
+        assert(~isempty(hif.A), 'hif.A must not be empty.');
+        [varargout{1:nargout}] = hifir4m_mex(HifEnum.M_SOLVE, hif.hdl, ...
+            x, hif.A.row_ptr, hif.A.col_ind, hif.A.val, nirs, ...
             op(end) == 'H' || op(end) == 'T' || op(end) == 'h' || op(end) == 't', rank);
     else
-        [y, varargout{1:nargout}] = hifir4m_mex(HifEnum.M_SOLVE, hif.hdl, x, ...
+        [varargout{1:nargout}] = hifir4m_mex(HifEnum.M_SOLVE, hif.hdl, x, ...
             op(end) == 'H' || op(end) == 'T' || op(end) == 'h' || op(end) == 't', rank);
     end
 else
     assert(op(1) == 'M' || op(1) == 'm',  ...
         'First character must be ''s'' or ''m''');
     assert(nargin <= 4, 'Multiply accepts up to four arguments.');
 
-    [y, varargout{1:nargout}] = hifir4m_mex(HifEnum.M_MULTIPLY, hif.hdl, x, ...
+    [varargout{1:nargout}] = hifir4m_mex(HifEnum.M_MULTIPLY, hif.hdl, x, ...
         op(end) == 'H' || op(end) == 'T' || op(end) == 'h' || op(end) == 't', rank);
 end
 

api/hifCreate.m

@@ -37,36 +37,37 @@
 if issparse(A)
     Astruct = hifir4m_sp2crs(A);
 else
+    Astruct = A;
     if hifir4m_isint64
         if ~isa(A.row_ptr, 'int64')
             Astruct.row_ptr = int64(A.row_ptr);
             Astruct.col_ind = int64(A.col_ind);
         end
-    elseif ~isint64
+    else
         if ~isa(A.row_ptr, 'int32')
             Astruct.row_ptr = int32(A.row_ptr);
             Astruct.col_ind = int32(A.col_ind);
         end
     end
-    Astruct.val = double(A.val);
 end
-Astruct.nrows = int32(numel(Astruct.row_ptr)-1);
 
 if nargin <= 1 || isempty(S)
     Sstruct = Astruct;
+elseif issparse(S)
+    Sstruct = hifir4m_sp2crs(S);
 else
-    if hifir4m_isint64
-        if ~isa(S.row_ptr, 'int64')
-            Sstruct.row_ptr = int64(S.row_ptr);
-            Sstruct.col_ind = int64(S.col_ind);
-        end
-    elseif ~isint64
-        if ~isa(S.row_ptr, 'int32')
-            Sstruct.row_ptr = int32(S.row_ptr);
-            Sstruct.col_ind = int32(S.col_ind);
-        end
+    Sstruct = S;
+end
+if hifir4m_isint64
+    if ~isa(Sstruct.row_ptr, 'int64')
+        Sstruct.row_ptr = int64(Sstruct.row_ptr);
+        Sstruct.col_ind = int64(Sstruct.col_ind);
+    end
+else
+    if ~isa(Sstruct.row_ptr, 'int32')
+        Sstruct.row_ptr = int32(Sstruct.row_ptr);
+        Sstruct.col_ind = int32(Sstruct.col_ind);
     end
-    Sstruct.val = double(S.val);
 end
 
 [varargout{1:nargout-1}] = hifir4m_mex(HifEnum.FACTORIZE, hif, ...

api/hifRefactorize.m

@@ -0,0 +1,38 @@
+function varargout = hifRefactorize(hif, S, varargin)
+% hifRefactorize refactorizes a preconditioner given a sparsifier S
+%
+%   hifRefactorize(hif, S)
+%   hifRefactorize(hif, S, 'alpha_L', 5, 'alpha_U', 5);
+%
+% See also hifUpdate hifCreate hifApply
+
+assert(~isempty(hif.hdl), 'the preconditioner cannot be empty.');
+
+if ~isempty(varargin) && isstruct(varargin{1})
+    params = varargin{1};
+else
+    params = HifParams(varargin{:});
+end
+
+% Setup sparsifier
+if issparse(S)
+    Sstruct = hifir4m_sp2crs(S);
+else
+    Sstruct = S;
+    if hifir4m_isint64
+        if ~isa(S.row_ptr, 'int64')
+            Sstruct.row_ptr = int64(S.row_ptr);
+            Sstruct.col_ind = int64(S.col_ind);
+        end
+    else
+        if ~isa(S.row_ptr, 'int32')
+            Sstruct.row_ptr = int32(S.row_ptr);
+            Sstruct.col_ind = int32(S.col_ind);
+        end
+    end
+end
+
+[varargout{1:nargout-1}] = hifir4m_mex(HifEnum.FACTORIZE, hif.hdl, ...
+    Sstruct.row_ptr, Sstruct.col_ind, Sstruct.val, params);
+
+end

api/hifUpdate.m

@@ -0,0 +1,29 @@
+function hif = hifUpdate(hif, A)
+% hifUpdate Updates a coefficient matrix in hif for iterative refinement
+%
+%   hif = hifUpdate(hif, A)
+%
+% See also hifApply hifRefactorize
+
+if issparse(A)
+    Astruct = hifir4m_sp2crs(A);
+else
+    Astruct = A;
+    if hifir4m_isint64
+        if ~isa(A.row_ptr, 'int64')
+            Astruct.row_ptr = int64(Astruct.row_ptr);
+            Astruct.col_ind = int64(Astruct.col_ind);
+        end
+    else
+        if ~isa(A.row_ptr, 'int32')
+            Astruct.row_ptr = int32(Astruct.row_ptr);
+            Astruct.col_ind = int32(Astruct.col_ind);
+        end
+    end
+end
+if ~isempty(hif.A)
+    assert(hif.A.nrows == Astruct.nrows, 'mismatched row sizes');
+    assert(hif.A.ncols == Astruct.ncols, 'mismatched column sizes');
+end
+hif.A = Astruct;
+end

examples/complex_eg.m

@@ -1,23 +1,23 @@
-clear;
-load('kim1.mat', 'Problem');
-is_mixed = false;
-is_complex = true;
-A = Problem.A;
-n = size(A,1);
-b = A*ones(n,1);
+%{
+    This example code uses HIF-preconditioned GMRES(30) to solve a
+    complex testing system "kim1" from SuiteSparse Matrix Collection.
+%}
 
-%% Initialize database
-h = HIF(is_mixed, is_complex);
+clear;
 
-%% Factorize A
-info = h.factorize(A);
-disp(info);
+load('kim1.mat', 'Problem');
 
-%% Solve for x=A\b
-[x, flag] = h.fgmres(A, b);
-if flag; fprintf(2, 'warning! solver failed with flag %d\n', flag); end
-res = norm(A*x-b)/norm(b);
-if res > 1e-6; fprintf(2, 'residual too large %.4g\n', res); end
+% RHS is A*1
+A = Problem.A;
+n = size(A, 1);
+b = A*ones(n, 1);
 
-%% Finalize
-clear h
+% solve with HIF-preconditioned GMRES
+[x, flag] = gmresHif(A, b, 'is_complex', true);
+if flag
+    fprintf(2, 'warning! solver failed with flag %d\n', flag);
+end
+res = norm(b - A * x) / norm(b);
+if res > 1e-6
+    fprintf(2, 'residual too large %.4g\n', res);
+end