Mercurial > matrix-functions

comparison toolbox/diagpiv.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 [L, D, P, rho] = diagpiv(A)
2 %DIAGPIV Diagonal pivoting factorization with partial pivoting.
3 % Given a Hermitian matrix A,
4 % [L, D, P, rho] = DIAGPIV(A) computes a permutation P,
5 % a unit lower triangular L, and a real block diagonal D
6 % with 1x1 and 2x2 diagonal blocks, such that
7 % P*A*P' = L*D*L'.
8 % The Bunch-Kaufman partial pivoting strategy is used.
9 % Rho is the growth factor.
10
11 % Reference:
12 % J.R. Bunch and L. Kaufman, Some stable methods for calculating
13 % inertia and solving symmetric linear systems, Math. Comp.,
14 % 31(137):163-179, 1977.
15
16 % This routine does not exploit symmetry and is not designed to be
17 % efficient.
18
19 if norm(triu(A,1)'-tril(A,-1),1), error('Must supply Hermitian matrix.'), end
20
21 n = max(size(A));
22 k = 1;
23 D = eye(n);
24 L = eye(n);
25 pp = 1:n;
26 normA = norm(A(:),inf);
27 rho = normA;
28
29 alpha = (1 + sqrt(17))/8;
30
31 while k < n
32 [lambda, r] = max( abs(A(k+1:n,k)) );
33 r = r(1) + k;
34
35 if lambda > 0
36 swap = 0;
37 if abs(A(k,k)) >= alpha*lambda
38 s = 1;
39 else
40 temp = A(k:n,r); temp(r-k+1) = 0;
41 sigma = norm(temp, inf);
42 if alpha*lambda^2 <= abs(A(k,k))*sigma
43 s = 1;
44 elseif abs(A(r,r)) >= alpha*sigma
45 swap = 1;
46 m1 = k; m2 = r;
47 s = 1;
48 else
49 swap = 1;
50 m1 = k+1; m2 = r;
51 s = 2;
52 end
53 end
54
55 if swap
56 A( [m1, m2], : ) = A( [m2, m1], : );
57 L( [m1, m2], : ) = L( [m2, m1], : );
58 A( :, [m1, m2] ) = A( :, [m2, m1] );
59 L( :, [m1, m2] ) = L( :, [m2, m1] );
60 pp( [m1, m2] ) = pp( [m2, m1] );
61 end
62
63 if s == 1
64
65 D(k,k) = A(k,k);
66 A(k+1:n,k) = A(k+1:n,k)/A(k,k);
67 L(k+1:n,k) = A(k+1:n,k);
68 i = k+1:n;
69 A(i,i) = A(i,i) - A(i,k) * A(k,i);
70
71 elseif s == 2
72
73 E = A(k:k+1,k:k+1);
74 D(k:k+1,k:k+1) = E;
75 C = A(k+2:n,k:k+1);
76 temp = C/E;
77 L(k+2:n,k:k+1) = temp;
78 A(k+2:n,k+2:n) = A(k+2:n,k+2:n) - temp*C';
79
80 end
81
82 % Make diagonal real (see LINPACK User's Guide, p. 5.17).
83 for i=k+s:n
84 A(i,i) = real(A(i,i));
85 end
86
87 if k+s <= n
88 rho = max(rho, max(max(abs(A(k+s:n,k+s:n)))) );
89 end
90
91 else % Nothing to do.
92
93 s = 1;
94 D(k,k) = A(k,k);
95
96 end
97
98 k = k + s;
99
100 if k == n
101 D(n,n) = A(n,n);
102 break
103 end
104
105 end
106
107 if nargout >= 3, P = eye(n); P = P(pp,:); end
108 rho = rho/normA;