//GMG last modified: 04/25/2000 
//////////////////////////////////////////////////////////////////////////////
version="$Id: linalg.lib,v 1.7 2000-12-22 14:12:32 greuel Exp $";
category="Linear Algebra";
info="
LIBRARY:  linalg.lib	Algorithmic Linear Algebra
AUTHOR:	  Ivor Saynisch (ivs@math.tu-cottbus.de)

PROCEDURES:
 inverse(A);	    matrix, the inverse of A (for constant matrices)
 inverse_B(A);      list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) )
 inverse_L(A);      list(matrix Inv,poly p),Inv*A=p*En ( using lift )
 sym_gauss(A);      symmetric gaussian algorithm
 orthogonalize(A);  Gram-Schmidt orthogonalization
 diag_test(A);      test whether A can be diagnolized
 busadj(A);         coefficients of Adj(E*t-A) and coefficients of det(E*t-A)
 charpoly(A,v);     characteristic polynomial of A ( using busadj(A) ) 
 adjoint(A);        adjoint of A ( using busadj(A) )
 det_B(A);          determinant of A. ( uses busadj(A) )
 gaussred(A);       gaussian reduction: P*A=U*S, S a row reduced form of A
 gaussred_pivot(A); gaussian reduction: P*A=U*S, uses row pivoting
 gauss_nf(A);       gaussian normal form of A
 mat_rk(A);         rank of constant matrix A
 U_D_O(A);          P*A=U*D*O, P,D,U,O = permutaion,diag,lower-,upper-triang
 pos_def(A,i);      test symmetric matrix for positive definiteness
";

LIB "matrix.lib";
LIB "ring.lib";

//////////////////////////////////////////////////////////////////////////////
// help functions
//////////////////////////////////////////////////////////////////////////////
static proc abs(poly c)
"RETURN: absolut value of c, c must be constants"
{
  if(c>=0){ return(c);}
  else{ return(-c);}
}

static 
proc const_mat(matrix A)
"RETURN:   1 (0) if A is (is not) a constant matrix"
{
  int i;
  int n=ncols(A);
  def BR=basering;
  changeord("@R","dp,c",BR);
  matrix A=fetch(BR,A);
  for(i=1;i<=n;i=i+1){
    if(deg(lead(A)[i])>=1){
      //"input is not a constant matrix";
      kill @R;
      setring BR;
      return(0);
    }
  }
  kill @R;
  setring BR;
  return(1);  
}
//////////////////////////////////////////////////////////////////////////////
static proc red(matrix A,int i,int j)
"USAGE:    red(A,i,j);  A = constant matrix
          reduces column j with respect to A[i,i] and column i
          reduces row j with respect to A[i,i] and row i
RETURN:   matrix
"
{	
  module m=module(A);
  
  if(A[i,i]==0){
    m[i]=m[i]+m[j];
    m=module(transpose(matrix(m)));
    m[i]=m[i]+m[j];
    m=module(transpose(matrix(m)));
  }
  
  A=matrix(m);
  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
  m=module(transpose(matrix(m)));
  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
  m=module(transpose(matrix(m)));
  
  return(matrix(m));
}

//////////////////////////////////////////////////////////////////////////////
proc inner_product(vector v1,vector v2)
"RETURN:   inner product <v1,v2> "
{
  int k;
  if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); } 
  return ((transpose(matrix(v1,k,1))*matrix(v2,k,1))[1,1]);
}

/////////////////////////////////////////////////////////////////////////////
// user functions
/////////////////////////////////////////////////////////////////////////////

proc inverse(matrix A)
"USAGE:    inverse(A);  A = constant, square matrix
RETURN:   matrix, the inverse matrix of A
NOTE:     parameters and minpoly are allowed, uses std applied to A enlarged
          bei unit matrix
