1######################################################################## 2## 3## Copyright (C) 2013-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {} ilu (@var{A}) 28## @deftypefnx {} {} ilu (@var{A}, @var{opts}) 29## @deftypefnx {} {[@var{L}, @var{U}] =} ilu (@dots{}) 30## @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} ilu (@dots{}) 31## 32## Compute the incomplete LU factorization of the sparse square matrix @var{A}. 33## 34## @code{ilu} returns a unit lower triangular matrix @var{L}, an upper 35## triangular matrix @var{U}, and optionally a permutation matrix @var{P}, such 36## that @code{@var{L}*@var{U}} approximates @code{@var{P}*@var{A}}. 37## 38## The factors given by this routine may be useful as preconditioners for a 39## system of linear equations being solved by iterative methods such as BICG 40## (BiConjugate Gradients) or GMRES (Generalized Minimum Residual Method). 41## 42## The factorization may be modified by passing options in a structure 43## @var{opts}. The option name is a field of the structure and the setting 44## is the value of field. Names and specifiers are case sensitive. 45## 46## @table @code 47## @item type 48## Type of factorization. 49## 50## @table @asis 51## @item @qcode{"nofill"} (default) 52## ILU factorization with no fill-in (ILU(0)). 53## 54## Additional supported options: @code{milu}. 55## 56## @item @qcode{"crout"} 57## Crout version of ILU factorization (@nospell{ILUC}). 58## 59## Additional supported options: @code{milu}, @code{droptol}. 60## 61## @item @qcode{"ilutp"} 62## ILU factorization with threshold and pivoting. 63## 64## Additional supported options: @code{milu}, @code{droptol}, @code{udiag}, 65## @code{thresh}. 66## @end table 67## 68## @item droptol 69## A non-negative scalar specifying the drop tolerance for factorization. The 70## default value is 0 which produces the complete LU factorization. 71## 72## Non-diagonal entries of @var{U} are set to 0 unless 73## 74## @code{abs (@var{U}(i,j)) >= droptol * norm (@var{A}(:,j))}. 75## 76## Non-diagonal entries of @var{L} are set to 0 unless 77## 78## @code{abs (@var{L}(i,j)) >= droptol * norm (@var{A}(:,j))/@var{U}(j,j)}. 79## 80## @item milu 81## Modified incomplete LU factorization: 82## 83## @table @asis 84## @item @qcode{"row"} 85## Row-sum modified incomplete LU factorization. 86## The factorization preserves row sums: 87## @code{@var{A} * e = @var{L} * @var{U} * e}, where e is a vector of ones. 88## 89## @item @qcode{"col"} 90## Column-sum modified incomplete LU factorization. 91## The factorization preserves column sums: 92## @code{e' * @var{A} = e' * @var{L} * @var{U}}. 93## 94## @item @qcode{"off"} (default) 95## Row and column sums are not necessarily preserved. 96## @end table 97## 98## @item udiag 99## If true, any zeros on the diagonal of the upper triangular factor are 100## replaced by the local drop tolerance 101## @code{droptol * norm (@var{A}(:,j))/@var{U}(j,j)}. The default is false. 102## 103## @item thresh 104## Pivot threshold for factorization. It can range between 0 (diagonal 105## pivoting) and 1 (default), where the maximum magnitude entry in the column 106## is chosen to be the pivot. 107## @end table 108## 109## If @code{ilu} is called with just one output, the returned matrix is 110## @code{@var{L} + @var{U} - speye (size (@var{A}))}, where @var{L} is unit 111## lower triangular and @var{U} is upper triangular. 112## 113## With two outputs, @code{ilu} returns a unit lower triangular matrix @var{L} 114## and an upper triangular matrix @var{U}. For @var{opts}.type == 115## @qcode{"ilutp"}, one of the factors is permuted based on the value of 116## @var{opts}.milu. When @var{opts}.milu == @qcode{"row"}, @var{U} is a 117## column permuted upper triangular factor. Otherwise, @var{L} is a 118## row-permuted unit lower triangular factor. 119## 120## If there are three named outputs and @var{opts}.milu != @qcode{"row"}, 121## @var{P} is returned such that @var{L} and @var{U} are incomplete factors 122## of @code{@var{P}*@var{A}}. When @var{opts}.milu == @qcode{"row"}, @var{P} 123## is returned such that @var{L} and @var{U} are incomplete factors of 124## @code{@var{A}*@var{P}}. 125## 126## EXAMPLES 127## 128## @example 129## @group 130## A = gallery ("neumann", 1600) + speye (1600); 131## opts.type = "nofill"; 132## nnz (A) 133## ans = 7840 134## 135## nnz (lu (A)) 136## ans = 126478 137## 138## nnz (ilu (A, opts)) 139## ans = 7840 140## @end group 141## @end example 142## 143## This shows that @var{A} has 7,840 nonzeros, the complete LU factorization 144## has 126,478 nonzeros, and the incomplete LU factorization, with 0 level of 145## fill-in, has 7,840 nonzeros, the same amount as @var{A}. Taken from: 146## @url{https://www.mathworks.com/help/matlab/ref/ilu.html} 147## 148## @example 149## @group 150## A = gallery ("wathen", 10, 10); 151## b = sum (A, 2); 152## tol = 1e-8; 153## maxit = 50; 154## opts.type = "crout"; 155## opts.droptol = 1e-4; 156## [L, U] = ilu (A, opts); 157## x = bicg (A, b, tol, maxit, L, U); 158## norm (A * x - b, inf) 159## @end group 160## @end example 161## 162## This example uses ILU as preconditioner for a random FEM-Matrix, which has a 163## large condition number. Without @var{L} and @var{U} BICG would not 164## converge. 165## 166## @seealso{lu, ichol, bicg, gmres} 167## @end deftypefn 168 169function [L, U, P] = ilu (A, opts = struct ()) 170 171 if (nargin < 1 || nargin > 2 || (nargout > 3)) 172 print_usage (); 173 endif 174 175 if (! (issparse (A) && issquare (A))) 176 error ("ilu: A must be a sparse square matrix"); 177 endif 178 179 if (! isstruct (opts)) 180 error ("ilu: OPTS must be a structure."); 181 endif 182 183 ## If A is empty then return empty L, U and P for Matlab compatibility 184 if (isempty (A)) 185 L = U = P = A; 186 return; 187 endif 188 189 ## Parse input options 190 if (! isfield (opts, "type")) 191 opts.type = "nofill"; # set default 192 else 193 type = tolower (getfield (opts, "type")); 194 if (! any (strcmp (type, {"nofill", "crout", "ilutp"}))) 195 error ("ilu: invalid TYPE specified"); 196 endif 197 opts.type = type; 198 endif 199 200 if (! isfield (opts, "droptol")) 201 opts.droptol = 0; # set default 202 else 203 if (! (isreal (opts.droptol) && isscalar (opts.droptol) 204 && opts.droptol >= 0)) 205 error ("ilu: DROPTOL must be a non-negative real scalar"); 206 endif 207 endif 208 209 if (! isfield (opts, "milu")) 210 opts.milu = "off"; # set default 211 else 212 milu = tolower (getfield (opts, "milu")); 213 if (! any (strcmp (milu, {"off", "col", "row"}))) 214 error ('ilu: MILU must be one of "off", "col", or "row"'); 215 endif 216 opts.milu = milu; 217 endif 218 219 if (! isfield (opts, "udiag")) 220 opts.udiag = 0; # set default 221 else 222 if (! isscalar (opts.udiag) || (opts.udiag != 0 && opts.udiag != 1)) 223 error ("ilu: UDIAG must be 0 or 1"); 224 endif 225 endif 226 227 if (! isfield (opts, "thresh")) 228 opts.thresh = 1; # set default 229 else 230 if (! (isreal (opts.thresh) && isscalar (opts.thresh)) 231 || opts.thresh < 0 || opts.thresh > 1) 232 error ("ilu: THRESH must be a scalar in the range [0, 1]"); 233 endif 234 endif 235 236 n = length (A); 237 238 ## Delegate to specialized ILU 239 switch (opts.type) 240 case "nofill" 241 [L, U] = __ilu0__ (A, opts.milu); 242 if (nargout == 3) 243 P = speye (length (A)); 244 endif 245 case "crout" 246 [L, U] = __iluc__ (A, opts.droptol, opts.milu); 247 if (nargout == 3) 248 P = speye (length (A)); 249 endif 250 case "ilutp" 251 if (nargout == 3) 252 [L, U, P] = __ilutp__ (A, opts.droptol, opts.thresh, 253 opts.milu, opts.udiag); 254 else 255 [L, U] = __ilutp__ (A, opts.droptol, opts.thresh, 256 opts.milu, opts.udiag); 257 endif 258 endswitch 259 260 if (nargout == 1) 261 L = L + U - speye (n); 262 endif 263 264endfunction 265 266 267%!shared n, dtol, A 268%! n = 1600; 269%! dtol = 0.1; 270%! A = gallery ("neumann", n) + speye (n); 271 272%!test 273%! opts.type = "nofill"; 274%! assert (nnz (ilu (A, opts)), 7840); 275 276## This test has been verified in both Matlab and Octave. 277%!test 278%! opts.type = "crout"; 279%! opts.milu = "row"; 280%! opts.droptol = dtol; 281%! [L, U] = ilu (A, opts); 282%! e = ones (size (A, 2),1); 283%! assert (norm (A*e - L*U*e), 1e-14, 1e-14); 284%!test 285%! opts.type = "crout"; 286%! opts.droptol = dtol; 287%! [L, U] = ilu (A, opts); 288%! assert (norm (A - L * U, "fro") / norm (A, "fro"), 0.05, 1e-2); 289 290## Check if the elements in U satisfy the non-dropping condition. 291%!test 292%! opts.type = "crout"; 293%! opts.droptol = dtol; 294%! [L, U] = ilu (A, opts); 295%! for j = 1:n 296%! cmp_value = dtol * norm (A(:, j)); 297%! non_zeros = nonzeros (U(:, j)); 298%! assert (abs (non_zeros) >= cmp_value); 299%! endfor 300%!test 301%! opts.type = "ilutp"; 302%! opts.droptol = dtol; 303%! [L, U] = ilu (A, opts); 304%! for j = 1:n 305%! cmp_value = dtol * norm (A(:, j)); 306%! non_zeros = nonzeros (U(:, j)); 307%! assert (abs (non_zeros) >= cmp_value); 308%! endfor 309 310## Check that the complete LU factorisation with crout and ilutp algorithms 311## produce the same result. 312%!test 313%! opts.type = "crout"; 314%! opts.droptol = 0; 315%! [L1, U1] = ilu (A, opts); 316%! opts.type = "ilutp"; 317%! opts.thresh = 0; 318%! [L2, U2] = ilu (A, opts); 319%! assert (norm (L1 - L2, "fro") / norm (L1, "fro"), 0, eps); 320%! assert (norm (U1 - U2, "fro") / norm (U1, "fro"), 0, eps); 321 322## Restore rand "state" value 323%!shared old_rand_state, restore_state 324%! ## Save and restore the state of the random number generator that is used by 325%! ## the unit tests in this file. 326%! old_rand_state = rand ("state"); 327%! restore_state = onCleanup (@() rand ("state", old_rand_state)); 328 329## Tests for real matrices of different sizes for ilu0, iluc and ilutp. 330## The difference A - L*U should be not greater than eps because with droptol 331## equal to 0, the LU complete factorization is performed. 332%!shared n_tiny, n_small, n_medium, n_large, A_tiny, A_small, A_medium, A_large 333%! n_tiny = 5; 334%! n_small = 40; 335%! n_medium = 600; 336%! n_large = 10000; 337%! 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]'); 338%! ## initialize generator to make behavior reproducible 339%! rand ("state", 42); 340%! A_small = sprand (n_small, n_small, 1/n_small) + speye (n_small); 341%! A_medium = sprand (n_medium, n_medium, 1/n_medium) + speye (n_medium); 342%! A_large = sprand (n_large, n_large, 1/n_large/10) + speye (n_large); 343 344%!test 345%! opts.type = "nofill"; 346%! [L, U] = ilu (A_tiny); 347%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), 0, n_tiny * eps); 348%!test 349%! opts.type = "nofill"; 350%! [L, U] = ilu (A_small); 351%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), 0, 1); 352%!test 353%! opts.type = "nofill"; 354%! [L, U] = ilu (A_medium); 355%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), 0, 1); 356%!test 357%! opts.type = "nofill"; 358%! [L, U] = ilu (A_large); 359%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), 0, 1); 360 361%!test 362%! opts.type = "crout"; 363%! [L, U] = ilu (A_tiny, opts); 364%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); 365%!test 366%! opts.type = "crout"; 367%! [L, U] = ilu (A_small, opts); 368%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); 369%!test 370%! opts.type = "crout"; 371%! [L, U] = ilu (A_medium, opts); 372%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); 373%!test 374%! opts.type = "crout"; 375%! [L, U] = ilu (A_large, opts); 376%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); 377 378%!test 379%! opts.type = "ilutp"; 380%! opts.droptol = 0; 381%! opts.thresh = 0; 382%! [L, U] = ilu (A_tiny, opts); 383%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); 384%!test 385%! opts.type = "ilutp"; 386%! opts.droptol = 0; 387%! opts.thresh = 0; 388%! [L, U] = ilu (A_small, opts); 389%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); 390%!test 391%! opts.type = "ilutp"; 392%! opts.droptol = 0; 393%! opts.thresh = 0; 394%! [L, U] = ilu (A_medium, opts); 395%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); 396%!test 397%! opts.type = "ilutp"; 398%! opts.droptol = 0; 399%! opts.thresh = 0; 400%! [L, U] = ilu (A_large, opts); 401%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); 402 403## Tests for complex matrices of different sizes for ilu0, iluc and ilutp. 404%!shared n_tiny, n_small, n_medium, n_large, A_tiny, A_small, A_medium, A_large 405%! n_tiny = 5; 406%! n_small = 40; 407%! n_medium = 600; 408%! n_large = 10000; 409%! 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]'); 410%! A_tiny(1,1) += 1i; 411%! ## initialize generator to make behavior reproducible 412%! rand ("state", 42); 413%! A_small = sprand (n_small, n_small, 1/n_small) + ... 414%! i * sprand (n_small, n_small, 1/n_small) + speye (n_small); 415%! A_medium = sprand (n_medium, n_medium, 1/n_medium) + ... 416%! i * sprand (n_medium, n_medium, 1/n_medium) + speye (n_medium); 417%! A_large = sprand (n_large, n_large, 1/n_large/10) + ... 418%! i * sprand (n_large, n_large, 1/n_large/10) + speye (n_large); 419 420%!test 421%! opts.type = "nofill"; 422%! [L, U] = ilu (A_tiny); 423%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), 0, n_tiny * eps); 424%!test 425%! opts.type = "nofill"; 426%! [L, U] = ilu (A_small); 427%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), 0, 1); 428%!test 429%! opts.type = "nofill"; 430%! [L, U] = ilu (A_medium); 431%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), 0, 1); 432%!test 433%! opts.type = "nofill"; 434%! [L, U] = ilu (A_large); 435%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), 0, 1); 436 437%!test 438%! opts.type = "crout"; 439%! [L, U] = ilu (A_tiny, opts); 440%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); 441%!test 442%! opts.type = "crout"; 443%! [L, U] = ilu (A_small, opts); 444%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); 445%!test 446%! opts.type = "crout"; 447%! [L, U] = ilu (A_medium, opts); 448%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); 449%!test 450%! opts.type = "crout"; 451%! [L, U] = ilu (A_large, opts); 452%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); 453 454%!test 455%! opts.type = "ilutp"; 456%! opts.droptol = 0; 457%! opts.thresh = 0; 458%! [L, U] = ilu (A_tiny, opts); 459%! assert (norm (A_tiny - L*U, "fro") / norm (A_tiny, "fro"), eps, eps); 460%!test 461%! opts.type = "ilutp"; 462%! opts.droptol = 0; 463%! opts.thresh = 0; 464%! [L, U] = ilu (A_small, opts); 465%! assert (norm (A_small - L*U, "fro") / norm (A_small, "fro"), eps, eps); 466%!test 467%! opts.type = "ilutp"; 468%! opts.droptol = 0; 469%! opts.thresh = 0; 470%! [L, U] = ilu (A_medium, opts); 471%! assert (norm (A_medium - L*U, "fro") / norm (A_medium, "fro"), eps, eps); 472%!test 473%! opts.type = "ilutp"; 474%! opts.droptol = 0; 475%! opts.thresh = 0; 476%! [L, U] = ilu (A_large, opts); 477%! assert (norm (A_large - L*U, "fro") / norm (A_large, "fro"), eps, eps); 478 479## Specific tests for ilutp 480%!shared A 481%! A = 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 ]); 482 483%!test 484%! opts.udiag = 1; 485%! opts.type = "ilutp"; 486%! opts.droptol = 0.2; 487%! [L, U, P] = ilu (A, opts); 488%! assert (norm (U, "fro"), 17.4577, 1e-4); 489%! assert (norm (L, "fro"), 2.4192, 1e-4); 490 491%!error <encountered a pivot equal to 0> 492%! opts.type = "ilutp"; 493%! opts.udiag = 0; 494%! opts.droptol = 0.2; 495%! ilu (A, opts); 496 497## Matlab R2017b doesn't error, but returns a singular L which isn't helpful. 498%!error <encountered a pivot equal to 0> 499%! A = 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]); 500%! opts.type = "ilutp"; 501%! opts.droptol = 0; 502%! opts.thresh = 0; 503%! opts.milu = "row"; 504%! [L, U, P] = ilu (A, opts); 505 506%!test <*53440> 507%! A = sparse (magic (4)); 508%! opts.type = "ilutp"; 509%! [L, U] = ilu (A, opts); 510%! assert (L * U, A, eps) 511 512## Tests for input validation 513%!shared A_tiny, opts 514%! 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]'); 515 516%!test 517%! [L, U] = ilu (sparse ([])); 518%! assert (isempty (L)); 519%! assert (isempty (U)); 520%! opts.type = "crout"; 521%! [L, U] = ilu (sparse ([]), opts); 522%! assert (isempty (L)); 523%! assert (isempty (U)); 524%! opts.type = "ilutp"; 525%! [L, U] = ilu (sparse ([]), opts); 526%! assert (isempty (L)); 527%! assert (isempty (U)); 528 529%!error <A must be a sparse square matrix> ilu (0) 530%!error <A must be a sparse square matrix> ilu ([]) 531%!error <zero on the diagonal> ilu (sparse (0)) 532 533%!error <invalid TYPE specified> 534%! opts.type = "foo"; 535%! ilu (A_tiny, opts); 536%!error <invalid TYPE specified> 537%! opts.type = 1; 538%! ilu (A_tiny, opts); 539%!error <invalid TYPE specified> 540%! opts.type = []; 541%! ilu (A_tiny, opts); 542 543%!error <DROPTOL must be a non-negative real scalar> 544%! clear opts; 545%! opts.droptol = -1; 546%! ilu (A_tiny, opts); 547%!error <DROPTOL must be a non-negative real scalar> 548%! opts.droptol = 0.5i; 549%! ilu (A_tiny, opts); 550%!error <DROPTOL must be a non-negative real scalar> 551%! opts.droptol = []; 552%! ilu (A_tiny, opts); 553 554%!error <MILU must be one of "off", "col", or "row"> 555%! clear opts; 556%! opts.milu = "foo"; 557%! ilu (A_tiny, opts); 558%!error <MILU must be one of "off", "col", or "row"> 559%! opts.milu = 1; 560%! ilu (A_tiny, opts); 561%!error <MILU must be one of "off", "col", or "row"> 562%! opts.milu = []; 563%! ilu (A_tiny, opts); 564 565%!error <UDIAG must be 0 or 1> 566%! clear opts; 567%! opts.udiag = -1; 568%! ilu (A_tiny, opts); 569%!error <UDIAG must be 0 or 1> 570%! opts.udiag = 0.5i; 571%! ilu (A_tiny, opts); 572%!error <UDIAG must be 0 or 1> 573%! opts.udiag = []; 574%! ilu (A_tiny, opts); 575 576%!error <THRESH must be a scalar in the range \[0, 1\]> 577%! clear opts; 578%! opts.thresh = -1; 579%! ilu (A_tiny, opts); 580%!error <THRESH must be a scalar in the range \[0, 1\]> 581%! opts.thresh = 0.5i; 582%! ilu (A_tiny, opts); 583%!error <THRESH must be a scalar in the range \[0, 1\]> 584%! opts.thresh = []; 585%! ilu (A_tiny, opts); 586