#
1;

function hlp

# usage : hlp
# description : supplies a short help of how to
# use the functions contained in this package.;

printf"The following functions compute the inverse of any invertible matrix:\n";
printf"\n";
printf"'inparpiv' , 'incompiv', 'invqrh' , 'inqrg'\n";
printf"\n";
printf"'invuptr' : computes the inverse of any upper triangular matrix:\n";
printf"The following functions solve the system Ax=b:\n";
printf"'linsyswp', 'linsyspp', 'linsyscp', 'linsyswf', 'linsysqrh',linsysqrg'\n";
printf"'forelm' : solves the system Ly=b; where L is a nonsingular lower triangular matrix\n";
printf"'backsub' : solves the system Ux=b; where U is a nonsingular upper triangular matrix\n";
printf" This space is for the other functions \n";
printf" For more help on how to use these function and to get a brief description on each of them type help and then the name of the function \n";

endfunction
#=======================

function y=forelm(l,b)
# description : (forward) solves the equation l*y=b;
# where l is a lower triangular matrix. Returns
# the solution vector y.
[n,m]=size(l);
y=zeros(n,1);
      for i=1:n
            s=0;
            for j=1:(i-1)
                  s=s+l(i,j)*y(j);
            endfor;
      y(i)=(1/l(i,i))*(b(i)-s);
      endfor;

endfunction;
 
function x=backsub(u,b)

# description : (backward) solves the system
# u*x=b; where u is upper triangular matrix.
# Returns the solution vector x.
[n,m]=size(u);
x=zeros(n,1);
x(n)=(b(n))/u(n,n);
for i=(n-1):-1:1
      s=0;
      for j=(i+1):n
            s=s+u(i,j)*x(j);
      endfor;
      x(i)=(1/u(i,i))*(b(i)-s);
endfor;

endfunction
 
#======================================

function pm=prmin(a)

# usage : pm=prmin(a)
# description : computes all the principal minors
# of the matrix (a) and returns them in a vector.
# Returns the vector of principal minors.

[n,m]=size(a);
if (n<>m)
      #pm=0;
      printf"ERROR : The matrix is not square. ";
      return;
endif;
if (n==m)
      pm=zeros(n,1);
      pm(1)=a(1,1);
      pm(n)=det(a);
      for i=2:(n-1)
            b=zeros(i,i);
            for j=1:i
                  for k=1:i
                        b(j,k)=a(j,k);
                  endfor;
            endfor;
            pm(i)=det(b);
      endfor;
endif;

endfunction


#========================


function id=id(n)

# usage : id=id(n)
# description : computes the n*n identity matrix.
# Returns In.

id=zeros(n,n);
for i=1:n
      id(i,i)=1;
endfor;

endfunction

#======================================

function mx=maxim(n,m)

# usage : mx=maxim(n,m)
# description : Returns the maximum of the two numbers n & m.

if (n>m)
      mx=n;
endif;
if (m>n)
      mx=m;
endif;
if (n==m)
      mx=n;
endif;

endfunction

#======================================

function elright=elr(a,prow)

# usage : elright=elr(a,prow)
# description : multiplies any matrix (a) with the elementery
# row matrix (prow). The elementery matrix (prow) is stored in a
# row vector and results from (ROW) exchanges of the identity.
# Returns a*prow.
[n,m]=size(prow);
if (m>n)
      temp=m;
      m=n;
      n=temp;
endif;
if (m<>1)
      printf"ERROR : the second input has to be a vector\n";
      return;
endif;
[n1,m1]=size(a);
if (m1 <> n)
      printf"ERROR : columns of (a) has to be the same as rows of (prow)\n";
      return;
endif;
b=zeros(n,n);
for i=1:n
      p=prow(i);
      for j=1:n
            b(j,p)=a(j,i);
      endfor;               
endfor;
elright=b;

endfunction

#======================================

function elcleft=ell(pcol,a)

# usage : elcleft=ell(pcol,a)
# description : multiplies any column permitation matrix (pcol)
# with the matrix (a). The elementery matrix (pcol) is stored in a
# row vector and results from (COLUMN) exchanges of the identity.
# Returns pcol*a.

[n,m]=size(pcol);
if (m>n)
      temp=m;
      m=n;
      n=temp;
endif;
if (m<>1)
      printf"ERROR : the first input has to be a vector\n";
      return;
endif;
[n1,m1]=size(a);
if (n1 <> n)
      printf"ERROR : columns of (pcol) has to be the same as rows of (a)\n";
      return;
endif;
b=zeros(n,n);
for i=1:n
      p=pcol(i);
      for j=1:n
            b(p,j)=a(i,j);
      endfor;               
endfor;
elcleft=b;

endfunction

#======================================

function l=lowdiag(l,d)

# usage : l=lowdiag(l,d)
# description : multiplies the lower triangular matrix (l)
# by the diagonal matrix(d). The diagonal matrix (d) is
# stored in a vector (i.e. (d) is a vector which contains
# the diagonal elements).
# Returns the product.

[n,m]=size(l);
for i=1:n
      for j=1:i
            l(i,j)=l(i,j)*d(j);
      endfor;
endfor;

endfunction

#=============

function a=elmul(m,a)


# usage : a=elmul(m,a)
# description : multiplies the elementery lower
# triangular matrix whose nondiagonal entries (multipliers)
# are stored in a vector (m). 
# Returns the product in (a).
[nn,mm]=size(a);
[w1,w2]=size(m);
if (w1<w2)
      m=m';
endif;
c=rows(a)-rows(m);
for i=(c+1):nn
      for j=1:mm
            a(i,j)=a(i,j)+(a(c,j))*(m(i-c));
      endfor;
endfor;

endfunction

#=======================

function pr=elrowmul(per,a)