EXAMPLE:  example inverse; shows an example"
{
  //erlaubt Parameterringe und minpoly

  int i;
  int n=nrows(A);
  string mp=string(minpoly);

  if (ncols(A)!=n){
    "input is not a square matrix";
    matrix invmat;
    return(invmat);
  }
  if(!const_mat(A)){
    "input is not a constant matrix";
    matrix invmat;
    return(invmat);
  }
	  
  def BR=basering;
  changeord("@R","c,dp",BR);
  execute("minpoly="+mp+";");
  matrix A=fetch(BR,A);
  
  matrix B=transpose(concat(A,unitmat(n)));
  matrix D=transpose(std(B));
  matrix D1=submat(D,1..n,1..n);
  matrix D2=submat(D,1..n,(n+1)..2*n);
  B=transpose(concat(D2,D1));

  changeord("@@R","C,dp",@R);
  execute("minpoly="+mp+";");
  matrix B=imap(@R,B);
 
  for(i=1;i<=n/2;i=i+1){
    B=permcol(B,i,n-i+1);
  }
  B=std(B);
  module C=module(B);
  
  for (i=1;i<=n;i=i+1){
    if (deg(B[n+i,i])<0){
      "matrix is not invertible";
      setring BR;
      kill @R;kill @@R;
      matrix invmat;
      return(invmat);
    }
    else{ C[i] = C[i]/B[n+i,i];}
  }
  
  B=transpose(matrix(C));	
  matrix invmat[n][n]=B[1..n,1..n];
  setring BR;
  matrix invmat=fetch(@@R,invmat);
  kill @R;kill @@R;	
  return(invmat);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x),lp;
  matrix A[4][4]=5,7,3,8,4,6,9,2,3,7,5,7,2,5,6,12;
  print(A);
  print(inverse(A));
  print(inverse(A)*A);         //test for correctness
}

//////////////////////////////////////////////////////////////////////////////
proc sym_gauss(matrix A)
"USAGE:    sym_gauss(A);  A = symmetric matrix
RETURN:   matrix, diagonalisation with symmetric gauss algorithm
EXAMPLE:  example sym_gauss; shows an example"
{
  int i,j;
  int n=nrows(A);

  if (ncols(A)!=n){
    "input is not a square matrix";;
    return(A);
  } 
  
  if(!const_mat(A)){
    "input is not a constant matrix";
    return(A);
  }
  
  if(deg(std(A-transpose(A))[1])!=-1){ 
    "input is not a symmetric matrix";
    return(A);
  }

  for(i=1; i<n; i++){
    for(j=i+1; j<=n; j++){
      if(A[i,j]!=0){ A=red(A,i,j); }
    }
  }
  
  return(A);
}
example
{"EXAMPLE:"; echo = 2;
  ring r=0,(x),lp;
  matrix A[2][2]=1,4,4,15;
  print(A);
  print(sym_gauss(A));
}

//////////////////////////////////////////////////////////////////////////////
proc orthogonalize(matrix A) 
"USAGE:    orthogonalize(A); A = constant matrix
RETURN:    matrix, orthogonal basis of the colum space of A
EXAMPLE:   example orthogonalize; shows an example "
{
  int i,j;
  int n=ncols(A);
  poly k;

  if(!const_mat(A)){
    "input is not a constant matrix";
    matrix B;
    return(B);
  }  

  module B=module(interred(A));
  
  for(i=1;i<=n;i=i+1) {
    for(j=1;j<i;j=j+1) {
      k=inner_product(B[j],B[j]);
      if (k==0) { "Error: vector with zero length"; return(matrix(B)); }
      B[i]=B[i]-(inner_product(B[i],B[j])/k)*B[j];
    }
  }
  
  return(matrix(B));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x),lp;
  matrix A[4][4]=5,6,12,4,7,3,2,6,12,1,1,2,6,4,2,10;
  print(A);
  print(orthogonalize(A));
}

////////////////////////////////////////////////////////////////////////////
proc diag_test(matrix A)
"USAGE:	  diag_test(A); A = const square matrix
NOTE:     Der Test ist nur fuer zerfallende Matrizen aussagefaehig.
          Parameterringe werden nicht unterstuetzt (benutzt factorize,gcd).
