///////////////////////////////////////////////////////////////////////////////
version="$Id: matrix.lib,v 1.26.2.2 2002/10/07 14:05:02 lossen Exp $";
category="Linear Algebra";
info="
LIBRARY:  matrix.lib    Elementary Matrix Operations

PROCEDURES:
 compress(A);           matrix, zero columns from A deleted
 concat(A1,A2,..);      matrix, concatenation of matrices A1,A2,...
 diag(p,n);             matrix, nxn diagonal matrix with entries poly p
 dsum(A1,A2,..);        matrix, direct sum of matrices A1,A2,...
 flatten(A);            ideal, generated by entries of matrix A
 genericmat(n,m[,id]);  generic nxm matrix [entries from id]
 is_complex(c);         1 if list c is a complex, 0 if not
 outer(A,B);            matrix, outer product of matrices A and B
 power(A,n);            matrix/intmat, n-th power of matrix/intmat A
 skewmat(n[,id]);       generic skew-symmetric nxn matrix [entries from id]
 submat(A,r,c);         submatrix of A with rows/cols specified by intvec r/c
 symmat(n[,id]);        generic symmetric nxn matrix [entries from id]
 tensor(A,B);           matrix, tensor product of matrices A nd B
 unitmat(n);            unit square matrix of size n
 gauss_col(A);          transform a matrix into col-reduced Gauss normal form
 gauss_row(A);          transform a matrix into row-reduced Gauss normal form
 addcol(A,c1,p,c2);     add p*(c1-th col) to c2-th column of matrix A, p poly
 addrow(A,r1,p,r2);     add p*(r1-th row) to r2-th row of matrix A, p poly
 multcol(A,c,p);        multiply c-th column of A with poly p
 multrow(A,r,p);        multiply r-th row of A with poly p
 permcol(A,i,j);        permute i-th and j-th columns
 permrow(A,i,j);        permute i-th and j-th rows
 rowred(A[,any]);       reduction of matrix A with elementary row-operations
 colred(A[,any]);       reduction of matrix A with elementary col-operations
 rm_unitrow(A);         remove unit rows and associated columns of A
 rm_unitcol(A);         remove unit columns and associated rows of A
          (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "ring.lib";
LIB "random.lib";
///////////////////////////////////////////////////////////////////////////////

proc compress (A)
"USAGE:   compress(A); A matrix/ideal/module/intmat/intvec
RETURN:  same type, zero columns/generators from A deleted
         (if A=intvec, zero elements are deleted)
EXAMPLE: example compress; shows an example
"
{
   if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); }
   if( typeof(A)=="intmat" or typeof(A)=="intvec" )
   {
      ring r=0,x,lp;
      if( typeof(A)=="intvec" ) { intmat C=transpose(A); kill A; intmat A=C; }
      module m = matrix(A);
      if ( size(m) == 0)
      { intmat B; }
      else
      { intmat B[nrows(A)][size(m)]; }
      int i,j;
      for( i=1; i<=ncols(A); i=i+1 )
      {
         if( m[i]!=[0] )
         {
            j=j+1;
            B[1..nrows(A),j]=A[1..nrows(A),i];
         }
      }
      if( defined(C) ) { return(intvec(B)); }
      return(B);
    }
   return(simplify(A,2));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   matrix A[3][4]=1,0,3,0,x,0,z,0,x2,0,z2,0;
   print(A);
   print(compress(A));
   module m=module(A); show(m);
   show(compress(m));
   intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0;
   compress(B);
   intvec C=0,0,1,2,0,3;
   compress(C);
}
///////////////////////////////////////////////////////////////////////////////
proc concat (list #)
"USAGE:   concat(A1,A2,..); A1,A2,... matrices
RETURN:  matrix, concatenation of A1,A2,.... Number of rows of result matrix
         is max(nrows(A1),nrows(A2),...)
EXAMPLE: example concat; shows an example
"
{
   int i;
   module B=module(#[1]);
   for( i=2; i<=size(#); i=i+1 ) { B=B,module(#[i]); }
   return(matrix(B));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   matrix A[3][3]=1,2,3,x,y,z,x2,y2,z2;
   matrix B[2][2]=1,0,2,0; matrix C[1][4]=4,5,x,y;
   print(A);
   print(B);
   print(C);
   print(concat(A,B,C));
}
///////////////////////////////////////////////////////////////////////////////

proc diag (list #)
"USAGE:   diag(p,n); p poly, n integer
         diag(A);   A matrix
RETURN:  diag(p,n): diagonal matrix, p times unit matrix of size n.
@*       diag(A)  : n*m x n*m diagonal matrix with entries all the entries of
                    the nxm matrix A, taken from the 1st row, 2nd row etc of A
EXAMPLE: example diag; shows an example
"
{
   if( size(#)==2 ) { return(matrix(#[1]*freemodule(#[2]))); }
   if( size(#)==1 )
   {
      int i; ideal id=#[1];
      int n=ncols(id); matrix A[n][n];
      for( i=1; i<=n; i=i+1 ) { A[i,i]=id[i]; }
   }
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),ds;
   print(diag(xy,4));
   matrix A[3][2] = 1,2,3,4,5,6;
   print(A);
   print(diag(A));
}
///////////////////////////////////////////////////////////////////////////////

proc dsum (list #)
"USAGE:   dsum(A1,A2,..); A1,A2,... matrices
RETURN:  matrix, direct sum of A1,A2,...
EXAMPLE: example dsum; shows an example
"
{
   int i,N,a;
   list L;
   for( i=1; i<=size(#); i=i+1 ) { N=N+nrows(#[i]); }
   for( i=1; i<=size(#); i=i+1 )
   {
      matrix B[N][ncols(#[i])];
      B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i];
      a=a+nrows(#[i]);
      L[i]=B; kill B;
   }
   return(concat(L));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),ds;
   matrix A[3][3] = 1,2,3,4,5,6,7,8,9;
   matrix B[2][2] = 1,x,y,z;
   print(A);
   print(B);
   print(dsum(A,B));
}
///////////////////////////////////////////////////////////////////////////////

proc flatten (matrix A)
"USAGE:   flatten(A); A matrix
RETURN:  ideal, generated by all entries from A
EXAMPLE: example flatten; shows an example
"
{
   return(ideal(A));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),ds;
   matrix A[2][3] = 1,2,x,y,z,7;
   print(A);
   flatten(A);
}
///////////////////////////////////////////////////////////////////////////////

proc genericmat (int n,int m,list #)
"USAGE:   genericmat(n,m[,id]);  n,m=integers, id=ideal
RETURN:  nxm matrix, with entries from id.
NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
         (default: id=maxideal(1)).
         genericmat(n,m); creates the generic nxm matrix
EXAMPLE: example genericmat; shows an example
"
{
   if( size(#)==0 ) { ideal id=maxideal(1); }
   if( size(#)==1 ) { ideal id=#[1]; }
   if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; }
   matrix B[n][m]=id;
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,x(1..16),lp;
   print(genericmat(3,3));      // the generic 3x3 matrix
   ring R1 = 0,(a,b,c,d),dp;
   matrix A = genericmat(3,4,maxideal(1)^3);
   print(A);
   int n,m = 3,2;
   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
   print(genericmat(n,m,i));    // matrix of generic linear forms
}
///////////////////////////////////////////////////////////////////////////////

proc is_complex (list c)
"USAGE:   is_complex(c); c = list of size-compatible modules or matrices
RETURN:  1 if c[i]*c[i+1]=0 for all i, 0 if not, hence checking whether the
         list of matrices forms a complex.
NOTE:    Ideals are treated internally as 1-line matrices.
         If printlevel > 0, the position where c is not a complex is shown.
EXAMPLE: example is_complex; shows an example
"
{
   int i;
   module @test;
   for( i=1; i<=size(c)-1; i=i+1 )
   {
      c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]);
      @test=c[i]*c[i+1];
      if (size(@test)!=0)
      {
        dbprint(printlevel-voice+2,"// not a complex at position " +string(i));
         return(0);
      }
   }
   return(1);
}
example
{ "EXAMPLE:";   echo = 2;
   ring r  = 32003,(x,y,z),ds;
   ideal i = x4+y5+z6,xyz,yx2+xz2+zy7;
   list L  = nres(i,0);
   is_complex(L);
   L[4]    = matrix(i);
   is_complex(L);
}
///////////////////////////////////////////////////////////////////////////////

proc outer (matrix A, matrix B)
"USAGE:   outer(A,B); A,B matrices
RETURN:  matrix, outer (tensor) product of A and B
EXAMPLE: example outer; shows an example
"
{
   int i,j; list L;
   int triv = nrows(B)*ncols(B);
   if( triv==1 )
   {
     return(B[1,1]*A);
   }
   else
   {
     int N = nrows(A)*nrows(B);
     matrix C[N][ncols(B)];
     for( i=1; i<=ncols(A); i=i+1 )
     {
       for( j=1; j<=nrows(A); j=j+1 )
       {
          C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B;
       }
       L[i]=C;
     }
     return(concat(L));
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),ds;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   matrix B[2][2]=x,y,0,z;
   print(A);
   print(B);
   print(outer(A,B));
}
///////////////////////////////////////////////////////////////////////////////

proc power ( A, int n)
"USAGE:   power(A,n);  A a square-matrix of type intmat or matrix, n=integer
RETURN:  intmat resp. matrix, the n-th power of A
NOTE:    for A=intmat and big n the result may be wrong because of int overflow
EXAMPLE: example power; shows an example
"
{
//---------------------------- type checking ----------------------------------
   if( typeof(A)!="matrix" and typeof(A)!="intmat" )
   {
      "// no matrix or intmat!";
      return (A);
   }
   if( ncols(A) != nrows(A) )
   {
      "// not a suare matrix!";
      return();
   }
//---------------------------- trivial cases ----------------------------------
   int ii;
   if( n <= 0 )
   {
      if( typeof(A)=="matrix" )
      {
         return (unitmat(nrows(A)));
      }
      if( typeof(A)=="intmat" )
      {
         intmat B[nrows(A)][nrows(A)];
         for( ii=1; ii<=nrows(A); ii++ )
         {
            B[ii,ii] = 1;
         }
         return (B);
      }
   }
   if( n == 1 ) { return (A); }
//---------------------------- sub procedure ----------------------------------
   proc matpow (A, int n)
   {
      def B = A*A;
      int ii= 2;
      int jj= 4;
      while( jj <= n )
      {
         B=B*B;
         ii=jj;
         jj=2*jj;
      }
      return(B,n-ii);
   }
//----------------------------- main program ----------------------------------
   list L = matpow(A,n);
   def B  = L[1];
   ii     = L[2];
   while( ii>=2 )
   {
      L = matpow(A,ii);
      B = B*L[1];
      ii= L[2];
   }
   if( ii == 0) { return(B); }
   if( ii == 1) { return(A*B); }
}
example
{ "EXAMPLE:"; echo = 2;
   intmat A[3][3]=1,2,3,4,5,6,7,8,9;
   print(power(A,3));"";
   ring r=0,(x,y,z),dp;
   matrix B[3][3]=0,x,y,z,0,0,y,z,0;
   print(power(B,3));"";
}
///////////////////////////////////////////////////////////////////////////////

proc skewmat (int n, list #)
"USAGE:   skewmat(n[,id]);  n integer, id ideal
RETURN:  skew-symmetric nxn matrix, with entries from id
         (default: id=maxideal(1))
         skewmat(n); creates the generic skew-symmetric matrix
NOTE:    if id has less than n*(n-1)/2 elements, the matrix is
         filled with 0's,
EXAMPLE: example skewmat; shows an example
"
{
   matrix B[n][n];
   if( size(#)==0 ) { ideal id=maxideal(1); }
   else { ideal id=#[1]; }
   id = id,B[1..n,1..n];
   int i,j;
   for( i=0; i<=n-2; i=i+1 )
   {
      B[i+1,i+2..n]=id[j+1..j+n-i-1];
      j=j+n-i-1;
   }
   matrix A=transpose(B);
   B=B-A;
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=0,x(1..5),lp;
   print(skewmat(4));    // the generic skew-symmetric matrix
   ring R1 = 0,(a,b,c),dp;
   matrix A=skewmat(4,maxideal(1)^2);
   print(A);
   int n=3;
   ideal i = ideal(randommat(1,n*(n-1) div 2,maxideal(1),9));
   print(skewmat(n,i));  // skew matrix of generic linear forms
   kill R1;
}
///////////////////////////////////////////////////////////////////////////////

proc submat (matrix A, intvec r, intvec c)
"USAGE:   submat(A,r,c);  A=matrix, r,c=intvec
RETURN:  matrix, submatrix of A with rows specified by intvec r
         and columns specified by intvec c.
EXAMPLE: example submat; shows an example
"
{
   matrix B[size(r)][size(c)]=A[r,c];
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=32003,(x,y,z),lp;
   matrix A[4][4]=x,y,z,0,1,2,3,4,5,6,7,8,9,x2,y2,z2;
   print(A);
   intvec v=1,3,4;
   matrix B=submat(A,v,1..3);
   print(B);
}
///////////////////////////////////////////////////////////////////////////////

proc symmat (int n, list #)
"USAGE:   symmat(n[,id]);  n integer, id ideal
RETURN:  symmetric nxn matrix, with entries from id (default: id=maxideal(1))
NOTE:    if id has less than n*(n+1)/2 elements, the matrix is filled with 0's,
         symmat(n); creates the generic symmetric matrix
EXAMPLE: example symmat; shows an example
"
{
   matrix B[n][n];
   if( size(#)==0 ) { ideal id=maxideal(1); }
   else { ideal id=#[1]; }
   id = id,B[1..n,1..n];
   int i,j;
   for( i=0; i<=n-1; i=i+1 )
   {
      B[i+1,i+1..n]=id[j+1..j+n-i];
      j=j+n-i;
   }
   matrix A=transpose(B);
   for( i=1; i<=n; i=i+1 ) {  A[i,i]=0; }
   B=A+B;
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=0,x(1..10),lp;
   print(symmat(4));    // the generic symmetric matrix
   ring R1 = 0,(a,b,c),dp;
   matrix A=symmat(4,maxideal(1)^3);
   print(A);
   int n=3;
   ideal i = ideal(randommat(1,n*(n+1) div 2,maxideal(1),9));
   print(symmat(n,i));  // symmetric matrix of generic linear forms
   kill R1;
}
///////////////////////////////////////////////////////////////////////////////

proc tensor (matrix A, matrix B)
"USAGE:   tensor(A,B); A,B matrices
RETURN:  matrix, tensor product of A and B
EXAMPLE: example tensor; shows an example
"
{
   if (ncols(A)==0)
   {
     int q=nrows(A)*nrows(B);
     matrix D[q][0];
     return(D);
   }
 
   int i,j;
   matrix C,D;
   for( i=1; i<=nrows(A); i++ )
   {
     C = A[i,1]*B;
     for( j=2; j<=ncols(A); j++ )
     {
       C = concat(C,A[i,j]*B);
     }
     D = concat(D,transpose(C));
   }
   D = transpose(D);
   return(submat(D,2..nrows(D),1..ncols(D)));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),(c,ds);
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   matrix B[2][2]=x,y,0,z;
   print(A);
   print(B);
   print(tensor(A,B));
}
///////////////////////////////////////////////////////////////////////////////

proc unitmat (int n)
"USAGE:   unitmat(n);  n integer >= 0
RETURN:  nxn unit matrix
NOTE:    needs a basering, diagonal entries are numbers (=1) in the basering
EXAMPLE: example unitmat; shows an example
"
{
   return(matrix(freemodule(n)));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   print(xyz*unitmat(4));
   print(unitmat(5));
}
///////////////////////////////////////////////////////////////////////////////

proc gauss_col (matrix A, list #)
"USAGE:   gauss_col(A[,e]); A a matrix, e any type
RETURN:  - a matrix B, if called with one argument; B is the complete column-
           reduced upper-triangular normal form of A if A is constant,
           (resp. as far as this is possible if A is a polynomial matrix;
           no division by polynomials).
@*       - a list L of two matrices, if called with two arguments;
           L satisfies L[1] = A * L[2] with L[1] the column-reduced form of A
           and L[2] the transformation matrix.
NOTE:    * The procedure just applies interred to A with ordering (C,dp).
           The transformation matrix is obtained by applying 'lift'.
           This should be faster than the procedure colred.
@*       * It should only be used with exact coefficient field (there is no
           pivoting and rounding error treatment).
@*       * Parameters are allowed. Hence, if the entries of A are parameters,
           B is the column-reduced form of A over the rational function field.
SEE ALSO:  colred
EXAMPLE: example gauss_col; shows an example
"
{
   def R=basering; int u;
   string mp = string(minpoly);
   int n = nrows(A);
   int m = ncols(A);
   module M = A;
   intvec v = option(get);
//------------------------ change ordering if necessary ----------------------
   if( ordstr(R) != ("C,dp("+string(nvars(R))+")") )
   {
     changeord("@R","C,dp",R); u=1;
     execute("minpoly="+mp+";");
     matrix A = imap(R,A);
     module M = A;
   }
//------------------------------ start computation ---------------------------
   option(redSB);
   M = simplify(interred(M),1);
   if(size(#) != 0)
   {
      module N = lift(A,M);
   }
//--------------- reset ring and options and return --------------------------
   if ( u==1 )
   {
      setring R;
      M=imap(@R,M);
      if (size(#) != 0)
      {
         module N = imap(@R,N);
      }
      kill @R;
   }
   option(set,v);
   // M = sort(M,size(M)..1)[1];
   A = matrix(M,n,m);
   if (size(#) != 0)
   {
     list L= A,matrix(N,m,m);
     return(L);
   }
   return(matrix(M,n,m));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[8][6]=
   0,    2*C, 0,    0,  0,   0,
   0,    -4*C,a*A,  0,  0,   0,
   b*B,  -A,  0,    0,  0,   0,
   -A,   B,   0,    0,  0,   0,
   -4*C, 0,   B,    2,  0,   0,
   2*A,  B,   0,    0,  0,   0,
   0,    3*B, 0,    0,  2b,  0,
   0,    AB,  0,    2*A,A,   2a;"";
   list L=gauss_col(m,1);
   print(L[1]);
   print(L[2]);

   ring S=0,x,(c,dp);
   matrix A[5][4] =
    3, 1, 1, 1,
   13, 8, 6,-7,
   14,10, 6,-7,
    7, 4, 3,-3,
    2, 1, 0, 3;
   print(gauss_col(A));
}
///////////////////////////////////////////////////////////////////////////////

proc gauss_row (matrix A, list #)
"USAGE:  gauss_row(A [,e]); A matrix, e any type
RETURN: - a matrix B, if called with one argument; B is the complete row-
          reduced lower-triangular normal form of A if A is constant,
          (resp. as far as this is possible if A is a polynomial matrix;
          no division by polynomials).
@*      - a list L of two matrices, if called with two arguments;
          L satisfies L[1] = L[2] * A with L[1] the row-reduced form of A
          and L[2] the transformation matrix.
NOTE:   * This procedure just applies gauss_col to the transposed matrix.
          The transformation matrix is obtained by applying lift.
          This should be faster than the procedure rowred.
@*      * It should only be used with exact coefficient field (there is no
          pivoting and rounding error treatment).
@*      * Parameters are allowed. Hence, if the entries of A are parameters,
          B is the row-reduced form of A over the rational function field.
SEE ALSO: rowred
EXAMPLE: example gauss_row; shows an example
"
{
   if(size(#) > 0)
   {
     list L = gauss_col(transpose(A),1);
     return(L);
   }
   A = gauss_col(transpose(A));
   return(transpose(A));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[6][8]=
   0, 0,  b*B, -A,-4C,2A,0, 0,
   2C,-4C,-A,B, 0,  B, 3B,AB,
   0,a*A,  0, 0, B,  0, 0, 0,
   0, 0,  0, 0, 2,  0, 0, 2A,
   0, 0,  0, 0, 0,  0, 2b, A,
   0, 0,  0, 0, 0,  0, 0, 2a;"";
   print(gauss_row(m));"";
   ring S=0,x,dp;
   matrix A[4][5] =  3, 1,1,-1,2,
                    13, 8,6,-7,1,
                    14,10,6,-7,1,
                     7, 4,3,-3,3;
   list L=gauss_row(A,1);
   print(L[1]);
   print(L[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc addcol (matrix A, int c1, poly p, int c2)
"USAGE:   addcol(A,c1,p,c2);  A matrix, p poly, c1, c2 positive integers
RETURN:  matrix,  A being modified by adding p times column c1 to column c2
EXAMPLE: example addcol; shows an example
"
{
   A[1..nrows(A),c2]=A[1..nrows(A),c2]+p*A[1..nrows(A),c1];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(addcol(A,1,xy,2));
}
///////////////////////////////////////////////////////////////////////////////

proc addrow (matrix A, int r1, poly p, int r2)
"USAGE:   addcol(A,r1,p,r2);  A matrix, p poly, r1, r2 positive integers
RETURN:  matrix,  A being modified by adding p times row r1 to row r2
EXAMPLE: example addrow; shows an example
"
{
   A[r2,1..ncols(A)]=A[r2,1..ncols(A)]+p*A[r1,1..ncols(A)];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(addrow(A,1,xy,3));
}
///////////////////////////////////////////////////////////////////////////////

proc multcol (matrix A, int c, poly p)
"USAGE:   addcol(A,c,p);  A matrix, p poly, c positive integer
RETURN:  matrix,  A being modified by multiplying column c with p
EXAMPLE: example multcol; shows an example
"
{
   A[1..nrows(A),c]=p*A[1..nrows(A),c];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(multcol(A,2,xy));
}
///////////////////////////////////////////////////////////////////////////////

proc multrow (matrix A, int r, poly p)
"USAGE:   multrow(A,r,p);  A matrix, p poly, r positive integer
RETURN:  matrix,  A being modified by multiplying row r with p
EXAMPLE: example multrow; shows an example
"
{
   A[r,1..ncols(A)]=p*A[r,1..ncols(A)];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(multrow(A,2,xy));
}
///////////////////////////////////////////////////////////////////////////////

proc permcol (matrix A, int c1, int c2)
"USAGE:   permcol(A,c1,c2);  A matrix, c1,c2 positive integers
RETURN:  matrix,  A being modified by permuting column c1 and c2
EXAMPLE: example permcol; shows an example
"
{
   matrix B=A;
   B[1..nrows(B),c1]=A[1..nrows(A),c2];
   B[1..nrows(B),c2]=A[1..nrows(A),c1];
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,x,3,4,y,6,7,z,9;
   print(A);
   print(permcol(A,2,3));
}
///////////////////////////////////////////////////////////////////////////////

proc permrow (matrix A, int r1, int r2)
"USAGE:   permrow(A,r1,r2);  A matrix, r1,r2 positive integers
RETURN:  matrix,  A being modified by permuting row r1 and r2
EXAMPLE: example permrow; shows an example
"
{
   matrix B=A;
   B[r1,1..ncols(B)]=A[r2,1..ncols(A)];
   B[r2,1..ncols(B)]=A[r1,1..ncols(A)];
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
   print(A);
   print(permrow(A,2,1));
}
///////////////////////////////////////////////////////////////////////////////

proc rowred (matrix A,list #)
"USAGE:   rowred(A[,e]);  A matrix, e any type
RETURN:  - a matrix B, being the row reduced form of A, if rowred is called
           with one argument.
           (as far as this is possible over the polynomial ring; no division
           by polynomials)
@*       - a list L of two matrices, such that L[1] = L[2] * A with L[1]
           the row-reduced form of A and L[2] the transformation matrix
           (if rowred is called with two arguments).
NOTE:    * This procedure is designed for teaching purposes mainly.
@*       * The straight forward Gaussian algorithm is implemented in the
           library (no standard basis computation).
           The transformation matrix is obtained by concatenating a unit
           matrix to A. proc gauss_row should be faster.
@*       * It should only be used with exact coefficient field (there is no
           pivoting) over the polynomial ring (ordering lp or dp).
@*       * Parameters are allowed. Hence, if the entries of A are parameters
           the computation takes place over the field of rational functions.
SEE ALSO:  gauss_row
EXAMPLE: example rowred; shows an example
"
{
   int m,n=nrows(A),ncols(A);
   int i,j,k,l,rk;
   poly p;
   matrix d[m][n];
   for (i=1;i<=m;i=i+1)
   {  for (j=1;j<=n;j=j+1)
      {  p = A[i,j];
         if (ord(p)==0)
         {  if (deg(p)==0) { d[i,j]=p; }
         }
      }
   }
   matrix b = A;
   if (size(#) != 0) { b = concat(b,unitmat(m)); }
      for (l=1;l<=n;l=l+1)
   {
      k  = findfirst(ideal(d[l]),rk+1);
      if (k)
      {  rk = rk+1;
         b  = permrow(b,rk,k);
         p  = b[rk,l];         p = 1/p;
         b  = multrow(b,rk,p);
         for (i=1;i<=m;i=i+1)
         {
            if (rk-i) { b = addrow(b,rk,-b[i,l],i);}
         }
         d = 0;
         for (i=rk+1;i<=m;i=i+1)
         {  for (j=l+1;j<=n;j=j+1)
            {  p = b[i,j];
               if (ord(p)==0)
               {  if (deg(p)==0) { d[i,j]=p; }
               }
            }
         }

      }
   }
   d = submat(b,1..m,1..n);
   if (size(#))
   {
      list L=d,submat(b,1..m,n+1..n+m);
      return(L);
   }
   return(d);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[6][8]=
   0, 0,  b*B, -A,-4C,2A,0, 0,
   2C,-4C,-A,B, 0,  B, 3B,AB,
   0,a*A,  0, 0, B,  0, 0, 0,
   0, 0,  0, 0, 2,  0, 0, 2A,
   0, 0,  0, 0, 0,  0, 2b, A,
   0, 0,  0, 0, 0,  0, 0, 2a;"";
   print(rowred(m));"";
   list L=rowred(m,1);
   print(L[1]);
   print(L[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc colred (matrix A,list #)
"USAGE:   colred(A[,e]);  A matrix, e any type
RETURN:  - a matrix B, being the column reduced form of A, if colred is
           called with one argument.
           (as far as this is possible over the polynomial ring;
           no division by polynomials)
@*       - a list L of two matrices, such that L[1] = A * L[2] with L[1]
           the column-reduced form of A and L[2] the transformation matrix
           (if colred is called with two arguments).
NOTE:    * This procedure is designed for teaching purposes mainly.
@*       * It applies rowred to the transposed matrix.
           proc gauss_col should be faster.
@*       * It should only be used with exact coefficient field (there is no
           pivoting) over the polynomial ring (ordering lp or dp).
@*       * Parameters are allowed. Hence, if the entries of A are parameters
           the computation takes place over the field of rational functions.
SEE ALSO:  gauss_col
EXAMPLE: example colred; shows an example
"
{
   A = transpose(A);
   if (size(#))
   { list L = rowred(A,1); return(transpose(L[1]),transpose(L[2]));}
   else
   { return(transpose(rowred(A)));}
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[8][6]=
   0,    2*C, 0,    0,  0,   0,
   0,    -4*C,a*A,  0,  0,   0,
   b*B,  -A,  0,    0,  0,   0,
   -A,   B,   0,    0,  0,   0,
   -4*C, 0,   B,    2,  0,   0,
   2*A,  B,   0,    0,  0,   0,
   0,    3*B, 0,    0,  2b,  0,
   0,    AB,  0,    2*A,A,   2a;"";
   print(colred(m));"";
   list L=colred(m,1);
   print(L[1]);
   print(L[2]);
}
//////////////////////////////////////////////////////////////////////////////

static proc findfirst (ideal i,int t)
{
   int n,k;
   for (n=t;n<=ncols(i);n=n+1)
   {
      if (i[n]!=0) { k=n;break;}
   }
   return(k);
}
//////////////////////////////////////////////////////////////////////////////

proc rm_unitcol(matrix A)
"USAGE:   rm_unitcol(A); A matrix (being row-reduced)
RETURN:  matrix, obtained from A by deleting unit columns (having just one 1
         and else 0 as entries) and associated rows
EXAMPLE: example rm_unitcol; shows an example
"
{
 int l,j;
 intvec v;
 for (j=1;j<=ncols(A);j=j+1)
 {
    if (gen(l+1)==module(A)[j]) {l=l+1;}
    else { v=v,j;}
 }
 if (size(v)>1)
    {  v = v[2..size(v)];
       return(submat(A,l+1..nrows(A),v));
    }
 else
    { return(0);}
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[6][8]=
   0,  0,    A,   0, 1,0,  0,0,
   0,  0,  -C2,   0, 0,0,  1,0,
   0,  0,    0,1/2B, 0,0,  0,1,
   0,  0,    B,  -A, 0,2A, 0,0,
   2C,-4C,  -A,   B, 0,B,  0,0,
   0,  A,    0,   0, 0,0,  0,0;
   print(rm_unitcol(m));
}
//////////////////////////////////////////////////////////////////////////////

proc rm_unitrow (matrix A)
"USAGE:   rm_unitrow(A); A matrix (being col-reduced)
RETURN:  matrix, obtained from A by deleting unit rows (having just one 1
         and else 0 as entries) and associated columns
EXAMPLE: example rm_unitrow; shows an example
"
{
 int l,j;
 intvec v;
 module M = transpose(A);
 for (j=1; j <= nrows(A); j=j+1)
 {
    if (gen(l+1) == M[j]) { l=l+1; }
    else { v=v,j; }
 }
 if (size(v) > 1)
    {  v = v[2..size(v)];
       return(submat(A,v,l+1..ncols(A)));
    }
 else
    { return(0);}
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[8][6]=
   0,0,  0,   0, 2C, 0,
   0,0,  0,   0, -4C,A,
   A,-C2,0,   B, -A, 0,
   0,0,  1/2B,-A,B,  0,
   1,0,  0,   0, 0,  0,
   0,0,  0,   2A,B,  0,
   0,1,  0,   0, 0,  0,
   0,0,  1,   0, 0,  0;
   print(rm_unitrow(m));
}
//////////////////////////////////////////////////////////////////////////////