Mercurial > octave-nkf

comparison scripts/sparse/ichol.m @ 19329:76a6ba7d65d0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
doc: Update documentation for ilu, ichol. * ichol.m, ilu.m: Clarify options available for factorization.
author Nir Krakauer <nkrakauer@ccny.cuny.edu> and Rik <rik@octave.org>
date Sat, 01 Nov 2014 17:24:23 -0700
parents 431dc1da050c
children cbce5d1bcaf9
comparison
equal deleted inserted replaced
19328:b80b396e7d54 19329:76a6ba7d65d0
1 ## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com> 1 ## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com>
2 ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> 2 ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com>
3 ## 3 ##
4 ## This file is part of Octave. 4 ## This file is part of Octave.
5 ## 5 ##
6 ## Octave is free software; you can redistribute it and/or modify it 6 ## Octave is free software; you can redistribute it and/or modify it
7 ## under the terms of the GNU General Public License as published by 7 ## under the terms of the GNU General Public License as published by
8 ## the Free Software Foundation; either version 3 of the License, or (at 8 ## the Free Software Foundation; either version 3 of the License, or (at
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {@var{L} =} ichol (@var{A}) 21 ## @deftypefn {Function File} {@var{L} =} ichol (@var{A})
22 ## @deftypefnx {Function File} {@var{L} =} ichol (@var{A}, @var{opts}) 22 ## @deftypefnx {Function File} {@var{L} =} ichol (@var{A}, @var{opts})
23 ## 23 ##
24 ## Compute the incomplete Cholesky factorization of the sparse square matrix 24 ## Compute the incomplete Cholesky factorization of the sparse square matrix
25 ## @var{A} with zero-fill. 25 ## @var{A}.
26 ## 26 ##
27 ## By default, ichol references the lower triangle of @var{A} and produces 27 ## By default, @code{ichol} uses only the lower triangle of @var{A} and
28 ## lower triangular factors. 28 ## produces a lower triangular factor @var{L} such that @tcode{L*L'}
29 ## approximates @var{A}.
29 ## 30 ##
30 ## The factor given by this routine may be useful as a preconditioner for a 31 ## The factor given by this routine may be useful as a preconditioner for a
31 ## system of linear equations being solved by iterative methods such as 32 ## system of linear equations being solved by iterative methods such as
32 ## PCG (Preconditioned Conjugate Gradient). 33 ## PCG (Preconditioned Conjugate Gradient).
33 ## 34 ##
34 ## The factorization may be modified by passing options in a structure 35 ## The factorization may be modified by passing options in a structure
35 ## @var{opts}. The option name is a field in the structure and the setting 36 ## @var{opts}. The option name is a field of the structure and the setting
36 ## is the value of field. Names and specifiers are case sensitive. 37 ## is the value of field. Names and specifiers are case sensitive.
37 ## 38 ##
38 ## @table @asis 39 ## @table @asis
39 ## @item type 40 ## @item type
40 ## Type of factorization. 41 ## Type of factorization.
41 ## String indicating which flavor of incomplete Cholesky to perform. Valid 42 ##
42 ## values of this field are @qcode{"nofill"} and @qcode{"ict"}. The 43 ## @table @asis
43 ## @qcode{"nofill"} variant performs incomplete Cholesky with zero-fill 44 ## @item @qcode{"nofill"} (default)
44 ## @nospell{[IC(0)]}. The @qcode{"ict"} variant performs incomplete Cholesky 45 ## Incomplete Cholesky factorization with no fill-in (@nospell{IC(0)}).
45 ## with threshold dropping @nospell{[ICT]}. The default value is 46 ##
46 ## @qcode{"nofill"}. 47 ## @item @qcode{"ict"}
48 ## Incomplete Cholesky factorization with threshold dropping (@nospell{ICT}).
49 ## @end table
50 ##
51 ## @item diagcomp
52 ## A non-negative scalar @var{alpha} for incomplete Cholesky factorization of
53 ## @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} instead of @var{A}.
54 ## This can be useful when @var{A} is not positive definite. The default value
55 ## is 0.
47 ## 56 ##
48 ## @item droptol 57 ## @item droptol
49 ## Drop tolerance when type is @qcode{"ict"}. 58 ## A non-negative scalar specifying the drop tolerance for factorization if
50 ## Non-negative scalar used as a drop tolerance when performing @nospell{ICT}@. 59 ## performing @nospell{ICT}@. The default value is 0 which produces the complete
51 ## Elements which are smaller in magnitude than a local drop tolerance are 60 ## Cholesky factorization.
52 ## dropped from the resulting factor except for the diagonal element which is 61 ##
53 ## never dropped. 62 ## Non-diagonal entries of @var{L} are set to 0 unless
54 ## The local drop tolerance at step j of the factorization is 63 ##
55 ## @code{norm (@var{A}(j:end, j), 1) * droptol}. @code{droptol} is ignored if 64 ## @code{abs (@var{L}(i,j)) >= droptol * norm (@var{A}(j:end, j), 1)}.
56 ## @code{type} is @qcode{"nofill"}. The default value is 0.
57 ## 65 ##
58 ## @item michol 66 ## @item michol
59 ## Indicate whether modified incomplete Cholesky [MIC] is performed. 67 ## Modified incomplete Cholesky factorization:
60 ## The field may be @qcode{"on"} or @qcode{"off"}. When performing MIC, the 68 ##
61 ## diagonal is compensated for dropped elements to enforce the relationship 69 ## @table @asis
62 ## @code{@var{A} * @var{e} = @var{L} * @var{L}' * @var{e}} where 70 ## @item @qcode{"off"} (default)
63 ## @code{@var{e} = ones (columns (@var{A}), 1)}. The default value is 71 ## Row and column sums are not necessarily preserved.
64 ## @qcode{"off"}. 72 ##
65 ## 73 ## @item @qcode{"on"}
66 ## @item diagcomp 74 ## The diagonal of @var{L} is modified so that row (and column) sums are
67 ## Perform compensated incomplete Cholesky with the specified coefficient. 75 ## preserved even when elements have been dropped during the factorization.
68 ## The coefficient is a real non-negative scalar used as a global diagonal 76 ## The relationship preserved is: @code{@var{A} * e = @var{L} * @var{L}' * e},
69 ## shift @code{@var{alpha}} in forming the incomplete Cholesky factor. That 77 ## where e is a vector of ones.
70 ## is, instead of performing incomplete Cholesky on @code{@var{A}}, the 78 ## @end table
71 ## factorization of @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is
72 ## formed. The default value is 0.
73 ## 79 ##
74 ## @item shape 80 ## @item shape
75 ## Determine which triangle is referenced and returned. 81 ##
76 ## Valid values are @qcode{"upper"} and @qcode{"lower"}. If @qcode{"upper"} 82 ## @table @asis
77 ## is specified, only the upper triangle of @code{@var{A}} is referenced and 83 ## @item @qcode{"lower"} (default)
78 ## @code{@var{R}} is constructed such that @code{@var{A}} is approximated by 84 ## Use only the lower triangle of @var{A} and return a lower triangular
79 ## @code{@var{R}' * @var{R}}. If @qcode{"lower"} is specified, only the 85 ## factor @var{L} such that @tcode{L*L'} approximates @var{A}.
80 ## lower triangle of @code{@var{A}} is referenced and @code{@var{L}} is 86 ##
81 ## constructed such that @code{@var{A}} is approximated by @code{@var{L} * 87 ## @item @qcode{"upper"}
82 ## @var{L}'}. The default value is @qcode{"lower"}. 88 ## Use only the upper triangle of @var{A} and return an upper triangular
89 ## factor @var{U} such that @code{U'*U} approximates @var{A}.
83 ## @end table 90 ## @end table
84 ## 91 ## @end table
85 ## Examples 92 ##
93 ## EXAMPLES
86 ## 94 ##
87 ## The following problem demonstrates how to factorize a sample symmetric 95 ## The following problem demonstrates how to factorize a sample symmetric
88 ## positive definite matrix with the full Cholesky decomposition and with the 96 ## positive definite matrix with the full Cholesky decomposition and with the
89 ## incomplete one. 97 ## incomplete one.
90 ## 98 ##
157 if (! isstruct (opts)) 165 if (! isstruct (opts))
158 error ("ichol: OPTS must be a structure."); 166 error ("ichol: OPTS must be a structure.");
159 endif 167 endif
160 168
161 ## If A is empty then return empty L for Matlab compatibility 169 ## If A is empty then return empty L for Matlab compatibility
162 if (isempty (A)) 170 if (isempty (A))
163 L = A; 171 L = A;
164 return; 172 return;
165 endif 173 endif
166 174
167 ## Parse input options 175 ## Parse input options
235 endswitch 243 endswitch
236 244
237 if (strcmp (opts.shape, "upper")) 245 if (strcmp (opts.shape, "upper"))
238 L = L'; 246 L = L';
239 endif 247 endif
240 248
241 endfunction 249 endfunction
242 250
243 251
244 %!shared A1, A2, A3, A4, A5, A6, A7 252 %!shared A1, A2, A3, A4, A5, A6, A7
245 %! A1 = [ 0.37, -0.05, -0.05, -0.07; 253 %! A1 = [ 0.37, -0.05, -0.05, -0.07;
340 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0195, 1e-4); 348 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0195, 1e-4);
341 %! opts.michol = "on"; 349 %! opts.michol = "on";
342 %! L = ichol (A5, opts); 350 %! L = ichol (A5, opts);
343 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4); 351 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4);
344 352
345 ## Negative pivot 353 ## Negative pivot
346 %!error <negative pivot> ichol (A6) 354 %!error <negative pivot> ichol (A6)
347 %!error ichol (A6) 355 %!error ichol (A6)
348 ## Complex entry in the diagonal 356 ## Complex entry in the diagonal
349 %!error <non-real pivot> ichol (A7) 357 %!error <non-real pivot> ichol (A7)
350 358
351 ## ICHOLT tests 359 ## ICHOLT tests
352 360
353 %%!test 361 %%!test
354 %! opts.type = "ict"; 362 %! opts.type = "ict";
355 %! opts.droptol = 1e-1; 363 %! opts.droptol = 1e-1;
356 %! opts.michol = "off"; 364 %! opts.michol = "off";
357 %! L = ichol (A1, opts); 365 %! L = ichol (A1, opts);