RETURN:   int,  1 falls A diagonalisierbar, 0 nicht diagonalisierbar
               -1 keine Aussage moeglich, da A nicht zerfallend.
EXAMPLE:  example diag_test; shows an example"
{
  int i,j;
  int n     = nrows(A);	
  string mp = string(minpoly);
  string cs = charstr(basering);
  int np=0;

  if(ncols(A) != n) { 
    "// input is not a square matrix";
    return(-1);
  }		

   if(!const_mat(A)){
    "// input is not a constant matrix";
    return(-1);
  }  

  //Parameterring wegen factorize nicht erlaubt
  for(i=1;i<size(cs);i=i+1){
    if(cs[i]==","){np=np+1;} //Anzahl der Parameter
  }
  if(np>0){
	"// rings with parameters not allowed";
	return(-1);
  }

  //speichern des aktuellen Rings  
  def BR=basering;
  //setze R[t]
  execute("ring rt=("+charstr(basering)+"),(@t,"+varstr(basering)+"),lp;");
  execute("minpoly="+mp+";");
  matrix A=imap(BR,A);	

  intvec z;
  intvec s;
  poly X;         //characteristisches Polynom 
  poly dXdt;      //Ableitung von X nach t
  ideal g;	      //ggT(X,dXdt)
  poly b;	      //Komponente der Busadjunkten-Matrix
  matrix E[n][n]; //Einheitsmatrix
	
  E=E+1;
  A=E*@t-A;
  X=det(A);

  matrix Xfactors=matrix(factorize(X,1));	//zerfaellt die Matrtix ?
  int nf=ncols(Xfactors);

  for(i=1;i<=nf;i++){
    if(lead(Xfactors[1,i])>=@t^2){
      //"Die Matrix ist nicht zerfallend";
      setring BR;
      return(-1);
    }
  }	

  dXdt=diff(X,@t);
  g=std(ideal(gcd(X,dXdt)));

  //Busadjunkte
  z=2..n;	
  for(i=1;i<=n;i++){
    s=2..n;  
    for(j=1;j<=n;j++){
      b=det(submat(A,z,s));
    
      if(0!=reduce(b,g)){
	//"Die Matrix ist nicht diagonalisierbar!";
	setring BR;
	return(0);
      }
      
      s[j]=j;
    }
    z[i]=i;
  }
  
  //"Die Matrix ist diagonalisierbar";  	
  setring BR;
  return(1);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x),dp;
  matrix A[4][4]=6,0,0,0,0,0,6,0,0,6,0,0,0,0,0,6;
  print(A);
  diag_test(A);
}

//////////////////////////////////////////////////////////////////////////////
proc busadj(matrix A)
"USAGE:   busadj(A);  A = square matrix
RETURN:  list L; 
         L[1] contains the (n+1) coefficients of the characteristic polynomial 
         X of A, i.e.
              X = L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1])*t^n
         L[2] contains the n (nxn)-matrices Hk which are the coefficients of
              the busadjoint bA=adjoint(E*t-A) of A, i.e.
              bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0,  ( Hk=L[2][k+1] ) 
EXAMPLE: example busadj; shows an example"
{
  int k;
  int n    = nrows(A);
  matrix E = unitmat(n);
  matrix H[n][n];
  matrix B[n][n];
  list bA, X, L;
  poly a;

  if(ncols(A) != n) { 
    "input is not a square matrix";
    return(L);
  }	

  bA   = E;
  X[1] = 1;
  for(k=1; k<n; k++){
    B  = A*bA[1];              //bA[1] is the last H
    a  = -trace(B)/k;
    H  = B+a*E;
    bA = insert(bA,H);
    X  = insert(X,a);
  }
  B = A*bA[1];
  a = -trace(B)/n;
  X = insert(X,a);

  L = insert(L,bA);
  L = insert(L,X);
  return(L);
}
example
{ "EXAMPLE"; echo = 2;
  ring r = 0,(t,x),lp;
  matrix A[2][2] = 1,x2,x,x2+3x;
  print(A);
  list L = busadj(A);
  poly X = L[1][1]+L[1][2]*t+L[1][3]*t2; X;
  matrix bA[2][2] = L[2][1]+L[2][2]*t; 
  print(bA);               //the busadjoint of A;
  print(bA*(t*unitmat(2)-A));  
}