# usage : pr=elrowmul(per,a)
# description : multiplies the elementery matrix
# (which is formed from row exchanges and whose
# row's indices are stored in the vector (per))
# with the matrix (a).
# Returns the product.
[t1,t2]=size(per);
if (t2>t1)
      per=per';
endif;
tt=rows(per);
orig=zeros(tt,1);
for q=1:tt
      orig(q)=q;
endfor;
if (per==orig)
      pr=a;
      return;
endif;
[n,m]=size(a);
pr=zeros(n,m);
for i=1:n
      k=per(i);
      for j=1:n
            pr(i,j)=a(k,j);
      endfor;
endfor;

endfunction

#================

function hu=hous(u)

# usage :  hu=hous(u)
# description : computes the Householder matrix
# (hu) of the vector u.

[n1,m1]=size(u);
n=maxim(n1,m1);
if (n1<n)
      u=u';
endif;
id=id(n);
hu=id-(2/(u'*u)*(u*u'));

endfunction

#======================================

function ha=housmul(u,a)

# usage :  ha=housmul(u,a)
# description : does a Householder multiply (multiplies
# h (the householder matrix of the vector u) by the
# matrix (a).
# Returns h*a.

[m,n]=size(a);
[n1,m1]=size(u);
mm=maxim(n1,m1);
if (n1<mm)
      u=u';
endif;
if (m<>mm)
      # ha=0;
      printf"ERROR : Check your dimensions \n"; #modified
      return;
