Mercurial > octave-antonio

comparison scripts/sparse/ichol.m @ 19055:38937efbee21

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 431dc1da050c
comparison
equal deleted inserted replaced
19054:df64071e538c 19055:38937efbee21
1 ## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com>
1 ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com> 2 ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com>
2 ## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@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 under the 6 ## Octave is free software; you can redistribute it and/or modify it
7 ## terms of the GNU General Public License as published by the Free Software 7 ## under the terms of the GNU General Public License as published by
8 ## Foundation; either version 3 of the License, or (at your option) any later 8 ## the Free Software Foundation; either version 3 of the License, or (at
9 ## version. 9 ## your option) any later version.
10 ## 10 ##
11 ## Octave is distributed in the hope that it will be useful, but WITHOUT ANY 11 ## Octave is distributed in the hope that it will be useful, but
12 ## WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
13 ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 ## details. 14 ## General Public License for more details.
15 ## 15 ##
16 ## You should have received a copy of the GNU General Public License along with 16 ## You should have received a copy of the GNU General Public License
17 ## Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. 17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>.
18 19
19 ## -*- texinfo -*- 20 ## -*- texinfo -*-
20 ## @deftypefn {Function File} ichol (@var{A}, @var{opts}) 21 ## @deftypefn {Function File} {@var{L} =} ichol (@var{A})
21 ## @deftypefnx {Function File} {@var{L} =} ichol (@var{A}, @var{opts}) 22 ## @deftypefnx {Function File} {@var{L} =} ichol (@var{A}, @var{opts})
22 ## 23 ##
23 ## @code{@var{L} = ichol (@var{A})} performs the incomplete Cholesky 24 ## Compute the incomplete Cholesky factorization of the sparse square matrix
24 ## factorization of A with zero-fill. 25 ## @var{A} with zero-fill.
25 ## 26 ##
26 ## @code{@var{L} = ichol (@var{A}, @var{opts})} performs the incomplete Cholesky 27 ## By default, ichol references the lower triangle of @var{A} and produces
27 ## factorization of A with options specified by opts. 28 ## lower triangular factors.
28 ##
29 ## By default, ichol references the lower triangle of A and produces lower
30 ## triangular factors.
31 ## 29 ##
32 ## The factor given by this routine may be useful as a preconditioner for a 30 ## The factor given by this routine may be useful as a preconditioner for a
33 ## system of linear equations being solved by iterative methods such as 31 ## system of linear equations being solved by iterative methods such as
34 ## PCG (Preconditioned conjugate gradient). 32 ## PCG (Preconditioned Conjugate Gradient).
35 ## 33 ##
36 ## ichol works only for sparse square matrices. 34 ## The factorization may be modified by passing options in a structure
37 ## 35 ## @var{opts}. The option name is a field in the structure and the setting
38 ## The fields of @var{opts} must be named exactly as shown below. You can 36 ## is the value of field. Names and specifiers are case sensitive.
39 ## include any number of these fields in the structure and define them in any
40 ## order. Any additional fields are ignored. Names and specifiers are case
41 ## sensitive.
42 ## 37 ##
43 ## @table @asis 38 ## @table @asis
44 ## @item type 39 ## @item type
45 ## Type of factorization. 40 ## Type of factorization.
46 ## String indicating which flavor of incomplete Cholesky to perform. Valid 41 ## String indicating which flavor of incomplete Cholesky to perform. Valid
47 ## values of this field are @samp{nofill} and @samp{ict}. The 42 ## values of this field are @qcode{"nofill"} and @qcode{"ict"}. The
48 ## @samp{nofill} variant performs incomplete Cholesky with zero-fill [IC(0)]. 43 ## @qcode{"nofill"} variant performs incomplete Cholesky with zero-fill
49 ## The @samp{ict} variant performs incomplete Cholesky with threshold dropping 44 ## [IC(0)]. The @qcode{"ict"} variant performs incomplete Cholesky with
50 ## [ICT]. The default value is @samp{nofill}. 45 ## threshold dropping [ICT]. The default value is @qcode{"nofill"}.
51 ## 46 ##
52 ## @item droptol 47 ## @item droptol
53 ## Drop tolerance when type is @samp{ict}. 48 ## Drop tolerance when type is @qcode{"ict"}.
54 ## Nonnegative scalar used as a drop tolerance when performing ICT. Elements 49 ## Non-negative scalar used as a drop tolerance when performing ICT@. Elements
55 ## which are smaller in magnitude than a local drop tolerance are dropped from 50 ## which are smaller in magnitude than a local drop tolerance are dropped from
56 ## the resulting factor except for the diagonal element which is never dropped. 51 ## the resulting factor except for the diagonal element which is never dropped.
57 ## The local drop tolerance at step j of the factorization is 52 ## The local drop tolerance at step j of the factorization is
58 ## @code{norm (@var{A}(j:end, j), 1) * droptol}. @samp{droptol} is ignored if 53 ## @code{norm (@var{A}(j:end, j), 1) * droptol}. @code{droptol} is ignored if
59 ## @samp{type} is @samp{nofill}. The default value is 0. 54 ## @code{type} is @qcode{"nofill"}. The default value is 0.
60 ## 55 ##
61 ## @item michol 56 ## @item michol
62 ## Indicates whether to perform modified incomplete Cholesky. 57 ## Indicate whether modified incomplete Cholesky [MIC] is performed.
63 ## Indicates whether or not modified incomplete Cholesky [MIC] is performed. 58 ## The field may be @qcode{"on"} or @qcode{"off"}. When performing MIC, the
64 ## The field may be @samp{on} or @samp{off}. When performing MIC, the diagonal 59 ## diagonal is compensated for dropped elements to enforce the relationship
65 ## is compensated for dropped elements to enforce the relationship
66 ## @code{@var{A} * @var{e} = @var{L} * @var{L}' * @var{e}} where 60 ## @code{@var{A} * @var{e} = @var{L} * @var{L}' * @var{e}} where
67 ## @code{@var{e} = ones (size (@var{A}, 2), 1))}. The default value is 61 ## @code{@var{e} = ones (columns (@var{A}), 1)}. The default value is
68 ## @samp{off}. 62 ## @qcode{"off"}.
69 ## 63 ##
70 ## @item diagcomp 64 ## @item diagcomp
71 ## Perform compensated incomplete Cholesky with the specified coefficient. 65 ## Perform compensated incomplete Cholesky with the specified coefficient.
72 ## Real nonnegative scalar used as a global diagonal shift @code{@var{alpha}} 66 ## The coefficient is a real non-negative scalar used as a global diagonal
73 ## in forming the incomplete Cholesky factor. That is, instead of performing 67 ## shift @code{@var{alpha}} in forming the incomplete Cholesky factor. That
74 ## incomplete Cholesky on @code{@var{A}}, the factorization of 68 ## is, instead of performing incomplete Cholesky on @code{@var{A}}, the
75 ## @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is formed. The default 69 ## factorization of @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is
76 ## value is 0. 70 ## formed. The default value is 0.
77 ## 71 ##
78 ## @item shape 72 ## @item shape
79 ## Determines which triangle is referenced and returned. 73 ## Determine which triangle is referenced and returned.
80 ## Valid values are @samp{upper} and @samp{lower}. If @samp{upper} is specified, 74 ## Valid values are @qcode{"upper"} and @qcode{"lower"}. If @qcode{"upper"}
81 ## only the upper triangle of @code{@var{A}} is referenced and @code{@var{R}} 75 ## is specified, only the upper triangle of @code{@var{A}} is referenced and
82 ## is constructed such that @code{@var{A}} is approximated by 76 ## @code{@var{R}} is constructed such that @code{@var{A}} is approximated by
83 ## @code{@var{R}' * @var{R}}. If @samp{lower} is specified, only the lower 77 ## @code{@var{R}' * @var{R}}. If @qcode{"lower"} is specified, only the
84 ## triangle of @code{@var{A}} is referenced and @code{@var{L}} is constructed 78 ## lower triangle of @code{@var{A}} is referenced and @code{@var{L}} is
85 ## such that @code{@var{A}} is approximated by @code{@var{L} * @var{L}'}. The 79 ## constructed such that @code{@var{A}} is approximated by @code{@var{L} *
86 ## default value is @samp{lower}. 80 ## @var{L}'}. The default value is @qcode{"lower"}.
87 ## @end table 81 ## @end table
88 ## 82 ##
89 ## EXAMPLES 83 ## Examples
90 ## 84 ##
91 ## The following problem demonstrates how to factorize a sample symmetric 85 ## The following problem demonstrates how to factorize a sample symmetric
92 ## positive definite matrix with the full Cholesky decomposition and with the 86 ## positive definite matrix with the full Cholesky decomposition and with the
93 ## incomplete one. 87 ## incomplete one.
94 ## 88 ##
95 ## @example 89 ## @example
90 ## @group
96 ## A = [ 0.37, -0.05, -0.05, -0.07; 91 ## A = [ 0.37, -0.05, -0.05, -0.07;
97 ## -0.05, 0.116, 0.0, -0.05; 92 ## -0.05, 0.116, 0.0, -0.05;
98 ## -0.05, 0.0, 0.116, -0.05; 93 ## -0.05, 0.0, 0.116, -0.05;
99 ## -0.07, -0.05, -0.05, 0.202]; 94 ## -0.07, -0.05, -0.05, 0.202];
100 ## A = sparse(A); 95 ## A = sparse (A);
101 ## nnz(tril (A)) 96 ## nnz (tril (A))
102 ## ans = 9 97 ## ans = 9
103 ## L = chol(A, "lower"); 98 ## L = chol (A, "lower");
104 ## nnz (L) 99 ## nnz (L)
105 ## ans = 10 100 ## ans = 10
106 ## norm (A - L * L', "fro") / norm (A, "fro") 101 ## norm (A - L * L', "fro") / norm (A, "fro")
107 ## ans = 1.1993e-16 102 ## ans = 1.1993e-16
108 ## opts.type = 'nofill'; 103 ## opts.type = "nofill";
109 ## L = ichol(A,opts); 104 ## L = ichol (A, opts);
110 ## nnz (L) 105 ## nnz (L)
111 ## ans = 9 106 ## ans = 9
112 ## norm (A - L * L', "fro") / norm (A, "fro") 107 ## norm (A - L * L', "fro") / norm (A, "fro")
113 ## ans = 0.019736 108 ## ans = 0.019736
109 ## @end group
114 ## @end example 110 ## @end example
115 ## 111 ##
116 ## Another example for decomposition is finite difference matrix to solve a 112 ## Another example for decomposition is a finite difference matrix used to
117 ## boundary value problem on the unit square. 113 ## solve a boundary value problem on the unit square.
118 ## 114 ##
119 ## @example 115 ## @example
116 ## @group
120 ## nx = 400; ny = 200; 117 ## nx = 400; ny = 200;
121 ## hx = 1 / (nx + 1); hy = 1 / (ny + 1); 118 ## hx = 1 / (nx + 1); hy = 1 / (ny + 1);
122 ## Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2); 119 ## Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
123 ## Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2); 120 ## [-1 0 1 ], nx, nx) / (hx ^ 2);
121 ## Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
122 ## [-1 0 1 ], ny, ny) / (hy ^ 2);
124 ## A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy); 123 ## A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
125 ## nnz (tril (A)) 124 ## nnz (tril (A))
126 ## ans = 239400 125 ## ans = 239400
127 ## opts.type = 'nofill'; 126 ## opts.type = "nofill";
128 ## L = ichol (A, opts); 127 ## L = ichol (A, opts);
129 ## nnz (tril (A)) 128 ## nnz (tril (A))
130 ## ans = 239400 129 ## ans = 239400
131 ## norm (A - L * L', "fro") / norm (A, "fro") 130 ## norm (A - L * L', "fro") / norm (A, "fro")
132 ## ans = 0.062327 131 ## ans = 0.062327
132 ## @end group
133 ## @end example 133 ## @end example
134 ## 134 ##
135 ## References for the implemented algorithms: 135 ## References for implemented algorithms:
136 ## 136 ##
137 ## [1] Saad, Yousef. "Preconditioning Techniques." Iterative Methods for Sparse Linear 137 ## [1] @nospell{Y. Saad}. "Preconditioning Techniques." @cite{Iterative
138 ## Systems. PWS Publishing Company, 1996. 138 ## Methods for Sparse Linear Systems}, PWS Publishing Company, 1996.
139 ## 139 ##
140 ## [2] Jones, Mark T. and Plassmann, Paul E.: An Improved Incomplete Cholesky 140 ## [2] @nospell{M. Jones, P. Plassmann}: @cite{An Improved Incomplete
141 ## Factorization (1992). 141 ## Cholesky Factorization}, 1992.
142 ## @seealso{chol, ilu, pcg}
142 ## @end deftypefn 143 ## @end deftypefn
143 144
144 function [L] = ichol (A, opts) 145 function L = ichol (A, opts = struct ())
145 146
146 if ((nargin > 2) || (nargin < 1) || (nargout > 1)) 147 if (nargin < 1 || nargin > 2 || nargout > 1)
147 print_usage (); 148 print_usage ();
148 endif 149 endif
149 150
150 % Check input matrix 151 if (! (issparse (A) && issquare (A)))
151 if (~issparse(A) || ~issquare (A)) 152 error ("ichol: A must be a sparse square matrix");
152 error ("ichol: Input A must be a non-empty sparse square matrix"); 153 endif
153 endif 154
154 155 if (! isstruct (opts))
155 % If A is empty and sparse then return empty L 156 error ("ichol: OPTS must be a structure.");
156 % Compatibility with Matlab 157 endif
158
159 ## If A is empty then return empty L for Matlab compatibility
157 if (isempty (A)) 160 if (isempty (A))
158 L = A; 161 L = A;
159 return; 162 return;
160 endif 163 endif
161 164
162 % Check input structure, otherwise set default values 165 ## Parse input options
163 if (nargin == 2) 166 if (! isfield (opts, "type"))
164 if (~isstruct (opts)) 167 opts.type = "nofill"; # set default
165 error ("ichol: Input \"opts\" must be a valid structure."); 168 else
169 type = tolower (getfield (opts, "type"));
170 if (! strcmp (type, "nofill") && ! strcmp (type, "ict"))
171 error ('ichol: TYPE must be "nofill" or "ict"');
166 endif 172 endif
167 else 173 opts.type = type;
168 opts = struct (); 174 endif
169 endif 175
170 176 if (! isfield (opts, "droptol"))
171 if (~isfield (opts, "type")) 177 opts.droptol = 0; # set default
172 opts.type = "nofill"; % set default 178 else
173 else 179 if (! (isreal (opts.droptol) && isscalar (opts.droptol)
174 type = tolower (getfield (opts, "type")); 180 && opts.droptol >= 0))
175 if ((strcmp (type, "nofill") == 0) 181 error ("ichol: DROPTOL must be a non-negative real scalar");
176 && (strcmp (type, "ict") == 0))
177 error ("ichol: Invalid field \"type\" in input structure.");
178 else
179 opts.type = type;
180 endif 182 endif
181 endif 183 endif
182 184
183 if (~isfield (opts, "droptol")) 185 michol = "";
184 opts.droptol = 0; % set default 186 if (! isfield (opts, "michol"))
185 else 187 opts.michol = "off"; # set default
186 if (~isscalar (opts.droptol) || (opts.droptol < 0)) 188 else
187 error ("ichol: Invalid field \"droptol\" in input structure."); 189 michol = tolower (getfield (opts, "michol"));
190 if (! strcmp (michol, "off") && ! strcmp (michol, "on"))
191 error ('ichol: MICHOL must be "on" or "off"');
188 endif 192 endif
189 endif 193 opts.michol = michol;
190 194 endif
191 michol = ""; 195
192 if (~isfield (opts, "michol")) 196 if (! isfield (opts, "diagcomp"))
193 opts.michol = "off"; % set default 197 opts.diagcomp = 0; # set default
194 else 198 else
195 michol = tolower (getfield (opts, "michol")); 199 if (! (isreal (opts.diagcomp) && isscalar (opts.diagcomp)
196 if ((strcmp (michol, "off") == 0) 200 && opts.diagcomp >= 0))
197 && (strcmp (michol, "on") == 0)) 201 error ("ichol: DIAGCOMP must be a non-negative real scalar");
198 error ("ichol: Invalid field \"michol\" in input structure.");
199 else
200 opts.michol = michol;
201 endif 202 endif
202 endif 203 endif
203 204
204 if (~isfield (opts, "diagcomp")) 205 if (! isfield (opts, "shape"))
205 opts.diagcomp = 0; % set default 206 opts.shape = "lower"; # set default
206 else 207 else
207 if (~isscalar (opts.diagcomp) || (opts.diagcomp < 0)) 208 shape = tolower (getfield (opts, "shape"));
208 error ("ichol: Invalid field \"diagcomp\" in input structure."); 209 if (! strcmp (shape, "lower") && ! strcmp (shape, "upper"))
210 error ('ichol: SHAPE must be "lower" or "upper"');
209 endif 211 endif
210 endif 212 opts.shape = shape;
211 213 endif
212 if (~isfield (opts, "shape")) 214
213 opts.shape = "lower"; % set default 215 ## Prepare input for specialized ICHOL
214 else
215 shape = tolower (getfield (opts, "shape"));
216 if ((strcmp (shape, "lower") == 0)
217 && (strcmp (shape, "upper") == 0))
218 error ("ichol: Invalid field \"shape\" in input structure.");
219 else
220 opts.shape = shape;
221 endif
222 endif
223
224 % Prepare input for specialized ICHOL
225 A_in = []; 216 A_in = [];
226 if (opts.diagcomp > 0) 217 if (opts.diagcomp > 0)
227 A += opts.diagcomp * diag (diag (A)); 218 A += opts.diagcomp * diag (diag (A));
228 endif 219 endif
229 if (strcmp (opts.shape, "upper") == 1) 220 if (strcmp (opts.shape, "upper"))
230 A_in = triu (A); 221 A_in = triu (A);
231 A_in = A_in'; 222 A_in = A_in';
232 else 223 else
233 A_in = tril (A); 224 A_in = tril (A);
234 endif 225 endif
235 226
236 % Delegate to specialized ICHOL 227 ## Delegate to specialized ICHOL
237 switch (opts.type) 228 switch (opts.type)
238 case "nofill" 229 case "nofill"
239 L = __ichol0__ (A_in, opts.michol); 230 L = __ichol0__ (A_in, opts.michol);
240 case "ict" 231 case "ict"
241 L = __icholt__ (A_in, opts.droptol, opts.michol); 232 L = __icholt__ (A_in, opts.droptol, opts.michol);
242 otherwise
243 printf ("The input structure is invalid.\n");
244 endswitch 233 endswitch
245 234
246 if (strcmp (opts.shape, "upper") == 1) 235 if (strcmp (opts.shape, "upper"))
247 L = L'; 236 L = L';
248 endif 237 endif
249 238
250
251 endfunction 239 endfunction
240
252 241
253 %!shared A1, A2, A3, A4, A5, A6, A7 242 %!shared A1, A2, A3, A4, A5, A6, A7
254 %! A1 = [ 0.37, -0.05, -0.05, -0.07; 243 %! A1 = [ 0.37, -0.05, -0.05, -0.07;
255 %! -0.05, 0.116, 0.0, -0.05; 244 %! -0.05, 0.116, 0.0, -0.05;
256 %! -0.05, 0.0, 0.116, -0.05; 245 %! -0.05, 0.0, 0.116, -0.05;
257 %! -0.07, -0.05, -0.05, 0.202]; 246 %! -0.07, -0.05, -0.05, 0.202];
258 %! A1 = sparse(A1); 247 %! A1 = sparse (A1);
259 %! A2 = gallery ('poisson', 30); 248 %! A2 = gallery ("poisson", 30);
260 %! A3 = gallery ('tridiag', 50); 249 %! A3 = gallery ("tridiag", 50);
261 %! nx = 400; ny = 200; 250 %! nx = 400; ny = 200;
262 %! hx = 1 / (nx + 1); hy = 1 / (ny + 1); 251 %! hx = 1 / (nx + 1); hy = 1 / (ny + 1);
263 %! Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2); 252 %! Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
264 %! Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2); 253 %! [-1 0 1 ], nx, nx) / (hx ^ 2);
254 %! Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
255 %! [-1 0 1 ], ny, ny) / (hy ^ 2);
265 %! A4 = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy); 256 %! A4 = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
266 %! A5 = [ 0.37, -0.05, -0.05, -0.07; 257 %! A5 = [ 0.37, -0.05, -0.05, -0.07;
267 %! -0.05, 0.116, 0.0, -0.05 + 0.05i; 258 %! -0.05, 0.116, 0.0, -0.05 + 0.05i;
268 %! -0.05, 0.0, 0.116, -0.05; 259 %! -0.05, 0.0, 0.116, -0.05;
269 %! -0.07, -0.05 - 0.05i, -0.05, 0.202]; 260 %! -0.07, -0.05 - 0.05i, -0.05, 0.202];
270 %! A5 = sparse(A5); 261 %! A5 = sparse (A5);
271 %! A6 = [ 0.37, -0.05 - i, -0.05, -0.07; 262 %! A6 = [ 0.37, -0.05 - i, -0.05, -0.07;
272 %! -0.05 + i, 0.116, 0.0, -0.05; 263 %! -0.05 + i, 0.116, 0.0, -0.05;
273 %! -0.05, 0.0, 0.116, -0.05; 264 %! -0.05, 0.0, 0.116, -0.05;
274 %! -0.07, -0.05, -0.05, 0.202]; 265 %! -0.07, -0.05, -0.05, 0.202];
275 %! A6 = sparse(A6); 266 %! A6 = sparse (A6);
276 %! A7 = A5; 267 %! A7 = A5;
277 %! A7(1) = 2i; 268 %! A7(1) = 2i;
278 %! 269
279 270 ## ICHOL0 tests
280 %!# Input validation tests
281
282 %!test
283 %!error ichol ([]);
284 %!error ichol (0);
285 %!error ichol (-0);
286 %!error ichol (1);
287 %!error ichol (-1);
288 %!error ichol (i);
289 %!error ichol (-i);
290 %!error ichol (1 + 1i);
291 %!error ichol (1 - 1i);
292 %!error ichol (sparse (0));
293 %!error ichol (sparse (-0));
294 %!error ichol (sparse (-1));
295 %!error ichol (sparse (-1));
296 %!test
297 %! opts.milu = 'foo';
298 %!error L = ichol (A1, opts);
299 %! opts.milu = 1;
300 %!error L = ichol (A1, opts);
301 %! opts.milu = [];
302 %!error L = ichol (A1, opts);
303 %!test
304 %! opts.droptol = -1;
305 %!error L = ichol (A1, opts);
306 %! opts.droptol = 0.5i;
307 %!error L = ichol (A1, opts);
308 %! opts.droptol = [];
309 %!error L = ichol (A1, opts);
310 %!test
311 %! opts.type = 'foo';
312 %!error L = ichol (A1, opts);
313 %! opts.type = 1;
314 %!error L = ichol (A1, opts);
315 %! opts.type = [];
316 %!error L = ichol (A1, opts);
317 %!test
318 %! opts.shape = 'foo';
319 %!error L = ichol (A1, opts);
320 %! opts.shape = 1;
321 %!error L = ichol (A1, opts);
322 %! opts.shape = [];
323 %!error L = ichol (A1, opts);
324 %!test
325 %! opts.diagcomp = -1;
326 %!error L = ichol (A1, opts);
327 %! opts.diagcomp = 0.5i;
328 %!error L = ichol (A1, opts);
329 %! opts.diagcomp = [];
330 %!error L = ichol (A1, opts);
331
332 %!# ICHOL0 tests
333 271
334 %!test 272 %!test
335 %! opts.type = "nofill"; 273 %! opts.type = "nofill";
336 %! opts.michol = "off"; 274 %! opts.michol = "off";
337 %! assert (nnz (tril (A1)), nnz (ichol (A1, opts))); 275 %! assert (nnz (tril (A1)), nnz (ichol (A1, opts)));
399 %! L = ichol (A5, opts); 337 %! L = ichol (A5, opts);
400 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0195, 1e-4); 338 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0195, 1e-4);
401 %! opts.michol = "on"; 339 %! opts.michol = "on";
402 %! L = ichol (A5, opts); 340 %! L = ichol (A5, opts);
403 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4); 341 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4);
404 %!test 342
405 %% Negative pivot 343 ## Negative pivot
406 %!error ichol (A6); 344 %!error <negative pivot> ichol (A6)
407 %% Complex entry in the diagonal 345 %!error ichol (A6)
408 %!error ichol (A7); 346 ## Complex entry in the diagonal
409 347 %!error <non-real pivot> ichol (A7)
410 %%!ICHOLT tests 348
349 ## ICHOLT tests
411 350
412 %%!test 351 %%!test
413 %! opts.type = "ict"; 352 %! opts.type = "ict";
414 %! opts.droptol = 1e-1; 353 %! opts.droptol = 1e-1;
415 %! opts.michol = "off"; 354 %! opts.michol = "off";
473 %! L = ichol (A5, opts); 412 %! L = ichol (A5, opts);
474 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.2044, 1e-4); 413 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.2044, 1e-4);
475 %! opts.michol = "on"; 414 %! opts.michol = "on";
476 %! L = ichol (A5, opts); 415 %! L = ichol (A5, opts);
477 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.3231, 1e-4); 416 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.3231, 1e-4);
478 %!test 417
479 %% Negative pivot 418 %% Input validation tests
480 %! opts.type = "ict"; 419
481 %!error ichol (A6, setup); 420 %!error <A must be a sparse square matrix> ichol ([])
482 %% Complex entry in the diagonal 421 %!error <A must be a sparse square matrix> ichol (0)
483 %!error ichol (A7, setup); 422 %!error <pivot equal to 0> ichol (sparse (0))
423 %!error <pivot equal to 0> ichol (sparse (-0))
424 %!error <negative pivot> ichol (sparse (-1))
425 %!test
426 %! opts.type = "foo";
427 %! fail ("ichol (A1, opts)", 'TYPE must be "nofill"');
428 %! opts.type = 1;
429 %! fail ("ichol (A1, opts)", 'TYPE must be "nofill"');
430 %! opts.type = [];
431 %! fail ("ichol (A1, opts)", 'TYPE must be "nofill"');
432 %!test
433 %! opts.droptol = -1;
434 %! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar");
435 %! opts.droptol = 0.5i;
436 %! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar");
437 %! opts.droptol = [];
438 %! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar");
439 %!test
440 %! opts.michol = "foo";
441 %! fail ("ichol (A1, opts)", 'MICHOL must be "on"');
442 %! opts.michol = 1;
443 %! fail ("ichol (A1, opts)", 'MICHOL must be "on"');
444 %! opts.michol = [];
445 %! fail ("ichol (A1, opts)", 'MICHOL must be "on"');
446 %!test
447 %! opts.diagcomp = -1;
448 %! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar");
449 %! opts.diagcomp = 0.5i;
450 %! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar");
451 %! opts.diagcomp = [];
452 %! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar");
453 %!test
454 %! opts.shape = "foo";
455 %! fail ("ichol (A1, opts)", 'SHAPE must be "lower"');
456 %! opts.shape = 1;
457 %! fail ("ichol (A1, opts)", 'SHAPE must be "lower"');
458 %! opts.shape = [];
459 %! fail ("ichol (A1, opts)", 'SHAPE must be "lower"');
460