//////////////////////////////////////////////////////////////////////////////
proc charpoly(matrix A, list #)
"USAGE:   charpoly(A,[v]); A = square matrix, v = name of a ringvariable
NOTE:    A must be constant in the variable v. The computation uses busadj(A).
RETURN:  poly, the characteristic polynomial det(E*v-A)
EXAMPLE: example charpoly; shows an example"
{
  int i,j;
  string s,v;
  list L;
  int n     = nrows(A);
  poly X    = 0;
  def BR    = basering;
  string mp = string(minpoly);

  if(ncols(A) != n){ 	
    "// input is not a square matrix";
    return(X);
  }	

  //test for correct variable
  if( size(#)==0 ){ 
    #[1] = varstr(1); 
  }
  if( typeof(#[1]) == "string"){
     v = #[1];
  }
  else{
    "// 2nd argument must be a name of a variable not contained in the matrix";
    return(X);
  }  
       
  j=-1;
  for(i=1; i<=nvars(BR); i++){
    if(varstr(i)==v){j=i;}
  }
  if(j==-1){
    "// "+v+" is not a variable in the basering";
    return(X);
  }
  
  //var can not be in A
  s="Wp(";
  for( i=1; i<=nvars(BR); i++ ){
    if(i!=j){ s=s+"0";}
    else{ s=s+"1";}
    if( i<nvars(BR)) {s=s+",";}
  }
  s=s+")";

  changeord("@R",s);
  execute("minpoly="+mp+";");
  matrix A = imap(BR,A);
  for(i=1; i<=n; i++){
    if(deg(lead(A)[i])>=1){
      "matrix must not contain the variable "+v;
      kill @R;
      setring BR;
      return(X);
    }
  }

  //get coefficients and build the char. poly
  kill @R;
  setring BR;
  L = busadj(A);
  for(i=1; i<=n+1; i++){
    execute("X=X+L[1][i]*"+v+"^"+string(i-1)+";");
  }
  
  return(X); 
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,t),dp;
  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
  print(A);
  charpoly(A,"t");
}

//////////////////////////////////////////////////////////////////////////////
proc adjoint(matrix A)
"USAGE:    adjoint(A);  A = square matrix
NOTE:     computation uses busadj(A)
RETURN:   adjoint 
EXAMPLE:  example adjoint; shows an example"
{
  int n=nrows(A);
  matrix Adj[n][n];
  list L;

  if(ncols(A) != n) { 	
    "// input is not a square matrix";
    return(Adj);
  }	
  
  L  = busadj(A);
  Adj= (-1)^(n-1)*L[2][1];
  return(Adj);
  
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(t,x),lp;
  matrix A[2][2]=1,x2,x,x2+3x;
  print(A);
  matrix Adj[2][2]=adjoint(A);
  print(Adj);                    //Adj*A=det(A)*E
  print(Adj*A);
}

//////////////////////////////////////////////////////////////////////////////
proc inverse_B(matrix A)
"USAGE:     inverse_B(A);  A = square matrix
RETURN:    list Inv with Inv[1]=matrix I and Inv[2]=poly p 
           and I*A = unitmat(n)*p; 
NOTE:      p=1 if 1/det(A) is computable and  p=det(A) if not
EXAMPLE:   example inverse_B; shows an example"
{
  int i;
  int n=nrows(A);
  matrix I[n][n];
  poly factor;
  list L;
  list Inv;
  
  if(ncols(A) != n) { 	
    "input is not a square matrix";
    return(I);
  }	
  
  L=busadj(A);
  I=module(-L[2][1]);        //+-Adj(A)   

  if(reduce(1,std(L[1][1]))==0){ 
    I=I*lift(L[1][1],1)[1][1];      
    factor=1;
  }
  else{ factor=L[1][1];}     //=+-det(A) or 1
  Inv=insert(Inv,factor);
  Inv=insert(Inv,matrix(I));

  return(Inv);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,y),lp;
  matrix A[3][3]=x,y,1,1,x2,y,x+y,6,7;
  print(A);
  list Inv=inverse_B(A);
  print(Inv[1]);
  print(Inv[2]);
  print(Inv[1]*A);
}

//////////////////////////////////////////////////////////////////////////////
proc det_B(matrix A)
"USAGE:     det_B(A);  A any matrix  
RETURN:    returns the determinant of A
NOTE:      the computation uses the busadj algorithm
EXAMPLE:   example det_B; shows an example"
{
  int n=nrows(A); 
  list L;

  if(ncols(A) != n){ return(0);}

  L=busadj(A);
  return((-1)^n*L[1][1]);
}
example
{ "EXAMPLE"; echo=2; 
  ring r=0,(x),dp;
  matrix A[10][10]=random(2,10,10)+unitmat(10)*x;
  print(A);
  det_B(A); 
}

//////////////////////////////////////////////////////////////////////////////
proc inverse_L(matrix A)
"USAGE:     inverse_L(A);  A = square matrix
RETURN:    list Inv representing a left inverse of A, i.e
           Inv[1] = matrix I and Inv[2]=poly p such that
           I*A = unitmat(n)*p; 
NOTE:      p=1 if 1/det(A) is computable and  p=det(A) if not
EXAMPLE:   example inverse_L; shows an example"
{
  int n=nrows(A);
  matrix I;
  matrix E[n][n]=unitmat(n);
  poly factor;
  poly d=1;
  list Inv;

  if (ncols(A)!=n){
    "// input is not a square matrix";
    return(I);
  }

  d=det(A);
  if(d==0){
    "// matrix is not invertible";
    return(Inv);
  }

  // test if 1/det(A) exists
  if(reduce(1,std(d))!=0){ E=E*d;}
  
  I=lift(A,E);
  if(I==unitmat(n)-unitmat(n)){ //catch error in lift 
    "// matrix is not invertible";
    return(Inv);
  }

  factor=d;      //=det(A) or 1
  Inv=insert(Inv,factor);
  Inv=insert(Inv,I);
  
  return(Inv);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,y),lp;
  matrix A[3][3]=x,y,1,1,x2,y,x+y,6,7;
  print(A);
  list Inv=inverse_L(A);
  print(Inv[1]);
  print(Inv[2]);
  print(Inv[1]*A);
}

