# Programmer: t. shores
# Last Rev: 4/27/94
# Description:
# Routines for generating real and complex matrix examples

# not a function file:
1;

function retval = tridiag(n,a,b,c)
# usage: tridiag(n,a,b,c)
# description: 
# this computes tridiagonal matrix with diags a,b,c of size n

# local variables:
# ii,jj: index variables

retval = zeros(n,n);
retval(1,1) = b;
retval(1,2) = c;
retval(n,n-1) = a;
retval(n,n) = b;
for ii = 2:(n-1)
  retval(ii,ii-1) = a;
  retval(ii,ii) = b;
  retval(ii,ii+1) = c;
endfor;

endfunction

#

function retval = golubreinsch(n,c)
# usage: golubreinsch(n,c)
# description: 
# this computes the golub-reinsch example of size n, cos c

# local variables:
# ii,jj: index variables
# s: sin theta, corresponding to input value c for cos
# tmpf: temporary float

retval = zeros(n,n);
s = sqrt(1-c*c);
for ii = 1:n
  tmpf = s^(ii-1);
  retval(ii,ii) = tmpf;
  for jj = (ii+1):n
    retval(ii,jj) = -c*tmpf;
    retval(jj,ii) = 0.0;
  endfor;
endfor;

endfunction

#

function retval = wilkplus(n)
# usage: wilkplus(n)
# description: 
# this computes an example of Wilkinson of size n

retval = zeros(n,n);
for ii = 1:n
  for jj = 1:n
    if (ii == jj-1)
      retval(ii,jj) = 1;
    elseif (ii == jj && ii < (n+1)/2)
      retval(ii,jj) = (n+1)/2 - ii;
    elseif (ii == jj && ii >= (n+1)/2)
      retval(ii,jj) = ii - (n+1)/2;
    elseif (ii == jj+1)
      retval(ii,jj) = 1;
    else
      retval(ii,jj) = 0.0;
    endif;
  endfor;
endfor;

endfunction

#

function retval = wilkminus(n)
# usage: wilkminus(n)
# description: 
# this computes a Wilkinson example of size n

retval = zeros(n,n);
for ii = 1:n
  for jj = 1:n
    if (ii == jj-1)
      retval(ii,jj) = 1;
    elseif (ii == jj)
      retval(ii,jj) = (n+1)/2-ii;
    elseif (ii == jj+1)
      retval(ii,jj) = 1;
    else
      retval(ii,jj) = 0.0;
    endif
  endfor
endfor

endfunction

#

function retval = moler(n)
# usage: moler(n)
# description: 
# this computes an example of Moler of size n

retval = zeros(n,n);
for ii = 1:n
  retval(ii,ii) = 1;
  for jj = 1:(ii-1)
    retval(ii,jj) = min(ii,jj)-2;
    retval(jj,ii) = retval(ii,jj);
  endfor;
endfor;

endfunction

#

function retval = cplxmat1(n)
# usage: cplxmat1(n)
# description: 
# this computes a complex matrix

# local variables:
# ii,jj: index variables

retval = zeros(n,n);
for ii = 1:n
  for jj = 1:n
    retval(ii,jj) = (3*ii-2*jj)+(2*ii+3*jj)*I;
  endfor;
endfor;

endfunction

#


function retval = cplxmat2a(v)
# usage: cplxmat2a(v)
# description: 
# this computes a complex matrix, similar to de rijk 2a
# the diagonal is multiplied by Housholder I-2vv^T/(v^Tv)

retval = diag(linspace(1,rows(v),rows(v)));
retval = retval*(eye(size(retval))- 2*v*v'/(v'*v));

endfunction

#

function retval = cplxmat2b(v)
# usage: cplxmat2b(v)
# description: 
# this computes a complex matrix, similar to de rijk 2b
# the diagonal is multiplied by Housholder I-2vv^T/(v^Tv)

n = rows(v);
retval = diag(linspace(n/4,-n/4+1/2,n));
retval = retval*(eye(n)- 2*v*v'/(v'*v));

endfunction

function retval = uperturb(v,tol)
# usage: uperturb(v,tol)
# description:
# this computes Householder I-2vv^T/(v^Tv), then perturbs it
# by a random complex matrix of entries at most tol.

n = rows(v);
retval = eye(n)- 2*v*v'/(v'*v) + (rand(n,n)+i*rand(n,n))*tol;

endfunction

function retval = rebar(n,lam)
#: usage: rebar(n,lam)
#description:
# this computes J'*J, where J is nxn Jordan block with lams on diag

jblock = tridiag(n,0,lam,1);
retval = jblock'*jblock;

endfunction

#
#[1 0 0 0 0 0; -2 1 0 0 0 0;3 -4  1 0 0 0; 5 -6 7 1 0 0; -8 9 -8 7 1 0;-6 5 -4 3 -2 1]
#[1 0 0 0 0 0; -1 1 0 0 0 0;1 -1  1 0 0 0; 1 -1 1 1 0 0; -1 1 -1 1 1 0;-1 1 -1 1 -1 1]
#[2 0 0 0 0 0 ;0 2 1 0 0 0 ;0 0 2 1 0 0;0 0 0 2 0 0;0 0 0 0 2 1;0 0 0 0 0 2]
#[1 0 0 0 0 0; -2 1 0 0 0 0;3 -4  1 0 0 0; 5 -4 3 1 0 0; -2 1 -2 3 1 0;-4 5 -4 3 -2 1]

function retval = rebar(n,lam)
#: usage: rebar(n,lam)
#description:
# this computes J'*J, where J is nxn Jordan block with lams on diag

jblock = tridiag(n,0,lam,1);
retval = jblock'*jblock;

endfunction