octave:2> A=rand(4); octave:3> b=A*ones(4,1); octave:4> x=A\b x = 1.00000 1.00000 1.00000 1.00000 octave:5> format long e octave:6> x=A\b x = 1.00000000000000e+00 1.00000000000000e+00 9.99999999999998e-01 9.99999999999998e-01 octave:7> [L,U]=lu(A); octave:8> U U = Columns 1 through 3: 6.77627188327638e-01 6.53274745327522e-01 3.61120573626780e-01 0.00000000000000e+00 -5.68175855259200e-01 -3.88062547892324e-02 0.00000000000000e+00 0.00000000000000e+00 2.71732762806760e-01 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 Column 4: 6.56378386216804e-01 -4.84580211516749e-01 -2.38787846856919e-01 9.13789837369361e-02 octave:9> format octave:10> U U = 0.67763 0.65327 0.36112 0.65638 0.00000 -0.56818 -0.03881 -0.48458 0.00000 0.00000 0.27173 -0.23879 0.00000 0.00000 0.00000 0.09138 octave:11> L L = 0.96046 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.32639 -0.95372 0.34265 1.00000 0.61181 -0.82252 1.00000 0.00000 octave:12> U\(L\b) ans = 1.00000 1.00000 1.00000 1.00000 octave:13> [L,U,P]=lu(A) L = 1.00000 0.00000 0.00000 0.00000 0.96046 1.00000 0.00000 0.00000 0.61181 -0.82252 1.00000 0.00000 0.32639 -0.95372 0.34265 1.00000 U = 0.67763 0.65327 0.36112 0.65638 0.00000 -0.56818 -0.03881 -0.48458 0.00000 0.00000 0.27173 -0.23879 0.00000 0.00000 0.00000 0.09138 P = Permutation Matrix 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 octave:14> octave:14> U\(L\(P*b)) ans = 1.00000 1.00000 1.00000 1.00000 octave:15> help lu 'lu' is a built-in function from the file libinterp/corefcn/lu.cc -- Built-in Function: [L, U] = lu (A) -- Built-in Function: [L, U, P] = lu (A) -- Built-in Function: [L, U, P, Q] = lu (S) -- Built-in Function: [L, U, P, Q, R] = lu (S) -- Built-in Function: [...] = lu (S, THRES) -- Built-in Function: Y = lu (...) -- Built-in Function: [...] = lu (..., "vector") Compute the LU decomposition of A. If A is full subroutines from LAPACK are used and if A is sparse then UMFPACK is used. The result is returned in a permuted form, according to the optional return value P. For example, given the matrix 'a = [1, 2; 3, 4]', [l, u, p] = lu (A) returns l = 1.00000 0.00000 0.33333 1.00000 u = 3.00000 4.00000 0.00000 0.66667 p = 0 1 1 0 The matrix is not required to be square. When called with two or three output arguments and a spare input matrix, 'lu' does not attempt to perform sparsity preserving column permutations. Called with a fourth output argument, the sparsity preserving column transformation Q is returned, such that 'P * A * Q = L * U'. Called with a fifth output argument and a sparse input matrix, 'lu' attempts to use a scaling factor R on the input matrix such that 'P * (R \ A) * Q = L * U'. This typically leads to a sparser and more stable factorization. An additional input argument THRES, that defines the pivoting threshold can be given. THRES can be a scalar, in which case it defines the UMFPACK pivoting tolerance for both symmetric and unsymmetric cases. If THRES is a 2-element vector, then the first element defines the pivoting tolerance for the unsymmetric UMFPACK pivoting strategy and the second for the symmetric strategy. By default, the values defined by 'spparms' are used ([0.1, 0.001]). Given the string argument "vector", 'lu' returns the values of P and Q as vector values, such that for full matrix, 'A (P,:) = L * U', and 'R(P,:) * A (:, Q) = L * U'. With two output arguments, returns the permuted forms of the upper and lower triangular matrices, such that 'A = L * U'. With one output argument Y, then the matrix returned by the LAPACK routines is returned. If the input matrix is sparse then the matrix L is embedded into U to give a return value similar to the full case. For both full and sparse matrices, 'lu' loses the permutation information. See also: luupdate, chol, hess, qr, qz, schur, svd. Additional help for built-in functions and operators is available in the online version of the manual. Use the command 'doc ' to search the manual index. Help and information about Octave is also available on the WWW at http://www.octave.org and via the help@octave.org mailing list. octave:16> [L,U,P]=lu(A,'vector') L = 1.00000 0.00000 0.00000 0.00000 0.96046 1.00000 0.00000 0.00000 0.61181 -0.82252 1.00000 0.00000 0.32639 -0.95372 0.34265 1.00000 U = 0.67763 0.65327 0.36112 0.65638 0.00000 -0.56818 -0.03881 -0.48458 0.00000 0.00000 0.27173 -0.23879 0.00000 0.00000 0.00000 0.09138 P = 2 1 4 3 octave:17> U\(L\(b(P))) ans = 1.00000 1.00000 1.00000 1.00000 octave:18> B=A*A'; octave:19> S=[1,2,3;4,5,6;7,8,9] S = 1 2 3 4 5 6 7 8 9 octave:20> det(S) ans = 6.6613e-16 octave:21> B B = 0.54326 0.68671 0.36513 0.56467 0.68671 1.44719 1.18296 1.40523 0.36513 1.18296 1.15113 1.26154 0.56467 1.40523 1.26154 1.51391 octave:22> R=chol(B) R = 0.73706 0.93169 0.49539 0.76612 0.00000 0.76101 0.94796 0.90859 0.00000 0.00000 0.08420 0.24589 0.00000 0.00000 0.00000 0.20242 octave:23> R'*R ans = 0.54326 0.68671 0.36513 0.56467 0.68671 1.44719 1.18296 1.40523 0.36513 1.18296 1.15113 1.26154 0.56467 1.40523 1.26154 1.51391 octave:24> format long e octave:25> det(S) ans = 6.66133814775094e-16 octave:26> ver ---------------------------------------------------------------------- GNU Octave Version 3.8.1 GNU Octave License: GNU General Public License Operating System: Linux 3.16.0-36-generic #48~14.04.1-Ubuntu SMP Wed Apr 15 13:11:28 UTC 2015 x86_64 ---------------------------------------------------------------------- no packages installed. octave:27> R' ans = Columns 1 through 3: 7.37059070914793e-01 0.00000000000000e+00 0.00000000000000e+00 9.31690341889998e-01 7.61012654731554e-01 0.00000000000000e+00 4.95392620593008e-01 9.47958200803291e-01 8.42000491054639e-02 7.66117672958429e-01 9.08592425374172e-01 2.45885852817047e-01 Column 4: 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 2.02420925683525e-01 octave:28> format octave:29> R' ans = 0.73706 0.00000 0.00000 0.00000 0.93169 0.76101 0.00000 0.00000 0.49539 0.94796 0.08420 0.00000 0.76612 0.90859 0.24589 0.20242 octave:30> R R = 0.73706 0.93169 0.49539 0.76612 0.00000 0.76101 0.94796 0.90859 0.00000 0.00000 0.08420 0.24589 0.00000 0.00000 0.00000 0.20242 octave:31> b=B*ones(4,1); octave:32> x=B\b x = 1.00000 1.00000 1.00000 1.00000 octave:33> R\(R'\b) ans = 1.00000 1.00000 1.00000 1.00000 octave:34> R'\(R\b) ans = -17.951 28.412 -468.655 625.508 octave:35> help pcg 'pcg' is a function from the file /usr/share/octave/3.8.1/m/sparse/pcg.m -- Function File: X = pcg (A, B, TOL, MAXIT, M1, M2, X0, ...) -- Function File: [X, FLAG, RELRES, ITER, RESVEC, EIGEST] = pcg (...) Solve the linear system of equations 'A * X = B' by means of the Preconditioned Conjugate Gradient iterative method. The input arguments are * A can be either a square (preferably sparse) matrix or a function handle, inline function or string containing the name of a function which computes 'A * X'. In principle, A should be symmetric and positive definite; if 'pcg' finds A not to be positive definite, a warning is printed and the FLAG output will be set. * B is the right-hand side vector. * TOL is the required relative tolerance for the residual error, 'B - A * X'. The iteration stops if 'norm (B - A * X)' <= TOL * norm (B). If TOL is omitted or empty then a tolerance of 1e-6 is used. * MAXIT is the maximum allowable number of iterations; if MAXIT is omitted or empty then a value of 20 is used. * M = M1 * M2 is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by 'pcg' 'P * X = M \ B', with 'P = M \ A'. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method. Instead of matrices M1 and M2, the user may pass two functions which return the results of applying the inverse of M1 and M2 to a vector (usually this is the preferred way of using the preconditioner). If M1 is omitted or empty '[]' then no preconditioning is applied. If M2 is omitted, M = M1 will be used as a preconditioner. * X0 is the initial guess. If X0 is omitted or empty then the function sets X0 to a zero vector by default. The arguments which follow X0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to 'pcg'. See the examples below for further details. The output arguments are * X is the computed approximation to the solution of 'A * X = B'. * FLAG reports on the convergence. A value of 0 means the solution converged and the tolerance criterion given by TOL is satisfied. A value of 1 means that the MAXIT limit for the iteration count was reached. A value of 3 indicates that the (preconditioned) matrix was found not to be positive definite. * RELRES is the ratio of the final residual to its initial value, measured in the Euclidean norm. * ITER is the actual number of iterations performed. * RESVEC describes the convergence history of the method. 'RESVEC(i,1)' is the Euclidean norm of the residual, and 'RESVEC(i,2)' is the preconditioned residual norm, after the (I-1)-th iteration, 'I = 1, 2, ..., ITER+1'. The preconditioned residual norm is defined as 'norm (R) ^ 2 = R' * (M \ R)' where 'R = B - A * X', see also the description of M. If EIGEST is not required, only 'RESVEC(:,1)' is returned. * EIGEST returns the estimate for the smallest 'EIGEST(1)' and largest 'EIGEST(2)' eigenvalues of the preconditioned matrix 'P = M \ A'. In particular, if no preconditioning is used, the estimates for the extreme eigenvalues of A are returned. 'EIGEST(1)' is an overestimate and 'EIGEST(2)' is an underestimate, so that 'EIGEST(2) / EIGEST(1)' is a lower bound for 'cond (P, 2)', which nevertheless in the limit should theoretically be equal to the actual value of the condition number. The method which computes EIGEST works only for symmetric positive definite A and M, and the user is responsible for verifying this assumption. Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A) n = 10; A = diag (sparse (1:n)); b = rand (n, 1); [l, u, p, q] = luinc (A, 1.e-3); EXAMPLE 1: Simplest use of 'pcg' x = pcg (A, b) EXAMPLE 2: 'pcg' with a function which computes 'A * X' function y = apply_a (x) y = [1:N]' .* x; endfunction x = pcg ("apply_a", b) EXAMPLE 3: 'pcg' with a preconditioner: L * U x = pcg (A, b, 1.e-6, 500, l*u) EXAMPLE 4: 'pcg' with a preconditioner: L * U. Faster than EXAMPLE 3 since lower and upper triangular matrices are easier to invert x = pcg (A, b, 1.e-6, 500, l, u) EXAMPLE 5: Preconditioned iteration, with full diagnostics. The preconditioner (quite strange, because even the original matrix A is trivial) is defined as a function function y = apply_m (x) k = floor (length (x) - 2); y = x; y(1:k) = x(1:k) ./ [1:k]'; endfunction [x, flag, relres, iter, resvec, eigest] = ... pcg (A, b, [], [], "apply_m"); semilogy (1:iter+1, resvec); EXAMPLE 6: Finally, a preconditioner which depends on a parameter K. function y = apply_M (x, varargin) K = varargin{1}; y = x; y(1:K) = x(1:K) ./ [1:K]'; endfunction [x, flag, relres, iter, resvec, eigest] = ... pcg (A, b, [], [], "apply_m", [], [], 3) References: 1. C.T. Kelley, 'Iterative Methods for Linear and Nonlinear Equations', SIAM, 1995. (the base PCG algorithm) 2. Y. Saad, 'Iterative Methods for Sparse Linear Systems', PWS 1996. (condition number estimate from PCG) Revised version of this book is available online at See also: sparse, pcr. Additional help for built-in functions and operators is available in the online version of the manual. Use the command 'doc ' to search the manual index. Help and information about Octave is also available on the WWW at http://www.octave.org and via the help@octave.org mailing list. octave:36> pcg(B,b,1e-4,10) ans = 0.82372 1.22087 0.78604 1.03900 octave:37> B B = 0.54326 0.68671 0.36513 0.56467 0.68671 1.44719 1.18296 1.40523 0.36513 1.18296 1.15113 1.26154 0.56467 1.40523 1.26154 1.51391 octave:38> pcg(B,b,1e-8,10) ans = 1.00000 1.00000 1.00000 1.00000 octave:39> pcg(A,b,1e-8,10) warning: pcg: maximum number of iterations (10) reached warning: the initial residual norm was reduced 4.82353 times. ans = 0.61822 1.82862 2.46870 1.60928 octave:40> F=@(x) [x(1)^2+x(2)^2-1;sin(pi*x(1)/2)+x(2)^3]; octave:41> F([1;1]) ans = 1 2 octave:42> J=@(x) [ > octave:42> J=@(x) [2*x(1),2*x(2);cos(pi*x(1)/2)*pi/2,3*x(2)^2]; octave:43> J([1;1]) ans = 2.0000e+00 2.0000e+00 9.6184e-17 3.0000e+00 octave:44> y=rand(2,1) y = 0.38967 0.51372 octave:45> epsilon=1e-4 epsilon = 1.0000e-04 octave:46> x=[1;1] x = 1 1 octave:47> (F(x+epsilon*v)-F(x))/epsilon error: 'v' undefined near line 1 column 14 error: evaluating argument list element number 1 octave:47> (F(x+epsilon*y)-F(x))/epsilon ans = 1.8068 1.5412 octave:48> J(x)*y ans = 1.8068 1.5412 octave:49> x x = 1 1 octave:50> tol=1e-6 tol = 1.0000e-06 octave:51> delta=-J(x)\F(x); octave:52> norm(delta,inf) ans = 0.66667 octave:53> while (norm(delta,inf)>tol) > x=x+delta; > delta=-J(x)\F(x); > norm(delta,inf) > end ans = 2.4144 ans = 1.1626 ans = 0.45504 ans = 0.20920 ans = 0.042336 ans = 0.0015990 ans = 2.7104e-06 ans = 7.0069e-12 octave:54> J=@(x) [2*x(1),2*x(2);cos(pi*x(1)/2),3*x(2)^2]; octave:55> octave:55> x=[1;1] x = 1 1 octave:56> delta=-J(x)\F(x); octave:57> while (norm(delta,inf)>tol) > x=x+delta; > delta=-J(x)\F(x); > norm(delta,inf) > end ans = 2.5912 ans = 1.1882 ans = 0.46950 ans = 0.37695 ans = 0.10931 ans = 0.017807 ans = 0.0037146 ans = 7.7675e-04 ans = 1.6358e-04 ans = 3.4406e-05 ans = 7.2386e-06 ans = 1.5228e-06 ans = 3.2037e-07 octave:58> a=1 a = 1 octave:59> f=@(x) x+a; octave:60> f(2) ans = 3 octave:61> f(3) ans = 4 octave:62> a=2 a = 2 octave:63> f93) parse error: syntax error >>> f93) ^ octave:63> f(3) ans = 4 octave:64> f=@(x) x+a; octave:65> f(3) ans = 5 octave:66> f=@(x,a) x+a; octave:67> A=[4,1,0;1,4,1;0,1,4] A = 4 1 0 1 4 1 0 1 4 octave:68> A=[4,1,0,0;1,4,1,0;0,1,4,1;0,0,1,4] A = 4 1 0 0 1 4 1 0 0 1 4 1 0 0 1 4 octave:69> toeplitz([4,1,0,0]) ans = 4 1 0 0 1 4 1 0 0 1 4 1 0 0 1 4 octave:70> toeplitz([4,1,zeros(1,8)]) ans = 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 1 4 octave:71> sparse(toeplitz([4,1,zeros(1,8)])) ans = Compressed Column Sparse (rows = 10, cols = 10, nnz = 28 [28%]) (1, 1) -> 4 (2, 1) -> 1 (1, 2) -> 1 (2, 2) -> 4 (3, 2) -> 1 (2, 3) -> 1 (3, 3) -> 4 (4, 3) -> 1 (3, 4) -> 1 (4, 4) -> 4 (5, 4) -> 1 (4, 5) -> 1 (5, 5) -> 4 (6, 5) -> 1 (5, 6) -> 1 (6, 6) -> 4 (7, 6) -> 1 (6, 7) -> 1 (7, 7) -> 4 (8, 7) -> 1 (7, 8) -> 1 (8, 8) -> 4 (9, 8) -> 1 (8, 9) -> 1 (9, 9) -> 4 (10, 9) -> 1 (9, 10) -> 1 (10, 10) -> 4 octave:72> B=sparse(toeplitz([4,1,zeros(1,8)])); octave:73> B(4,3) ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 1 [100%]) (1, 1) -> 1 octave:74> det(B) ans = 5.6472e+05 octave:75> B=sparse(toeplitz([4,1,zeros(1,8)])); octave:76> toeplitz([4,1,sparse(zeros(1,8))]) ans = Compressed Column Sparse (rows = 10, cols = 10, nnz = 28 [28%]) (1, 1) -> 4 (2, 1) -> 1 (1, 2) -> 1 (2, 2) -> 4 (3, 2) -> 1 (2, 3) -> 1 (3, 3) -> 4 (4, 3) -> 1 (3, 4) -> 1 (4, 4) -> 4 (5, 4) -> 1 (4, 5) -> 1 (5, 5) -> 4 (6, 5) -> 1 (5, 6) -> 1 (6, 6) -> 4 (7, 6) -> 1 (6, 7) -> 1 (7, 7) -> 4 (8, 7) -> 1 (7, 8) -> 1 (8, 8) -> 4 (9, 8) -> 1 (8, 9) -> 1 (9, 9) -> 4 (10, 9) -> 1 (9, 10) -> 1 (10, 10) -> 4 octave:77> octave:77> sparse([1,1],[1,2],[4,1],1,8) ans = Compressed Column Sparse (rows = 1, cols = 8, nnz = 2 [25%]) (1, 1) -> 4 (1, 2) -> 1 octave:78> toeplitz(sparse([1,1],[1,2],[4,1],1,8)) ans = Compressed Column Sparse (rows = 8, cols = 8, nnz = 22 [34%]) (1, 1) -> 4 (2, 1) -> 1 (1, 2) -> 1 (2, 2) -> 4 (3, 2) -> 1 (2, 3) -> 1 (3, 3) -> 4 (4, 3) -> 1 (3, 4) -> 1 (4, 4) -> 4 (5, 4) -> 1 (4, 5) -> 1 (5, 5) -> 4 (6, 5) -> 1 (5, 6) -> 1 (6, 6) -> 4 (7, 6) -> 1 (6, 7) -> 1 (7, 7) -> 4 (8, 7) -> 1 (7, 8) -> 1 (8, 8) -> 4 octave:79> diary off