//////////////////////////////////////////////////////////////////////////////
proc gaussred(matrix A)
"USAGE:     gaussred(A);   A any constant matrix 
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
           gives a row reduced matrix S, a permutation matrix P and a
           normalized lower triangular matrix U, with P*A=U*S            
EXAMPLE:   example gaussred; shows an example"
{
  int i,j,l,k,jp,rang;
  poly c,pivo;
  list Z;
  int n = nrows(A);
  int m = ncols(A);
  int mr= n; //max. rang
  matrix P[n][n] = unitmat(n);
  matrix U[n][n] = P;

  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(Z);
  }

  if(n>m){mr=m;} //max. rang

  for(i=1;i<=mr;i=i+1){
    if((i+k)>m){break}; 
    
    //Test: Diagonalelement=0
    if(A[i,i+k]==0){
      jp=i;pivo=0;
      for(j=i+1;j<=n;j=j+1){
	c=abs(A[j,i+k]);
	if(pivo<c){ pivo=c;jp=j;}
      }
      if(jp != i){       //Zeilentausch
	for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
	  c=A[i,j];         
	  A[i,j]=A[jp,j];
	  A[jp,j]=c;
	}
	for(j=1;j<=n;j=j+1){ //Zeilentausch in P
	  c=P[i,j];        
	  P[i,j]=P[jp,j];
	  P[jp,j]=c;  
	}
      }
      if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
    }                          //i sollte im naechsten Lauf nicht erhoeht sein
    
    //Eliminationsschritt   
    for(j=i+1;j<=n;j=j+1){ 
      c=A[j,i+k]/A[i,i+k];
      for(l=i+k+1;l<=m;l=l+1){
	A[j,l]=A[j,l]-A[i,l]*c;
      }
      A[j,i+k]=0;  // nur wichtig falls k>0 ist
      A[j,i]=c;    // bildet U
    }
  rang=i;    
  }
 
  for(i=1;i<=mr;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      U[j,i]=A[j,i];
      A[j,i]=0;
    }
  }
  
  Z=insert(Z,rang);
  Z=insert(Z,A);
  Z=insert(Z,U);
  Z=insert(Z,P);
  
  return(Z);
}
example
{ "EXAMPLE";echo=2;
  ring r=0,(x),dp;
  matrix A[5][4]=1,3,-1,4,2,5,-1,3,1,3,-1,4,0,4,-3,1,-3,1,-5,-2;
  print(A);
  list Z=gaussred(A);   //construct P,U,S s.t. P*A=U*S
  print(Z[1]);          //P
  print(Z[2]);          //U
  print(Z[3]);          //S
  print(Z[4]);          //rank
  print(Z[1]*A);        //P*A
  print(Z[2]*Z[3]);     //U*S
}

