# source:git/Singular/LIB/matrix.lib@18dd47

spielwiese
Last change on this file since 18dd47 was 18dd47, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: changed int / into div in HNPuiseux.lib finvar.lib general.lib homolog.lib matrix.lib poly.lib prim_dec.lib primdec.lib random.lib git-svn-id: file:///usr/local/Singular/svn/trunk@612 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 16.3 KB
Line
1// \$Id: matrix.lib,v 1.3 1997-08-12 14:01:08 Singular Exp \$
3///////////////////////////////////////////////////////////////////////////////
4LIBRARY:  matrix.lib    PROCEDURES FOR MATRIX OPERATIONS
5
6 compress(A);           matrix, zero columns from A deleted
7 concat(A1,A2,..);      matrix, concatenation of matrices A1,A2,...
8 diag(p,n);             matrix, nxn diagonal matrix with entries poly p
9 dsum(A1,A2,..);        matrix, direct sum of matrices A1,A2,...
10 flatten(A);            ideal, generated by entries of matrix A
11 genericmat(n,m[,id]);  generic nxm matrix [entries from id]
12 is_complex(c);         1 if list c is a complex, 0 if not
13 outer(A,B);            matrix, outer product of matrices A and B
14 skewmat(n[,id]);       generic skew-symmetric nxn matrix [entries from id]
15 submat(A,r,c);         submatrix of A with rows/cols specified by intvec r/c
16 symmat(n[,id]);        generic symmetric nxn matrix [entries from id]
17 tensor(A,B);           matrix, tensor product of matrices A nd B
18 unitmat(n);            unit square matrix of size n
19 gauss_col(A);          transform constant matrix A into col-reduced nf
20 gauss_row(A);          transform constant matrix A into row-reduced nf
21 addcol(A,c1,p,c2);     add p*(c1-th col) to c2-th column of matrix A, p poly
22 addrow(A,r1,p,r2);     add p*(r1-th row) to r2-th row of matrix A, p poly
23 multcol(A,c,p);        multiply c-th column of A with poly p
24 multrow(A,r,p);        multiply r-th row of A with poly p
25 permcol(A,i,j);        permute i-th and j-th columns
26 permrow(A,i,j);        permute i-th and j-th rows
27           (parameters in square brackets [] are optional)
28
29LIB "inout.lib";
30LIB "ring.lib";
31LIB "random.lib";
32///////////////////////////////////////////////////////////////////////////////
33
34proc compress (A)
35USAGE:   compress(A); A matrix/ideal/module/intmat/intvec
36RETURN:  same type, zero columns/generators from A deleted
37         (in an intvec zero elements are deleted)
38EXAMPLE: example compress; shows an example
39{
40   if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); }
41   if( typeof(A)=="intmat" or typeof(A)=="intvec" )
42   {
43      ring r=0,x,lp;
44      if( typeof(A)=="intvec" ) { intmat C=transpose(A); kill A; intmat A=C; }
45      module m = matrix(A);
46      intmat B[nrows(A)][size(m)];
47      int i,j;
48      for( i=1; i<=ncols(A); i=i+1 )
49      {
50         if( m[i]!=[0] )
51         {
52            j=j+1;
53            B[1..nrows(A),j]=A[1..nrows(A),i];
54         }
55      }
56      if( defined(C) ) { return(intvec(B)); }
57      return(B);
58    }
59   return(simplify(A,2));
60}
61example
62{ "EXAMPLE:"; echo = 2;
63   ring r=0,(x,y,z),ds;
64   matrix A[3][4]=1,0,3,0,x,0,z,0,x2,0,z2,0;
65   print(A);
66   print(compress(A));
67   module m=module(A); show(m);
68   show(compress(m));
69   intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0;
70   compress(B);
71   intvec C=0,0,1,2,0,3;
72   compress(C);
73}
74////////////////////////////////////////////////////////////////////////////////
75proc concat (list #)
76USAGE:   concat(A1,A2,..); A1,A2,... matrices
77RETURN:  matrix, concatenation of A1,A2,... . Number of rows of result matrix is
78         max(nrows(A1),nrows(A2),...)
79EXAMPLE: example concat; shows an example
80{
81   int i;
82   module B=module(#[1]);
83   for( i=2; i<=size(#); i=i+1 ) { B=B,module(#[i]); }
84   return(matrix(B));
85}
86example
87{ "EXAMPLE:"; echo = 2;
88   ring r=0,(x,y,z),ds;
89   matrix A[3][3]=1,2,3,x,y,z,x2,y2,z2;
90   matrix B[2][2]=1,0,2,0; matrix C[1][4]=4,5,x,y;
91   print(A);
92   print(B);
93   print(C);
94   print(concat(A,B,C));
95}
96////////////////////////////////////////////////////////////////////////////////
97
98proc diag (list #)
99USAGE:   diag(p,n); p poly, n integer
100         diag(A);   A matrix
101RETURN:  diag(p,n): diagonal matrix, p times unitmatrix of size n
102         diag(A)  : n*mxn*m diagonal matrix with entries all the entries of the
103                    nxm matrix A, taken from the 1st row, 2nd row etc of A
104EXAMPLE: example diag; shows an example
105{
106   if( size(#)==2 ) { return(matrix(#[1]*freemodule(#[2]))); }
107   if( size(#)==1 )
108   {
109      int i; ideal id=#[1];
110      int n=ncols(id); matrix A[n][n];
111      for( i=1; i<=n; i=i+1 ) { A[i,i]=id[i]; }
112   }
113   return(A);
114}
115example
116{ "EXAMPLE:"; echo = 2;
117   ring r=0,(x,y,z),ds;
118   print(diag(xy,4));
119   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
120   print(A);
121   print(diag(A));
122}
123////////////////////////////////////////////////////////////////////////////////
124
125proc dsum (list #)
126USAGE:   dsum(A1,A2,..); A1,A2,... matrices
127RETURN:  matrix, direct sum of A1,A2,...
128EXAMPLE: example dsum; shows an example
129{
130   int i,N,a;
131   list L;
132   for( i=1; i<=size(#); i=i+1 ) { N=N+nrows(#[i]); }
133   for( i=1; i<=size(#); i=i+1 )
134   {
135      matrix B[N][ncols(#[i])];
136      B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i];
137      a=a+nrows(#[i]);
138      L[i]=B; kill B;
139   }
140   return(concat(L));
141}
142example
143{ "EXAMPLE:"; echo = 2;
144   ring r=0,(x,y,z),ds;
145   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
146   matrix B[2][2]=1,x,y,z;
147   matrix C[1][4]=4,5,x,y;
148   print(A);
149   print(B);
150   print(C);
151   print(dsum(A,B,C));
152}
153////////////////////////////////////////////////////////////////////////////////
154
155proc flatten (matrix A)
156USAGE:   flatten(A); A matrix
157RETURN:  ideal, generated by all entries from A
158EXAMPLE: example flatten; shows an example
159{
160   return(ideal(A));
161}
162example
163{ "EXAMPLE:"; echo = 2;
164   ring r=0,(x,y,z),ds;
165   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
166   print(A);
167   flatten(A);
168}
169////////////////////////////////////////////////////////////////////////////////
170
171proc genericmat (int n,int m,list #)
172USAGE:   genericmat(n,m[,id]);  n,m=integers, id=ideal
173RETURN:  nxm matrix, with entries from id (default: id=maxideal(1))
174NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
175         genericmat(n,m); creates the generic nxm matrix
176EXAMPLE: example genericmat; shows an example
177{
178   if( size(#)==0 ) { ideal id=maxideal(1); }
179   if( size(#)==1 ) { ideal id=#[1]; }
180   if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; }
181   matrix B[n][m]=id;
182   return(B);
183}
184example
185{ "EXAMPLE:"; echo = 2;
186   ring R=0,x(1..16),lp;
187   print(genericmat(4,4));    // the generic 4x4 matrix
188   changevar("R1",A_Z("a",4),R);
189   matrix A=genericmat(4,5,maxideal(1)^3);
190   print(A);
191   int n,m=4,3;
192   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
193   print(genericmat(n,m,i));  // matrix of generic linear forms
194   kill R1;
195}
196///////////////////////////////////////////////////////////////////////////////
197
198proc is_complex (list c)
199USAGE:   is_complex(c); c = list of size-compatible modules or matrices
200RETURN:  1 if c[i]*c[i+1]=0 for all i, 0 if not.
201NOTE:    Ideals are treated internally as 1-line matrices
202EXAMPLE: example is_complex; shows an example
203{
204   int i;
205   module @test;
206   for( i=1; i<=size(c)-1; i=i+1 )
207   {
208      c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]);
209      @test=c[i]*c[i+1];
210      if (size(@test)!=0)
211      {
212         if( voice==2 ) { "// argument is not a complex at position",i; }
213         return(0);
214      }
215   }
216   if( voice==2 ) { "// argument is a complex"; }
217   return(1);
218}
219example
220{ "EXAMPLE:";   echo = 2;
221   ring r=32003,(x,y,z),ds;
222   ideal i=x4+y5+z6,xyz,yx2+xz2+zy7;
223   list L=res(i,0);
224   is_complex(L);
225   L[4]=matrix(i);
226   is_complex(L);
227}
228////////////////////////////////////////////////////////////////////////////////
229
230proc outer (matrix A, matrix B)
231USAGE:   outer(A,B); A,B matrices
232RETURN:  matrix, outer product of A and B
233EXAMPLE: example outer; shows an example
234{
235   int i,j; list L;
236   int triv = nrows(B)*ncols(B);
237   if( triv==1 )
238   {
239     return(B[1,1]*A);
240   }
241   else
242   {
243     int N = nrows(A)*nrows(B);
244     matrix C[N][ncols(B)];
245     for( i=1; i<=ncols(A); i=i+1 )
246     {
247       for( j=1; j<=nrows(A); j=j+1 )
248       {
249          C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B;
250       }
251       L[i]=C;
252     }
253     return(concat(L));
254   }
255}
256example
257{ "EXAMPLE:"; echo = 2;
258   ring r=32003,(x,y,z),ds;
259   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
260   matrix B[2][2]=x,y,0,z;
261   print(A);
262   print(B);
263   print(outer(A,B));
264}
265////////////////////////////////////////////////////////////////////////////////
266
267proc skewmat (int n, list #)
268USAGE:   skewmat(n[,id]);  n integer, id ideal
269RETURN:  skew-symmetric nxn matrix, with entries from id
270         (default: id=maxideal(1))
271NOTE:    if id has less than n*(n-1)/2 elements, the matrix is filled with 0's,
272         skewmat(n); creates the generic skew-symmetric matrix
273EXAMPLE: example skewmat; shows an example
274{
275   matrix B[n][n];
276   if( size(#)==0 ) { ideal id=maxideal(1); }
277   else { ideal id=#[1]; }
278   id = id,B[1..n,1..n];
279   int i,j;
280   for( i=0; i<=n-2; i=i+1 )
281   {
282      B[i+1,i+2..n]=id[j+1..j+n-i-1];
283      j=j+n-i-1;
284   }
285   matrix A=transpose(B);
286   B=B-A;
287   return(B);
288}
289example
290{ "EXAMPLE:"; echo = 2;
291   ring R=0,x(1..5),lp;
292   print(skewmat(4));    // the generic skew-symmetric matrix
293   changevar("R1",A_Z("a",5),R);
294   matrix A=skewmat(6,maxideal(1)^2);
295   print(A);
296   int n=4;
297   ideal i = ideal(randommat(1,n*(n-1) div 2,maxideal(1),9));
298   print(skewmat(n,i));  // skew matrix of generic linear forms
299   kill R1;
300}
301////////////////////////////////////////////////////////////////////////////////
302
303proc submat (matrix A, intvec r, intvec c)
304USAGE:   submat(A,r,c);  A=matrix, r,c=intvec
305RETURN:  matrix, submatrix of A with rows specified by intvec r and columns
306         specified by intvec c
307EXAMPLE: example submat; shows an example
308{
309   matrix B[size(r)][size(c)]=A[r,c];
310   return(B);
311}
312example
313{ "EXAMPLE:"; echo = 2;
314   ring R=32003,(x,y,z),lp;
315   matrix A[4][4]=x,y,z,0,1,2,3,4,5,6,7,8,9,x2,y2,z2;
316   print(A);
317   intvec v=1,3,4;
318   matrix B=submat(A,v,1..3);
319   print(B);
320}
321////////////////////////////////////////////////////////////////////////////////
322
323proc symmat (int n, list #)
324USAGE:   symmat(n[,id]);  n integer, id ideal
325RETURN:  symmetric nxn matrix, with entries from id (default: id=maxideal(1))
326NOTE:    if id has less than n*(n+1)/2 elements, the matrix is filled with 0's,
327         symmat(n); creates the generic symmetric matrix
328EXAMPLE: example symmat; shows an example
329{
330   matrix B[n][n];
331   if( size(#)==0 ) { ideal id=maxideal(1); }
332   else { ideal id=#[1]; }
333   id = id,B[1..n,1..n];
334   int i,j;
335   for( i=0; i<=n-1; i=i+1 )
336   {
337      B[i+1,i+1..n]=id[j+1..j+n-i];
338      j=j+n-i;
339   }
340   matrix A=transpose(B);
341   for( i=1; i<=n; i=i+1 ) {  A[i,i]=0; }
342   B=A+B;
343   return(B);
344}
345example
346{ "EXAMPLE:"; echo = 2;
347   ring R=0,x(1..10),lp;
348   print(symmat(4));    // the generic symmetric matrix
349   changevar("R1",A_Z("a",5),R);
350   matrix A=symmat(5,maxideal(1)^2);
351   print(A);
352   int n=3;
353   ideal i = ideal(randommat(1,n*(n+1) div 2,maxideal(1),9));
354   print(symmat(n,i));  // symmetric matrix of generic linear forms
355   kill R1;
356}
357////////////////////////////////////////////////////////////////////////////////
358
359proc tensor (matrix A, matrix B)
360USAGE:   tensor(A,B); A,B matrices
361RETURN:  matrix, tensor product of A and B
362EXAMPLE: example tensor; shows an example
363{
364   int i,j;
365   matrix C=B;
366   for( i=2; i<=nrows(A); i=i+1 ) { C=dsum(C,B); }
367   matrix D[nrows(C)][ncols(A)*nrows(B)];
368   for( j=1; j<=nrows(B); j=j+1 )
369   {
370      for( i=1; i<=nrows(A); i=i+1 )
371      {
372         D[(i-1)*nrows(B)+j,(j-1)*ncols(A)+1..j*ncols(A)]=A[i,1..ncols(A)];
373      }
374   }
375   return(concat(C,D));
376}
377example
378{ "EXAMPLE:"; echo = 2;
379   ring r=32003,(x,y,z),(c,ds);
380   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
381   matrix B[2][2]=x,y,0,z;
382   print(A);
383   print(B);
384   print(tensor(A,B));
385}
386////////////////////////////////////////////////////////////////////////////////
387
388proc unitmat (int n)
389USAGE:   unitmat(n);  n integer >= 0
390RETURN:  nxn unit matrix
391NOTE:    needs a basering, diagonal entries are numbers (=1) in the basering
392EXAMPLE: example unitmat; shows an example
393{
394   return(matrix(freemodule(n)));
395}
396example
397{ "EXAMPLE:"; echo = 2;
398   ring r=32003,(x,y,z),lp;
399   print(xyz*unitmat(4));
400   print(unitmat(5));
401}
402///////////////////////////////////////////////////////////////////////////////
403
404proc gauss_col (matrix A)
405USAGE:   gauss_col(A); A=matrix with constant coefficients
406RETURN:  matrix = col-reduced lower-triagonal normal form of A
407NOTE:    the procedure sets the global option-command: option(noredSB);
408EXAMPLE: example gauss_col; shows an example
409{
410   def R=basering;
411   changeord("@R","ds,c",R);
412   option(redSB); option(nointStrategy);
413   matrix A = imap(R,A);
414   A = matrix(std(A),nrows(A),ncols(A));
415   setring R;
416   A=imap(@R,A);
417   option(noredSB);
418   kill @R;
419   return(A);
420}
421example
422{ "EXAMPLE:"; echo = 2;
423   ring S=0,x,dp;
424   matrix A[5][4] =  3, 1,1,-1,
425                    13, 8,6,-7,
426                    14,10,6,-7,
427                     7, 4,3,-3,
428                     2, 1,0, 3;
429   print(gauss_col(A));
430}
431///////////////////////////////////////////////////////////////////////////////
432
433proc gauss_row (matrix A)
434USAGE:   gauss_row(A); A=matrix with constant coefficients
435RETURN:  matrix = row-reduced upper-triangular normal form of A
436NOTE:    may be used to solve a system of linear equations
437         see proc 'linearsolve' from 'solve.lib' for a different method
438         the procedure sets the global option-command: option(noredSB);
439EXAMPLE: example gauss_row; shows an example
440{
441   A = gauss_col(transpose(A));
442   return(transpose(A));
443}
444example
445{ "EXAMPLE:"; echo = 2;
446   ring S=0,x,dp;
447   matrix A[4][5] =  3, 1,1,-1,2,
448                    13, 8,6,-7,1,
449                    14,10,6,-7,1,
450                     7, 4,3,-3,3;
451   print(gauss_row(A));
452}
453////////////////////////////////////////////////////////////////////////////////
454
455proc addcol (matrix A, int c1, poly p, int c2)
456USAGE:   addcol(A,c1,p,c2);  A matrix, p poly, c1, c2 positive integers
457RETURN:  matrix,  A being modified by adding p times column c1 to column c2
458EXAMPLE: example addcol; shows an example
459{
460   A[1..nrows(A),c2]=A[1..nrows(A),c2]+p*A[1..nrows(A),c1];
461   return(A);
462}
463example
464{ "EXAMPLE:"; echo = 2;
465   ring r=32003,(x,y,z),lp;
466   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
467   print(A);
469}
470////////////////////////////////////////////////////////////////////////////////
471
472proc addrow (matrix A, int r1, poly p, int r2)
473USAGE:   addcol(A,r1,p,r2);  A matrix, p poly, r1, r2 positive integers
474RETURN:  matrix,  A being modified by adding p times row r1 to row r2
475EXAMPLE: example addrow; shows an example
476{
477   A[r2,1..ncols(A)]=A[r2,1..ncols(A)]+p*A[r1,1..ncols(A)];
478   return(A);
479}
480example
481{ "EXAMPLE:"; echo = 2;
482   ring r=32003,(x,y,z),lp;
483   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
484   print(A);
486}
487////////////////////////////////////////////////////////////////////////////////
488
489proc multcol (matrix A, int c, poly p)
490USAGE:   addcol(A,c,p);  A matrix, p poly, c positive integer
491RETURN:  matrix,  A being modified by multiplying column c with p
492EXAMPLE: example multcol; shows an example
493{
494   A[1..nrows(A),c]=p*A[1..nrows(A),c];
495   return(A);
496}
497example
498{ "EXAMPLE:"; echo = 2;
499   ring r=32003,(x,y,z),lp;
500   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
501   print(A);
502   print(multcol(A,2,xy));
503}
504////////////////////////////////////////////////////////////////////////////////
505
506proc multrow (matrix A, int r, poly p)
507USAGE:   addcol(A,r,p);  A matrix, p poly, r positive integer
508RETURN:  matrix,  A being modified by multiplying row r with p
509EXAMPLE: example multrow; shows an example
510{
511   A[r,1..ncols(A)]=p*A[r,1..ncols(A)];
512   return(A);
513}
514example
515{ "EXAMPLE:"; echo = 2;
516   ring r=32003,(x,y,z),lp;
517   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
518   print(A);
519   print(multrow(A,2,xy));
520}
521////////////////////////////////////////////////////////////////////////////////
522
523proc permcol (matrix A, int c1, int c2)
524USAGE:   permcol(A,c1,c2);  A matrix, c1,c2 positive integers
525RETURN:  matrix,  A being modified by permuting column c1 and c2
526EXAMPLE: example permcol; shows an example
527{
528   matrix B=A;
529   B[1..nrows(B),c1]=A[1..nrows(A),c2];
530   B[1..nrows(B),c2]=A[1..nrows(A),c1];
531   return(B);
532}
533example
534{ "EXAMPLE:"; echo = 2;
535   ring r=32003,(x,y,z),lp;
536   matrix A[3][3]=1,x,3,4,y,6,7,z,9;
537   print(A);
538   print(permcol(A,2,3));
539}
540////////////////////////////////////////////////////////////////////////////////
541
542proc permrow (matrix A, int r1, int r2)
543USAGE:   permrow(A,r1,r2);  A matrix, r1,r2 positive integers
544RETURN:  matrix,  A being modified by permuting row r1 and r2
545EXAMPLE: example permrow; shows an example
546{
547   matrix B=A;
548   B[r1,1..ncols(B)]=A[r2,1..ncols(A)];
549   B[r2,1..ncols(B)]=A[r1,1..ncols(A)];
550   return(B);
551}
552example
553{ "EXAMPLE:"; echo = 2;
554   ring r=32003,(x,y,z),lp;
555   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
556   print(A);
557   print(permrow(A,2,1));
558}
559////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.