endif;
if (m==mm)
      b=2/(u'*u);
      for j=1:n
            aj=zeros(m,1);
            for k=1:m
                  aj(k)=a(k,j);
            endfor;
            c=u'*aj;
            c=b*c;
            for i=1:m
                  a(i,j)=a(i,j)-c*u(i);
            endfor;
      endfor;
endif;
if (m==mm)
      ha=a;
endif;

endfunction


#======================================
function [l,u]=lugsel(a)

# description : Factorizes the matrix a into LU
# using Guassian Elimination without pivoting.
# Returns L and U.
[n,m]=size(a);
for j=1:(n-1)
      for i=(j+1):n
            a(i,j)=a(i,j)/a(j,j);
            for k=(j+1):n
                  a(i,k)=a(i,k)-a(i,j)*a(j,k);
            endfor;
      endfor;
endfor;
for i=1:n
      for j=1:(i-1)
            l(i,j)=a(i,j);
      endfor;
      l(i,i)=1;
endfor;
for i=1:n
      for j=i:n
            u(i,j)=a(i,j);
      endfor;
endfor;

endfunction

#==============

function x=linsyswp(a,b)

# description : solves the system a*x=b, by factorizing
# a into it's LU without pivoting.
# Returns the solution vector x.
#condition=(det(a))*(det(inv(a))) #mod
            
[l,u]=lugsel(a);
y=forelm(l,b);
x=backsub(u,y);
                                       

endfunction

#=================


function [per,L,U] = parpiv1(a)

# usage : [per,L,U]=parpiv1(a)
# description : Finds LU factorization of the
# matrix a, by using G.E. with partial pivoting
# and implicit row scaling.
# Returns the permitation vector "per" and the
# LU fact.

[n,m]=size(a);
s=zeros(n,1);
for i=1:n
      max=a(i,1);
      for j=1:n
            if (a(i,j)>max)
                  max=a(i,j);
            endif;
      endfor;
      s(i)=max;
endfor;
p=zeros(n,1);
temp=zeros(n,1);
for i=1:n
      p(i)=i;
endfor;
for j=1:(n-1)
      flag=j;
      max=abs(a(flag,flag))/s(flag);
      for k=j:n
            if ((abs(a(k,j))/s(j))>max)
                  max=(abs(a(k,j))/s(j));
                  tp=p(k);
                  p(k)=p(flag);
                  p(flag)=tp;
                  for q=1:n
                        temp(q)=a(k,q);
                  endfor;
                  for q=1:n
                        a(k,q)=a(flag,q);
                        a(flag,q)=temp(q);
                  endfor;
                  tp=s(k);
                  s(k)=s(flag);
                  s(flag)=tp;
            endif;
      endfor;
            for i=(j+1):n
                  a(i,j)=(a(i,j)/a(j,j));
                  for k=(j+1):n
                        a(i,k)=a(i,k)-a(i,j)*a(j,k);
                  endfor;
            endfor;
endfor;
per=p;
L=zeros(n,n);
U=zeros(n,n);
for i=1:n
      L(i,i)=1;
      for j=1:(i-1)
            L(i,j)=a(i,j);
      endfor;
endfor;
for i=1:n
      for j=i:n
            U(i,j)=a(i,j);
      endfor;
endfor;

endfunction

#=======================

function x=linsyspp1(a,b)

# description : solves the system a*x=b, by factorizing
# a into it's LU which is computed with partial pivoting
# and implicit row scaling.
# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))

            
[per,l,u]=parpiv1(a);
[n,m]=size(a);
temp=zeros(n,1);
for i=1:n
      tp=per(i);
      temp(i)=b(tp);
endfor;
b=temp;
y=forelm(l,b);

x=backsub(u,y);
                              

endfunction

#==========================

function a=invuptr(a)

# usage : a=invuptr(a)
# description : finds the inverse of the upper triangular
# matrix (a), with a flop count of (n^3)/3.
# Returns the inverse in (a) itself.
n=rows(a);
for k=n:-1:1
      a(k,k)=1/a(k,k);
      for i=(k-1):-1:1
            s=0;
            for j=(i+1):k # modified
                  s=s+a(i,j)*a(j,k);
            endfor;
            a(i,k)=(-1/a(i,i))*s;
       endfor;
endfor;

endfunction
#=======================

function a=invlowtr(a)

# usage : a=invlowtr(a)
# description : finds the inverse of the lower triangular
# matrix (a), with a flop count of (n^3)/3.
# Returns the inverse in (a) itself.

a=a';
a=invuptr(a);
a=a';

endfunction


#==============================

function is=inlu(a)

# usage : is=inlu(a)
# description : computes the inverse of the nonsingular
# square matrix (a), using LU factorization without pivoting.
# Returns the inverse if (a) is invertible and (0) if not.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
            #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
       [L,U] = lu(a);
       is = invuptr(U) * invlowtr(L);
       #is = inv(U) * inv(L);
endif;
                 

endfunction

#==========================
#========================

function [q,r] = housqr(a)

# usage : [q,r] = housqr(a)
# description : Finds the QR factorization of
# the matrix a using Householder's method.
# Returns Q and R.
[n,m]=size(a);
b=a;
#h=zeros(n-1,n,n);
ii=zeros(n,n);
for i=1:n
      ii(i,i)=1;
endfor;
q=ii;
r=ii;
for i=1:(n-1)
      ai=zeros(n,1);
      for j=i:n
            ai(j)=a(j,i);
      endfor;
      ei=zeros(n,1);
      ei(i)=1; 
      yi=ei*(-1*norm(ai,2));
      ui=ai-yi;
      hi=ii-(2/((ui)'*ui))*(ui*(ui)');
      a=hi*a;
      q=q*hi;
      r=hi*r;
endfor;
r=r*b;

endfunction

#============================

function x=linsysqrh(a,b)

# usage : x=linsysqrh(a,b)
# description : solves the system a*x=b, by factorizing
# a into it's Householder QR factorization. 
# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))

                     
[n,m]=size(a);
[q,r]=housqr(a);
y=(q')*b;
x=backsub(r,y);
                 

endfunction

#============================
function [prow,pcol,L,U] = compiv1(a)

# usage :  [prow,pcol,L,U] = compiv1(a)
# description : this algorithm uses Guassian's Elimination
# with complete pivoting to factorize the matrix a into
# LU, where L is a lower triangular matrix and U is an
# upper triangular one. This function returns L and U
# together with the row & column permitation vectors
# prow and pcol respectively.

[n,m]=size(a);
s=zeros(n,1);
#for i=1:n
#        max=a(i,1);
#        for j=1:n
#                if (a(i,j)>max)
#                        max=a(i,j);
#                endif;
#        endfor;
#        s(i)=max;
#endfor;
prow=zeros(n,1);
pcol=zeros(n,1);
tempr=zeros(n,1);
tempc=zeros(n,1);
for i=1:n
      prow(i)=i;
      pcol(i)=i;
endfor;
for j=1:(n-1)
      flag=j;
      max=abs(a(flag,flag));
      for k=j:n
            for l=j:n
                  if ((abs(a(k,l)))>max)
                        max=(abs(a(k,l)));
                        tpr=prow(k);
                        tpc=pcol(l);
                        prow(k)=prow(flag);
                        pcol(l)=pcol(flag);
                        prow(flag)=tpr;
                        pcol(flag)=tpc;
                        for q=1:n
                              tempr(q)=a(k,q);
                        endfor;
                        for q=1:n
                              a(k,q)=a(flag,q);
                        endfor;
                        for q=1:n
                              a(flag,q)=tempr(q);
                        endfor;
                       
                        for q=1:n
                              tempc(q)=a(q,l);
                        endfor;

                         for q=1:n
                              a(q,l)=a(q,flag);
                         endfor;
                         for q=1:n
                              a(q,flag)=tempc(q);
                         endfor;

                        # tp=s(k);
                        # s(k)=s(flag);
                        # s(flag)=tp;
                  endif;
            endfor;
      endfor;
            for i=(j+1):n
                  a(i,j)=(a(i,j)/a(j,j));
                  for k=(j+1):n
                        a(i,k)=a(i,k)-a(i,j)*a(j,k);
                  endfor;
            endfor;
endfor;
#perrow=prow;
#percol=pcol;
L=zeros(n,n);
U=zeros(n,n);
for i=1:n
      L(i,i)=1;
      for j=1:(i-1)
            L(i,j)=a(i,j);
      endfor;
endfor;
for i=1:n
      for j=i:n
            U(i,j)=a(i,j);
      endfor;
endfor;

endfunction

#========================
function x=linsyscp1(a,b)

# usage : x=linsyscp1(a,b)
# usage : solves the system ax=b, by factorizing the matrix
# (a) into its LU (using complete pivoting).
# Returns the solution vector (x).

[n,m]=size(a);
[prow,pcol,L,U] = compiv1(a);
temp=zeros(n,1);
for i=1:n
       
      tp=prow(i);
      temp(i)=b(tp);
endfor;
b=temp;
z=forelm(L,b);
y=backsub(U,z);
q=zeros(n,n);
for i=1:n
      q(pcol(i),i)=1;
endfor;
x=q*y;

endfunction

#========================

function x=linsyscp(a,b)

# usage : x=linsyscp(a,b)
# usage : solves the system ax=b, by factorizing the matrix
# (a) into its LU (using complete pivoting).
# Returns the solution vector (x).

                     
[n,m]=size(a);
[prow,pcol,L,U] = compiv1(a);
temp=zeros(n,1);
for i=1:n
       
      tp=prow(i);
      temp(i)=b(tp);
endfor;
b=temp;
z=forelm(L,b);
y=backsub(U,z);
q=zeros(n,n);
for i=1:n
      q(pcol(i),i)=1;
endfor;
x=q*y;
                 

endfunction

#========================


function is=incompiv1(a)

# usage : is=incompiv1(a)
# description : computes the inverse of the nonsingular
# square matrix (a), using LU factorization with complete pivoting.
# Returns the inverse if (a) is invertible and (0) if not.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
             #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
      if (det(a)<>0)
      [prow,pcol,L,U]=compiv1(a);
      pcol1=zeros(n,n);
       # for qq=1:n
       #         pos=pcol(qq);
       #         pcol1(pos,qq)=1;
       # endfor;
      inL=invlowtr(L);
      inU=invuptr(U);
      prod1=ell(pcol,inU);
      prod2=elr(inL,prow);       
      is=prod1*prod2;
      #is=(pcol1')*inv(L)*inv(U)*prow1';
      # [C,U,M,Q] = compiv1(a);   # Data
      #T =  prow * (inv(U)) * M;  #Data    
      endif;
endif;
                 

endfunction

#==========================

function [y,z] = inter(y,z);

#INTER   Interchange two vectors
#[y,z] = inter(y,z) interchanges two vectors y and z
#so that on output y contains z and z contains y.
#input  : Vectors y and z
#output : Vectors y and z

      c = y ;
      y = z;
      z = c;
     

     
endfunction

#===============

function [A,U,M,Q] = compiv(A)

#COMPIV Triangularization using Gaussian elimination with complete pivoting
#[A,U,M,Q] = compiv(A) produces an upper triangular matrix U,
#a permuted lower triangular matrix M,  and a permutation matrix Q,
#using complete pivoting so that MAQ = U.  The lower triangular part
#of the output matrix A contains the multipliers and the upper
#triangular part of A contains U.
#This program implements Algorithm 5.2.3 of the book.
#input  : Matrix A
#output : Matrices A, U, M and Q. 

      [m,n] = size(A);
      if (m~=n)
            printf"matrix A  is not square\n"; 
            return;
      endif;
      M = eye(n,n);
      Q = eye(n,n);
      U = zeros(n,n);
      for k = 1:n-1
      M1 = eye(n,n);
        d = k ;
        e = k;
        s = A(k,k);
        for i = k:n
         for j =  k : n
           if abs(A(i,j)) > abs(s)
           s = A(i,j);
           d = i;
           e = j;
        endif;
        endfor;
        endfor;
       if  (s == 0) 
        printf"the algorithm has encountered a zero pivot";
        A=[];
        U=[];
        M=[];
        Q=[];
        return;
       endif;
       p(k) = d;
        if (d ~= k | e ~=k)
         [A(k,k:n),A(d,k:n)] = inter(A(k,k:n),A(d,k:n));
         [Q(:,k),Q(:,e)] = inter(Q(:,k),Q(:,e));
         [A(:,k),A(:,e)] = inter(A(:,k),A(:,e));
          [M(k,:),M(d,:)] = inter(M(k,:),M(d,:));
        endif;
        M1(k+1:n,k) = -A(k+1:n,k)/A(k,k);
        A(k+1:n,k) = M1(k+1:n,k);
          A(k+1:n,k+1:n) = A(k+1:n,k+1:n) +M1(k+1:n,k) * A(k,k+1:n);
#-updating M---
          M(k+1:n,1:n) = M(k+1:n,1:n) +M1(k+1:n,k) * M(k,1:n);
#----
      endfor;
      U = triu(A);
      end;



function T = incompiv(A)

#INCOMPIV Inverse by complete pivoting
#T = incompiv(A) computes the inverse of a matrix A using complete
#pivoting : inv(A) = Q*inv(U)*M. 
#input  : Matrix A
#output : Matrix T
                     
      [m,n] = size(A);
      if (m~=n)
            printf"matrix A  is not square\n"  ;
            return;
      endif;
      [C,U,M,Q] = compiv(A);
      T =  Q * (inv(U)) * M;

                 

endfunction


#=========================

function [per,M,U]=parpiv(a)


# usage : [per,M,U]=parpiv(a)
# description : Finds MA=U factorization of the
# matrix a, by using G.E. with partial pivoting
# and implicit row scaling.
# Returns the permitation vector "per" and the
# M,U ; where M*A=U.
b=a;

[n,m]=size(a);
s=zeros(n,1);
for i=1:n
      max=a(i,1);
      for j=1:n
            if (a(i,j)>max)
                  max=a(i,j);
            endif;
      endfor;
      s(i)=max;
endfor;
p=zeros(n,1);
per=zeros(n,1);
for i=1:n
      per(i)=i;
endfor;
temp=zeros(n,1);
aa=id(n);
ident=id(n);
for j=1:(n-1)
      flag=j;
      id=zeros(n,1);
      for q=1:n
            id(q)=q;
      endfor;
      idd=id;
      max=abs(a(flag,flag));
      for k=j:n
            if ((abs(a(k,j)))>max)
                  max=(abs(a(k,j)));
                  tp=per(k); #mod
                  p(flag)=k;
                  per(k)=per(flag); #mod
                  per(flag)=tp; #mod
                  tem=id(k);
                  id(k)=id(flag);
                  id(flag)=tem;
                 
                  for q=1:n
                        temp(q)=a(k,q);
                  endfor;
                  
                  for q=1:n
                        a(k,q)=a(flag,q);
                        a(flag,q)=temp(q);
                  endfor;
                  tp=s(k);
                  s(k)=s(flag);
                  s(flag)=tp;
            endif;
             
      endfor;
              
                   
                  aa=elrowmul(id,aa);
       
            for i=(j+1):n
                  a(i,j)=-1*(a(i,j)/a(j,j)); #mod
                  for k=(j+1):n
                        a(i,k)=a(i,k)+a(i,j)*a(j,k); #mod
                  endfor;
             endfor;

M=zeros(n-j,1);

for q=1:n-j
      M(q)=a(j+q,j);
endfor;     
aa=elmul(M,aa);
#U=aa;
#M=U\b;
M=aa;
U=M*b;
#M*a
p(n)=n; # modified
#per=p; #mod

#per=p; #mod
endfor;
 
endfunction

#========================

function is=inparpiv(a)

# usage : is=inparpiv(a)
# description : computes the inverse of the nonsingular
# square matrix (a), using LU factorization with partial pivoting.
# Returns the inverse if (a) is invertible and (0) if not.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
            #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
      if (det(a)<>0)
            [C,M,U] = parpiv(a); #@@@@@@@@@@@@@ change them
             is = U \ M; #@@@@@@@@@@@@@@@@@@@@ change them
       endif;
endif;
                 

endfunction

#==========================
function is=inparpiv2(a)

# usage : is=inparpiv2(a)
# description : computes the inverse of the nonsingular
# square matrix (a), using LU factorization with partial pivoting.
# Returns the inverse if (a) is invertible and (0) if not.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
            #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
      if (det(a)<>0)
            [per,L,U] = parpiv1(a); #@@@@@@@@@@@@@ change them
             inL=invlowtr(L);
             inU=invuptr(U);
             prod=elr(inL,per);
             is = inU*prod; #@@@@@@@@@@@@@@@@@@@@ change them
       endif;
endif;
                 

endfunction

#==========================

function is=inparpiv1(a)

# usage : is=inparpiv1(a)
# description : computes the inverse of the nonsingular
# square matrix (a).
# Returns the inverse if (a) is invertible and (0) if not.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
            #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
      if (det(a)<>0)
            is=zeros(0,0);
            for i=1:n
                  ei=zeros(n,1);
                  ei(i)=1;
                  x=linsyspp1(a,ei);
                  is=[is,x];
            endfor;
       endif;
endif;
                 

endfunction

#==========================
 


#========================
function x=linsyspp(a,b)

# usage : x=linsyspp(a,b)
# description : solves the system a*x=b, by finding two
# matrices M & U (upper triangular) such that
# MA=U ( computed with partial pivoting)

# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))

                     
[per,M,U]=parpiv(a);
[n,mm]=size(a);
#z=M*b;           #modified
#b=z;
#temp=zeros(n,1);
#for i=1:n
#        tp=per(i);
#        temp(i)=b(tp);
#endfor;
#b=temp;
y=M*b;        #modified
x=backsub(U,y);
                 

endfunction

#========================
function cholesky = choles(a)

# usage: ch=choles(a). a is a square symmetric positive definite matrix.
# n is the dimension of the matrix a .
# description: returns the cholesk factorization of a.
[n,m]=size(a);
#Dif=a'-a
#DetA=det(a)
pm=prmin(a);
flag=0;
for i=1:n
      if (pm(i)<=0)
            flag=1;
      endif;
endfor;
if (flag==1)
      #cholesky=0;
      printf"ERROR : The matrix is not symmetric positive definite." ;
      return;
endif;
# The body of Cholesk's Algorithm.
if (flag==0);
      for i=1:n
            for j=1:(i-1)
                  s=0;
                  for k=1:(j-1)
                        s=s+a(i,k)*a(j,k);
                  endfor;
                  a(i,j)=(1/a(j,j))*(a(i,j)-s);
            endfor;
            s=0;
            for k=1:(i-1)
                  s=s+(a(i,k))^2;
            endfor;

            a(i,i)=(a(i,i)-s)^.5;
      endfor;
      cholesky=zeros(n,n);

      for i=1:n
            for j=1:i
                  cholesky(i,j)=a(i,j);
            endfor;
      endfor;
#cholesky=a;    
endif;

endfunction

#==============================

function x=linsysch(a,b)

# description : solves the system a*x=b, by factorizing
# (a) into it's Cholesky factorization (i.e. finding a
# lower triangular matrix (h) such that a=h*h') using
# Cholesky's method.
# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))
                     
h=choles(a);
y=forelm(h,b);
x=backsub(h',y);
                 

endfunction


#==============================
function h=chlgaus(a)

# usage : h=chlgaus(a)
# description : computes the Cholesky factorization
# of the symmetric positive definite matrix (a)
# (i.e. find a matrix (h) such that a=h*h') by
# first finding the lu factorization (without pivoting).
# Returns the matrix (h).
[l,u]=lugsel(a);
[n,mm]=size(a);
d=zeros(n,1);
for i=1:n
      if (u(i,i)>0)
            d(i)=(u(i,i))^(.5);
      endif;
endfor;
h=lowdiag(l,d);

endfunction

#==========================

function x=linsyscg(a,b)

# description : solves the system a*x=b, by factorizing
# (a) into it's Cholesky factorization (i.e. finding a
# lower triangular matrix (h) such that a=h*h') using LU
# method.
# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))

                     
h=chlgaus(a);
y=forelm(h,b);
x=backsub(h',y);
                 

endfunction

#==============================


function  a=invqrh(a)

# usage : a=invqrh(a)
# description : This procedure finds the inverse of
# the matrix (a), by first factorizing it into its QR
# (using Householder's method) and then calling invuptr(r).
# This is done by using the fact that inv(a) = inv(R) * Q'. 

                     
n=rows(a);
[q,r]=housqr(a);
r=invuptr(r);
a=r*q';
                 

endfunction

#========================


function x=linsyswf(a,b)

# usage : x=linsyswf(a,b)
# description : solves the system ax=b, by finding two
# matrices (m) and (u) such that  ma=u, where u is an
# upper triangular matrix.

                     
[n,m]=size(a);
temp=zeros(n,1);
for k=1:(n-1)
      flag=k;
      max=abs(a(flag,flag));
      for i=k:n
            if ((abs(a(i,k)))>max)
                  max=(abs(a(i,k)));
                  for q=1:n
                        temp(q)=a(i,q);
                  endfor;
                  
                  for q=1:n
                        a(i,q)=a(flag,q);
                  endfor;
                  for q=1:n
                        a(flag,q)=temp(q);
                  endfor;
                  tp=b(i);
                  b(i)=b(flag);
                  b(flag)=tp;
            endif;

       endfor;
       if (max==0)
            printf"ERROR : all pivots are zeros \n" ;
            return;
       endif;

            for i=(k+1):n
                  a(i,k)=(-1)*(a(i,k)/a(k,k));
                  for j=(k+1):n
                        a(i,j)=a(i,j)+a(i,k)*a(k,j);
                  endfor;

             endfor;
             for i=(k+1):n
                  b(i)=b(i)+a(i,k)*b(k);
             endfor;
endfor;
x=backsub(a,b);
                 

endfunction

#========================

function [per,M,U]=MAU(a)

# usage : [per,M,U]=MAU(a)
# description : Finds MA=U factorization of the
# matrix a, by using G.E. with partial pivoting
# and implicit row scaling.
# Returns the permitation vector "per" and the
# M,U ; where M*A=U.
[n,m]=size(a);
s=zeros(n,1);
for i=1:n
      max=a(i,1);
      for j=1:n
            if (a(i,j)>max)
                  max=a(i,j);
            endif;
      endfor;
      s(i)=max;
endfor;
p=zeros(n,1);
per=zeros(n,1);
for i=1:n
      per(i)=i;
endfor;
temp=zeros(n,1);
ident=zeros(n,n);
#aa=a;
b=a;
for q=1:n
      ident(q,q)=1;
endfor;
aa=ident;
for j=1:(n-1)
      flag=j;
      id=ident;
      max=abs(a(flag,flag));
      for k=j:n
            if ((abs(a(k,j)))>max)
                  max=(abs(a(k,j)));
                  tp=per(k); #mod
                  p(flag)=k;
                  per(k)=per(flag); #mod
                  per(flag)=tp; #mod
                 
                  for ii=1:n
                        id(k,ii)=0;
                  endfor;
                  for ii=1:n
                        id(flag,ii)=0;
                  endfor;

                  id(flag,k)=1;
                  id(k,flag)=1;
                  for q=1:n
                        temp(q)=a(k,q);
                  endfor;
                  
                  for q=1:n
                        a(k,q)=a(flag,q);
                        a(flag,q)=temp(q);
                  endfor;
                  tp=s(k);
                  s(k)=s(flag);
                  s(flag)=tp;
            endif;
             
      endfor;
            aa=id*aa;
            for i=(j+1):n
                  a(i,j)=-1*(a(i,j)/a(j,j)); #mod
                  for k=(j+1):n
                        a(i,k)=a(i,k)+a(i,j)*a(j,k); #mod
                  endfor;
             endfor;
M=ident;
for q=(j+1):n
      M(q,j)=-1*a(i,j);
endfor;     
aa=M*aa;
#U=aa;
#M=U\b;
M=aa;
U=M*b;
#M*a
p(n)=n; # modified
#per=p; #mod
endfor;


endfunction

#====================

function x=linsysmau(a,b)

# usage : x=linsysmau(a,b)
# description : solves the system a*x=b, by finding two
# matrices M & U (upper triangular) such that
# MA=U ( computed with partial pivoting)

# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))
[per,M,U]=MAU(a);
[n,mm]=size(a);
#z=M*b;           #modified
#b=z;
#temp=zeros(n,1);
#for i=1:n
#        tp=per(i);
#        temp(i)=b(tp);
#endfor;
#b=temp;
y=M*b;        #modified
x=backsub(U,y);

endfunction

#========================

function is=invermau(a)

# usage : b=invermau(a)
# description : computes the inverse of the matrix (a),
# by first finding two matrices m,u such that ma=u;
# where u is upper triangular.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
            #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
      if (det(a)<>0)
            is=zeros(0,0);
            for i=1:n
                  ei=zeros(n,1);
                  ei(i)=1;
                  x=linsyswf(a,ei);
                  is=[is,x];
            endfor;
       endif;
endif;
                 

endfunction

#==========================

function is=invermau1(a)

# usage : is=invermau1(a)
# description : computes the inverse of the matrix (a),
# by first finding two matrices m,u such that ma=u;
# where u is upper triangular.

                     
[n,m]=size(a);
if (n<>m)
      #is=0;
      printf"ERROR : The matrix has to be a square one \n";
      return;
endif;
if (n==m)
      if (det(a)==0)
            #is=0;
            printf"ERROR : The matrix has zero determinant \n"; #mod
            return;
      endif;
endif;
if (n==m)
      if (det(a)<>0)
            [per,M,U]=MAU(a);
            inU=invuptr(U);
            is=inv(elrowmul(per,U))*(elrowmul(per,M));
      endif;
endif;
                 

endfunction

#==========================


function [c,s] = givzero(x)

# GIVZERO Givens zeroing in a vector x.
# [c,s] = givzero(x) produces the Givens parameters c and s
# for a 2-vector x such that J(1,2,theta)x has a zero
# in the second place. c = cos(theta), s = sin(theta).
# This program implements Algorithm 5.5.1 of the book.
# input   : Vector x
# output  : Scalars c and s


      if abs(x(2)) > abs(x(1))
         t = x(1)/x(2);
          s = 1/((1+t*t)^.5);
          c = s*t ;
      else
         t = x(2)/x(1);
          c = 1/((1+t*t)^.5);
          s = c*t ;
      endif;
       
     
endfunction
#=========================
function [Q,R] = givqr(A)

# GIVQR QR factorization by Givens rotation
# [Q,R] = givqr(A) produces  an orthogonal matrix Q
# and a matrix R of the same size as A
# with zeros below the diagonal, such that A = QR,
# using Givens rotations.
# This program calls MATCOM program GIVZERO.
# This program implements Algorithm 5.5.4 of the book and
# also explicitly computes Q and R.
# input  : Matrix A
# output : Matrices Q and R

      [m,n] = size(A);
      Q = eye(m,m);
      for k= 1:min(n,m-1)
       for l = k+1:m
         x=[A(k,k) A(l,k)]';
         [c,s] = givzero(x);
         A1 = c * A(k,:) + s*A(l,:);
         A2 =-s * A(k,:) + c *A(l,:);
         A(k,:) = A1;
         A(l,:) = A2;
         A3 = c * Q(:,k) + s*Q(:,l);
         A4 = -s * Q(:,k) + c *Q(:,l);
         Q(:,k) = A3;
         Q(:,l) = A4;
       endfor;
      endfor;
      R = A;
 
endfunction

#=======================

function x=linsysqrg(a,b)

# usage : x=linsysqrg(a,b)
# description : solves the system a*x=b, by factorizing
# a into it's Given's QR factorization. 
# Returns the solution vector x.
#condition=(det(a))*(det(inv(a)))

                     
[n,m]=size(a);
[q,r]=givqr(a);
y=(q')*b;
x=backsub(r,y);
                 

endfunction

#======================

function  a=invqrg(a)

# usage : a=invqrg(a)
# description : This procedure finds the inverse of
# the matrix (a), by first factorizing it into its QR
# (using Given's method) and then calling invuptr(r).
# This is done by using the fact that inv(a) = inv(R) * Q'. 

                                                                 
n=rows(a);
[q,r]=givqr(a);
r=invuptr(r);
a=r*q';
                 

endfunction
#==========================

function g=gen(n)

# usage : g=gen(n)
# description : computes the matrix used in excercise 6.6 of order n.
# Returns the matrix.

# Computing the tridiagonal ((n-1)*(n-1)) matrix whose diagonal 4 and #subdiagonals -1.
# This procedure has been copied from your samples package. I can create this #using a  different method, but I didn't do this in order to save time.
a1=-1;
b1=4;
c1=-1;
m=n-1;
m1=m;

retval = zeros(m1,m1);
retval(1,1) = b1;
retval(1,2) = c1;
retval(m1,m1-1) = a1;
retval(m1,m1) = b1;
for ii = 2:(m1-1)
  retval(ii,ii-1) = a1;
  retval(ii,ii) = b1;
  retval(ii,ii+1) = c1;
endfor;
a=retval;

#----------------
# Computing the (n-1)*(n-1) identity matrix.
#
In=zeros(m1,m1);
for i=1:m1
      In(i,i)=1;             
endfor;
#---------------
#Computing the m*m zero matrix.
zer=zeros(m,m);
#---------------
# Computing - In.
#
p=-1*In;
#---------------
# Computing the system's coefficient matrix.
#
row1=[a,p,zeros(m,m*m-2*m)];
mid=[p,a,p];
row2=[mid,zeros(m,m*m-3*m)];
union=[row1;row2];
for i=1:m-4
      row=[zeros(m,i*m),mid,zeros(m,m*m-(i+3)*m)];
      union=[union;row];
endfor;
rowblast=[zeros(m,m*m-3*m),mid];
rowlast=[zeros(m,m*m-2*m),p,a];
union=[union;rowblast;rowlast];
coef=union;
if (m==3)
      coef=[a,p,zer;p,a,p;zer,p,a];
endif;
g=coef;

endfunction



#==========================

function sol=axb(a,b)
     
# usage : sol=axb(a,b).
# description : computes the solution to the system ax=b, using diffrent
# methods and outputs the solutions in three matrices. The first column
# of the first matrix contains the solution using 'linsyswp'. The second
# column of the first matrix contains the solution using inv(a)*b, where
# inv(a) is found using 'inlu'. The first column of the second matrix
# contains the solution using 'linsyspp'. The second column of the second        # matrix contains the solution using inv(a)*b, where inv(a) is found using       # 'inparpiv'. The first column of the third matrix contains the solution using      # 'linsyscp'. The second column of the third matrix contains the solution using      # inv(a)*b, where inv(a) is found using 'incompiv'. Then it computes the          # solution using 'linsyswf'. After that inv(a) using MA=U factorization (i.e.    # using the function 'invermau'). The next step is computing the solution using # 'linsysqrh' (i.e. Householder's QR factorization). Then using 'linsysqrg
#(i.e
# using the Given's method of QR-factorization
# Returns table 6.3 page 310.
#
#
printf"The first column in each matrix is for the approximated solution\n";
printf"and the second is for the exact\n";
len=rows(b);
one=zeros(len,1);
for i=1:len
        one(i)=1;
endfor;

nor2=norm(one,2);
norinf=norm(one,inf);

printf"Solution to Q #2 parts a and b respectively :\n";
printf"For linsyswp :\n";
tic();
x1=linsyswp(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"-----------\n";
printf"for inlu(a)*b :\n";
tic();
h=inlu(a);
x2=h*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"===========\n";
#--------------------
printf"Solution to Q #3 parts a and b respectively \n";
printf"(using linsyspp,inparpiv) \n";
printf"For linsyspp :\n";
tic();
x1=linsyspp(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"---------\n";
printf"(For inparpiv(a)*b) \n";
tic();
h=inparpiv(a);
x2=h*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"---------\n";
#-------------------
printf"Solution to Q #3 parts a and b respectively \n";
printf"(using linsyspp1,inparpiv1)\n";
printf"For linsyspp1 :\n";
tic();
x1=linsyspp1(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"--------\n";
printf"For inparpiv1*b :\n";
tic();
h=inparpiv1(a);
x2=h*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"-------\n";
printf"For inparpiv2*b :\n";
tic();
x3=inparpiv2(a)*b;
time=toc()
Resid2=norm((b-a*x3),2)
RelErr2=norm((one-x3),2)/nor2
RelErrInf=norm((one-x3),inf)/norinf
printf"=============\n";
#-----------------
printf"Solution to Q #4 parts a and b respectively (using Matcom's functions) :\n";
printf"For linsyscp :\n";
tic();
x1=linsyscp(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"-------\n";
printf"For incompiv(a)*b :\n";
tic();
h=incompiv(a);
x2=h*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf

printf"===============\n";
#---------------
printf"Solution to Q #4 parts a and b respectively (using my functions) :\n";
printf"For linsyscp1 :\n";
tic();
x1=linsyscp1(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"---------\n";
printf"For incompiv1(a)*b :\n";
tic();
h=incompiv1(a);
x2=h*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"===============\n";
#--------------

printf"Solution to Q #5 part a  :\n";
Resid2=norm((b-a*x1),2)
printf"For linsyswf :\n";
tic();
x1=linsyswf(a,b);
time=toc()
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"==============\n";
#---------------
printf"Solution to Q #6 : Using linsysqrh :\n";
printf"For linsysqrh :\n";
tic();
x1=linsysqrh(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"---------------\n";
printf"For invqrh(a)*b :\n";
tic();
x2=(invqrh(a))*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"================\n";
#-----------------
printf"Solution to Q #6 : Using linsysqrg :\n";
printf"For linsysqrg :\n";
tic();
x1=linsysqrg(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"------------\n";
printf"For invqrg(a)*b :\n";
tic();
x2=(invqrg(a))*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"===============\n";
#---------------

printf"Solution to Q #6 : Using linsysmau :\n";

printf"For linsysmau :\n";
tic();
x1=linsysmau(a,b);
time=toc()
Resid2=norm((b-a*x1),2)
RelErr2=norm((one-x1),2)/nor2
RelErrInf=norm((one-x1),inf)/norinf
printf"---------\n";
printf"For invermau(a)*b :\n";
tic();
x2=(invermau(a))*b;
time=toc()
Resid2=norm((b-a*x2),2)
RelErr2=norm((one-x2),2)/nor2
RelErrInf=norm((one-x2),inf)/norinf
printf"-------------\n";
printf"For invermau1(a)*b :\n";
tic();
x3=(invermau1(a))*b;
time=toc()
Resid2=norm((b-a*x3),2)
RelErr2=norm((one-x3),2)/nor2
RelErrInf=norm((one-x3),inf)/norinf
printf"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n";

endfunction

#==========================

function invers(a)

# usage : invers(a).
# description : computes the inverse of the invertble matrix a,
# using diffrent methods and outputs the solutions in order.
# The first inverse is found by 'linsyswf'(i.e. by 'invermau').
# The second by 'invqrh'. The third by 'invqrg'. The fourth by
# 'inlu'. The fifth by 'inparpiv'. The sixth by 'inparpiv1'.
# The seventh by incompiv.

B=inv(a);
norB2=norm(B,2);
norBinf=norm(B,inf);
#---------------
printf"Inverse by invermau :\n";
tic();
C=invermau(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"-------------\n";
#-------------
#printf"Inverse by invermau1 :\n";
#tic();
#C=invermau1(a);
#time=toc()
#RelErr2=(norm((B-C),2))/norB2
#RelErrInf=(norm((B-C),inf))/norBinf
#printf"--------\n";
#---------------

printf"Inverse by invqrh :\n";
tic();
C=invqrh(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"--------\n";
#------------
printf"Inverse by invqrg :\n";
tic();
C=invqrg(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"-----------";
#-------------
printf"Inverse by inlu :\n";
tic();
C=inlu(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"------------\n";
#----------------
printf"Inverse by inparpiv :\n";
tic();
C=inparpiv(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"------------\n";
#-------------------
printf"Inverse by inparpiv1 :\n";
tic();
C=inparpiv1(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"-----------\n";
#--------------
printf"Inverse by inparpiv2 :\n";
tic();
C=inparpiv2(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"----------\n";
#---------------

printf"Inverse by incompiv :\n";
tic();
C=incompiv(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf
printf"---------\n";
#----------------
printf"Inverse by incompiv1 :\n";
tic();
C=incompiv1(a);
time=toc()
RelErr2=(norm((B-C),2))/norB2
RelErrInf=(norm((B-C),inf))/norBinf


endfunction

#==========================

function linsys(a)

# usage : linsys(a)
# description : computes all the solutions to the system ax=b using
# the different methods described in the project and the matrix (a).


n=rows(a);
#printf"Here";
b=ones(n,1);
b=a*b;
axb(a,b);

endfunction

#========================

function linsystem

# usage : linsystem
# description : provides the above function by the required matrix and
# computes the solutions.

#printf"Hilbert matrix\n";
#a=hilb(10)
#linsys(a);
printf"Matrix of problem 6.6\n";
a=gen(5)
linsys(a);
printf"Matrix # 5 page 308\n";
a=[.00001,1,1;0,.00001,1;0,0,.00001]
linsys(a);
vect=zeros(10,1);
for i=1:10
      vect(i)=i;
endfor;
printf"Hankel matrix\n";
a=hankel(vect)
linsys(a);
printf"Vandermonde matrix :\n";
a=vander(vect)
linsys(a);
# printf"Random matrix :\n";
#a=rand(10)
#linsys(a);

endfunction

#==========================


function inverses

# usage : inverses
# description : computes the inverse for all the matrices in the
# project, using all the required methods.
 
#printf"Hilbert matrix\n";
#a=hilb(10)
#invers(a);
printf"Matrix of problem 6.6\n";
a=gen(5)
invers(a);
printf"Matrix # 5 page 308\n";
a=[.00001,1,1;0,.00001,1;0,0,.00001]
invers(a);
vect=zeros(10,1);
for i=1:10
      vect(i)=i;
endfor;
printf"Hankel matrix\n";
a=hankel(vect)
invers(a);
printf"Vandermonde matrix :\n";
a=vander(vect)
invers(a);
# printf"Random matrix :\n";
#a=rand(10);
#invers(a);

endfunction


#==========================
printf"If you need a help on the functions contained in this package, type 'hlp'\n";