//////////////////////////////////////////////////////////////////////////////
proc gaussred_pivot(matrix A)
"USAGE:     gaussred_pivot(A);   A any constant matrix
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
           gives n row reduced matrix S, a permutation matrix P and a
           normalized lower triangular matrix U, with P*A=U*S
NOTE:      with row pivoting            
EXAMPLE:   example gaussred_pivot; shows an example"
{
  int i,j,l,k,jp,rang;
  poly c,pivo;
  list Z;
  int n=nrows(A);
  int m=ncols(A);
  int mr=n; //max. rang
  matrix P[n][n]=unitmat(n);
  matrix U[n][n]=P;

  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(Z);
  }

  if(n>m){mr=m;} //max. rang

  for(i=1;i<=mr;i=i+1){
    if((i+k)>m){break}; 
    
    //Pivotisierung
    pivo=abs(A[i,i+k]);jp=i;
    for(j=i+1;j<=n;j=j+1){
      c=abs(A[j,i+k]);
      if(pivo<c){ pivo=c;jp=j;}
    }
    if(jp != i){ //Zeilentausch
      for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
	c=A[i,j];         
	A[i,j]=A[jp,j];
	A[jp,j]=c;
      }
      for(j=1;j<=n;j=j+1){ //Zeilentausch in P
	c=P[i,j];        
	P[i,j]=P[jp,j];
	P[jp,j]=c;  
      }
    }
    if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
                               //i sollte im naechsten Lauf nicht erhoeht sein
    //Eliminationsschritt   
    for(j=i+1;j<=n;j=j+1){ 
      c=A[j,i+k]/A[i,i+k];
      for(l=i+k+1;l<=m;l=l+1){
	A[j,l]=A[j,l]-A[i,l]*c;
      }
      A[j,i+k]=0;  // nur wichtig falls k>0 ist
      A[j,i]=c;    // bildet U
    }
  rang=i;    
  }
 
  for(i=1;i<=mr;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      U[j,i]=A[j,i];
      A[j,i]=0;
    }
  }
  
  Z=insert(Z,rang);
  Z=insert(Z,A);
  Z=insert(Z,U);
  Z=insert(Z,P);
  
  return(Z);
}
example
{ "EXAMPLE";echo=2;
  ring r=0,(x),dp;
  matrix A[5][4] = 1, 3,-1,4,
                   2, 5,-1,3,
                   1, 3,-1,4,
                   0, 4,-3,1,
                  -3,1,-5,-2;
  list Z=gaussred_pivot(A);  //construct P,U,S s.t. P*A=U*S
  print(Z[1]);               //P
  print(Z[2]);               //U
  print(Z[3]);               //S
  print(Z[4]);               //rank
  print(Z[1]*A);             //P*A
  print(Z[2]*Z[3]);          //U*S
}

