A3 Scripts

(See ShowP at the end--it shows how one can do permutation matrix operations effectively in Matlab)

P1

     function [d,c,e] = HermFactor(a,g,h)
   % a is a real column n-vector.
   % g and h are real column (n-1)-vectors.
   % Let A be the n-by-n Hermitian tridiagonal matrix with diagonal a and subdiagonal 
   % c + i*e,  where i^2 = -1. Assume that A is positive definite so that it has
   % the factorization A = L*D*L' where D is diagonal and L is unit lower bidiagonal.
   % 
   % d is a real column n-vector so D = diag(d).
   % c and e are real column (n-1)-vectors so that the subdiagonal
   % of L is  c + i*e.




% Script P1: Tests HermFactor

clc
% Generate a random example...
randn('seed',0)
n = 6;
B = triu(tril(randn(n,n) + sqrt(-1)*randn(n,n)),-1);
A = B'*B;

% Extract the diagonal and the real and imaginary part of the
% subdiagonal...
a = diag(A); g = real(diag(A,-1)); h = imag(diag(A,-1));

% Get the LDL' factorization
[d,c,e] = HermFactor(a,g,h);
L = eye(n,n) + diag( c + sqrt(-1)*e,-1)
D = diag(d)
error = A - L*D*L'

P2

     function u = NegVec(A)
   % A is an n-by-n symmetric matrix with at least one nonpositive eigenvalue.
   % v is a unit 2-norm vector with the property that v'*A*v <= 0.
   

   
% Script P2: Tests NegVec

clc
rand('seed',0)
n = 6;
% Generate a random symmetric matrix with two negative eigenvalues
D = diag([ -.1 rand(1,n-2) -.3]);
[Q,R] = qr(rand(n,n));
A = 10*Q*D*Q'
Eigenvalues = eig(A)

% Get a vector of "negative curvature" and verify...
v = NegVec(A)
alfa = v'*A*v


     function [L,D] = LDLT(A)
   % A is an n-by-n symmetric matrix
   % L is unit lower triangular, and D is diagonal
   % so A = LDL'

     [n,n] = size(A);
     D = zeros(n,n);
     L = eye(n,n);
     for k=1:n
        D(k,k) = A(k,k);
        v = A(k+1:n,k);
        L(k+1:n,k) = v/D(k,k);
        A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - L(k+1:n,k)*v';
     end

P3

     function [L,mu,pVec] = SkewBunchParlett(A)
   % A is an n-by-n symmetric matrix with n = 2m.
   % mu is a column m-vector, L is unit lower triangular, and pVec is a permutation of 1:n 
   % so that if P = I(:,pVec) then PAP' = LDL' is the Skew Bunch-Parlett factorization
   % where D = diag(D1,...,Dm), with Dk = [0 mu(k); -mu(k) 0] for k=1:m.




% Script P3: Tests SkewBunchParlett

clc
randn('seed',0)
disp('A small example..')

% Generate a random skew symmetric matrix...
n = 6;
A = randn(n,n); A = A- A'

% Compute the SkewBunchParlett factorization

[L,mu,pVec] = SkewBunchParlett(A);

% Display P*A*P' - L*D*L' ...

v = zeros(n-1,1);
v(1:2:n) = mu; 
D = -diag(v,-1) + diag(v,+1)
L = L
error = A(pVec,pVec) - L*D*L'


% A bigger problem..
disp(' ')
disp('An n = 500 example...')
n = 500;
A = randn(n,n); A = A- A';
[L,mu,pVec] = SkewBunchParlett(A);
v = zeros(n-1,1); v(1:2:n) = mu; 
D = -diag(v,-1) + diag(v,+1);
Error = norm(A(pVec,pVec) - L*D*L',1)/norm(A,1)

P4

     function x = SolveSkew(L,mu,pVec,b)
   % b is a column n-vector with n = 2m. Solves Ax = b where A is a nonsingular n-by-n 
   % skew-symmetric matrix and L, mu, and pVec are output from SkewBunchParlett(A).




% Script P4: Tests SolveSkew

clc
randn('seed',0);
n = 20;
A = randn(n,n); A = A - A';
[L,mu,pVec] = SkewBunchParlett(A);
b = randn(n,1);
x = SolveSkew(L,mu,pVec,b)

% Permutation matrix representation and use

  n = 5;
  I = eye(n,n);
  A = magic(5); A = A + A';
  
  clc
  disp('Representation:')
  pVec = [ 2 4 1 5 3]
  P = I(pVec,:)
  pause
  
  clc
  disp('P times a vector')
  x = 10*(1:n)'
  disp('y = P*x'); disp(' ' ); disp(P*x)
  disp('z = x(pVec)') ; disp(' ');  disp(x(pVec))
  pause
  
  clc
  disp('P'' times a vector')
  P = P
  P_transpose = P'
  x = 10*(1:n)'
  disp('y = P''*x'); disp(' ' ); disp(P'*x)
  disp('z(pVec) = x') ; disp(' ');  z(pVec) = x; disp(z)
  pause
  
  clc
  disp('P times a matix')
  A = magic(5)
  P = P
  disp('B = P*A'); disp(' '); disp(P*A)
  disp('C = A(pVec,:)'); disp(' '); disp(A(pVec,:))
  pause
  
  clc
  disp('P times a matrix times P''')
  A = magic(5); A = A+A'
  P = P
  disp('B = P*A*P'''); disp(' '); disp(P*A*P')
  disp('C = A(pVec,pVec)'); disp(' '); disp(A(pVec,pVec))
  pause

  clc
  disp('P = P2*P1')
  disp(' ')
  disp('pVec1 = [2 4 1 5 3]'); pVec1 = [2 4 1 5 3];
  disp('P1 = I(pVec1,:)'); P1 = I(pVec1,:); disp(' '); disp(P1)
  disp(' ')
  disp('pVec2 = [4 5 3 1 2]'); pVec2 = [4 5 3 1 2];
  disp('P2 = I(pVec2,:)'); P2 = I(pVec2,:); disp(' '); disp(P2)
  disp(' ')
  disp('P = P2*P1'); P = P2*P1; disp(' '); disp(P)
  disp('pVec = pVec1(pVec2)'); pVec = pVec1(pVec2); disp(' '); disp(pVec)