source: git/Singular/LIB/matrix.lib @ 5dc4ea

// $Id: matrix.lib,v 1.4 1997-09-18 09:58:24 Singular Exp $
// (GMG/BM, last modified 22.06.96)
///////////////////////////////////////////////////////////////////////////////
LIBRARY:  matrix.lib    PROCEDURES FOR MATRIX OPERATIONS

 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 constant matrix A into col-reduced nf
 gauss_row(A);          transform constant matrix A into row-reduced nf
 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
           (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
         (in an 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);
      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 unitmatrix of size n
         diag(A)  : n*mxn*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][3]=1,2,3,4,5,6,7,8,9;
   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;
   matrix C[1][4]=4,5,x,y;
   print(A);
   print(B);
   print(C);
   print(dsum(A,B,C));
}
////////////////////////////////////////////////////////////////////////////////

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[3][3]=1,2,3,x,y,z,7,8,9;
   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 (default: id=maxideal(1))
NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
         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(4,4));    // the generic 4x4 matrix
   changevar("R1",A_Z("a",4),R);
   matrix A=genericmat(4,5,maxideal(1)^3);
   print(A);
   int n,m=4,3;
   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
   print(genericmat(n,m,i));  // matrix of generic linear forms
   kill R1;
}
///////////////////////////////////////////////////////////////////////////////

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.
NOTE:    Ideals are treated internally as 1-line matrices
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)
      {
         if( voice==2 ) { "// argument is not a complex at position",i; }
         return(0);
      }
   }
   if( voice==2 ) { "// argument is a complex"; }
   return(1);
}
example
{ "EXAMPLE:";   echo = 2;
   ring r=32003,(x,y,z),ds;
   ideal i=x4+y5+z6,xyz,yx2+xz2+zy7;
   list L=res(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 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:  inmat resp. matrix, the n-th power of A
NOTE:    for intamt 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[4][4]=0,x,y,z,0,0,y,z,0,0,0,z,x,y,z,0;
   print(power(B,3));"";
   matrix C[3][3]=1,2,3,4,5,6,7,8,9;
   power(C,50);
}
////////////////////////////////////////////////////////////////////////////////

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))
NOTE:    if id has less than n*(n-1)/2 elements, the matrix is filled with 0's,
         skewmat(n); creates the generic skew-symmetric matrix
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
   changevar("R1",A_Z("a",5),R);
   matrix A=skewmat(6,maxideal(1)^2);
   print(A);
   int n=4;
   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
   changevar("R1",A_Z("a",5),R);
   matrix A=symmat(5,maxideal(1)^2);
   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
{
   int i,j;
   matrix C=B;
   for( i=2; i<=nrows(A); i=i+1 ) { C=dsum(C,B); }
   matrix D[nrows(C)][ncols(A)*nrows(B)];
   for( j=1; j<=nrows(B); j=j+1 )
   {
      for( i=1; i<=nrows(A); i=i+1 )
      {
         D[(i-1)*nrows(B)+j,(j-1)*ncols(A)+1..j*ncols(A)]=A[i,1..ncols(A)];
      }
   }
   return(concat(C,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)
USAGE:   gauss_col(A); A=matrix with constant coefficients
RETURN:  matrix = col-reduced lower-triagonal normal form of A
NOTE:    the procedure sets the global option-command: option(noredSB);
EXAMPLE: example gauss_col; shows an example
{
   def R=basering;
   changeord("@R","ds,c",R);
   option(redSB); option(nointStrategy);
   matrix A = imap(R,A);
   A = matrix(std(A),nrows(A),ncols(A));
   setring R;
   A=imap(@R,A);
   option(noredSB);
   kill @R;
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring S=0,x,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)
USAGE:   gauss_row(A); A=matrix with constant coefficients
RETURN:  matrix = row-reduced upper-triangular normal form of A
NOTE:    may be used to solve a system of linear equations
         see proc 'linearsolve' from 'solve.lib' for a different method
         the procedure sets the global option-command: option(noredSB);
EXAMPLE: example gauss_row; shows an example
{
   A = gauss_col(transpose(A));
   return(transpose(A));
}
example
{ "EXAMPLE:"; echo = 2;
   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;
   print(gauss_row(A));
}
////////////////////////////////////////////////////////////////////////////////

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:   addcol(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));
}
////////////////////////////////////////////////////////////////////////////////