//////////////////////////////////////////////////////////////////////////////
proc gauss_nf(matrix A)
"USAGE:    gauss_nf(A); A any constant matrix
RETURN:   matrix; gauss normal form of A
EXAMPLE:  example gauss_nf; shows an example"
{
  list Z;
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(A);
  }
  Z = gaussred(A);
  return(Z[3]);
}
example
{ "EXAMPLE";echo=2;
  ring r = 0,(x),dp;
  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
  print(gauss_nf(A));
}

//////////////////////////////////////////////////////////////////////////////
proc mat_rk(matrix A)
"USAGE:    mat_rk(A); A any constant matrix
RETURN:   int, rank of A 
EXAMPLE:  example mat_rk; shows an example"
{
  list Z;
  if(!const_mat(A)){
    "//input is not a constant matrix";
    return(-1);
  }
  Z = gaussred(A);
  return(Z[4]);
}
example
{ "EXAMPLE";echo=2;
  ring r = 0,(x),dp;
  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
  mat_rk(A);
}

//////////////////////////////////////////////////////////////////////////////
proc U_D_O(matrix A)
"USAGE:     U_D_O(A);   constant invertible matrix A
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=D , Z[4]=O
           gives a permutation matrix P, 
           a normalized lower triangular matrix U ,
           a diagonal matrix D, and a normalized upper triangular matrix O
           with P*A=U*D*O
NOTE:      Z[1]=-1 means that A is not regular
EXAMPLE:   example U_D_O; shows an example"
{
  int i,j;
  list Z,L;
  int n=nrows(A);
  matrix O[n][n]=unitmat(n);
  matrix D[n][n];

  if (ncols(A)!=n){
    "// input is not a square matrix";
    return(Z);
  }
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(Z);
  }
  
  L=gaussred(A);

  if(L[4]!=n){
    "// input is not an invertible matrix";
    Z=insert(Z,-1);  //hint for calling procedures 
    return(Z);
  }

  D=L[3];

  for(i=1; i<=n; i++){
    for(j=i+1; j<=n; j++){
      O[i,j] = D[i,j]/D[i,i];
      D[i,j] = 0;
    }
  }

  Z=insert(Z,O);
  Z=insert(Z,D);
  Z=insert(Z,L[2]);
  Z=insert(Z,L[1]);
  return(Z);
}
example
{ "EXAMPLE";echo=2;
  ring r = 0,(x),dp;
  matrix A[5][5] = 10, 4,  0, -9,  8, 
                   -3, 6, -6, -4,  9, 
                    0, 3, -1, -9, -8,
                   -4,-2, -6, -10,10,
                   -9, 5, -1, -6,  5;
  list Z = U_D_O(A);              //construct P,U,D,O s.t. P*A=U*D*O 
  print(Z[1]);                    //P
  print(Z[2]);                    //U
  print(Z[3]);                    //D
  print(Z[4]);                    //O
  print(Z[1]*A);                  //P*A
  print(Z[2]*Z[3]*Z[4]);          //U*D*O
}

