Mercurial > matrix-functions

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