Mercurial > hg > octave-nkf

comparison scripts/sparse/ichol.m @ 19246:df64071e538c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Removed ichol0.cc, icholt.cc, ilu0.cc, iluc.cc, ilutp.cc. Created __ichol__.cc and __ilu__.cc. Minor bugs fixed. * libinterp/dldfcn/__ichol__.cc: File created now contains __ichol0__ and __icholt__ functions. * libinterp/dldfcn/__ilu__.cc: File created now contains __ilu0__ __iluc__ and __ilutp__ functions. * scripts/sparse/ichol.m: Tests have been moved from .cc files to this one. * changed scripts/sparse/ilu.m: Tests have been moved from .cc files to this one.
author Eduardo Ramos (edu159) <eduradical951@gmail.com>
date Mon, 18 Aug 2014 12:32:16 +0100
parents 168c0aa9bb05
children 38937efbee21
comparison
equal deleted inserted replaced
19245:168c0aa9bb05 19246:df64071e538c
101 ## nnz(tril (A)) 101 ## nnz(tril (A))
102 ## ans = 9 102 ## ans = 9
103 ## L = chol(A, "lower"); 103 ## L = chol(A, "lower");
104 ## nnz (L) 104 ## nnz (L)
105 ## ans = 10 105 ## ans = 10
106 ## norm (A - L * L', 'fro') / norm (A, 'fro') 106 ## norm (A - L * L', "fro") / norm (A, "fro")
107 ## ans = 1.1993e-16 107 ## ans = 1.1993e-16
108 ## opts.type = 'nofill'; 108 ## opts.type = 'nofill';
109 ## L = ichol(A,opts); 109 ## L = ichol(A,opts);
110 ## nnz (L) 110 ## nnz (L)
111 ## ans = 9 111 ## ans = 9
112 ## norm (A - L * L', 'fro') / norm (A, 'fro') 112 ## norm (A - L * L', "fro") / norm (A, "fro")
113 ## ans = 0.019736 113 ## ans = 0.019736
114 ## @end example 114 ## @end example
115 ## 115 ##
116 ## Another example for decomposition is finite difference matrix to solve a 116 ## Another example for decomposition is finite difference matrix to solve a
117 ## boundary value problem on the unit square. 117 ## boundary value problem on the unit square.
126 ## ans = 239400 126 ## ans = 239400
127 ## opts.type = 'nofill'; 127 ## opts.type = 'nofill';
128 ## L = ichol (A, opts); 128 ## L = ichol (A, opts);
129 ## nnz (tril (A)) 129 ## nnz (tril (A))
130 ## ans = 239400 130 ## ans = 239400
131 ## norm (A - L * L', 'fro') / norm (A, 'fro') 131 ## norm (A - L * L', "fro") / norm (A, "fro")
132 ## ans = 0.062327 132 ## ans = 0.062327
133 ## @end example 133 ## @end example
134 ## 134 ##
135 ## References for the implemented algorithms: 135 ## References for the implemented algorithms:
136 ## 136 ##
146 if ((nargin > 2) || (nargin < 1) || (nargout > 1)) 146 if ((nargin > 2) || (nargin < 1) || (nargout > 1))
147 print_usage (); 147 print_usage ();
148 endif 148 endif
149 149
150 % Check input matrix 150 % Check input matrix
151 if (isempty (A) || ~issparse(A) || ~issquare (A)) 151 if (~issparse(A) || ~issquare (A))
152 error ("ichol: Input A must be a non-empty sparse square matrix"); 152 error ("ichol: Input A must be a non-empty sparse square matrix");
153 endif
154
155 % If A is empty and sparse then return empty L
156 % Compatibility with Matlab
157 if (isempty (A))
158 L = A;
159 return;
153 endif 160 endif
154 161
155 % Check input structure, otherwise set default values 162 % Check input structure, otherwise set default values
156 if (nargin == 2) 163 if (nargin == 2)
157 if (~isstruct (opts)) 164 if (~isstruct (opts))
218 A_in = []; 225 A_in = [];
219 if (opts.diagcomp > 0) 226 if (opts.diagcomp > 0)
220 A += opts.diagcomp * diag (diag (A)); 227 A += opts.diagcomp * diag (diag (A));
221 endif 228 endif
222 if (strcmp (opts.shape, "upper") == 1) 229 if (strcmp (opts.shape, "upper") == 1)
223 disp("entro");
224 A_in = triu (A); 230 A_in = triu (A);
225 A_in = A_in'; 231 A_in = A_in';
226
227 else 232 else
228 A_in = tril (A); 233 A_in = tril (A);
229 endif 234 endif
230 235
231 % Delegate to specialized ICHOL 236 % Delegate to specialized ICHOL
232 switch (opts.type) 237 switch (opts.type)
233 case "nofill" 238 case "nofill"
234 L = ichol0 (A_in, opts.michol); 239 L = __ichol0__ (A_in, opts.michol);
235 case "ict" 240 case "ict"
236 L = icholt (A_in, opts.droptol, opts.michol); 241 L = __icholt__ (A_in, opts.droptol, opts.michol);
237 otherwise 242 otherwise
238 printf ("The input structure is invalid.\n"); 243 printf ("The input structure is invalid.\n");
239 endswitch 244 endswitch
240 245
241 if (strcmp (opts.shape, "upper") == 1) 246 if (strcmp (opts.shape, "upper") == 1)
243 endif 248 endif
244 249
245 250
246 endfunction 251 endfunction
247 252
248 %!shared A1, A2 253 %!shared A1, A2, A3, A4, A5, A6, A7
249 %! A1 = [ 0.37, -0.05, -0.05, -0.07; 254 %! A1 = [ 0.37, -0.05, -0.05, -0.07;
250 %! -0.05, 0.116, 0.0, -0.05; 255 %! -0.05, 0.116, 0.0, -0.05;
251 %! -0.05, 0.0, 0.116, -0.05; 256 %! -0.05, 0.0, 0.116, -0.05;
252 %! -0.07, -0.05, -0.05, 0.202]; 257 %! -0.07, -0.05, -0.05, 0.202];
253 %! A1 = sparse(A1); 258 %! A1 = sparse(A1);
259 %! A2 = gallery ('poisson', 30);
260 %! A3 = gallery ('tridiag', 50);
254 %! nx = 400; ny = 200; 261 %! nx = 400; ny = 200;
255 %! hx = 1 / (nx + 1); hy = 1 / (ny + 1); 262 %! hx = 1 / (nx + 1); hy = 1 / (ny + 1);
256 %! Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2); 263 %! Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2);
257 %! Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2); 264 %! Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2);
258 %! A2 = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy); 265 %! A4 = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
259 %! 266 %! A5 = [ 0.37, -0.05, -0.05, -0.07;
267 %! -0.05, 0.116, 0.0, -0.05 + 0.05i;
268 %! -0.05, 0.0, 0.116, -0.05;
269 %! -0.07, -0.05 - 0.05i, -0.05, 0.202];
270 %! A5 = sparse(A5);
271 %! A6 = [ 0.37, -0.05 - i, -0.05, -0.07;
272 %! -0.05 + i, 0.116, 0.0, -0.05;
273 %! -0.05, 0.0, 0.116, -0.05;
274 %! -0.07, -0.05, -0.05, 0.202];
275 %! A6 = sparse(A6);
276 %! A7 = A5;
277 %! A7(1) = 2i;
278 %!
279
280 %!# Input validation tests
281
260 %!test 282 %!test
261 %!error ichol ([]); 283 %!error ichol ([]);
262 %!error ichol (0); 284 %!error ichol (0);
263 %!error ichol (-0); 285 %!error ichol (-0);
264 %!error ichol (1); 286 %!error ichol (1);
269 %!error ichol (1 - 1i); 291 %!error ichol (1 - 1i);
270 %!error ichol (sparse (0)); 292 %!error ichol (sparse (0));
271 %!error ichol (sparse (-0)); 293 %!error ichol (sparse (-0));
272 %!error ichol (sparse (-1)); 294 %!error ichol (sparse (-1));
273 %!error ichol (sparse (-1)); 295 %!error ichol (sparse (-1));
274 %! 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
275 %!test 334 %!test
276 %! opts.type = "nofill"; 335 %! opts.type = "nofill";
277 %! opts.michol = "off"; 336 %! opts.michol = "off";
278 %! assert (nnz (tril (A1)), nnz (ichol (A1, opts))); 337 %! assert (nnz (tril (A1)), nnz (ichol (A1, opts)));
279 %! assert (nnz (tril (A2)), nnz (ichol (A2, opts))); 338 %! assert (nnz (tril (A2)), nnz (ichol (A2, opts)));
280 %! 339 %! assert (nnz (tril (A3)), nnz (ichol (A3, opts)));
281 %!test 340 %! assert (nnz (tril (A4)), nnz (ichol (A4, opts)));
282 %! opts.type = "nofill"; 341 %! assert (nnz (tril (A5)), nnz (ichol (A5, opts)));
283 %! opts.michol = "off"; 342 %!
284 %! L = ichol (A1, opts); 343 %!test
285 %! assert (norm (A1 - L * L', 'fro') / norm (A1, 'fro'), 0.01, 0.01); 344 %! opts.type = "nofill";
345 %! opts.michol = "off";
346 %! L = ichol (A1, opts);
347 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.0197, 1e-4);
348 %! opts.shape = "upper";
349 %! U = ichol (A1, opts);
350 %! assert (norm (A1 - U' * U, "fro") / norm (A1, "fro"), 0.0197, 1e-4);
351 %! opts.shape = "lower";
352 %! L = ichol (A1, opts);
353 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.0197, 1e-4);
354 %!
355 %!test
356 %! opts.michol = "on";
357 %! opts.shape = "lower";
358 %! opts.type = "nofill";
359 %! L = ichol (A1, opts);
360 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.0279, 1e-4);
361 %! opts.shape = "upper";
362 %! U = ichol (A1, opts);
363 %! assert (norm (A1 - U' * U, "fro") / norm (A1, "fro"), 0.0279, 1e-4);
364 %! opts.shape = "lower";
365 %! opts.diagcomp = 3e-3;
366 %! L = ichol (A1, opts);
367 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.0272, 1e-4);
368 %!
369 %!test
370 %! opts.type = "nofill";
371 %! opts.michol = "off";
286 %! L = ichol (A2, opts); 372 %! L = ichol (A2, opts);
287 %! assert (norm (A2 - L * L', 'fro') / norm (A2, 'fro'), 0.06, 0.01); 373 %! assert (norm (A2 - L*L', "fro") / norm (A2, "fro"), 0.0893, 1e-4)
288 %! 374 %! opts.michol = "on";
375 %! L = ichol (A2, opts);
376 %! assert (norm (A2 - L*L', "fro") / norm (A2, "fro"), 0.2377, 1e-4)
377 %!
378 %!test
379 %! opts.type = "nofill";
380 %! opts.michol = "off";
381 %! L = ichol (A3, opts);
382 %! assert (norm (A3 - L*L', "fro") / norm (A3, "fro"), eps, eps);
383 %! opts.michol = "on";
384 %! L = ichol (A3, opts);
385 %! assert (norm (A3 - L*L', "fro") / norm (A3, "fro"), eps, eps);
386 %!
387 %!test
388 %! opts.type = "nofill";
389 %! opts.michol = "off";
390 %! L = ichol (A4, opts);
391 %! assert (norm (A4 - L*L', "fro") / norm (A4, "fro"), 0.0623, 1e-4);
392 %! opts.michol = "on";
393 %! L = ichol (A4, opts);
394 %! assert (norm (A4 - L*L', "fro") / norm (A4, "fro"), 0.1664, 1e-4);
395 %!
396 %!test
397 %! opts.type = "nofill";
398 %! opts.michol = "off";
399 %! L = ichol (A5, opts);
400 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0195, 1e-4);
401 %! opts.michol = "on";
402 %! L = ichol (A5, opts);
403 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4);
404 %!test
405 %% Negative pivot
406 %!error ichol (A6);
407 %% Complex entry in the diagonal
408 %!error ichol (A7);
409
410 %%!ICHOLT tests
411
289 %%!test 412 %%!test
290 %%! opts.type = "nofill"; 413 %! opts.type = "ict";
291 %%! opts.michol = "off"; 414 %! opts.droptol = 1e-1;
292 %%! opts.shape = "upper"; 415 %! opts.michol = "off";
293 %%! U = ichol (A1, opts); 416 %! L = ichol (A1, opts);
294 %%! assert (norm (A1 - U' * U, 'fro') / norm (A1, 'fro'), 0.01, 0.01); 417 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.2065, 1e-4);
295 %! 418 %! opts.shape = "upper";
296 %!test 419 %! U = ichol (A1, opts);
297 %! opts.type = "nofill"; 420 %! assert (norm (A1 - U' * U, "fro") / norm (A1, "fro"), 0.2065, 1e-4);
298 %! opts.michol = "off";
299 %! opts.shape = "lower"; 421 %! opts.shape = "lower";
300 %! L = ichol (A1, opts); 422 %! L = ichol (A1, opts);
301 %! assert (norm (A1 - L * L', 'fro') / norm (A1, 'fro'), 0.01, 0.01); 423 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.2065, 1e-4);
302 %! 424 %!
303 %!test 425 %%!test
304 %! opts.type = "nofill"; 426 %! opts.type = "ict";
305 %! opts.michol = "on"; 427 %! opts.droptol = 1e-1;
306 %! L = ichol (A1, opts); 428 %! opts.michol = "on";
307 %! assert (norm (A1 - L * L', 'fro') / norm (A1, 'fro'), 0.02, 0.01); 429 %! L = ichol (A1, opts);
308 %! 430 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.3266, 1e-4);
309 %!test 431 %! opts.shape = "upper";
310 %! opts.type = "nofill"; 432 %! U = ichol (A1, opts);
311 %! opts.michol = "on"; 433 %! assert (norm (A1 - U' * U, "fro") / norm (A1, "fro"), 0.3266, 1e-4);
434 %! opts.shape = "lower";
312 %! opts.diagcomp = 3e-3; 435 %! opts.diagcomp = 3e-3;
313 %! L = ichol (A1, opts); 436 %! L = ichol (A1, opts);
314 %! assert (norm (A1 - L * L', 'fro') / norm (A1, 'fro'), 0.02, 0.01); 437 %! assert (norm (A1 - L * L', "fro") / norm (A1, "fro"), 0.3266, 1e-4);
315 %! 438 %!
316 %!test 439 %!test
317 %! opts.type = "ict"; 440 %! opts.type = "ict";
318 %! opts.michol = "off"; 441 %! opts.droptol = 1e-1;
319 %! opts.droptol = 1e-4; 442 %! opts.michol = "off";
320 %! L = ichol (A1, opts);
321 %! assert (norm (A1 - L * L', 'fro') / norm (A1, 'fro'), eps, eps);
322 %!
323 %!test
324 %! opts.type = "ict";
325 %! opts.michol = "off";
326 %! opts.droptol = 1e-4;
327 %! L = ichol (A2, opts); 443 %! L = ichol (A2, opts);
328 %! assert (norm (A2 - L * L', 'fro') / norm (A2, 'fro'), 5e-4, 5e-4); 444 %! assert (norm (A2 - L*L', "fro") / norm (A2, "fro"), 0.0893, 1e-4)
445 %! opts.michol = "on";
446 %! L = ichol (A2, opts);
447 %! assert (norm (A2 - L*L', "fro") / norm (A2, "fro"), 0.2377, 1e-4)
448 %!
449 %!test
450 %! opts.type = "ict";
451 %! opts.droptol = 1e-1;
452 %! opts.michol = "off";
453 %! L = ichol (A3, opts);
454 %! assert (norm (A3 - L*L', "fro") / norm (A3, "fro"), eps, eps);
455 %! opts.michol = "on";
456 %! L = ichol (A3, opts);
457 %! assert (norm (A3 - L*L', "fro") / norm (A3, "fro"), eps, eps);
458 %!
459 %!test
460 %! opts.type = "ict";
461 %! opts.droptol = 1e-1;
462 %! opts.michol = "off";
463 %! L = ichol (A4, opts);
464 %! assert (norm (A4 - L*L', "fro") / norm (A4, "fro"), 0.1224, 1e-4);
465 %! opts.michol = "on";
466 %! L = ichol (A4, opts);
467 %! assert (norm (A4 - L*L', "fro") / norm (A4, "fro"), 0.2118, 1e-4);
468 %!
469 %!test
470 %! opts.type = "ict";
471 %! opts.droptol = 1e-1;
472 %! opts.michol = "off";
473 %! L = ichol (A5, opts);
474 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.2044, 1e-4);
475 %! opts.michol = "on";
476 %! L = ichol (A5, opts);
477 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.3231, 1e-4);
478 %!test
479 %% Negative pivot
480 %! opts.type = "ict";
481 %!error ichol (A6, setup);
482 %% Complex entry in the diagonal
483 %!error ichol (A7, setup);