//////////////////////////////////////////////////////////////////////////////
proc pos_def(matrix A)
"USAGE:     pos_def(A); A = constant, symmetric square matrix       
RETURN:    int; 1 if A is positive definit , 0 if not, -1 if unknown 
EXAMPLE:   example pos_def; shows an example"
{
  int j;
  list Z;
  int n = nrows(A);
  matrix H[n][n];

  if (ncols(A)!=n){
    "// input is not a square matrix";
    return(0);
  }
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(-1);
  }	
  if(deg(std(A-transpose(A))[1])!=-1){ 
    "// input is not a hermitian (symmetric) matrix";
    return(-1);
  }
  
  Z=U_D_O(A);

  if(Z[1]==-1){
    return(0);
  }  //A not regular, therefore not pos. definit

  H=Z[1];
  //es fand Zeilentausch statt: also nicht positiv definit 
  if(deg(std(H-unitmat(n))[1])!=-1){
    return(0);
  }
  
  H=Z[3];
  
  for(j=1;j<=n;j=j+1){
    if(H[j,j]<=0){ 
      return(0);
    } //eigenvalue<=0, not pos.definit
  }

  return(1); //positiv definit;
}
example
{ "EXAMPLE"; echo=2;
  ring r = 0,(x),dp;
  matrix A[5][5] = 20,  4,  0, -9,   8,
                    4, 12, -6, -4,   9,
                    0, -6, -2, -9,  -8, 
                   -9, -4, -9, -20, 10, 
                    8,  9, -8,  10, 10;
  pos_def(A);
  matrix B[3][3] =  3,  2,  0,
                    2, 12,  4,
                    0,  4,  2;
  pos_def(B);
}

//////////////////////////////////////////////////////////////////////////////
proc linsolve(matrix A, matrix b)
"USAGE:     linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix
RETURN:    a solution of inhomogeneous linear system A*X = b 
EXAMPLE:   example linsolve; shows an example"
{
  int i,j,k,rc,r;
  poly c;
  list Z;
  int n  = nrows(A);
  int m  = ncols(A);
  int n_b= nrows(b);
  matrix Ab[n][m+1];
  matrix X[m][1];
  	
  if(ncols(b)!=1){
    "// right hand side b is not a nx1 matrix";
    return(X);
  }

  if(!const_mat(A)){
    "// input hand is not a constant matrix";
    return(X);
  }	
  
  if(n_b>n){
    for(i=n; i<=n_b; i++){
      if(b[i,1]!=0){
	"// right hand side b not in Image(A)";
	return X;
      }
    }    
  }
  
  if(n_b<n){ 
    matrix copy[n_b][1]=b;
    matrix b[n][1]=0;
    for(i=1;i<=n_b;i=i+1){
      b[i,1]=copy[i,1];
    }
  }
  
  r=mat_rk(A);
	
  //1. b constant vector
  if(const_mat(b)){	
    //extend A with b
    for(i=1; i<=n; i++){
      for(j=1; j<=m; j++){
        Ab[i,j]=A[i,j];
      }
      Ab[i,m+1]=b[i,1];
    }
      
    //Gauss reduction
    Z  = gaussred(Ab);
    Ab = Z[3];  //normal form
    rc = Z[4];  //rank(Ab)  
    //print(Ab);

    if(r<rc){
	"// no solution";
        return(X); 
    }
    k=m;	
    for(i=r;i>=1;i=i-1){
      
      j=1;	
      while(Ab[i,j]==0){j=j+1;}// suche Ecke 
   
      for(;k>j;k=k-1){ X[k]=0;}//springe zur Ecke
     

      c=Ab[i,m+1]; //i-te Komponene von b
      for(j=m;j>k;j=j-1){
	c=c-X[j,1]*Ab[i,j];
      }
      if(Ab[i,k]==0){
	X[k,1]=1; //willkuerlich
      }
      else{ 
  	X[k,1]=c/Ab[i,k];
      }
      k=k-1;
      if(k==0){break;}
    }
    
    
  }//endif (const b)
  else{  //b not constant
    "// !not implemented!";
 
  }

  return(X);
}
example
{ "EXAMPLE";echo=2;
  ring r=0,(x),dp;
  matrix A[3][2] = -4,-6,
                    2, 3, 
                   -5, 7;
  matrix b[3][1] = 10,
                   -5,
                    2;
  matrix X = linsolve(A,b);
  print(X);
  print(A*X);
}
//////////////////////////////////////////////////////////////////////////////