Mercurial > octave-nkf
changeset 19089:38937efbee21
Incorporate new functions ichol and ilu into Octave.
* NEWS: Announce new functions.
* aspell-octave.en.pws: Add new functions names to custom Octave dictionary.
* sparse.txi: Add functions to Octave manual.
* __unimplemented__.m: Remove functions from unimplemented list.
* lu.cc (Flu), luinc.cc (Fluinc), chol.cc (Fchol): Add seealso links in
docstrings.
* __ichol__.cc: Wrap long lines to less than 80 chars. Remove trailing
whitespace. Don't repeat input validation done in ichol.m for internal
functions. Avoid resizing retval vector.
* __ilu__.cc: Wrap long lines to less than 80 chars. Remove trailing
whitespace. Don't repeat input validation done in ichol.m for internal
functions. Avoid resizing retval vector.
* ichol.m: Rewrite docstring. Use Octave coding conventions (double quotes
hash for comments, ! instead of ~). Replace %!error tests not being run
with fail().
* ilu.m: Rewrite docstring. Use Octave coding conventions (double quotes
hash for comments, ! instead of ~). Replace %!error tests not being run
with fail().
| author | Rik <rik@octave.org> |
|---|---|
| date | Tue, 26 Aug 2014 15:27:21 -0700 |
| parents | df64071e538c |
| children | 4630a18757b3 |
| files | NEWS doc/interpreter/doccheck/aspell-octave.en.pws doc/interpreter/sparse.txi libinterp/corefcn/lu.cc libinterp/corefcn/luinc.cc libinterp/dldfcn/__ichol__.cc libinterp/dldfcn/__ilu__.cc libinterp/dldfcn/chol.cc scripts/help/__unimplemented__.m scripts/sparse/ichol.m scripts/sparse/ilu.m |
| diffstat | 11 files changed, 980 insertions(+), 1109 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Mon Aug 18 12:32:16 2014 +0100 +++ b/NEWS Tue Aug 26 15:27:21 2014 -0700 @@ -66,10 +66,11 @@ ** Other new functions added in 4.2: - bandwidth isbanded javachk - dir_in_loadpath isdiag linkaxes - hgload istril numfields - hgsave istriu + bandwidth ilu javachk + dir_in_loadpath isbanded linkaxes + hgload isdiag numfields + hgsave istril + ichol istriu ** Deprecated functions.
--- a/doc/interpreter/doccheck/aspell-octave.en.pws Mon Aug 18 12:32:16 2014 +0100 +++ b/doc/interpreter/doccheck/aspell-octave.en.pws Tue Aug 26 15:27:21 2014 -0700 @@ -416,6 +416,7 @@ hygernd Hypergeometric hypergeometric +ichol IEC ieee IEEE @@ -424,6 +425,7 @@ ifftn ignorecase ij +ilu Im imag ImageMagick
--- a/doc/interpreter/sparse.txi Mon Aug 18 12:32:16 2014 +0100 +++ b/doc/interpreter/sparse.txi Tue Aug 26 15:27:21 2014 -0700 @@ -485,11 +485,11 @@ @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm} @item Linear algebra: - @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank}, - @dfn{spaugment}, @dfn{svds} + @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, + @dfn{normest}, @dfn{sprank}, @dfn{spaugment}, @dfn{svds} @item Iterative techniques: - @dfn{luinc}, @dfn{pcg}, @dfn{pcr} + @dfn{ichol}, @dfn{ilu}, @dfn{luinc}, @dfn{pcg}, @dfn{pcr} @c @dfn{bicg}, @dfn{bicgstab}, @dfn{cholinc}, @dfn{cgs}, @dfn{gmres}, @c @dfn{lsqr}, @dfn{minres}, @dfn{qmr}, @dfn{symmlq} @@ -860,7 +860,7 @@ The left division @code{\} and right division @code{/} operators, discussed in the previous section, use direct solvers to resolve a linear equation of the form @code{@var{x} = @var{A} \ @var{b}} or -@code{@var{x} = @var{b} / @var{A}}. Octave equally includes a number of +@code{@var{x} = @var{b} / @var{A}}. Octave also includes a number of functions to solve sparse linear equations using iterative techniques. @DOCSTRING(pcg) @@ -873,6 +873,10 @@ @var{A} \ @var{b}} is solved instead. Typical pre-conditioning matrices are partial factorizations of the original matrix. +@DOCSTRING(ichol) + +@DOCSTRING(ilu) + @DOCSTRING(luinc) @node Real Life Example
--- a/libinterp/corefcn/lu.cc Mon Aug 18 12:32:16 2014 +0100 +++ b/libinterp/corefcn/lu.cc Tue Aug 26 15:27:21 2014 -0700 @@ -137,7 +137,7 @@ is embedded into @var{U} to give a return value similar to the full case.\n\ For both full and sparse matrices, @code{lu} loses the permutation\n\ information.\n\ -@seealso{luupdate, chol, hess, qr, qz, schur, svd}\n\ +@seealso{luupdate, ilu, chol, hess, qr, qz, schur, svd}\n\ @end deftypefn") { octave_value_list retval;
--- a/libinterp/corefcn/luinc.cc Mon Aug 18 12:32:16 2014 +0100 +++ b/libinterp/corefcn/luinc.cc Tue Aug 26 15:27:21 2014 -0700 @@ -93,7 +93,7 @@ \n\ Given the string argument @qcode{\"vector\"}, @code{luinc} returns the\n\ values of @var{p} @var{q} as vector values.\n\ -@seealso{sparse, lu}\n\ +@seealso{sparse, lu, ilu, ichol}\n\ @end deftypefn") { int nargin = args.length ();
--- a/libinterp/dldfcn/__ichol__.cc Mon Aug 18 12:32:16 2014 +0100 +++ b/libinterp/dldfcn/__ichol__.cc Tue Aug 26 15:27:21 2014 -0700 @@ -1,22 +1,25 @@ -/** - * Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> - * Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> - * - * This file is part of Octave. - * - * Octave is free software; you can redistribute it and/or modify it under the - * terms of the GNU General Public License as published by the Free Software - * Foundation; either version 3 of the License, or (at your option) any later - * version. - * - * Octave is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - * details. - * - * You should have received a copy of the GNU General Public License along with - * Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. - */ +/* + +Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> +Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ #ifdef HAVE_CONFIG_H #include <config.h> @@ -25,49 +28,46 @@ #include "defun-dld.h" #include "parse.h" -// Secondary functions for complex and real case used -// in ichol algorithms. +// Secondary functions for complex and real case used in ichol algorithms. Complex ichol_mult_complex (Complex a, Complex b) { b.imag (-std::imag (b)); return a * b; } +double ichol_mult_real (double a, double b) +{ + return a * b; +} + bool ichol_checkpivot_complex (Complex pivot) { if (pivot.imag () != 0) { - error ("ichol: Non-real pivot encountered. \ - The matrix must be hermitian."); + error ("ichol: non-real pivot encountered. The matrix must be hermitian."); return false; } else if (pivot.real () < 0) { - error ("ichol: Non-positive pivot encountered."); + error ("ichol: negative pivot encountered"); return false; } return true; - } bool ichol_checkpivot_real (double pivot) { if (pivot < 0) { - error ("ichol: Non-positive pivot encountered."); + error ("ichol: negative pivot encountered"); return false; } return true; } -double ichol_mult_real (double a, double b) -{ - return a * b; -} - -template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T), +template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T), bool (*ichol_checkpivot) (T)> -void ichol_0 (octave_matrix_t& sm, const std::string michol = "off") +void ichol_0 (octave_matrix_t& sm, const std::string michol = "off") { const octave_idx_type n = sm.cols (); @@ -100,7 +100,7 @@ dropsums[i] = 0; } - // Main loop + // Main loop for (k = 0; k < n; k++) { j1 = cidx[k]; @@ -110,7 +110,7 @@ jrow = Llist [k]; // Iterate over each non-zero element in the actual row. - while (jrow != -1) + while (jrow != -1) { jjrow = Lfirst[jrow]; jend = cidx[jrow+1]; @@ -122,16 +122,15 @@ if (jw != -1) data[jw] -= tl; else - // Because the simetry of the matrix we know the drops - // in the column r are also in the column k. + // Because of the symmetry of the matrix, we know + // the drops in the column r are also in the column k. if (opt == ON) { dropsums[r] -= tl; dropsums[k] -= tl; } } - // Update the linked list and the first entry of the - // actual column. + // Update the linked list and the first entry of the actual column. if ((jjrow + 1) < jend) { Lfirst[jrow]++; @@ -149,7 +148,7 @@ if (ridx[j1] != k) { - error ("ichol: There is a pivot equal to zero."); + error ("ichol: encountered a pivot equal to 0"); break; } @@ -158,13 +157,12 @@ data[cidx[k]] = std::sqrt (data[j1]); - // Update Llist and Lfirst with the k-column information. - // Also scale the column elements by the pivot and reset - // the working array iw. - if (k < (n - 1)) + // Update Llist and Lfirst with the k-column information. Also, + // scale the column elements by the pivot and reset the working array iw. + if (k < (n - 1)) { iw[ridx[j1]] = -1; - for(i = j1 + 1; i < j2; i++) + for (i = j1 + 1; i < j2; i++) { iw[ridx[i]] = -1; data[i] /= data[j1]; @@ -182,8 +180,8 @@ } DEFUN_DLD (__ichol0__, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{L} =} __ichol0__ (@var{A})\n\ -@deftypefnx {Loadable Function} {@var{L} =} __ichol0__ (@var{A}, @var{michol})\n\ +@deftypefn {Loadable Function} {@var{L} =} __ichol0__ (@var{A})\n\ +@deftypefnx {Loadable Function} {@var{L} =} __ichol0__ (@var{A}, @var{michol})\n\ Undocumented internal function.\n\ @end deftypefn") @@ -192,7 +190,6 @@ int nargin = args.length (); std::string michol = "off"; - if (nargout > 1 || nargin < 1 || nargin > 2) { @@ -200,65 +197,47 @@ return retval; } - if (args(0).is_scalar_type () || ! args(0).is_sparse_type ()) - error ("__ichol0__: 1. parameter must be a sparse square matrix."); - - if (args(0).is_empty ()) - { - retval(0) = octave_value (SparseMatrix ()); - return retval; - } - + if (nargin == 2) + michol = args(1).string_value (); - if (nargin == 2) - { - michol = args(1).string_value (); - if (error_state || ! (michol == "on" || michol == "off")) - error ("__ichol0__: 2. parameter must be 'on' or 'off' character string."); - } - - - if (!error_state) + // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved + // so it's structure does not change during the algorithm. The same input + // matrix is used to build the output matrix due to that fact. + octave_value_list param_list; + if (!args(0).is_complex_type ()) { - // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved so - // it's structure does not change during the algorithm. The same input - // matrix is used to build the output matrix due to that fact. - octave_value_list param_list; - if (!args(0).is_complex_type ()) - { - SparseMatrix sm = args(0).sparse_matrix_value (); - param_list.append (sm); - sm = feval ("tril", param_list)(0).sparse_matrix_value (); - ichol_0 <SparseMatrix, double, ichol_mult_real, ichol_checkpivot_real> (sm, michol); - if (! error_state) - retval(0) = octave_value (sm); - } - else - { - SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); - param_list.append (sm); - sm = feval ("tril", param_list)(0).sparse_complex_matrix_value (); - ichol_0 <SparseComplexMatrix, Complex, ichol_mult_complex, ichol_checkpivot_complex> (sm, michol); - if (! error_state) - retval(0) = octave_value (sm); - } - + SparseMatrix sm = args(0).sparse_matrix_value (); + param_list.append (sm); + sm = feval ("tril", param_list)(0).sparse_matrix_value (); + ichol_0 <SparseMatrix, double, ichol_mult_real, + ichol_checkpivot_real> (sm, michol); + if (! error_state) + retval(0) = sm; + } + else + { + SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); + param_list.append (sm); + sm = feval ("tril", param_list)(0).sparse_complex_matrix_value (); + ichol_0 <SparseComplexMatrix, Complex, ichol_mult_complex, + ichol_checkpivot_complex> (sm, michol); + if (! error_state) + retval(0) = sm; } return retval; } -template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T), +template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T), bool (*ichol_checkpivot) (T)> void ichol_t (const octave_matrix_t& sm, octave_matrix_t& L, const T* cols_norm, const T droptol, const std::string michol = "off") - + { const octave_idx_type n = sm.cols (); - octave_idx_type j, jrow, jend, jjrow, jw, i, k, jj, Llist_len, total_len, w_len, - max_len, ind; - + octave_idx_type j, jrow, jend, jjrow, jw, i, k, jj, Llist_len, total_len, + w_len, max_len, ind; char opt; enum {OFF, ON}; if (michol == "on") @@ -271,13 +250,13 @@ octave_idx_type* ridx = sm.ridx (); T* data = sm.data (); - // Output matrix data structures. Because it is not known the - // final zero pattern of the output matrix due to fill-in elements, - // an heuristic approach has been adopted for memory allocation. The - // size of ridx_out_l and data_out_l is incremented 10% of their actual - // size (nnz(A) in the beginning). If that amount is less than n, their - // size is just incremented in n elements. This way the number of - // reallocations decrease throughout the process, obtaining a good performance. + // Output matrix data structures. Because the final zero pattern pattern of + // the output matrix is not known due to fill-in elements, a heuristic + // approach has been adopted for memory allocation. The size of ridx_out_l + // and data_out_l is incremented 10% of their actual size (nnz (A) in the + // beginning). If that amount is less than n, their size is just incremented + // in n elements. This way the number of reallocations decreases throughout + // the process, obtaining a good performance. max_len = sm.nnz (); max_len += (0.1 * max_len) > n ? 0.1 * max_len : n; Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1)); @@ -295,7 +274,6 @@ std::vector <octave_idx_type> vec; vec.resize (n); - T zero = T (0); cidx_l[0] = cidx[0]; for (i = 0; i < n; i++) @@ -321,7 +299,7 @@ } } jrow = Llist[k]; - while (jrow != -1) + while (jrow != -1) { jjrow = Lfirst[jrow]; jend = cidx_l[jrow+1]; @@ -329,18 +307,17 @@ { j = ridx_l[jj]; // If the element in the j position of the row is zero, - // then it will become non-zero, so we add it to the - // vector that keeps track of non-zero elements in the working row. + // then it will become non-zero, so we add it to the + // vector that tracks non-zero elements in the working row. if (w_data[j] == zero) { - vec[ind] = j; + vec[ind] = j; ind++; } w_data[j] -= ichol_mult (data_l[jj], data_l[jjrow]); - } - // Update the actual column first element and update the - // linked list of the jrow row. + // Update the actual column first element and + // update the linked list of the jrow row. if ((jjrow + 1) < jend) { Lfirst[jrow]++; @@ -362,21 +339,20 @@ ridx_out_l.resize (dim_vector (max_len, 1)); ridx_l = ridx_out_l.fortran_vec (); } - + // The sorting of the non-zero elements of the working column can be - // handled in a couple of ways. The most efficient two I found, are - // keeping the elements in an ordered binary search tree dinamically - // or keep them unsorted in a vector and at the end of the outer - // iteration order them. The last approach exhibit lower execution - // times. + // handled in a couple of ways. The most efficient two I found, are + // keeping the elements in an ordered binary search tree dynamically or + // keep them unsorted in a vector and at the end of the outer iteration + // order them. The last approach exhibits lower execution times. std::sort (vec.begin (), vec.begin () + ind); data_l[total_len] = w_data[k]; ridx_l[total_len] = k; w_len = 1; - // Extract then non-zero elements of working column and drop the - // elements that are lower than droptol * cols_norm[k]. + // Extract the non-zero elements of working column and + // drop the elements that are lower than droptol * cols_norm[k]. for (i = 0; i < ind ; i++) { jrow = vec[i]; @@ -407,29 +383,29 @@ if (data_l[total_len] == zero) { - error ("ichol: There is a pivot equal to zero."); + error ("ichol: encountered a pivot equal to 0"); break; } else if (! ichol_checkpivot (data_l[total_len])) break; - // Once the elements are dropped and compensation of columns - // sums are done, scale the elements by the pivot. + // Once elements are dropped and compensation of column sums are done, + // scale the elements by the pivot. data_l[total_len] = std::sqrt (data_l[total_len]); for (jj = total_len + 1; jj < (total_len + w_len); jj++) data_l[jj] /= data_l[total_len]; total_len += w_len; - // Check if there are too many elements to be indexed with octave_idx_type - // type due to fill-in during the process. + // Check if there are too many elements to be indexed with + // octave_idx_type type due to fill-in during the process. if (total_len < 0) { - error ("ichol: Integer overflow. Too many fill-in elements in L"); + error ("ichol: integer overflow. Too many fill-in elements in L"); break; } cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len; // Update Llist and Lfirst with the k-column information. - if (k < (n - 1)) + if (k < (n - 1)) { Lfirst[k] = cidx_l[k]; if ((Lfirst[k] + 1) < cidx_l[k+1]) @@ -440,8 +416,7 @@ Llist[jjrow] = k; } } - - } + } if (! error_state) { @@ -455,13 +430,12 @@ L.data (i) = data_l[i]; } } - } DEFUN_DLD (__icholt__, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{L} =} __icholt__ (@var{A})\n\ -@deftypefnx {Loadable Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol})\n\ -@deftypefnx {Loadable Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})\n\ +@deftypefn {Loadable Function} {@var{L} =} __icholt__ (@var{A})\n\ +@deftypefnx {Loadable Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol})\n\ +@deftypefnx {Loadable Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})\n\ Undocumented internal function.\n\ @end deftypefn") { @@ -470,7 +444,6 @@ // Default values of parameters std::string michol = "off"; double droptol = 0; - if (nargout > 1 || nargin < 1 || nargin > 3) { @@ -478,69 +451,51 @@ return retval; } - if (args(0).is_scalar_type () || ! args(0).is_sparse_type ()) - error ("__icholt__: 1. parameter must be a sparse square matrix."); + // Don't repeat input validation of arguments done in ichol.m - if (args(0).is_empty ()) - { - retval(0) = octave_value (SparseMatrix ()); - return retval; - } + if (nargin >= 2) + droptol = args(1).double_value (); - if (! error_state && (nargin >= 2)) - { - droptol = args(1).double_value (); - if (error_state || (droptol < 0) || ! args(1).is_real_scalar ()) - error ("__icholt__: 2. parameter must be a positive real scalar."); - } + if (nargin == 3) + michol = args(2).string_value (); - if (! error_state && (nargin == 3)) - { - michol = args(2).string_value (); - if (error_state || !(michol == "on" || michol == "off")) - error ("__icholt__: 3. parameter must be 'on' or 'off' character string."); - } - - if (!error_state) + octave_value_list param_list; + if (! args(0).is_complex_type ()) { - octave_value_list param_list; - if (! args(0).is_complex_type ()) - { - Array <double> cols_norm; - SparseMatrix L; - param_list.append (args(0).sparse_matrix_value ()); - SparseMatrix sm_l = feval ("tril", - param_list)(0).sparse_matrix_value (); - param_list(0) = sm_l; - param_list(1) = 1; - param_list(2) = "cols"; - cols_norm = feval ("norm", param_list)(0).vector_value (); - param_list.clear (); - ichol_t <SparseMatrix, - double, ichol_mult_real, ichol_checkpivot_real> - (sm_l, L, cols_norm.fortran_vec (), droptol, michol); - if (! error_state) - retval(0) = octave_value (L); - } - else - { - Array <Complex> cols_norm; - SparseComplexMatrix L; - param_list.append (args(0).sparse_complex_matrix_value ()); - SparseComplexMatrix sm_l = feval ("tril", - param_list)(0).sparse_complex_matrix_value (); - param_list(0) = sm_l; - param_list(1) = 1; - param_list(2) = "cols"; - cols_norm = feval ("norm", param_list)(0).complex_vector_value (); - param_list.clear (); - ichol_t <SparseComplexMatrix, - Complex, ichol_mult_complex, ichol_checkpivot_complex> - (sm_l, L, cols_norm.fortran_vec (), Complex (droptol), michol); - if (! error_state) - retval(0) = octave_value (L); - } - + Array <double> cols_norm; + SparseMatrix L; + param_list.append (args(0).sparse_matrix_value ()); + SparseMatrix sm_l = + feval ("tril", param_list)(0).sparse_matrix_value (); + param_list(0) = sm_l; + param_list(1) = 1; + param_list(2) = "cols"; + cols_norm = feval ("norm", param_list)(0).vector_value (); + param_list.clear (); + ichol_t <SparseMatrix, + double, ichol_mult_real, ichol_checkpivot_real> + (sm_l, L, cols_norm.fortran_vec (), droptol, michol); + if (! error_state) + retval(0) = L; + } + else + { + Array <Complex> cols_norm; + SparseComplexMatrix L; + param_list.append (args(0).sparse_complex_matrix_value ()); + SparseComplexMatrix sm_l = + feval ("tril", param_list)(0).sparse_complex_matrix_value (); + param_list(0) = sm_l; + param_list(1) = 1; + param_list(2) = "cols"; + cols_norm = feval ("norm", param_list)(0).complex_vector_value (); + param_list.clear (); + ichol_t <SparseComplexMatrix, + Complex, ichol_mult_complex, ichol_checkpivot_complex> + (sm_l, L, cols_norm.fortran_vec (), + Complex (droptol), michol); + if (! error_state) + retval(0) = L; } return retval; @@ -550,3 +505,4 @@ ## No test needed for internal helper function. %!assert (1) */ +
--- a/libinterp/dldfcn/__ilu__.cc Mon Aug 18 12:32:16 2014 +0100 +++ b/libinterp/dldfcn/__ilu__.cc Tue Aug 26 15:27:21 2014 -0700 @@ -1,22 +1,25 @@ -/** - * Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> - * Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> - * - * This file is part of Octave. - * - * Octave is free software; you can redistribute it and/or modify it under the - * terms of the GNU General Public License as published by the Free Software - * Foundation; either version 3 of the License, or (at your option) any later - * version. - * - * Octave is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - * details. - * - * You should have received a copy of the GNU General Public License along with - * Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. - */ +/* + +Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> +Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ #ifdef HAVE_CONFIG_H #include <config.h> @@ -25,14 +28,14 @@ #include "defun-dld.h" #include "parse.h" -// That function implements the IKJ and JKI variants of gaussian elimination to -// perform the ILUTP decomposition. The behaviour is controlled by milu -// parameter. If milu = ['off'|'col'] the JKI version is performed taking -// advantage of CCS format of the input matrix. If milu = 'row' the input matrix -// has to be transposed to obtain the equivalent CRS structure so we can work -// efficiently with rows. In this case IKJ version is used. +// That function implements the IKJ and JKI variants of Gaussian elimination to +// perform the ILUTP decomposition. The behaviour is controlled by milu +// parameter. If milu = ['off'|'col'] the JKI version is performed taking +// advantage of CCS format of the input matrix. If milu = 'row' the input +// matrix has to be transposed to obtain the equivalent CRS structure so we can +// work efficiently with rows. In this case IKJ version is used. template <typename octave_matrix_t, typename T> -void ilu_0 (octave_matrix_t& sm, const std::string milu = "off") +void ilu_0 (octave_matrix_t& sm, const std::string milu = "off") { const octave_idx_type n = sm.cols (); @@ -70,7 +73,7 @@ r = 0; j = j1; jrow = ridx[j]; - while ((jrow < k) && (j <= j2)) + while ((jrow < k) && (j <= j2)) { if (opt == ROW) { @@ -97,23 +100,23 @@ jrow = ridx[j]; } uptr[k] = j; - if(opt != OFF) + if (opt != OFF) data[uptr[k]] -= r; if (opt != ROW) for (jj = uptr[k] + 1; jj < cidx[k+1]; jj++) data[jj] /= data[uptr[k]]; if (k != jrow) { - error ("ilu: Your input matrix has a zero in the diagonal."); + error ("ilu: A has a zero on the diagonal"); break; } if (data[j] == T(0)) { - error ("ilu: There is a pivot equal to zero."); + error ("ilu: encountered a pivot equal to 0"); break; } - for(i = j1; i <= j2; i++) + for (i = j1; i <= j2; i++) iw[ridx[i]] = -1; } if (opt == ROW) @@ -131,7 +134,6 @@ int nargin = args.length (); std::string milu; - if (nargout > 2 || nargin < 1 || nargin > 2) { @@ -139,67 +141,42 @@ return retval; } - if (args (0).is_empty ()) - { - retval(0) = octave_value (SparseMatrix()); - retval(1) = octave_value (SparseMatrix()); - return retval; - } - - if (args(0).is_scalar_type () || ! args(0).is_sparse_type ()) - error ("__ilu0__: 1. parameter must be a sparse square matrix."); - - if (nargin == 2) + // In ILU0 algorithm the zero-pattern of the input matrix is preserved so + // it's structure does not change during the algorithm. The same input + // matrix is used to build the output matrix due to that fact. + octave_value_list param_list; + if (! args(0).is_complex_type ()) { - milu = args(1).string_value (); - if (error_state || !(milu == "row" || milu == "col" || milu == "off")) - error ( - "__ilu0__: 2. parameter must be 'row', 'col' or 'off' character string."); - } - - - if (! error_state) - { - // In ILU0 algorithm the zero-pattern of the input matrix is preserved so - // it's structure does not change during the algorithm. The same input - // matrix is used to build the output matrix due to that fact. - octave_value_list param_list; - if (! args(0).is_complex_type ()) + SparseMatrix sm = args(0).sparse_matrix_value (); + ilu_0 <SparseMatrix, double> (sm, milu); + if (!error_state) { - SparseMatrix sm = args(0).sparse_matrix_value (); - ilu_0 <SparseMatrix, double> (sm, milu); - if (!error_state) - { - param_list.append (sm); - retval(1) = octave_value ( - feval ("triu", param_list)(0).sparse_matrix_value ()); - SparseMatrix eye = feval ("speye", - octave_value_list ( - octave_value (sm.cols ())))(0).sparse_matrix_value (); - param_list.append (-1); - retval(0) = octave_value ( - eye + feval ("tril", param_list)(0).sparse_matrix_value ()); - - } + param_list.append (sm); + retval(1) = feval ("triu", param_list)(0).sparse_matrix_value (); + SparseMatrix eye = + feval ("speye", octave_value_list ( + octave_value (sm.cols ())))(0).sparse_matrix_value (); + param_list.append (-1); + retval(0) = eye + + feval ("tril", param_list)(0).sparse_matrix_value (); } - else + } + else + { + SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); + ilu_0 <SparseComplexMatrix, Complex> (sm, milu); + if (! error_state) { - SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); - ilu_0 <SparseComplexMatrix, Complex> (sm, milu); - if (! error_state) - { - param_list.append (sm); - retval(1) = octave_value ( - feval ("triu", param_list)(0).sparse_complex_matrix_value ()); - SparseComplexMatrix eye = feval ("speye", - octave_value_list ( - octave_value (sm.cols ())))(0).sparse_complex_matrix_value (); - param_list.append (-1); - retval(0) = octave_value (eye + - feval ("tril", param_list)(0).sparse_complex_matrix_value ()); - } + param_list.append (sm); + retval(1) = + feval ("triu", param_list)(0).sparse_complex_matrix_value (); + SparseComplexMatrix eye = + feval ("speye", octave_value_list ( + octave_value (sm.cols ())))(0).sparse_complex_matrix_value (); + param_list.append (-1); + retval(0) = + eye + feval ("tril", param_list)(0).sparse_complex_matrix_value (); } - } return retval; @@ -212,7 +189,7 @@ const std::string milu = "off") { - // Map the strings into chars to faster comparation inside loops + // Map the strings into chars for faster comparing inside loops char opt; enum {OFF, ROW, COL}; if (milu == "row") @@ -264,7 +241,7 @@ OCTAVE_LOCAL_BUFFER (T, cr_sum, n); T zero = T (0); - + cidx_u[0] = cidx_in_u[0]; cidx_l[0] = cidx_in_l[0]; for (i = 0; i < n; i++) @@ -281,8 +258,7 @@ for (k = 0; k < n; k++) { - - // Load the working column and working row + // Load the working column and working row for (i = cidx_in_l[k]; i < cidx_in_l[k+1]; i++) w_data_l[ridx_in_l[i]] = data_in_l[i]; @@ -380,10 +356,10 @@ // Check if the pivot is zero if (data_u[total_len_u] == zero) { - error ("ilu: There is a pivot equal to zero."); - break; + error ("ilu: encountered a pivot equal to 0"); + break; } - + // Scale the elements in L by the pivot for (i = total_len_l ; i < (total_len_l + w_len_l); i++) data_l[i] /= data_u[total_len_u]; @@ -391,24 +367,24 @@ total_len_u += w_len_u; total_len_l += w_len_l; - // Check if there are too many elements to be indexed with octave_idx_type - // type due to fill-in during the process. + // Check if there are too many elements to be indexed with + // octave_idx_type type due to fill-in during the process. if (total_len_l < 0 || total_len_u < 0) { - error ("ilu: Integer overflow. Too many fill-in elements in L or U"); + error ("ilu: integer overflow. Too many fill-in elements in L or U"); break; } cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u; cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l; - // The tricky part of the algorithm. The arrays pointing to the first + // The tricky part of the algorithm. The arrays pointing to the first // working element of each column in the next iteration (Lfirst) or - // the first working element of each row (Ufirst) are updated. Also the - // arrays working as lists cols_list and rows_list are filled with indexes - // pointing to Ufirst and Lfirst respectively. + // the first working element of each row (Ufirst) are updated. Also the + // arrays working as lists cols_list and rows_list are filled with + // indices pointing to Ufirst and Lfirst respectively. // TODO: Maybe the -1 indicating in Ufirst and Lfirst, that no elements - // have to be considered in a certain column or row in next iteration, can - // be removed. It feels safer to me using such an indicator. + // have to be considered in a certain column or row in next iteration, + // can be removed. It feels safer to me using such an indicator. if (k < (n - 1)) { if (w_len_u > 0) @@ -437,7 +413,7 @@ jj = ridx_u[Ufirst[i]]; } } - if (jj == (k + 1)) + if (jj == (k + 1)) { cols_list[cols_list_len] = i; cols_list_len++; @@ -448,7 +424,7 @@ { jj = ridx_l[Lfirst[i]]; if (jj < (k + 1)) - if(Lfirst[i] < (cidx_l[i+1])) + if (Lfirst[i] < (cidx_l[i+1])) { Lfirst[i]++; if (Lfirst[i] == cidx_l[i+1]) @@ -456,7 +432,7 @@ else jj = ridx_l[Lfirst[i]]; } - if (jj == (k + 1)) + if (jj == (k + 1)) { rows_list[rows_list_len] = i; rows_list_len++; @@ -497,7 +473,6 @@ Undocumented internal function.\n\ @end deftypefn") { - octave_value_list retval; int nargin = args.length (); std::string milu = "off"; @@ -509,103 +484,86 @@ return retval; } - // To be matlab compatible - if (args(0).is_empty ()) - { - retval(0) = octave_value (SparseMatrix()); - retval(1) = octave_value (SparseMatrix()); - return retval; - } + // Don't repeat input validation of arguments done in ilu.m + if (nargin >= 2) + droptol = args(1).double_value (); - if (args(0).is_scalar_type () || ! args(0).is_sparse_type ()) - error ("__iluc__: 1. parameter must be a sparse square matrix."); + if (nargin == 3) + milu = args(2).string_value (); - if (! error_state && (nargin >= 2)) + octave_value_list param_list; + if (! args(0).is_complex_type ()) { - droptol = args(1).double_value (); - if (error_state || (droptol < 0) || ! args(1).is_real_scalar ()) - error ("__iluc__: 2. parameter must be a positive real scalar."); + Array<double> cols_norm, rows_norm; + param_list.append (args(0).sparse_matrix_value ()); + SparseMatrix sm_u = feval ("triu", param_list)(0).sparse_matrix_value (); + param_list.append (-1); + SparseMatrix sm_l = feval ("tril", param_list)(0).sparse_matrix_value (); + param_list(1) = "rows"; + rows_norm = feval ("norm", param_list)(0).vector_value (); + param_list(1) = "cols"; + cols_norm = feval ("norm", param_list)(0).vector_value (); + param_list.clear (); + SparseMatrix U; + SparseMatrix L; + ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U, + cols_norm.fortran_vec (), + rows_norm.fortran_vec (), + droptol, milu); + if (! error_state) + { + param_list.append (octave_value (L.cols ())); + SparseMatrix eye = + feval ("speye", param_list)(0).sparse_matrix_value (); + retval(1) = U; + retval(0) = L + eye; + } } - - if (! error_state && (nargin == 3)) + else { - milu = args(2).string_value (); - if (error_state || !(milu == "row" || milu == "col" || milu == "off")) - error ("__iluc__: 3. parameter must be 'row', 'col' or 'off' character string."); + Array<Complex> cols_norm, rows_norm; + param_list.append (args(0).sparse_complex_matrix_value ()); + SparseComplexMatrix sm_u = + feval("triu", param_list)(0).sparse_complex_matrix_value (); + param_list.append (-1); + SparseComplexMatrix sm_l = + feval("tril", param_list)(0).sparse_complex_matrix_value (); + param_list(1) = "rows"; + rows_norm = feval ("norm", param_list)(0).complex_vector_value (); + param_list(1) = "cols"; + cols_norm = feval ("norm", param_list)(0).complex_vector_value (); + param_list.clear (); + SparseComplexMatrix U; + SparseComplexMatrix L; + ilu_crout < SparseComplexMatrix, Complex > + (sm_l, sm_u, L, U, cols_norm.fortran_vec () , + rows_norm.fortran_vec (), Complex (droptol), milu); + if (! error_state) + { + param_list.append (octave_value (L.cols ())); + SparseComplexMatrix eye = + feval ("speye", param_list)(0).sparse_complex_matrix_value (); + retval(1) = U; + retval(0) = L + eye; + } } - if (! error_state) - { - octave_value_list param_list; - if (! args(0).is_complex_type ()) - { - Array<double> cols_norm, rows_norm; - param_list.append (args(0).sparse_matrix_value ()); - SparseMatrix sm_u = feval ("triu", param_list)(0).sparse_matrix_value (); - param_list.append (-1); - SparseMatrix sm_l = feval ("tril", param_list)(0).sparse_matrix_value (); - param_list(1) = "rows"; - rows_norm = feval ("norm", param_list)(0).vector_value (); - param_list(1) = "cols"; - cols_norm = feval ("norm", param_list)(0).vector_value (); - param_list.clear (); - SparseMatrix U; - SparseMatrix L; - ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U, cols_norm.fortran_vec (), - rows_norm.fortran_vec (), droptol, milu); - if (! error_state) - { - param_list.append (octave_value (L.cols ())); - SparseMatrix eye = feval ("speye", param_list)(0).sparse_matrix_value (); - retval(0) = octave_value (L + eye); - retval(1) = octave_value (U); - } - } - else - { - Array<Complex> cols_norm, rows_norm; - param_list.append (args(0).sparse_complex_matrix_value ()); - SparseComplexMatrix sm_u = feval("triu", - param_list)(0).sparse_complex_matrix_value (); - param_list.append (-1); - SparseComplexMatrix sm_l = feval("tril", - param_list)(0).sparse_complex_matrix_value (); - param_list(1) = "rows"; - rows_norm = feval ("norm", param_list)(0).complex_vector_value (); - param_list(1) = "cols"; - cols_norm = feval ("norm", param_list)(0).complex_vector_value (); - param_list.clear (); - SparseComplexMatrix U; - SparseComplexMatrix L; - ilu_crout < SparseComplexMatrix, Complex > - (sm_l, sm_u, L, U, cols_norm.fortran_vec () , - rows_norm.fortran_vec (), Complex (droptol), milu); - if (! error_state) - { - param_list.append (octave_value (L.cols ())); - SparseComplexMatrix eye = feval ("speye", - param_list)(0).sparse_complex_matrix_value (); - retval(0) = octave_value (L + eye); - retval(1) = octave_value (U); - } - } - } return retval; } -// That function implements the IKJ and JKI variants of gaussian elimination -// to perform the ILUTP decomposition. The behaviour is controlled by milu -// parameter. If milu = ['off'|'col'] the JKI version is performed taking -// advantage of CCS format of the input matrix. Row pivoting is performed. -// If milu = 'row' the input matrix has to be transposed to obtain the -// equivalent CRS structure so we can work efficiently with rows. In that +// That function implements the IKJ and JKI variants of gaussian elimination +// to perform the ILUTP decomposition. The behaviour is controlled by milu +// parameter. If milu = ['off'|'col'] the JKI version is performed taking +// advantage of CCS format of the input matrix. Row pivoting is performed. +// If milu = 'row' the input matrix has to be transposed to obtain the +// equivalent CRS structure so we can work efficiently with rows. In that // case IKJ version is used and column pivoting is performed. template <typename octave_matrix_t, typename T> -void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U, - octave_idx_type nnz_u, octave_idx_type nnz_l, T* cols_norm, +void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U, + octave_idx_type nnz_u, octave_idx_type nnz_l, T* cols_norm, Array <octave_idx_type>& perm_vec, const T droptol = T(0), - const T thresh = T(0), const std::string milu = "off", + const T thresh = T(0), const std::string milu = "off", const double udiag = 0) { char opt; @@ -616,7 +574,7 @@ opt = COL; else opt = OFF; - + const octave_idx_type n = sm.cols (); // That is necessary for the JKI (milu = "row") variant. @@ -627,7 +585,7 @@ octave_idx_type* cidx_in = sm.cidx (); octave_idx_type* ridx_in = sm.ridx (); T* data_in = sm.data (); - octave_idx_type jrow, i, j, k, jj, c, total_len_l, total_len_u, p_perm, + octave_idx_type jrow, i, j, k, jj, c, total_len_l, total_len_u, p_perm, max_ind, max_len_l, max_len_u; T tl, aux, maximum; @@ -690,7 +648,7 @@ it = iw_u.begin (); jrow = *it; total_sum = zero; - while ((jrow < k) && (it != iw_u.end ())) + while ((jrow < k) && (it != iw_u.end ())) { if (opt == COL) partial_col_sum = w_data[jrow]; @@ -712,7 +670,7 @@ } else { - tl = data_l[jj] * w_data[jrow]; + tl = data_l[jj] * w_data[jrow]; w_data[p_perm] -= tl; if (opt == COL) partial_col_sum += tl; @@ -727,26 +685,27 @@ } } - // Drop element from the U part in IKJ and L part in JKI - // variant (milu = [col|off]) - if ((std::abs (w_data[jrow]) < (droptol * cols_norm[k])) - && (w_data[jrow] != zero)) - { - if (opt == COL) - total_sum += partial_col_sum; - else if (opt == ROW) - total_sum += partial_row_sum; - w_data[jrow] = zero; - it2 = it; - it++; - iw_u.erase (it2); - jrow = *it; - continue; - } - else - // This is the element scaled by the pivot in the actual iteration - if (opt == ROW) - w_data[jrow] = tl; + // Drop element from the U part in IKJ and L part in JKI + // variant (milu = [col|off]) + if ((std::abs (w_data[jrow]) < (droptol * cols_norm[k])) + && (w_data[jrow] != zero)) + { + if (opt == COL) + total_sum += partial_col_sum; + else if (opt == ROW) + total_sum += partial_row_sum; + w_data[jrow] = zero; + it2 = it; + it++; + iw_u.erase (it2); + jrow = *it; + continue; + } + else + // This is the element scaled by the pivot + // in the actual iteration + if (opt == ROW) + w_data[jrow] = tl; } jrow = *(++it); } @@ -757,14 +716,14 @@ { maximum = std::abs (w_data[k]) / thresh; max_ind = perm[k]; - for (it = iw_l.begin (); it != iw_l.end (); ++it) + for (it = iw_l.begin (); it != iw_l.end (); ++it) { p_perm = iperm[*it]; if (std::abs (w_data[p_perm]) > maximum) { maximum = std::abs (w_data[p_perm]); max_ind = *it; - it2 = it; + it2 = it; } } // If the pivot is not the diagonal element update all. @@ -775,7 +734,7 @@ if (w_data[k] != zero) iw_l.insert (perm[k]); else - iw_u.insert (k); + iw_u.insert (k); // Swap data and update permutation vectors aux = w_data[k]; iperm[perm[p_perm]] = k; @@ -786,13 +745,13 @@ w_data[k] = w_data[p_perm]; w_data[p_perm] = aux; } - - } + + } // Drop elements in the L part in the IKJ and from the U part in the JKI // version. it = iw_l.begin (); - while (it != iw_l.end ()) + while (it != iw_l.end ()) { p_perm = iperm[*it]; if (droptol > zero) @@ -810,14 +769,15 @@ it++; } - // If milu =[row|col] sumation is preserved --> Compensate diagonal element. + // If milu == [row|col] summation is preserved. + // Compensate diagonal element. if (opt != OFF) { if ((total_sum > zero) && (w_data[k] == zero)) iw_u.insert (k); w_data[k] += total_sum; } - + // Check if the pivot is zero and if udiag is active. @@ -832,19 +792,19 @@ } else { - error ("ilu: There is a pivot equal to zero."); + error ("ilu: encountered a pivot equal to 0"); break; } } - // Scale the elements on the L part for IKJ version (milu = [col|off]) + // Scale the elements on the L part for IKJ version (milu = [col|off]) if (opt != ROW) - for (it = iw_l.begin (); it != iw_l.end (); ++it) + for (it = iw_l.begin (); it != iw_l.end (); ++it) { - p_perm = iperm[*it]; - w_data[p_perm] = w_data[p_perm] / w_data[k]; + p_perm = iperm[*it]; + w_data[p_perm] = w_data[p_perm] / w_data[k]; } - + if ((max_len_u - total_len_u) < n) { @@ -891,11 +851,11 @@ } total_len_u += w_len_u; total_len_l += w_len_l; - // Check if there are too many elements to be indexed with octave_idx_type - // type due to fill-in during the process. + // Check if there are too many elements to be indexed with + // octave_idx_type type due to fill-in during the process. if (total_len_l < 0 || total_len_u < 0) { - error ("ilu: Integer overflow. Too many fill-in elements in L or U"); + error ("ilu: Integer overflow. Too many fill-in elements in L or U"); break; } if (opt == ROW) @@ -909,11 +869,11 @@ if (! error_state) { - octave_matrix_t *L_ptr; + octave_matrix_t *L_ptr; octave_matrix_t *U_ptr; octave_matrix_t diag (n, n, n); - - // L and U are interchanged if milu = 'row'. It is a matter + + // L and U are interchanged if milu = 'row'. It is a matter // of nomenclature to re-use code with both IKJ and JKI // versions of the algorithm. if (opt == ROW) @@ -939,7 +899,7 @@ U_ptr->cidx (i) -= i; } - for (i = 0; i < n; i++) + for (i = 0; i < n; i++) { if (opt == ROW) diag.elem (i,i) = data_u[uptr[i]]; @@ -959,7 +919,7 @@ if (opt == ROW) { // The diagonal is removed from L if milu = 'row'. - // That is because is convenient to have it inside + // That is because is convenient to have it inside // the L part to carry out the process. if (ridx_u[j] == i) { @@ -975,7 +935,7 @@ } } - if (opt == ROW) + if (opt == ROW) { U = U.transpose (); // The diagonal, conveniently permuted is added to U @@ -1002,157 +962,133 @@ double droptol, thresh; double udiag; - if (nargout < 2 || nargout > 3 || nargin < 1 || nargin > 5) { print_usage (); return retval; } - // To be matlab compatible - if (args(0).is_empty ()) - { - retval(0) = octave_value (SparseMatrix ()); - retval(1) = octave_value (SparseMatrix ()); - if (nargout == 3) - retval(2) = octave_value (SparseMatrix ()); - return retval; - } + // Don't repeat input validation of arguments done in ilu.m + if (nargin >= 2) + droptol = args(1).double_value (); - if (args(0).is_scalar_type () || ! args(0).is_sparse_type () ) - error ("__ilutp__: 1. parameter must be a sparse square matrix."); + if (nargin >= 3) + thresh = args(2).double_value (); - if (! error_state && (nargin >= 2)) - { - droptol = args(1).double_value (); - if (error_state || (droptol < 0) || ! args(1).is_real_scalar ()) - error ("__ilutp__: 2. parameter must be a positive scalar."); - } + if (nargin >= 4) + milu = args(3).string_value (); - if (! error_state && (nargin >= 3)) - { - thresh = args(2).double_value (); - if (error_state || ! args(2).is_real_scalar () || (thresh < 0) || thresh > 1) - error ("__ilutp__: 3. parameter must be a scalar 0 <= thresh <= 1."); - } + if (nargin == 5) + udiag = args(4).double_value (); - if (! error_state && (nargin >= 4)) - { - milu = args(3).string_value (); - if (error_state || !(milu == "row" || milu == "col" || milu == "off")) - error ("__ilutp__: 4. parameter must be 'row', 'col' or 'off' character string."); - } - - if (! error_state && (nargin == 5)) + octave_value_list param_list; + octave_idx_type nnz_u, nnz_l; + if (! args(0).is_complex_type ()) { - udiag = args(4).double_value (); - if (error_state || ! args(4).is_real_scalar () || ((udiag != 0) - && (udiag != 1))) - error ("__ilutp__: 5. parameter must be a scalar with value 1 or 0."); - } - - if (! error_state) - { - octave_value_list param_list; - octave_idx_type nnz_u, nnz_l; - if (! args(0).is_complex_type ()) + Array <double> rc_norm; + SparseMatrix sm = args(0).sparse_matrix_value (); + param_list.append (sm); + nnz_u = (feval ("triu", param_list)(0).sparse_matrix_value ()).nnz (); + param_list.append (-1); + nnz_l = (feval ("tril", param_list)(0).sparse_matrix_value ()).nnz (); + if (milu == "row") + param_list (1) = "rows"; + else + param_list (1) = "cols"; + rc_norm = feval ("norm", param_list)(0).vector_value (); + param_list.clear (); + Array <octave_idx_type> perm (dim_vector (sm.cols (), 1)); + SparseMatrix U; + SparseMatrix L; + ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l, + rc_norm.fortran_vec (), + perm, droptol, thresh, milu, udiag); + if (! error_state) { - Array <double> rc_norm; - SparseMatrix sm = args(0).sparse_matrix_value (); - param_list.append (sm); - nnz_u = (feval ("triu", param_list)(0).sparse_matrix_value ()).nnz (); - param_list.append (-1); - nnz_l = (feval ("tril", param_list)(0).sparse_matrix_value ()).nnz (); + param_list.append (octave_value (L.cols ())); + SparseMatrix eye = + feval ("speye", param_list)(0).sparse_matrix_value (); if (milu == "row") - param_list (1) = "rows"; + { + if (nargout == 3) + { + retval(2) = eye.index (idx_vector::colon, perm); + retval(1) = U.index (idx_vector::colon, perm); + } + else if (nargout == 2) + retval(1) = U; + retval(0) = L + eye; + } else - param_list (1) = "cols"; - rc_norm = feval ("norm", param_list)(0).vector_value (); - param_list.clear (); - Array <octave_idx_type> perm (dim_vector (sm.cols (), 1)); - SparseMatrix U; - SparseMatrix L; - ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), - perm, droptol, thresh, milu, udiag); - if (! error_state) { - param_list.append (octave_value (L.cols ())); - SparseMatrix eye = feval ("speye", param_list)(0).sparse_matrix_value (); - if (milu == "row") + if (nargout == 3) { - retval(0) = octave_value (L + eye); - if (nargout == 2) - retval(1) = octave_value (U); - else if (nargout == 3) - { - retval(1) = octave_value (U.index (idx_vector::colon, perm)); - retval(2) = octave_value (eye.index (idx_vector::colon, perm)); - } + retval(2) = eye.index (perm, idx_vector::colon); + retval(1) = U; + retval(0) = L.index (perm, idx_vector::colon) + eye; } else { - retval(1) = octave_value (U); - if (nargout == 2) - retval(0) = octave_value (L + eye.index (perm, idx_vector::colon)); - else if (nargout == 3) - { - retval(0) = octave_value (L.index (perm, idx_vector::colon) + eye); - retval(2) = octave_value (eye.index (perm, idx_vector::colon)); - } + retval(1) = U; + retval(0) = L + eye.index (perm, idx_vector::colon); } } } + } + else + { + Array <Complex> rc_norm; + SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); + param_list.append (sm); + nnz_u = + feval ("triu", param_list)(0).sparse_complex_matrix_value ().nnz (); + param_list.append (-1); + nnz_l = + feval ("tril", param_list)(0).sparse_complex_matrix_value ().nnz (); + if (milu == "row") + param_list(1) = "rows"; else + param_list(1) = "cols"; + rc_norm = feval ("norm", param_list)(0).complex_vector_value (); + Array <octave_idx_type> perm (dim_vector (sm.cols (), 1)); + param_list.clear (); + SparseComplexMatrix U; + SparseComplexMatrix L; + ilu_tp < SparseComplexMatrix, Complex> + (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm, + Complex (droptol), Complex (thresh), milu, udiag); + + if (! error_state) { - Array <Complex> rc_norm; - SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); - param_list.append (sm); - nnz_u = feval ("triu", param_list)(0).sparse_complex_matrix_value ().nnz (); - param_list.append (-1); - nnz_l = feval ("tril", param_list)(0).sparse_complex_matrix_value ().nnz (); + param_list.append (octave_value (L.cols ())); + SparseComplexMatrix eye = + feval ("speye", param_list)(0).sparse_complex_matrix_value (); if (milu == "row") - param_list (1) = "rows"; + { + if (nargout == 3) + { + retval(2) = eye.index (idx_vector::colon, perm); + retval(1) = U.index (idx_vector::colon, perm); + } + else if (nargout == 2) + retval(1) = U; + retval(0) = L + eye; + } else - param_list (1) = "cols"; - rc_norm = feval ("norm", param_list)(0).complex_vector_value (); - Array <octave_idx_type> perm (dim_vector (sm.cols (), 1)); - param_list.clear (); - SparseComplexMatrix U; - SparseComplexMatrix L; - ilu_tp < SparseComplexMatrix, Complex> - (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm, - Complex (droptol), Complex (thresh), milu, udiag); - - if (! error_state) { - param_list.append (octave_value (L.cols ())); - SparseComplexMatrix eye = feval ("speye", - param_list)(0).sparse_complex_matrix_value (); - if (milu == "row") + if (nargout == 3) { - retval(0) = octave_value (L + eye); - if (nargout == 2) - retval(1) = octave_value (U); - else if (nargout == 3) - { - retval(1) = octave_value (U.index (idx_vector::colon, perm)); - retval(2) = octave_value (eye.index (idx_vector::colon, perm)); - } + retval(2) = eye.index (perm, idx_vector::colon); + retval(1) = U; + retval(0) = L.index (perm, idx_vector::colon) + eye; } else { - retval(1) = octave_value (U); - if (nargout == 2) - retval(0) = octave_value (L + eye.index (perm, idx_vector::colon)) ; - else if (nargout == 3) - { - retval(0) = octave_value (L.index (perm, idx_vector::colon) + eye); - retval(2) = octave_value (eye.index (perm, idx_vector::colon)); - } + retval(1) = U; + retval(0) = L + eye.index (perm, idx_vector::colon); } } } - } return retval; @@ -1162,3 +1098,4 @@ ## No test needed for internal helper function. %!assert (1) */ +
--- a/libinterp/dldfcn/chol.cc Mon Aug 18 12:32:16 2014 +0100 +++ b/libinterp/dldfcn/chol.cc Tue Aug 26 15:27:21 2014 -0700 @@ -135,7 +135,7 @@ \n\ In general the lower triangular factorization is significantly faster for\n\ sparse matrices.\n\ -@seealso{hess, lu, qr, qz, schur, svd, cholinv, chol2inv, cholupdate, cholinsert, choldelete, cholshift}\n\ +@seealso{hess, lu, qr, qz, schur, svd, ichol, cholinv, chol2inv, cholupdate, cholinsert, choldelete, cholshift}\n\ @end deftypefn") { octave_value_list retval;
--- a/scripts/help/__unimplemented__.m Mon Aug 18 12:32:16 2014 +0100 +++ b/scripts/help/__unimplemented__.m Tue Aug 26 15:27:21 2014 -0700 @@ -650,8 +650,6 @@ "hgexport", "hgsetget", "hgtransform", - "ichol", - "ilu", "im2frame", "im2java", "imapprox",
--- a/scripts/sparse/ichol.m Mon Aug 18 12:32:16 2014 +0100 +++ b/scripts/sparse/ichol.m Tue Aug 26 15:27:21 2014 -0700 @@ -1,335 +1,273 @@ +## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> -## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> ## ## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it under the -## terms of the GNU General Public License as published by the Free Software -## Foundation; either version 3 of the License, or (at your option) any later -## version. -## -## Octave is distributed in the hope that it will be useful, but WITHOUT ANY -## WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -## details. -## -## You should have received a copy of the GNU General Public License along with -## Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} ichol (@var{A}, @var{opts}) +## @deftypefn {Function File} {@var{L} =} ichol (@var{A}) ## @deftypefnx {Function File} {@var{L} =} ichol (@var{A}, @var{opts}) ## -## @code{@var{L} = ichol (@var{A})} performs the incomplete Cholesky -## factorization of A with zero-fill. +## Compute the incomplete Cholesky factorization of the sparse square matrix +## @var{A} with zero-fill. ## -## @code{@var{L} = ichol (@var{A}, @var{opts})} performs the incomplete Cholesky -## factorization of A with options specified by opts. -## -## By default, ichol references the lower triangle of A and produces lower -## triangular factors. +## By default, ichol references the lower triangle of @var{A} and produces +## lower triangular factors. ## ## The factor given by this routine may be useful as a preconditioner for a ## system of linear equations being solved by iterative methods such as -## PCG (Preconditioned conjugate gradient). -## -## ichol works only for sparse square matrices. +## PCG (Preconditioned Conjugate Gradient). ## -## The fields of @var{opts} must be named exactly as shown below. You can -## include any number of these fields in the structure and define them in any -## order. Any additional fields are ignored. Names and specifiers are case -## sensitive. +## The factorization may be modified by passing options in a structure +## @var{opts}. The option name is a field in the structure and the setting +## is the value of field. Names and specifiers are case sensitive. ## ## @table @asis ## @item type ## Type of factorization. -## String indicating which flavor of incomplete Cholesky to perform. Valid -## values of this field are @samp{nofill} and @samp{ict}. The -## @samp{nofill} variant performs incomplete Cholesky with zero-fill [IC(0)]. -## The @samp{ict} variant performs incomplete Cholesky with threshold dropping -## [ICT]. The default value is @samp{nofill}. +## String indicating which flavor of incomplete Cholesky to perform. Valid +## values of this field are @qcode{"nofill"} and @qcode{"ict"}. The +## @qcode{"nofill"} variant performs incomplete Cholesky with zero-fill +## [IC(0)]. The @qcode{"ict"} variant performs incomplete Cholesky with +## threshold dropping [ICT]. The default value is @qcode{"nofill"}. ## ## @item droptol -## Drop tolerance when type is @samp{ict}. -## Nonnegative scalar used as a drop tolerance when performing ICT. Elements +## Drop tolerance when type is @qcode{"ict"}. +## Non-negative scalar used as a drop tolerance when performing ICT@. Elements ## which are smaller in magnitude than a local drop tolerance are dropped from ## the resulting factor except for the diagonal element which is never dropped. ## The local drop tolerance at step j of the factorization is -## @code{norm (@var{A}(j:end, j), 1) * droptol}. @samp{droptol} is ignored if -## @samp{type} is @samp{nofill}. The default value is 0. +## @code{norm (@var{A}(j:end, j), 1) * droptol}. @code{droptol} is ignored if +## @code{type} is @qcode{"nofill"}. The default value is 0. ## ## @item michol -## Indicates whether to perform modified incomplete Cholesky. -## Indicates whether or not modified incomplete Cholesky [MIC] is performed. -## The field may be @samp{on} or @samp{off}. When performing MIC, the diagonal -## is compensated for dropped elements to enforce the relationship +## Indicate whether modified incomplete Cholesky [MIC] is performed. +## The field may be @qcode{"on"} or @qcode{"off"}. When performing MIC, the +## diagonal is compensated for dropped elements to enforce the relationship ## @code{@var{A} * @var{e} = @var{L} * @var{L}' * @var{e}} where -## @code{@var{e} = ones (size (@var{A}, 2), 1))}. The default value is -## @samp{off}. +## @code{@var{e} = ones (columns (@var{A}), 1)}. The default value is +## @qcode{"off"}. ## ## @item diagcomp ## Perform compensated incomplete Cholesky with the specified coefficient. -## Real nonnegative scalar used as a global diagonal shift @code{@var{alpha}} -## in forming the incomplete Cholesky factor. That is, instead of performing -## incomplete Cholesky on @code{@var{A}}, the factorization of -## @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is formed. The default -## value is 0. +## The coefficient is a real non-negative scalar used as a global diagonal +## shift @code{@var{alpha}} in forming the incomplete Cholesky factor. That +## is, instead of performing incomplete Cholesky on @code{@var{A}}, the +## factorization of @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is +## formed. The default value is 0. ## ## @item shape -## Determines which triangle is referenced and returned. -## Valid values are @samp{upper} and @samp{lower}. If @samp{upper} is specified, -## only the upper triangle of @code{@var{A}} is referenced and @code{@var{R}} -## is constructed such that @code{@var{A}} is approximated by -## @code{@var{R}' * @var{R}}. If @samp{lower} is specified, only the lower -## triangle of @code{@var{A}} is referenced and @code{@var{L}} is constructed -## such that @code{@var{A}} is approximated by @code{@var{L} * @var{L}'}. The -## default value is @samp{lower}. +## Determine which triangle is referenced and returned. +## Valid values are @qcode{"upper"} and @qcode{"lower"}. If @qcode{"upper"} +## is specified, only the upper triangle of @code{@var{A}} is referenced and +## @code{@var{R}} is constructed such that @code{@var{A}} is approximated by +## @code{@var{R}' * @var{R}}. If @qcode{"lower"} is specified, only the +## lower triangle of @code{@var{A}} is referenced and @code{@var{L}} is +## constructed such that @code{@var{A}} is approximated by @code{@var{L} * +## @var{L}'}. The default value is @qcode{"lower"}. ## @end table ## -## EXAMPLES +## Examples ## ## The following problem demonstrates how to factorize a sample symmetric ## positive definite matrix with the full Cholesky decomposition and with the ## incomplete one. ## ## @example +## @group ## A = [ 0.37, -0.05, -0.05, -0.07; ## -0.05, 0.116, 0.0, -0.05; ## -0.05, 0.0, 0.116, -0.05; ## -0.07, -0.05, -0.05, 0.202]; -## A = sparse(A); -## nnz(tril (A)) +## A = sparse (A); +## nnz (tril (A)) ## ans = 9 -## L = chol(A, "lower"); +## L = chol (A, "lower"); ## nnz (L) ## ans = 10 ## norm (A - L * L', "fro") / norm (A, "fro") ## ans = 1.1993e-16 -## opts.type = 'nofill'; -## L = ichol(A,opts); +## opts.type = "nofill"; +## L = ichol (A, opts); ## nnz (L) ## ans = 9 ## norm (A - L * L', "fro") / norm (A, "fro") ## ans = 0.019736 +## @end group ## @end example ## -## Another example for decomposition is finite difference matrix to solve a -## boundary value problem on the unit square. +## Another example for decomposition is a finite difference matrix used to +## solve a boundary value problem on the unit square. ## ## @example +## @group ## nx = 400; ny = 200; ## hx = 1 / (nx + 1); hy = 1 / (ny + 1); -## Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2); -## Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2); +## Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)], +## [-1 0 1 ], nx, nx) / (hx ^ 2); +## Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)], +## [-1 0 1 ], ny, ny) / (hy ^ 2); ## A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy); ## nnz (tril (A)) ## ans = 239400 -## opts.type = 'nofill'; +## opts.type = "nofill"; ## L = ichol (A, opts); ## nnz (tril (A)) ## ans = 239400 ## norm (A - L * L', "fro") / norm (A, "fro") ## ans = 0.062327 +## @end group ## @end example ## -## References for the implemented algorithms: +## References for implemented algorithms: ## -## [1] Saad, Yousef. "Preconditioning Techniques." Iterative Methods for Sparse Linear -## Systems. PWS Publishing Company, 1996. +## [1] @nospell{Y. Saad}. "Preconditioning Techniques." @cite{Iterative +## Methods for Sparse Linear Systems}, PWS Publishing Company, 1996. ## -## [2] Jones, Mark T. and Plassmann, Paul E.: An Improved Incomplete Cholesky -## Factorization (1992). +## [2] @nospell{M. Jones, P. Plassmann}: @cite{An Improved Incomplete +## Cholesky Factorization}, 1992. +## @seealso{chol, ilu, pcg} ## @end deftypefn -function [L] = ichol (A, opts) +function L = ichol (A, opts = struct ()) - if ((nargin > 2) || (nargin < 1) || (nargout > 1)) + if (nargin < 1 || nargin > 2 || nargout > 1) print_usage (); endif - % Check input matrix - if (~issparse(A) || ~issquare (A)) - error ("ichol: Input A must be a non-empty sparse square matrix"); + if (! (issparse (A) && issquare (A))) + error ("ichol: A must be a sparse square matrix"); endif - % If A is empty and sparse then return empty L - % Compatibility with Matlab + if (! isstruct (opts)) + error ("ichol: OPTS must be a structure."); + endif + + ## If A is empty then return empty L for Matlab compatibility if (isempty (A)) L = A; return; endif - % Check input structure, otherwise set default values - if (nargin == 2) - if (~isstruct (opts)) - error ("ichol: Input \"opts\" must be a valid structure."); - endif - else - opts = struct (); - endif - - if (~isfield (opts, "type")) - opts.type = "nofill"; % set default + ## Parse input options + if (! isfield (opts, "type")) + opts.type = "nofill"; # set default else type = tolower (getfield (opts, "type")); - if ((strcmp (type, "nofill") == 0) - && (strcmp (type, "ict") == 0)) - error ("ichol: Invalid field \"type\" in input structure."); - else - opts.type = type; + if (! strcmp (type, "nofill") && ! strcmp (type, "ict")) + error ('ichol: TYPE must be "nofill" or "ict"'); endif + opts.type = type; endif - if (~isfield (opts, "droptol")) - opts.droptol = 0; % set default + if (! isfield (opts, "droptol")) + opts.droptol = 0; # set default else - if (~isscalar (opts.droptol) || (opts.droptol < 0)) - error ("ichol: Invalid field \"droptol\" in input structure."); + if (! (isreal (opts.droptol) && isscalar (opts.droptol) + && opts.droptol >= 0)) + error ("ichol: DROPTOL must be a non-negative real scalar"); endif endif michol = ""; - if (~isfield (opts, "michol")) - opts.michol = "off"; % set default + if (! isfield (opts, "michol")) + opts.michol = "off"; # set default else michol = tolower (getfield (opts, "michol")); - if ((strcmp (michol, "off") == 0) - && (strcmp (michol, "on") == 0)) - error ("ichol: Invalid field \"michol\" in input structure."); - else - opts.michol = michol; + if (! strcmp (michol, "off") && ! strcmp (michol, "on")) + error ('ichol: MICHOL must be "on" or "off"'); + endif + opts.michol = michol; + endif + + if (! isfield (opts, "diagcomp")) + opts.diagcomp = 0; # set default + else + if (! (isreal (opts.diagcomp) && isscalar (opts.diagcomp) + && opts.diagcomp >= 0)) + error ("ichol: DIAGCOMP must be a non-negative real scalar"); endif endif - if (~isfield (opts, "diagcomp")) - opts.diagcomp = 0; % set default + if (! isfield (opts, "shape")) + opts.shape = "lower"; # set default else - if (~isscalar (opts.diagcomp) || (opts.diagcomp < 0)) - error ("ichol: Invalid field \"diagcomp\" in input structure."); + shape = tolower (getfield (opts, "shape")); + if (! strcmp (shape, "lower") && ! strcmp (shape, "upper")) + error ('ichol: SHAPE must be "lower" or "upper"'); endif + opts.shape = shape; endif - if (~isfield (opts, "shape")) - opts.shape = "lower"; % set default - else - shape = tolower (getfield (opts, "shape")); - if ((strcmp (shape, "lower") == 0) - && (strcmp (shape, "upper") == 0)) - error ("ichol: Invalid field \"shape\" in input structure."); - else - opts.shape = shape; - endif - endif - - % Prepare input for specialized ICHOL + ## Prepare input for specialized ICHOL A_in = []; if (opts.diagcomp > 0) A += opts.diagcomp * diag (diag (A)); endif - if (strcmp (opts.shape, "upper") == 1) + if (strcmp (opts.shape, "upper")) A_in = triu (A); A_in = A_in'; else A_in = tril (A); endif - % Delegate to specialized ICHOL + ## Delegate to specialized ICHOL switch (opts.type) case "nofill" L = __ichol0__ (A_in, opts.michol); case "ict" L = __icholt__ (A_in, opts.droptol, opts.michol); - otherwise - printf ("The input structure is invalid.\n"); endswitch - if (strcmp (opts.shape, "upper") == 1) + if (strcmp (opts.shape, "upper")) L = L'; endif +endfunction -endfunction %!shared A1, A2, A3, A4, A5, A6, A7 %! A1 = [ 0.37, -0.05, -0.05, -0.07; -%! -0.05, 0.116, 0.0, -0.05; -%! -0.05, 0.0, 0.116, -0.05; -%! -0.07, -0.05, -0.05, 0.202]; -%! A1 = sparse(A1); -%! A2 = gallery ('poisson', 30); -%! A3 = gallery ('tridiag', 50); +%! -0.05, 0.116, 0.0, -0.05; +%! -0.05, 0.0, 0.116, -0.05; +%! -0.07, -0.05, -0.05, 0.202]; +%! A1 = sparse (A1); +%! A2 = gallery ("poisson", 30); +%! A3 = gallery ("tridiag", 50); %! nx = 400; ny = 200; %! hx = 1 / (nx + 1); hy = 1 / (ny + 1); -%! Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2); -%! Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2); +%! Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)], +%! [-1 0 1 ], nx, nx) / (hx ^ 2); +%! Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)], +%! [-1 0 1 ], ny, ny) / (hy ^ 2); %! A4 = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy); -%! A5 = [ 0.37, -0.05, -0.05, -0.07; -%! -0.05, 0.116, 0.0, -0.05 + 0.05i; -%! -0.05, 0.0, 0.116, -0.05; -%! -0.07, -0.05 - 0.05i, -0.05, 0.202]; -%! A5 = sparse(A5); -%! A6 = [ 0.37, -0.05 - i, -0.05, -0.07; -%! -0.05 + i, 0.116, 0.0, -0.05; -%! -0.05, 0.0, 0.116, -0.05; -%! -0.07, -0.05, -0.05, 0.202]; -%! A6 = sparse(A6); +%! A5 = [ 0.37, -0.05, -0.05, -0.07; +%! -0.05, 0.116, 0.0, -0.05 + 0.05i; +%! -0.05, 0.0, 0.116, -0.05; +%! -0.07, -0.05 - 0.05i, -0.05, 0.202]; +%! A5 = sparse (A5); +%! A6 = [ 0.37, -0.05 - i, -0.05, -0.07; +%! -0.05 + i, 0.116, 0.0, -0.05; +%! -0.05, 0.0, 0.116, -0.05; +%! -0.07, -0.05, -0.05, 0.202]; +%! A6 = sparse (A6); %! A7 = A5; %! A7(1) = 2i; -%! -%!# Input validation tests - -%!test -%!error ichol ([]); -%!error ichol (0); -%!error ichol (-0); -%!error ichol (1); -%!error ichol (-1); -%!error ichol (i); -%!error ichol (-i); -%!error ichol (1 + 1i); -%!error ichol (1 - 1i); -%!error ichol (sparse (0)); -%!error ichol (sparse (-0)); -%!error ichol (sparse (-1)); -%!error ichol (sparse (-1)); -%!test -%! opts.milu = 'foo'; -%!error L = ichol (A1, opts); -%! opts.milu = 1; -%!error L = ichol (A1, opts); -%! opts.milu = []; -%!error L = ichol (A1, opts); -%!test -%! opts.droptol = -1; -%!error L = ichol (A1, opts); -%! opts.droptol = 0.5i; -%!error L = ichol (A1, opts); -%! opts.droptol = []; -%!error L = ichol (A1, opts); -%!test -%! opts.type = 'foo'; -%!error L = ichol (A1, opts); -%! opts.type = 1; -%!error L = ichol (A1, opts); -%! opts.type = []; -%!error L = ichol (A1, opts); -%!test -%! opts.shape = 'foo'; -%!error L = ichol (A1, opts); -%! opts.shape = 1; -%!error L = ichol (A1, opts); -%! opts.shape = []; -%!error L = ichol (A1, opts); -%!test -%! opts.diagcomp = -1; -%!error L = ichol (A1, opts); -%! opts.diagcomp = 0.5i; -%!error L = ichol (A1, opts); -%! opts.diagcomp = []; -%!error L = ichol (A1, opts); - -%!# ICHOL0 tests +## ICHOL0 tests %!test %! opts.type = "nofill"; @@ -401,13 +339,14 @@ %! opts.michol = "on"; %! L = ichol (A5, opts); %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4); -%!test -%% Negative pivot -%!error ichol (A6); -%% Complex entry in the diagonal -%!error ichol (A7); -%%!ICHOLT tests +## Negative pivot +%!error <negative pivot> ichol (A6) +%!error ichol (A6) +## Complex entry in the diagonal +%!error <non-real pivot> ichol (A7) + +## ICHOLT tests %%!test %! opts.type = "ict"; @@ -475,9 +414,47 @@ %! opts.michol = "on"; %! L = ichol (A5, opts); %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.3231, 1e-4); + +%% Input validation tests + +%!error <A must be a sparse square matrix> ichol ([]) +%!error <A must be a sparse square matrix> ichol (0) +%!error <pivot equal to 0> ichol (sparse (0)) +%!error <pivot equal to 0> ichol (sparse (-0)) +%!error <negative pivot> ichol (sparse (-1)) %!test -%% Negative pivot -%! opts.type = "ict"; -%!error ichol (A6, setup); -%% Complex entry in the diagonal -%!error ichol (A7, setup); +%! opts.type = "foo"; +%! fail ("ichol (A1, opts)", 'TYPE must be "nofill"'); +%! opts.type = 1; +%! fail ("ichol (A1, opts)", 'TYPE must be "nofill"'); +%! opts.type = []; +%! fail ("ichol (A1, opts)", 'TYPE must be "nofill"'); +%!test +%! opts.droptol = -1; +%! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar"); +%! opts.droptol = 0.5i; +%! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar"); +%! opts.droptol = []; +%! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar"); +%!test +%! opts.michol = "foo"; +%! fail ("ichol (A1, opts)", 'MICHOL must be "on"'); +%! opts.michol = 1; +%! fail ("ichol (A1, opts)", 'MICHOL must be "on"'); +%! opts.michol = []; +%! fail ("ichol (A1, opts)", 'MICHOL must be "on"'); +%!test +%! opts.diagcomp = -1; +%! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar"); +%! opts.diagcomp = 0.5i; +%! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar"); +%! opts.diagcomp = []; +%! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar"); +%!test +%! opts.shape = "foo"; +%! fail ("ichol (A1, opts)", 'SHAPE must be "lower"'); +%! opts.shape = 1; +%! fail ("ichol (A1, opts)", 'SHAPE must be "lower"'); +%! opts.shape = []; +%! fail ("ichol (A1, opts)", 'SHAPE must be "lower"'); +
--- a/scripts/sparse/ilu.m Mon Aug 18 12:32:16 2014 +0100 +++ b/scripts/sparse/ilu.m Tue Aug 26 15:27:21 2014 -0700 @@ -1,96 +1,103 @@ +## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> -## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> -## ## ## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it under the -## terms of the GNU General Public License as published by the Free Software -## Foundation; either version 3 of the License, or (at your option) any later -## version. -## -## Octave is distributed in the hope that it will be useful, but WITHOUT ANY -## WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -## details. -## -## You should have received a copy of the GNU General Public License along with -## Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. - +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. +## ## -*- texinfo -*- -## @deftypefn {Function File} ilu (@var{A}, @var{setup}) -## @deftypefnx {Function File} {[@var{L}, @var{U}] =} ilu (@var{A}, @var{setup}) -## @deftypefnx {Function File} {[@var{L}, @var{U}, @var{P}] =} ilu (@var{A}, @var{setup}) -## ilu produces a unit lower triangular matrix, an upper triangular matrix, and -## a permutation matrix. +## @deftypefn {Function File} {@var{L}, @var{U}] =} ilu (@var{A}) +## @deftypefnx {Function File} {@var{L}, @var{U}] =} ilu (@var{A}, @var{opts}) +## @deftypefnx {Function File} {[@var{L}, @var{U}, @var{P}] =} ilu (@dots{}) +## +## Compute the incomplete LU factorization of the sparse square matrix @var{A} +## into a unit lower triangular matrix (@var{L}), an upper triangular matrix +## (@var{U}), and a permutation matrix (@var{P}). ## -## These incomplete factorizations may be useful as preconditioners for a system -## of linear equations being solved by iterative methods such as BICG -## (BiConjugate Gradients), GMRES (Generalized Minimum Residual Method). +## These incomplete factorizations may be useful as preconditioners for a +## system of linear equations being solved by iterative methods such as BICG +## (BiConjugate Gradients) or GMRES (Generalized Minimum Residual Method). ## -## @code{ilu (@var{A}, @var{setup})} computes the incomplete LU factorization -## of @var{A}. @var{setup} is an input structure with up to five setup options. -## The fields must be named exactly as shown below. You can include any number -## of these fields in the structure and define them in any order. Any -## additional fields are ignored. +## The factorization may be modified by passing options in a structure +## @var{opts}. The option name is a field in the structure and the setting +## is the value of field. Names and specifiers are case sensitive. +## +## @table @code +## @item type +## Type of factorization. Values for type include: ## ## @table @asis -## @item type -## Type of factorization. Values for type include: +## @item @qcode{"nofill"} +## Perform ILU factorization with 0 level of fill in, known as ILU(0). With +## type set to @qcode{"nofill"}, only the @code{milu} option is used; all other +## fields are ignored. ## -## @table @asis -## @item @samp{nofill} -## Performs ILU factorization with 0 level of fill in, known as ILU(0). With -## type set to @samp{nofill}, only the milu setup option is used; all other -## fields are ignored. -## @item @samp{crout} -## Performs the Crout version of ILU factorization, known as ILUC. With type -## set to @samp{crout}, only the droptol and milu setup options are used; all -## other fields are ignored. -## @item @samp{ilutp} +## @item @qcode{"crout"} +## Perform the Crout version of ILU factorization, known as ILUC@. With type +## set to @qcode{crout}, only the @code{droptol} and @code{milu} options are +## used; all other fields are ignored. +## +## @item @qcode{"ilutp"} ## (default) Performs ILU factorization with threshold and pivoting. ## @end table ## ## If type is not specified, the ILU factorization with pivoting ILUTP is -## performed. Pivoting is never performed with type set to @samp{nofill} or -## @samp{crout}. +## performed. Pivoting is never performed with type set to @qcode{"nofill"} or +## @qcode{"crout"}. ## ## @item droptol -## Drop tolerance of the incomplete LU factorization. droptol is a non-negative -## scalar. The default value is 0, which produces the complete LU factorization. +## Drop tolerance of the incomplete LU factorization. @code{droptol} is a +## non-negative scalar. The default value is 0, which produces the complete +## LU factorization. ## -## The nonzero entries of U satisfy +## The nonzero entries of @var{U} satisfy ## ## @code{abs (@var{U}(i,j)) >= droptol * norm ((@var{A}:,j))} ## ## with the exception of the diagonal entries, which are retained regardless of -## satisfying the criterion. The entries of @var{L} are tested against the +## satisfying the criterion. The entries of @var{L} are tested against the ## local drop tolerance before being scaled by the pivot, so for nonzeros in ## @var{L} ## -## @code{abs(@var{L}(i,j)) >= droptol * norm(@var{A}(:,j))/@var{U}(j,j)}. +## @code{abs (@var{L}(i,j)) >= droptol * norm (@var{A}(:,j))/@var{U}(j,j)}. ## ## @item milu -## Modified incomplete LU factorization. Values for milu +## Modified incomplete LU factorization. Values for @code{milu} ## include: +## ## @table @asis -## @item @samp{row} -## Produces the row-sum modified incomplete LU factorization. Entries from the +## @item @qcode{"row"} +## Produces the row-sum modified incomplete LU factorization. Entries from the ## newly-formed column of the factors are subtracted from the diagonal of the -## upper triangular factor, @var{U}, preserving column sums. That is, +## upper triangular factor, @var{U}, preserving column sums. That is, ## @code{@var{A} * e = @var{L} * @var{U} * e}, where e is the vector of ones. -## @item @samp{col} -## Produces the column-sum modified incomplete LU factorization. Entries from +## +## @item @qcode{"col"} +## Produces the column-sum modified incomplete LU factorization. Entries from ## the newly-formed column of the factors are subtracted from the diagonal of -## the upper triangular factor, @var{U}, preserving column sums. That is, +## the upper triangular factor, @var{U}, preserving column sums. That is, ## @code{e'*@var{A} = e'*@var{L}*@var{U}}. -## @item @samp{off} +## +## @item @qcode{"off"} ## (default) No modified incomplete LU factorization is produced. ## @end table ## ## @item udiag -## If udiag is 1, any zeros on the diagonal of the upper -## triangular factor are replaced by the local drop tolerance. The default is 0. +## If @code{udiag} is 1, any zeros on the diagonal of the upper +## triangular factor are replaced by the local drop tolerance. The default +## is 0. ## ## @item thresh ## Pivot threshold between 0 (forces diagonal pivoting) and 1, @@ -103,161 +110,151 @@ ## lower triangular matrix and @var{U} is an upper triangular matrix. ## ## @code{[@var{L}, @var{U}] = ilu (@var{A},@var{setup})} returns a unit lower -## triangular matrix in @var{L} and an upper triangular matrix in @var{U}. When -## SETUP.type = 'ilutp', the role of @var{P} is determined by the value of -## SETUP.milu. For SETUP.type == 'ilutp', one of the factors is permuted -## based on the value of SETUP.milu. When SETUP.milu == 'row', U is a column -## permuted upper triangular factor. Otherwise, L is a row-permuted unit lower -## triangular factor. +## triangular matrix in @var{L} and an upper triangular matrix in @var{U}. When +## SETUP.type = @qcode{"ilutp"}, the role of @var{P} is determined by the +## value of SETUP.milu. For SETUP.type == @qcode{"ilutp"}, one of the +## factors is permuted based on the value of SETUP.milu. When SETUP.milu == +## @qcode{"row"}, U is a column permuted upper triangular factor. Otherwise, +## L is a row-permuted unit lower triangular factor. ## ## @code{[@var{L}, @var{U}, @var{P}] = ilu (@var{A},@var{setup})} returns a ## unit lower triangular matrix in @var{L}, an upper triangular matrix in -## @var{U}, and a permutation matrix in @var{P}. When SETUP.milu ~= 'row', @var{P} -## is returned such that @var{L} and @var{U} are incomplete factors of @var{P}*@var{A}. -## When SETUP.milu == 'row', @var{P} is returned such that and @var{U} are -## incomplete factors of A*P. +## @var{U}, and a permutation matrix in @var{P}. When @code{milu} ! = +## @qcode{"row"}, @var{P} is returned such that @var{L} and @var{U} are +## incomplete factors of @var{P}*@var{A}. When @code{milu} == @qcode{"row"}, +## @var{P} is returned such that and @var{U} are incomplete factors of A*P. ## -## @strong{NOTE}: ilu works on sparse square matrices only. -## -## EXAMPLES +## Examples ## ## @example -## A = gallery('neumann', 1600) + speye(1600); -## setup.type = 'nofill'; -## nnz(A) +## @group +## A = gallery ("neumann", 1600) + speye (1600); +## opts.type = "nofill"; +## nnz (A) ## ans = 7840 ## -## nnz(lu(A)) +## nnz (lu (A)) ## ans = 126478 ## -## nnz(ilu(A,setup)) +## nnz (ilu (A, opts)) ## ans = 7840 +## @end group ## @end example ## ## This shows that @var{A} has 7840 nonzeros, the complete LU factorization has ## 126478 nonzeros, and the incomplete LU factorization, with 0 level of -## fill-in, has 7840 nonzeros, the same amount as @var{A}. Taken from: +## fill-in, has 7840 nonzeros, the same amount as @var{A}. Taken from: ## http://www.mathworks.com/help/matlab/ref/ilu.html ## ## @example -## A = gallery ('wathen', 10, 10); -## b = sum (A,2); +## @group +## A = gallery ("wathen", 10, 10); +## b = sum (A, 2); ## tol = 1e-8; ## maxit = 50; -## opts.type = 'crout'; +## opts.type = "crout"; ## opts.droptol = 1e-4; ## [L, U] = ilu (A, opts); ## x = bicg (A, b, tol, maxit, L, U); -## norm(A * x - b, inf) +## norm (A * x - b, inf) +## @end group ## @end example ## ## This example uses ILU as preconditioner for a random FEM-Matrix, which has a -## bad condition. Without @var{L} and @var{U} BICG would not converge. +## large condition number. Without @var{L} and @var{U} BICG would not converge. ## +## @seealso{lu, ichol, bicg, gmres} ## @end deftypefn -function [L, U, P] = ilu (A, setup) +function [L, U, P] = ilu (A, opts = struct ()) - if ((nargin > 2) || (nargin < 1) || (nargout > 3)) + if (nargin < 1 || nargin > 2 || (nargout > 3)) print_usage (); endif - - % Check input matrix - if (~issparse (A) || ~issquare (A)) - error ("ilu: Input A must be a sparse square matrix."); + if (! (issparse (A) && issquare (A))) + error ("ichol: A must be a sparse square matrix"); endif - % If A is empty and sparse then return empty L, U and P - % Compatibility with Matlab + if (! isstruct (opts)) + error ("ichol: OPTS must be a structure."); + endif + + ## If A is empty then return empty L, U and P for Matlab compatibility if (isempty (A)) - L = A; - U = A; - P = A; + L = U = P = A; return; endif - % Check input structure, otherwise set default values - if (nargin == 2) - if (~isstruct (setup)) - error ("ilu: Input 'setup' must be a valid structure."); + ## Parse input options + if (! isfield (opts, "type")) + opts.type = "nofill"; # set default + else + type = tolower (getfield (opts, "type")); + if (! any (strcmp (type, {"nofill", "crout", "ilutp"}))) + error ("ilu: invalid TYPE specified"); endif - else - setup = struct (); + opts.type = type; endif - if (~isfield (setup, "type")) - setup.type = "nofill"; % set default + if (! isfield (opts, "droptol")) + opts.droptol = 0; # set default else - type = tolower (getfield (setup, "type")); - if ((strcmp (type, "nofill") == 0) - && (strcmp (type, "crout") == 0) - && (strcmp (type, "ilutp") == 0)) - error ("ilu: Invalid field \"type\" in input structure."); - else - setup.type = type; + if (! (isreal (opts.droptol) && isscalar (opts.droptol) + && opts.droptol >= 0)) + error ("ilu: DROPTOL must be a non-negative real scalar"); endif endif - if (~isfield (setup, "droptol")) - setup.droptol = 0; % set default + if (! isfield (opts, "milu")) + opts.milu = "off"; # set default else - if (~isscalar (setup.droptol) || (setup.droptol < 0)) - error ("ilu: Invalid field \"droptol\" in input structure."); + milu = tolower (getfield (opts, "milu")); + if (! any (strcmp (milu, {"off", "col", "row"}))) + error ('ilu: MILU must be one of "off", "col", or "row"'); + endif + opts.milu = milu; + endif + + if (! isfield (opts, "udiag")) + opts.udiag = 0; # set default + else + if (! isscalar (opts.udiag) || (opts.udiag != 0 && opts.udiag != 1)) + error ("ilu: UDIAG must be 0 or 1"); endif endif - if (~isfield (setup, "milu")) - setup.milu = "off"; % set default + if (! isfield (opts, "thresh")) + opts.thresh = 1; # set default else - milu = tolower (getfield (setup, "milu")); - if ((strcmp (milu, "off") == 0) - && (strcmp (milu, "col") == 0) - && (strcmp (milu, "row") == 0)) - error ("ilu: Invalid field \"milu\" in input structure."); - else - setup.milu = milu; - endif - endif - - if (~isfield (setup, "udiag")) - setup.udiag = 0; % set default - else - if (~isscalar (setup.udiag) || ((setup.udiag ~= 0) && (setup.udiag ~= 1))) - error ("ilu: Invalid field \"udiag\" in input structure."); - endif - endif - - if (~isfield (setup, "thresh")) - setup.thresh = 1; % set default - else - if (~isscalar (setup.thresh) || (setup.thresh < 0) || (setup.thresh > 1)) - error ("ilu: Invalid field \"thresh\" in input structure."); + if (! (isreal (opts.thresh) && isscalar (opts.thresh)) + || opts.thresh < 0 || opts.thresh > 1) + error ("ilu: THRESH must be a scalar in the range [0, 1]"); endif endif n = length (A); - % Delegate to specialized ILU - switch (setup.type) + ## Delegate to specialized ILU + switch (opts.type) case "nofill" - [L, U] = __ilu0__ (A, setup.milu); + [L, U] = __ilu0__ (A, opts.milu); if (nargout == 3) P = speye (length (A)); endif case "crout" - [L, U] = __iluc__ (A, setup.droptol, setup.milu); + [L, U] = __iluc__ (A, opts.droptol, opts.milu); if (nargout == 3) P = speye (length (A)); endif case "ilutp" if (nargout == 2) - [L, U] = __ilutp__ (A, setup.droptol, setup.thresh, setup.milu, setup.udiag); + [L, U] = __ilutp__ (A, opts.droptol, opts.thresh, + opts.milu, opts.udiag); elseif (nargout == 3) - [L, U, P] = __ilutp__ (A, setup.droptol, setup.thresh, setup.milu, setup.udiag); + [L, U, P] = __ilutp__ (A, opts.droptol, opts.thresh, + opts.milu, opts.udiag); endif - otherwise - printf ("The input structure is invalid.\n"); endswitch if (nargout == 1) @@ -266,120 +263,63 @@ endfunction + %!shared n, dtol, A %! n = 1600; %! dtol = 0.1; -%! A = gallery ('neumann', n) + speye (n); +%! A = gallery ("neumann", n) + speye (n); %!test -%! setup.type = 'nofill'; -%! assert (nnz (ilu (A, setup)), 7840); -%! # This test is taken from the mathworks and should work for full support. +%! opts.type = "nofill"; +%! assert (nnz (ilu (A, opts)), 7840); +## This test has been verified in both Matlab and Octave. %!test -%! setup.type = 'crout'; -%! setup.milu = 'row'; -%! setup.droptol = dtol; -%! [L, U] = ilu (A, setup); +%! opts.type = "crout"; +%! opts.milu = "row"; +%! opts.droptol = dtol; +%! [L, U] = ilu (A, opts); %! e = ones (size (A, 2),1); %! assert (norm (A*e - L*U*e), 1e-14, 1e-14); %!test -%! setup.type = 'crout'; -%! setup.droptol = dtol; -%! [L, U] = ilu(A, setup); -%! assert (norm (A - L * U, 'fro') / norm (A, 'fro'), 0.05, 1e-2); +%! opts.type = "crout"; +%! opts.droptol = dtol; +%! [L, U] = ilu (A, opts); +%! assert (norm (A - L * U, "fro") / norm (A, "fro"), 0.05, 1e-2); -%! # Tests for input validation -%!test -%! [L, U] = ilu (sparse ([])); -%! assert (isempty (L), logical (1)); -%! assert (isempty (U), logical (1)); -%! setup.type = 'crout'; -%! [L, U] = ilu (sparse ([]), setup); -%! assert (isempty (L), logical (1)); -%! assert (isempty (U), logical (1)); -%! setup.type = 'ilutp'; -%! [L, U] = ilu (sparse ([]), setup); -%! assert (isempty (L), logical (1)); -%! assert (isempty (U), logical (1)); -%!error [L, U] = ilu (0); -%!error [L, U] = ilu ([]); -%!error [L, U] = ilu (sparse (0)); -%!test -%! setup.type = 'foo'; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.type = 1; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.type = []; -%!error [L, U] = ilu (A_tiny, setup); +## Check if the elements in U satisfy the non-dropping condition. %!test -%! setup.droptol = -1; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.droptol = 0.5i; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.droptol = []; -%!error [L, U] = ilu (A_tiny, setup); -%!test -%! setup.thresh= -1; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.thresh = 0.5i; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.thresh = []; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.thresh = 2; -%!error [L, U] = ilu (A_tiny, setup); -%!test -%! setup.diag = 0.5; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.diag = []; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.diag = -1; -%!error [L, U] = ilu (A_tiny, setup); -%!test -%! setup.milu = 'foo'; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.milu = 1; -%!error [L, U] = ilu (A_tiny, setup); -%! setup.milu = []; -%!error [L, U] = ilu (A_tiny, setup); - -%! # Check if the elements in U satisfy the non-dropping condition. -%!test -%! setup.type = 'crout'; -%! setup.droptol = dtol; -%! [L, U] = ilu (A, setup); +%! opts.type = "crout"; +%! opts.droptol = dtol; +%! [L, U] = ilu (A, opts); %! for j = 1:n %! cmp_value = dtol * norm (A(:, j)); %! non_zeros = nonzeros (U(:, j)); -%! for i = 1:length (non_zeros); -%! assert (abs (non_zeros (i)) >= cmp_value, logical (1)); -%! endfor +%! assert (abs (non_zeros) >= cmp_value); %! endfor %!test -%! setup.type = 'ilutp'; -%! setup.droptol = dtol; -%! [L, U] = ilu (A, setup); +%! opts.type = "ilutp"; +%! opts.droptol = dtol; +%! [L, U] = ilu (A, opts); %! for j = 1:n %! cmp_value = dtol * norm (A(:, j)); %! non_zeros = nonzeros (U(:, j)); -%! for i = 1:length (non_zeros); -%! assert (abs (non_zeros (i)) >= cmp_value, logical (1)); -%! endfor +%! assert (abs (non_zeros) >= cmp_value); %! endfor -%! # Check that the complete LU factorisation with crout and ilutp algorithms -%! # output the same result. +## Check that the complete LU factorisation with crout and ilutp algorithms +## produce the same result. %!test -%! setup.type = 'crout'; -%! setup.droptol = 0; -%! [L1, U1] = ilu (A, setup); -%! setup.type = 'ilutp'; -%! setup.thresh = 0; -%! [L2, U2] = ilu (A, setup); -%! assert (norm (L1 - L2, 'fro') / norm (L1, 'fro'), 0, eps); -%! assert (norm (U1 - U2, 'fro') / norm (U1, 'fro'), 0, eps); +%! opts.type = "crout"; +%! opts.droptol = 0; +%! [L1, U1] = ilu (A, opts); +%! opts.type = "ilutp"; +%! opts.thresh = 0; +%! [L2, U2] = ilu (A, opts); +%! assert (norm (L1 - L2, "fro") / norm (L1, "fro"), 0, eps); +%! assert (norm (U1 - U2, "fro") / norm (U1, "fro"), 0, eps); -%! # Tests for real matrices of different sizes for ilu0, iluc and ilutp. -%! # The difference A - L*U should be not greater than eps because with droptol -%! # equaling 0, the LU complete factorization is performed. +## Tests for real matrices of different sizes for ilu0, iluc and ilutp. +## The difference A - L*U should be not greater than eps because with droptol +## equaling 0, the LU complete factorization is performed. %!shared n_tiny, n_small, n_medium, n_large, A_tiny, A_small, A_medium, A_large %! n_tiny = 5; %! n_small = 40; @@ -391,66 +331,65 @@ %! A_large = sprand (n_large, n_large, 1/n_large/10) + speye (n_large); %! %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_tiny); -%! assert (norm (A_tiny - L * U, "fro") / norm (A_tiny, "fro"), 0, n_tiny * eps); +%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), 0, n_tiny * eps); %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_small); -%! assert (norm (A_small - L * U, "fro") / norm (A_small, "fro"), 0, 1); +%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), 0, 1); %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_medium); -%! assert (norm (A_medium - L * U, "fro") / norm (A_medium, "fro"), 0, 1); +%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), 0, 1); %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_large); -%! assert (norm (A_large - L * U, "fro") / norm (A_large, "fro"), 0, 1); +%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), 0, 1); %! %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_tiny, setup); -%! assert (norm (A_tiny - L * U, "fro") / norm (A_tiny, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_tiny, opts); +%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_small, setup); -%! assert (norm (A_small - L * U, "fro") / norm (A_small, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_small, opts); +%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_medium, setup); -%! assert (norm (A_medium - L * U, "fro") / norm (A_medium, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_medium, opts); +%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_large, setup); -%! assert (norm (A_large - L * U, "fro") / norm (A_large, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_large, opts); +%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); %! %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_tiny, setup); -%! assert (norm (A_tiny - L * U, "fro") / norm (A_tiny, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_tiny, opts); +%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_small, setup); -%! assert (norm (A_small - L * U, "fro") / norm (A_small, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_small, opts); +%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_medium, setup); -%! assert (norm (A_medium - L * U, "fro") / norm (A_medium, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_medium, opts); +%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_large, setup); -%! assert (norm (A_large - L * U, "fro") / norm (A_large, "fro"), eps, eps); -%! +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_large, opts); +%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); -%! # Tests for complex matrices of different sizes for ilu0, iluc and ilutp. +## Tests for complex matrices of different sizes for ilu0, iluc and ilutp. %!shared n_tiny, n_small, n_medium, n_large, A_tiny, A_small, A_medium, A_large %! n_tiny = 5; %! n_small = 40; @@ -458,90 +397,147 @@ %! n_large = 10000; %! A_tiny = spconvert ([1 4 2 3 3 4 2 5; 1 1 2 3 4 4 5 5; 1 2 3 4 5 6 7 8]'); %! A_tiny(1,1) += 1i; -%! A_small = sprand(n_small, n_small, 1/n_small) + ... -%! i * sprand(n_small, n_small, 1/n_small) + speye (n_small); -%! A_medium = sprand(n_medium, n_medium, 1/n_medium) + ... -%! i * sprand(n_medium, n_medium, 1/n_medium) + speye (n_medium); -%! A_large = sprand(n_large, n_large, 1/n_large/10) + ... -%! i * sprand(n_large, n_large, 1/n_large/10) + speye (n_large); +%! A_small = sprand (n_small, n_small, 1/n_small) + ... +%! i * sprand (n_small, n_small, 1/n_small) + speye (n_small); +%! A_medium = sprand (n_medium, n_medium, 1/n_medium) + ... +%! i * sprand (n_medium, n_medium, 1/n_medium) + speye (n_medium); +%! A_large = sprand (n_large, n_large, 1/n_large/10) + ... +%! i * sprand (n_large, n_large, 1/n_large/10) + speye (n_large); %! %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_tiny); -%! assert (norm (A_tiny - L * U, "fro") / norm (A_tiny, "fro"), 0, n_tiny * eps); +%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), 0, n_tiny * eps); %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_small); -%! assert (norm (A_small - L * U, "fro") / norm (A_small, "fro"), 0, 1); +%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), 0, 1); %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_medium); -%! assert (norm (A_medium - L * U, "fro") / norm (A_medium, "fro"), 0, 1); +%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), 0, 1); %!test -%! setup.type = "nofill"; +%! opts.type = "nofill"; %! [L, U] = ilu (A_large); -%! assert (norm (A_large - L * U, "fro") / norm (A_large, "fro"), 0, 1); +%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), 0, 1); %! %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_tiny, setup); -%! assert (norm (A_tiny - L * U, "fro") / norm (A_tiny, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_tiny, opts); +%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_small, setup); -%! assert (norm (A_small - L * U, "fro") / norm (A_small, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_small, opts); +%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_medium, setup); -%! assert (norm (A_medium - L * U, "fro") / norm (A_medium, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_medium, opts); +%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); %!test -%! setup.type = "crout"; -%! [L, U] = ilu (A_large, setup); -%! assert (norm (A_large - L * U, "fro") / norm (A_large, "fro"), eps, eps); +%! opts.type = "crout"; +%! [L, U] = ilu (A_large, opts); +%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); %! %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_tiny, setup); -%! assert (norm (A_tiny - L * U, "fro") / norm (A_tiny, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_tiny, opts); +%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_small, setup); -%! assert (norm (A_small - L * U, "fro") / norm (A_small, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_small, opts); +%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_medium, setup); -%! assert (norm (A_medium - L * U, "fro") / norm (A_medium, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_medium, opts); +%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! [L, U] = ilu (A_large, setup); -%! assert (norm (A_large - L * U, "fro") / norm (A_large, "fro"), eps, eps); +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! [L, U] = ilu (A_large, opts); +%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); -%! #Specific tests for ilutp +## Specific tests for ilutp + %!shared a1, a2 %! a1 = sparse ([0 0 4 3 1; 5 1 2.3 2 4.5; 0 0 0 2 1;0 0 8 0 2.2; 0 0 9 9 1 ]); %! a2 = sparse ([3 1 0 0 4; 3 1 0 0 -2;0 0 8 0 0; 0 4 0 4 -4.5; 0 -1 0 0 1]); %!test -%! setup.udiag = 1; -%! setup.type = "ilutp"; -%! setup.droptol = 0.2; -%! [L, U, P] = ilu (a1, setup); +%! opts.udiag = 1; +%! opts.type = "ilutp"; +%! opts.droptol = 0.2; +%! [L, U, P] = ilu (a1, opts); %! assert (norm (U, "fro"), 17.4577, 1e-4); %! assert (norm (L, "fro"), 2.4192, 1e-4); -%! setup.udiag = 0; -%!error [L, U, P] = ilu (a1, setup); +%! opts.udiag = 0; +%! #fail ("ilu (a1, opts)"); %! %!test -%! setup.type = "ilutp"; -%! setup.droptol = 0; -%! setup.thresh = 0; -%! setup.milu = "row"; -%!error [L, U] = ilu (a2, setup); -%! +%! opts.type = "ilutp"; +%! opts.droptol = 0; +%! opts.thresh = 0; +%! opts.milu = "row"; +%! #fail ("ilu (a2, opts)"); + +%% Tests for input validation +%!shared A_tiny +%! A_tiny = spconvert ([1 4 2 3 3 4 2 5; 1 1 2 3 4 4 5 5; 1 2 3 4 5 6 7 8]'); + +%!test +%! [L, U] = ilu (sparse ([])); +%! assert (isempty (L)); +%! assert (isempty (U)); +%! opts.type = "crout"; +%! [L, U] = ilu (sparse ([]), opts); +%! assert (isempty (L)); +%! assert (isempty (U)); +%! opts.type = "ilutp"; +%! [L, U] = ilu (sparse ([]), opts); +%! assert (isempty (L)); +%! assert (isempty (U)); +%!error <A must be a sparse square matrix> ilu (0) +%!error <A must be a sparse square matrix> ilu ([]) +%!error <zero on the diagonal> ilu (sparse (0)) + +%!test +%! opts.type = "foo"; +%! fail ("ilu (A_tiny, opts)", "invalid TYPE specified"); +%! opts.type = 1; +%! fail ("ilu (A_tiny, opts)", "invalid TYPE specified"); +%! opts.type = []; +%! fail ("ilu (A_tiny, opts)", "invalid TYPE specified"); +%!test +%! opts.droptol = -1; +%! fail ("ilu (A_tiny, opts)", "DROPTOL must be a non-negative real scalar"); +%! opts.droptol = 0.5i; +%! fail ("ilu (A_tiny, opts)", "DROPTOL must be a non-negative real scalar"); +%! opts.droptol = []; +%! fail ("ilu (A_tiny, opts)", "DROPTOL must be a non-negative real scalar"); +%!test +%! opts.milu = "foo"; +%! fail ("ilu (A_tiny, opts)", 'MILU must be one of "off"'); +%! opts.milu = 1; +%! fail ("ilu (A_tiny, opts)", 'MILU must be one of "off"'); +%! opts.milu = []; +%! fail ("ilu (A_tiny, opts)", 'MILU must be one of "off"'); +%!test +%! opts.udiag = -1; +%! fail ("ilu (A_tiny, opts)", "UDIAG must be 0 or 1"); +%! opts.udiag = 0.5i; +%! fail ("ilu (A_tiny, opts)", "UDIAG must be 0 or 1"); +%! opts.udiag = []; +%! fail ("ilu (A_tiny, opts)", "UDIAG must be 0 or 1"); +%!test +%! opts.thresh = -1; +%! fail ("ilu (A_tiny, opts)", "THRESH must be a scalar"); +%! opts.thresh = 0.5i; +%! fail ("ilu (A_tiny, opts)", "THRESH must be a scalar"); +%! opts.thresh = []; +%! fail ("ilu (A_tiny, opts)", "THRESH must be a scalar"); +