# HG changeset patch # User John W. Eaton # Date 1523128460 14400 # Node ID baadbf803b10ec8619170489be9e364d8c1f2ee4 # Parent dd774017866bcb74e8bd9f345c6a683109f2d656# Parent dedc0128645a2ca68b9c91fa2144468296a5976f maint: Merge stable to default. diff -r dd774017866b -r baadbf803b10 scripts/sparse/pcg.m --- a/scripts/sparse/pcg.m Sat Apr 07 15:06:20 2018 -0400 +++ b/scripts/sparse/pcg.m Sat Apr 07 15:14:20 2018 -0400 @@ -33,7 +33,7 @@ ## function @code{Afun} such that @code{Afun(x) = A * x}. Additional ## parameters to @code{Afun} are passed after @var{x0}. ## -## @var{A} has to be Hermitian and Positive Definite (@nospell{HPD}). If +## @var{A} has to be Hermitian and Positive Definite (@nospell{HPD})@. If ## @code{pcg} detects @var{A} not to be positive definite, a warning ## is printed and the @var{flag} output is set. ## @@ -456,190 +456,12 @@ endif endfunction -%!test -%! ## Check that all type of inputs work -%! A = toeplitz (sparse ([2, 1 ,0, 0, 0])); -%! b = A * ones (5, 1); -%! M1 = diag (sqrt (diag (A))); -%! M2 = M1; # M1 * M2 is the Jacobi preconditioner -%! Afun = @(z) A*z; -%! M1_fun = @(z) M1 \ z; -%! M2_fun = @(z) M2 \ z; -%! [x, flag, ~, iter] = pcg (A,b); -%! assert (flag, 0) -%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2); -%! assert (flag, 0) -%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2); -%! assert (flag, 0) -%! [x, flag] = pcg (A, b, [], [], M1_fun, M2_fun); -%! assert (flag, 0) -%! [x, flag] = pcg (A, b,[],[], M1_fun, M2); -%! assert (flag, 0) -%! [x, flag] = pcg (A, b,[],[], M1, M2_fun); -%! assert (flag, 0) -%! [x, flag] = pcg (Afun, b); -%! assert (flag, 0) -%! [x, flag] = pcg (Afun, b,[],[], M1 * M2); -%! assert (flag, 0) -%! [x, flag] = pcg (Afun, b,[],[], M1, M2); -%! assert (flag, 0) -%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2); -%! assert (flag, 0) -%! [x, flag] = pcg (Afun, b,[],[], M1, M2_fun); -%! assert (flag, 0) -%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2_fun); -%! assert (flag, 0) - -%!test -%! ## solve small diagonal system -%! -%! N = 10; -%! A = diag ([1:N]); b = rand (N, 1); -%! X = A \ b; # X is the true solution -%! [x, flag] = pcg (A, b, [], N+1); -%! assert (norm (x - X) / norm (X), 0, 1e-10); -%! assert (flag, 0); - -%!test -%! ## A not positive definite -%! ## The indefiniteness of A is detected. -%! -%! N = 10; -%! A = -diag ([1:N]); b = sum (A, 2); -%! [x, flag] = pcg (A, b, [], N + 1); -%! assert (flag, 4) - - -%!test -%! ## solve tridiagonal system, do not converge in default 20 iterations -%! N = 100; -%! A = zeros (N, N); -%! for i = 1 : N - 1 # form 1-D Laplacian matrix -%! A(i:i+1, i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = ones (N, 1); -%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12); -%! assert (flag); -%! assert (relres >= 1.0); - -%!warning -%! ## solve tridiagonal system with "perfect" preconditioner which converges -%! ## in one iteration, so the eigest does not work and issues a warning. -%! -%! N = 100; -%! A = zeros (N, N); -%! for i = 1 : N - 1 # form 1-D Laplacian matrix -%! A(i:i+1, i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = ones (N, 1); -%! X = A \ b; # X is the true solution -%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b); -%! assert (norm (x - X) / norm (X), 0, 1e-6); -%! assert (flag, 0); -%! assert (iter, 1); # should converge in one iteration -%! assert (isnan (eigest), isnan ([NaN, NaN])); - -%!test -%! ## pcg detect a non-Hermitian matrix, with a considerable imaginary part -%! ## With this example, Matlab doesn't recognize the wrong type of matrix and -%! ## makes iterations until it reaches maxit -%! N = 10; -%! A = diag (1:N) + 1i * 1e-04; -%! b = ones (N, 1); -%! [x,flag] = pcg (A, b, []); -%! assert (flag, 4) - -%!test -%! ## The imaginary part is not influent (it is too small), so pcg doesn't stop -%! N = 10; -%! A = diag (1:N) + 1i * 1e-10; -%! b = ones (N, 1); -%! [x,flag] = pcg (A, b, [], N+1); -%! assert (flag, 0) -%! assert (x, A \ b, -1e-6) - -%!test -%! ## pcg solves linear system with A Hermitian positive definite -%! N = 20; -%! A = toeplitz (sparse ([4, 1, zeros(1, 18)])) + ... -%! 1i * toeplitz (sparse ([0, 1, zeros(1, 18)]), ... -%! sparse ([0, -1, zeros(1,18)])); -%! b = A * ones (N, 1); -%! Hermitian_A = ishermitian (A); -%! [x,flag] = pcg (A, b, [], 2*N); -%! assert (Hermitian_A, true) -%! assert (flag, 0) -%! assert (x, ones (N, 1), -1e-4) - -%!testif HAVE_CHOLMOD -%! ## pcg solves preconditioned linear system with A HPD -%! N = 20; -%! A = toeplitz (sparse ([4, 1, zeros(1, 18)])) + ... -%! 1i * toeplitz (sparse ([0, 1, zeros(1, 18)]), ... -%! sparse ([0, -1, zeros(1,18)])); -%! b = A * ones (N, 1); -%! M2 = chol (A + 0.1 * eye (N)); # factor of a perturbed matrix -%! M = M2' * M2; -%! Hermitian_A = ishermitian (A); -%! Hermitian_M = ishermitian (M); -%! [x,flag] = pcg (A, b, [], 2*N, M); -%! assert (Hermitian_A, true) -%! assert (Hermitian_M, true) -%! assert (flag, 0) -%! assert (x, ones (N, 1), -1e-4) - -%!test -%! ## pcg recognizes that the preconditioner matrix is singular -%! N = 3; -%! A = toeplitz ([2, 1, 0]); -%! M = [1 0 0; 0 1 0; 0 0 0]; # the last rows is zero -%! [x, flag] = pcg (A, ones(3, 1), [], [], M); -%! assert (flag, 2) - -%!test -%! A = rand (4); -%! A = A' * A; -%! [x, flag] = pcg (A, zeros (4, 1), [], [], [], [], ones (4, 1)); -%! assert (x, zeros (4, 1)) - -%!test -%! A = single (1); -%! b = 1; -%! [x, flag] = pcg (A, b); -%! assert (class (x), "single") - -%!test -%! A = 1; -%! b = single (1); -%! [x, flag] = pcg (A, b); -%! assert (class (x), "single") - -%!test -%! A = single (1); -%! b = single (1); -%! [x, flag] = pcg (A, b); -%! assert (class (x), "single") - -%!test -%!function y = Afun (x) -%! A = toeplitz ([2, 1, 0, 0]); -%! y = A * x; -%!endfunction -%! [x, flag] = pcg ("Afun", [3; 4; 4; 3]); -%! assert (x, ones(4, 1), 1e-6) - -%!test # unpreconditioned residual -%! A = toeplitz (sparse ([4, 1, 0, 0, 0])); -%! b = sum (A, 2); -%! M = toeplitz (sparse ([2, 1, 0, 0, 0])); -%! [x, flag, relres] = pcg (A, b, [], 2, M); -%! assert (relres, norm ((b - A * x)) / norm (b), 8 * eps) %!demo # simplest use %! n = 10; %! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n)); %! b = A * ones (n, 1); -%! M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)' +%! M1 = ichol (A); # for this tridiagonal case it corresponds to chol (A)' %! M2 = M1'; %! M = M1 * M2; %! x = pcg (A, b); @@ -652,27 +474,201 @@ %! M1fun = @(x) M1 \ x; %! M2fun = @(x) M2 \ x; %! x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun); -%! function y = Ap (A, x, p) # compute A^p * x +%! function y = Ap (A, x, p) # compute A^p * x %! y = x; %! for i = 1:p %! y = A * y; %! endfor %! endfunction %! Afun = @(x, p) Ap (A, x, p); -%! x = pcg (Afun, b, [], [], [], [], [], 2); # solution of A^2 * x = b +%! ## solution of A^2 * x = b +%! x = pcg (Afun, b, [], [], [], [], [], 2); %!demo %! n = 10; %! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n)); %! b = A * ones (n, 1); -%! M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed +%! M1 = ichol (A + 0.1 * eye (n)); # Perturb the factorization of A %! M2 = M1'; %! M = M1 * M2; %! -%! ## reference solution computed by pcg after two iterations +%! ## Reference solution computed by pcg after two iterations %! [x_ref, fl] = pcg (A, b, [], 2, M); %! x_ref %! -%! ## split preconditioning +%! ## Split preconditioning %! [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2); -%! x = M2 \ y # compare x and x_ref +%! x = M2 \ y # compare x and x_ref +%!test +%! ## Check that all type of inputs work +%! A = toeplitz (sparse ([2, 1 ,0, 0, 0])); +%! b = A * ones (5, 1); +%! M1 = diag (sqrt (diag (A))); +%! M2 = M1; # M1 * M2 is the Jacobi preconditioner +%! Afun = @(z) A*z; +%! M1_fun = @(z) M1 \ z; +%! M2_fun = @(z) M2 \ z; +%! [x, flag, ~, iter] = pcg (A,b); +%! assert (flag, 0); +%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2); +%! assert (flag, 0); +%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2); +%! assert (flag, 0); +%! [x, flag] = pcg (A, b, [], [], M1_fun, M2_fun); +%! assert (flag, 0); +%! [x, flag] = pcg (A, b,[],[], M1_fun, M2); +%! assert (flag, 0); +%! [x, flag] = pcg (A, b,[],[], M1, M2_fun); +%! assert (flag, 0); +%! [x, flag] = pcg (Afun, b); +%! assert (flag, 0); +%! [x, flag] = pcg (Afun, b,[],[], M1 * M2); +%! assert (flag, 0); +%! [x, flag] = pcg (Afun, b,[],[], M1, M2); +%! assert (flag, 0); +%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2); +%! assert (flag, 0); +%! [x, flag] = pcg (Afun, b,[],[], M1, M2_fun); +%! assert (flag, 0); +%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2_fun); +%! assert (flag, 0); + +%!test +%! ## solve a small diagonal system +%! N = 10; +%! A = diag ([1:N]); b = rand (N, 1); +%! [x, flag] = pcg (A, b, [], N+1); +%! assert (flag, 0); +%! assert (norm (b - A*x) / norm (b), 0, 1e-6); + +%!test +%! ## A is not positive definite +%! ## The indefiniteness of A is detected. +%! N = 10; +%! A = -diag ([1:N]); b = sum (A, 2); +%! [x, flag] = pcg (A, b, [], N + 1); +%! assert (flag, 4); + +%!test +%! ## solve tridiagonal system, do not converge in default 20 iterations +%! N = 100; +%! A = zeros (N, N); +%! for i = 1 : N - 1 # form 1-D Laplacian matrix +%! A(i:i+1, i:i+1) = [2 -1; -1 2]; +%! endfor +%! b = ones (N, 1); +%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12); +%! assert (flag); +%! assert (relres >= 1.0); + +%!warning +%! ## solve tridiagonal system with "perfect" preconditioner which converges +%! ## in one iteration, so the eigest does not work and issues a warning. +%! N = 100; +%! A = zeros (N, N); +%! for i = 1 : N - 1 # form 1-D Laplacian matrix +%! A(i:i+1, i:i+1) = [2 -1; -1 2]; +%! endfor +%! b = ones (N, 1); +%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b); +%! assert (flag, 0); +%! assert (norm (b - A*x) / norm (b), 0, 1e-6); +%! +%! assert (isnan (eigest), isnan ([NaN, NaN])); + +%!test +%! ## pcg detect a non-Hermitian matrix, with a considerable imaginary part. +%! ## In this example, Matlab does not recognize the wrong type of matrix and +%! ## makes iterations until it reaches maxit. +%! N = 10; +%! A = diag (1:N) + 1e-4*i; +%! b = ones (N, 1); +%! [x, flag] = pcg (A, b, []); +%! assert (flag, 4); + +%!test +%! ## The imaginary part is not influent (it is too small), so pcg doesn't stop +%! N = 10; +%! A = diag (1:N) + 1e-10*i; +%! b = ones (N, 1); +%! [x, flag] = pcg (A, b, [], N+1); +%! assert (flag, 0); +%! assert (norm (b - A*x) / norm (b), 0, 1e-6); + +%!test +%! ## pcg solves linear system with A Hermitian positive definite +%! N = 20; +%! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ... +%! i * sparse (toeplitz ([0, 1, zeros(1, 18)], [0, -1, zeros(1,18)])); +%! b = A * ones (N, 1); +%! Hermitian_A = ishermitian (A); +%! [x, flag] = pcg (A, b, [], 2*N); +%! assert (Hermitian_A, true); +%! assert (flag, 0); +%! assert (x, ones (N, 1), -1e-4); + +%!testif HAVE_CHOLMOD +%! ## pcg solves preconditioned linear system with A HPD +%! N = 20; +%! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ... +%! i * sparse (toeplitz ([0, 1, zeros(1, 18)], [0, -1, zeros(1,18)])); +%! b = A * ones (N, 1); +%! M2 = chol (A + 0.1 * eye (N)); # Factor of a perturbed matrix +%! M = M2' * M2; +%! Hermitian_A = ishermitian (A); +%! Hermitian_M = ishermitian (M); +%! [x, flag] = pcg (A, b, [], 2*N, M); +%! assert (Hermitian_A, true); +%! assert (Hermitian_M, true); +%! assert (flag, 0); +%! assert (x, ones (N, 1), -1e-4); + +%!test +%! ## pcg recognizes that the preconditioner matrix is singular +%! N = 3; +%! A = toeplitz ([2, 1, 0]); +%! M = [1 0 0; 0 1 0; 0 0 0]; # the last row is zero +%! [x, flag] = pcg (A, ones (3, 1), [], [], M); +%! assert (flag, 2); + +%!test +%! A = rand (4); +%! A = A' * A; +%! [x, flag] = pcg (A, zeros (4, 1), [], [], [], [], ones (4, 1)); +%! assert (x, zeros (4, 1)); + +## Test return types +%!test +%! A = single (1); +%! b = 1; +%! [x, flag] = pcg (A, b); +%! assert (class (x), "single"); + +%!test +%! A = 1; +%! b = single (1); +%! [x, flag] = pcg (A, b); +%! assert (class (x), "single"); + +%!test +%! A = single (1); +%! b = single (1); +%! [x, flag] = pcg (A, b); +%! assert (class (x), "single"); + +%!test +%!function y = Afun (x) +%! A = toeplitz ([2, 1, 0, 0]); +%! y = A * x; +%!endfunction +%! [x, flag] = pcg ("Afun", [3; 4; 4; 3]); +%! assert (x, ones (4, 1), 1e-6); + +%!test +%! ## unpreconditioned residual +%! A = toeplitz (sparse ([4, 1, 0, 0, 0])); +%! b = sum (A, 2); +%! M = toeplitz (sparse ([2, 1, 0, 0, 0])); +%! [x, flag, relres] = pcg (A, b, [], 2, M); +%! assert (norm (b - A * x) / norm (b), relres, 8 * eps); +