function L=Laplace(n)
  L=zeros(n,n)
  for i=2:n-1
      L(i,i)=2; L(i,i-1)=-1; L(i,i+1)=-1;
  end
  L(1,1)=2; L(n,n)=2; L(1,2)=-1; L(n,n-1)=-1;
endfunction

function H=Hilbert(n)
  J = 1:n;
  J = J(ones(n,1),:);
  I = J';
  E = ones(n,n);
  H = E./(I+J-1);
endfunction

function H = Inv_Hilbert(n)
  p = n;
  H = zeros(n,n);
  for i = 1:n
      if i > 1, p = ((n-i+1)*p*(n+i-1))/(i-1)^2; end
      r = p*p;
      H(i,i) = r/(2*i-1);
      for j = i+1:n
          r = -((n-j+1)*r*(n+j-1))/(j-1)^2;
          H(i,j) = r/(i+j-1);
          H(j,i) = r/(i+j-1);
      end
  end
endfunction

function M=vanderMonde(x)
  n = length(x);
  if n == 1  &  x == int(x)
    n= x;  x= 1.2 ^(0:n-1);
  end
  x = x(:);
  M = ones(n,n);
  for j = 1:n
      M(:,j) = x.^(n-j);
  end
endfunction

function P=Pascal(n)
  P = diag((-1).^(0:n-1));
  P(:,1) = ones(n,1);
  // Generate the Pascal Cholesky factor (up to signs).
  for j=2:n-1
      for i=j+1:n
          P(i,j) = P(i-1,j) - P(i-1,j-1);
      end
  end
  P= P*P';
endfunction

function P=mPascal(n)
  n1= n+1;
  P = diag((-1).^(0:n-1));
  P(:,1) = ones(n,1);
  // Generate the Pascal Cholesky factor (up to signs).
  for j=2:n-1
      for i=j+1:n
          P(i,j) = P(i-1,j) - P(i-1,j-1);
      end
  end
  P= P*P';
  
  nh= int(n/2);
  for i=1:nh
    for k=1:nh
      t= P(n1-i,n1-k); P(n1-i,n1-k)= P(i,k); P(i,k)= t;
    end
  end
endfunction

function W=Wilkinson(n,c)
  if argn(2) < 2, c=8/7; end
  W=-tril(ones(n,n)) + 2*eye(n,n);
  W(:,n)= c*ones(n,1);
  
endfunction

function L=Lotkin(n)
  // first like Hilbert
    J = 1:n;
  J = J(ones(n,1),:);
  I = J';
  E = ones(n,n);
  L = E./(I+J-1);
  // first column = ones
  L(:,1)= ones(n,1);
endfunction

function L=Lehmer(n)
  L= eye(n,n);
  for i=1:n-1
    L(i,i+1:n)= i./(i+1:n);
    L(i+1:n,i)= i./(i+1:n)';
  end
endfunction

function R=illRand(n,m,d)
//  R = n x n random matrix
//  with diagonal = const d
//  set strict upper triangular part to 0 except the last m columns

  if argn(2) < 3, d=1e-3;
    if argn(2) < 2, m= 2; end
  end
  R=rand(n,n);
  for i=1:n, R(i,i)= d; end
  for i=1:n-m-1, for j=i+1:n-m, R(i,j)= 0; end, end
endfunction

function R=illMatrix(n,loff,roff,rho,Lfac,Rfac)
  np= argn(2);
  if np < 6, Rfac= 0.01;
    if np < 5, Lfac= 2;
      if np < 4, rho= 0.2;
        if np < 3, roff= 2;
          if np < 2, loff= 4;
          end
        end
      end
    end
  end
  U= diag(rho^(0:n-1));
  
  L= eye(n,n);
  for i=2:n
    for j=max(i-loff,1):i-1, L(i,j)=(i+j)*Lfac; end
  end
  
  for i=1:n-1
    for j=max(n-roff,i+1):n, U(i,j)=(i+j)*Rfac; end
  end
  
  R= L*U;
endfunction

function G=graded_Matrix(n,p)
  if argn(2) < 2, p= 1.2; end
  G=zeros(n,n);
  for i=1:n
    G(i,:)= 1.0 ./((n:-1:1).^p + i);
  end
endfunction

function [LS,RS,minSing,maxSing]=extreme_Perturbation(A)
  [m,n]= size(A);
  [U,S,V]= svd(A,0);   //  A =  U * S * V'
  LS=[U(:,n),U(:,1)];  RS=[V(:,n),V(:,1)];
  minSing=S(min(m,n),min(m,n));
  maxSing=S(1,1);
endfunction

function [xe,dA,db,maxPert,Kond]=pert_Matrix(A,option,rhs_only)
  [m,n]= size(A);
  if  argn(2) < 3, rhs_only= %f; end

  [LS,RS,minSing,maxSing]=extreme_Perturbation(A);
  Kond=maxSing/minSing;
  
  select option
    case 'w' then
      maxPert= minSing;
      dA= -LS(:,1)*RS(:,1)';
      db= LS(:,1);
      if  rhs_only,  xe= RS(:,2);
      else           xe= RS(:,1);
      end
    
    case 'b' then
      maxPert= maxSing;
      dA= LS(:,2)*RS(:,2)';
      db= LS(:,2);
      if  rhs_only,  xe= RS(:,1);
      else           xe= RS(:,2);
      end
    
    case 'r' then
      maxPert= minSing;
      dA= rand(A);
      db= rand(m,1);  db= db/norm(db);
      xe= rand(m,1);  xe= xe/norm(xe);
    else
      warning 'illegal option for pert_Matrix';
  end
  dA= dA/norm(dA);
endfunction

function [P,cnd]=_pert_RHS(A,option,f)
  [m,n]= size(A);
  [LS,RS,minSing,maxSing]=extreme_Perturbation(A);
  cnd=maxSing/minSing;
  if argn(2) < 3, f=1.5; end
  select option
    case 'w' then, P= LS(:,1);
    case 'b' then, P= LS(:,2);
    case 'r' then, P= rand(m);
    case 'g' then, P= (1/f).^(0:m-1)';
    case 'z' then, P= RS(:,1);
    else
      warning 'illegal option for pert_RHS';
  end
  P= P/norm(P);

endfunction


function x=svd_solve(A,b)
[U,S,V]=svd(A);
x=U'*b;
x= x ./ diag(S);
x= V*x;
// one iterative refinement
  r= A*x-b;
  dx= U'*r;
  dx= dx ./ diag(S);
  dx= V*dx;
  x= x - dx;
endfunction

function x = backsub(R, b)
  n = length(b);
  x = b;
  x(n) = x(n) / R(n,n);
  for i = (n-1):-1:1
    x(i) = ( x(i) - R(i,(i+1):n)*x((i+1):n) ) / R(i,i);
  end
endfunction

function x = lr_solve(A, b)
//  B = sparse(A);
//  [h, rk] = lufact(B);
//  x = lusolve(h, b);
//  ludel(h);

  n = length(b);

  [L, R, P] = lu(A);

  y = P * b;
  for i = 2:n
    y(i)= y(i) - L(i,1:(i-1))*y(1:(i-1))
  end

  x = backsub(R, y);    
endfunction


function x = qr_solve(A, b)
  [Q, R] = qr(A);
  x = backsub(R, Q'*b);
endfunction