Mercurial > matrix-functions

comparison toolbox/cholp.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function [R, P, I] = cholp(A, piv)
2 %CHOLP Cholesky factorization with pivoting of a pos. semidefinite matrix.
3 % [R, P] = CHOLP(A) returns R and a permutation matrix P such that
4 % R'*R = P'*A*P. Only the upper triangular part of A is used.
5 % If A is not positive definite, an error message is printed.
6 %
7 % [R, P, I] = CHOLP(A) never produces an error message.
8 % If A is positive definite then I = 0 and R is the Cholesky factor.
9 % If A is not positive definite then I is positive and
10 % R is (I-1)-by-N with P'*A*P - R'*R zeros in columns 1:I-1 and
11 % rows 1:I-1.
12 % [R, I] = CHOLP(A, 0) forces P = EYE(SIZE(A)), and therefore behaves
13 % like [R, I] = CHOL(A).
14
15 % This routine is based on the LINPACK routine CCHDC. It works
16 % for both real and complex matrices.
17 %
18 % Reference:
19 % G.H. Golub and C.F. Van Loan, Matrix Computations, Second
20 % Edition, Johns Hopkins University Press, Baltimore, Maryland,
21 % 1989, sec. 4.2.9.
22
23 if nargin == 1, piv = 1; end
24
25 n = length(A);
26 pp = 1:n;
27 I = 0;
28
29 for k = 1:n
30
31 if piv
32 d = diag(A);
33 [big, m] = max( d(k:n) );
34 m = m+k-1;
35 else
36 big = A(k,k); m = k;
37 end
38 if big <= 0, I = k; break, end
39
40 % Symmetric row/column permutations.
41 if m ~= k
42 A(:, [k m]) = A(:, [m k]);
43 A([k m], :) = A([m k], :);
44 pp( [k m] ) = pp( [m k] );
45 end
46
47 A(k,k) = sqrt( A(k,k) );
48 if k == n, break, end
49 A(k, k+1:n) = A(k, k+1:n) / A(k,k);
50
51 % For simplicity update the whole of the remaining submatrix (rather
52 % than just the upper triangle).
53
54 j = k+1:n;
55 A(j,j) = A(j,j) - A(k,j)'*A(k,j);
56
57 end
58
59 R = triu(A);
60 if I > 0
61 if nargout < 3, error('Matrix must be positive definite.'), end
62 R = R(1:I-1,:);
63 end
64
65 if piv == 0
66 P = I;
67 else
68 P = eye(n); P = P(:,pp);
69 end