source: git/gfanlib/gfanlib_matrix.h @ c4d065

spielwiese
Last change on this file since c4d065 was c4d065, checked in by Frank Seelisch <seelisch@…>, 13 years ago
coding at Goettingen (cones&fans) git-svn-id: file:///usr/local/Singular/svn/trunk@13677 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.1 KB
Line 
1/*
2 * lib_zmatrix.h
3 *
4 *  Created on: Sep 28, 2010
5 *      Author: anders
6 */
7
8#ifndef LIB_ZMATRIX_H_
9#define LIB_ZMATRIX_H_
10
11#include <vector>
12#include <algorithm>
13#include "gfanlib_vector.h"
14
15namespace gfan{
16
17template <class typ> class Matrix{
18  int width,height;
19  std::vector<Vector<typ> > rows;
20public:
21  inline int getHeight()const{return height;};
22  inline int getWidth()const{return width;};
23  Matrix(const Matrix &a):rows(a.rows),width(a.getWidth()),height(a.getHeight()){
24  }
25  Matrix(int height_, int width_):rows(height_),height(height_),width(width_){
26    assert(height>=0);
27    assert(width>=0);
28    for(int i=0;i<getHeight();i++)rows[i]=Vector<typ>(width);
29  };
30  ~Matrix(){
31  };
32  Matrix():width(0),height(0){
33  };
34  static Matrix rowVectorMatrix(Vector<typ> const &v)
35  {
36    Matrix ret(1,v.size());
37    for(int i=0;i<v.size();i++)ret[0][i]=v[i];
38    return ret;
39  }
40  Vector<typ> column(int i)const
41    {
42      assert(i>=0);
43      assert(i<getWidth());
44      Vector<typ> ret(getHeight());
45      for(int j=0;j<getHeight();j++)ret[j]=rows[j][i];
46      return ret;
47    }
48  Matrix transposed()const
49    {
50      Matrix ret(getWidth(),getHeight());
51      for(int i=0;i<getWidth();i++)
52        ret.rows[i]=column(i);
53      return ret;
54    }
55  static Matrix identity(int n)
56    {
57      Matrix m(n,n);
58      for(int i=0;i<n;i++)m.rows[i]=Vector<typ>::standardVector(n,i);
59      return m;
60    }
61  void append(Matrix const &m)
62    {
63      for(int i=0;i<m.height;i++)
64        {
65          rows.push_back(m[i]);
66        }
67      height+=m.height;
68    }
69  void appendRow(Vector<typ> const &v)
70    {
71      assert(v.size()==width);
72      rows.push_back(v);
73      height++;
74    }
75  void eraseLastRow()
76  {
77    assert(rows.size()>0);
78    rows.resize(rows.size()-1);
79    height--;
80  }
81  /*IntegerVector vectormultiply(IntegerVector const &v)const
82    {
83      assert(v.size()==width);
84      IntegerVector ret(height);
85      for(int i=0;i<height;i++)
86        ret[i]=dot(rows[i],v);
87      return ret;
88    }*/
89  /**
90   * Decides if v is in the kernel of the matrix.
91   */
92/*  bool inKernel(IntegerVector const &v)const
93    {
94      assert(v.size()==width);
95      for(int i=0;i<height;i++)
96          if(dotLong(rows[i],v)!=0)return false;
97      return true;
98    }
99*/
100  friend Matrix operator*(const typ &s, const Matrix& q)
101    {
102      Matrix p=q;
103      for(int i=0;i<q.height;i++)p[i]=s*(q[i]);
104      return p;
105    }
106/*  friend Matrix operator*(const Matrix& a, const Matrix& b)
107    {
108      assert(a.width==b.height);
109      Matrix ret(b.width,a.height);
110      for(int i=0;i<b.width;i++)
111        ret[i]=a.vectormultiply(b.column(i));
112      return ret.transposed();
113    }*/
114  /*  template<class T>
115    Matrix<T>(const Matrix<T>& c):v(c.size()){
116    for(int i=0;i<size();i++)v[i]=typ(c[i]);}
117  */
118  friend Matrix operator-(const Matrix &b)
119  {
120    Matrix ret(b.height,b.width);
121    for(int i=0;i<b.height;i++)ret[i]=-b[i];
122    return ret;
123  }
124
125  /**
126     Returns the specified submatrix. The endRow and endColumn are not included.
127   */
128  Matrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const
129  {
130    assert(startRow>=0);
131    assert(startColumn>=0);
132    assert(endRow>=startRow);
133    assert(endColumn>=startColumn);
134    assert(endRow<=height);
135    assert(endColumn<=width);
136    Matrix ret(endRow-startRow,endColumn-startColumn);
137    for(int i=startRow;i<endRow;i++)
138      for(int j=startColumn;j<endColumn;j++)
139        ret[i-startRow][j-startColumn]=rows[i][j];
140    return ret;
141  }
142  const Vector<typ>& operator[](int n)const{assert(n>=0 && n<getHeight());return (rows[n]);}
143// Bugfix for gcc4.5 (passing assertion to the above operator):
144  Vector<typ>& operator[](int n){if(!(n>=0 && n<getHeight())){(*(const Matrix<typ>*)(this))[n];}return (rows[n]);}
145
146
147  bool operator<(const Matrix & b)const
148  {
149    if(getWidth()<b.getWidth())return true;
150    if(b.getWidth()<getWidth())return false;
151    if(getHeight()<b.getHeight())return true;
152    if(b.getHeight()<getHeight())return false;
153
154    for(int i=0;i<getHeight();i++)
155      {
156        if((*this)[i]<b[i])return true;
157        if(b[i]<(*this)[i])return false;
158      }
159    return false;
160  }
161  /**
162     Adds a times the i th row to the j th row.
163  */
164  void madd(int i, typ a, int j)
165  {
166    assert(i!=j);
167    assert(i>=0 && i<height);
168    assert(j>=0 && j<height);
169
170    if(!a.isZero())
171    for(int k=0;k<width;k++)
172      if(!rows[i][k].isZero())
173        rows[j][k].madd(rows[i][k],a);
174  }
175
176  friend std::ostream &operator<<(std::ostream &f, Matrix const &a){
177    f<<"{";
178    for(int i=0;i<a.getHeight();i++)
179      {
180        if(i)f<<","<<std::endl;
181        f<<a.rows[i];
182      }
183    f<<"}"<<std::endl;
184    return f;
185  }
186  /**
187     Swaps the i th and the j th row.
188   */
189  void swapRows(int i, int j)
190  {
191    std::swap(rows[i],rows[j]);
192  }
193  /**
194     This method is used for iterating through the pivots in a matrix
195     in row echelon form. To find the first pivot put i=-1 and
196     j=-1 and call this routine. The (i,j) th entry of the matrix is a
197     pivot. Call the routine again to get the next pivot. When no more
198     pivots are found the routine returns false.
199  */
200  bool nextPivot(int &i, int &j)const
201  {
202    i++;
203    if(i>=height)return false;
204    while(++j<width)
205      {
206        if(!rows[i][j].isZero()) return true;
207      }
208    return false;
209  }
210  /**
211     Returns the indices of the columns containing a pivot.
212     The returned list is sorted.
213     The matrix must be in row echelon form.
214   */
215  std::vector<int> pivotColumns()const
216  {
217    std::vector<int> ret;
218    int pivotI=-1;
219    int pivotJ=-1;
220    while(nextPivot(pivotI,pivotJ))ret.push_back(pivotJ);
221    return ret;
222  }
223  /**
224     Returns the indices of the columns not containing a pivot.
225     The returned list is sorted.
226     The matrix must be in row echelon form.
227   */
228  std::vector<int> nonPivotColumns()const
229  {
230    std::vector<int> ret;
231    int pivotI=-1;
232    int pivotJ=-1;
233    int firstPossiblePivot=0;
234    while(nextPivot(pivotI,pivotJ))
235      {
236        for(int j=firstPossiblePivot;j<pivotJ;j++)
237          ret.push_back(j);
238        firstPossiblePivot=pivotJ+1;
239      }
240    for(int j=firstPossiblePivot;j<getWidth();j++)
241      ret.push_back(j);
242
243    return ret;
244  }
245  /**
246     This routine removes the zero rows of the matrix.
247   */
248  void removeZeroRows()
249  {
250    int nonZeros=0;
251    for(int i=0;i<height;i++)if(!(*this)[i].isZero())nonZeros++;
252    if(nonZeros==height)return;
253
254    Matrix b(nonZeros,width);
255
256    int j=0;
257    for(int i=0;i<height;i++)
258      {
259        if(!(*this)[i].isZero())
260          {
261            b[j]=(*this)[i];
262            j++;
263          }
264      }
265    *this=b;
266  }
267  /**
268     Returns the index of a row whose index is at least currentRow and
269     which has a non-zero element on the column th entry. If no such
270     row exists then -1 is returned. This routine is used in the Gauss
271     reduction. To make the reduction more efficient the routine
272     chooses its row with as few non-zero entries as possible.
273   */
274  int findRowIndex(int column, int currentRow)const
275  {
276    int best=-1;
277    int bestNumberOfNonZero=0;
278    for(int i=currentRow;i<height;i++)
279      if(!rows[i][column].isZero())
280        {
281          int nz=0;
282          for(int k=column+1;k<width;k++)
283            if(!rows[i][k].isZero())nz++;
284          if(best==-1)
285            {
286              best=i;
287              bestNumberOfNonZero=nz;
288            }
289          else if(nz<bestNumberOfNonZero)
290            {
291              best=i;
292              bestNumberOfNonZero=nz;
293            }
294        }
295    return best;
296  }
297  /**
298     Performs a Gauss reduction and returns the number of row swaps
299     done.  The result is a matrix in row echelon form. The pivots may
300     not be all 1.  In terms of Groebner bases, what is computed is a
301     minimal (not necessarily reduced) Groebner basis of the linear
302     ideal generated by the rows.  The number of swaps is need if one
303     wants to compute the determinant afterwards. In this case it is
304     also a good idea to set the flag returnIfZeroDeterminant which
305     make the routine terminate before completion if it discovers that
306     the determinant is zero.
307  */
308  int reduce(bool returnIfZeroDeterminant=false, bool integral=false)
309  {
310    assert(integral || typ::isField());
311    int retSwaps=0;
312    int currentRow=0;
313
314    for(int i=0;i<width;i++)
315      {
316        int s=findRowIndex(i,currentRow);
317
318        if(s!=-1)
319          {
320            if(s!=currentRow)
321              {
322                swapRows(currentRow,s);
323                retSwaps++;
324              }
325            for(int j=currentRow+1;j<height;j++)
326              if(integral)
327                {
328                  if(!rows[j][i].isZero())
329                    {
330                      typ s;
331                      typ t;
332
333                      typ g=typ::gcd(rows[currentRow][i],rows[j][i],s,t);
334                      typ u=-rows[j][i]/g;
335                      typ v=rows[currentRow][i]/g;
336                        /* We want the (s,t) vector to be as small as possible.
337                         * We are allowed to adjust by multiples of (u,v).
338                         * The following computes the correct multiplier (in most cases).
339                         */
340/*                        {
341                          FieldElement multiplier=(s*u+t*v)*((u*u+v*v).inverse());
342                          double d=mpq_get_d(*(multiplier.getGmpRationalTemporaryPointer()));
343                          multiplier=multiplier.getField()->zHomomorphism(-(((int)(d+10000.5))-10000));
344                          s.madd(multiplier,u);
345                          t.madd(multiplier,v);
346                        }*/
347                        for(int k=0;k<width;k++)
348                          {
349                            typ A=rows[currentRow][k];
350                            typ B=rows[j][k];
351
352                            rows[currentRow][k]=s*A+t*B;
353                            rows[j][k]=u*A+v*B;
354                          }
355                      }
356                  }
357                else
358                  {
359                    if(!rows[j][i].isZero())
360                      madd(currentRow,-rows[j][i]/rows[currentRow][i],j);
361                  }
362              currentRow++;
363            }
364          else
365            if(returnIfZeroDeterminant)return -1;
366        }
367
368      return retSwaps;
369    }
370  /**
371     Computes a reduced row echelon form from a row echelon form. In
372     Groebner basis terms this is the same as tranforming a minimal
373     Groebner basis to a reduced one except that we do not force
374     pivots to be 1 unless the scalePivotsToOne parameter is set.
375   */
376  int REformToRREform(bool scalePivotsToOne=false)
377  {
378    int ret=0;
379    int pivotI=-1;
380    int pivotJ=-1;
381    while(nextPivot(pivotI,pivotJ))
382      {
383        if(scalePivotsToOne)
384          rows[pivotI]=rows[pivotI]/rows[pivotI][pivotJ];
385        for(int i=0;i<pivotI;i++)
386          if(!rows[i][pivotJ].isZero())
387            madd(pivotI,-rows[i][pivotJ]/rows[pivotI][pivotJ],i);
388      }
389    return ret;
390  }
391  /**
392     This function may be called if the FieldMatrix is in Row Echelon
393     Form. The input is a FieldVector which is rewritten modulo the
394     rows of the matrix. The result is unique and is the same as the
395     normal form of the input vector modulo the Groebner basis of the
396     linear ideal generated by the rows of the matrix.
397  */
398  Vector<typ> canonicalize(Vector<typ> v)const
399  {
400    assert(typ::isField());
401    assert(v.size()==getWidth());
402
403    int pivotI=-1;
404    int pivotJ=-1;
405
406    while(nextPivot(pivotI,pivotJ))
407      if(!v[pivotJ].isZero())
408      {
409        typ s=-v[pivotJ]/rows[pivotI][pivotJ];
410
411        for(int k=0;k<width;k++)
412          if(!rows[pivotI][k].isZero())
413            v[k].madd(rows[pivotI][k],s);
414      }
415    return v;
416  }
417  /**
418     Calls reduce() and constructs matrix whose rows forms a basis of
419     the kernel of the linear map defined by the original matrix. The
420     return value is the new matrix.
421   */
422  Matrix reduceAndComputeKernel()
423  {
424    Matrix ret(width-reduceAndComputeRank(),width);
425
426    REformToRREform();
427
428    int k=0;
429    int pivotI=-1;
430    int pivotJ=-1;
431    bool pivotExists=nextPivot(pivotI,pivotJ);
432    for(int j=0;j<width;j++)
433      {
434        if(pivotExists && (pivotJ==j))
435          {
436            pivotExists=nextPivot(pivotI,pivotJ);
437            continue;
438          }
439        int pivot2I=-1;
440        int pivot2J=-1;
441        while(nextPivot(pivot2I,pivot2J))
442          {
443            ret[k][pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J];
444          }
445        ret[k][j]=typ(-1);
446        k++;
447      }
448    return ret;
449  }
450  /**
451     Assumes that the matrix has a kernel of dimension 1.
452     Calls reduce() and returns a non-zero vector in the kernel.
453     If the matrix is an (n-1)x(n) matrix then the returned vector has
454     the property that if it was appended as a row to the original matrix
455     then the determinant of that matrix would be positive. Of course
456     this property, as described here, only makes sense for ordered fields.
457     Only allowed for fields at the moment.
458   */
459  Vector<typ> reduceAndComputeVectorInKernel()
460  {
461    assert(typ::isField());
462    // TODO: (optimization) if the field is ordered, then it is better to just keep track of signs rather than
463    // multiplying by sign*diagonalProduct*lastEntry at the end.
464    int nswaps=this->reduce();
465    typ sign=typ(1-2*(nswaps&1));
466    int rank=reduceAndComputeRank();
467    assert(rank+1==width);
468
469    REformToRREform();
470
471    Vector<typ> ret(width);
472
473    typ diagonalProduct(1);
474    {
475      int pivot2I=-1;
476      int pivot2J=-1;
477      while(nextPivot(pivot2I,pivot2J))
478        {
479          diagonalProduct*=rows[pivot2I][pivot2J];
480        }
481    }
482    {
483      int j=nonPivotColumns().front();
484      int pivot2I=-1;
485      int pivot2J=-1;
486      ret[j]=typ(-1);
487      // Pretend that we are appending ret to the matrix, and reducing this
488      // new row by the previous ones. The last entry of the resulting matrix
489      // is lastEntry.
490      typ lastEntry=ret[j];
491      while(nextPivot(pivot2I,pivot2J))
492        {
493          ret[pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J];
494          lastEntry-=ret[pivot2J]*ret[pivot2J];
495        }
496      ret=(sign*(diagonalProduct*lastEntry))*ret;
497    }
498
499    return ret;
500  }
501
502  /**
503     Calls reduce() and returns the rank of the matrix.
504   */
505  int reduceAndComputeRank()
506  {
507    reduce();
508    int ret=0;
509    int pivotI=-1;
510    int pivotJ=-1;
511    while(nextPivot(pivotI,pivotJ))ret++;
512    return ret;
513  }
514  /**
515   * Sort the rows of the matrix.
516   */
517  void sortRows()
518  {
519    std::sort(rows.begin(),rows.end());
520  }
521  /**
522   * Sort the rows of the matrix and remove duplicate rows.
523   */
524  void sortAndRemoveDuplicateRows()
525  {
526    sortRows();
527    if(getHeight()==0)return;
528    Matrix B(0,getWidth());
529    B.appendRow((*this)[0]);
530    for(int i=1;i<getHeight();i++)
531      if(rows[i]!=rows[i-1])B.appendRow((*this)[i]);
532    *this=B;
533  }
534  /**
535     Takes two matrices with the same number of columns and construct
536     a new matrix which has the rows of the matrix top on the top and
537     the rows of the matrix bottom at the bottom. The return value is
538     the constructed matrix.
539   */
540  friend Matrix combineOnTop(Matrix const &top, Matrix const &bottom)
541  {
542    assert(bottom.getWidth()==top.getWidth());
543    Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth());
544    for(int i=0;i<top.getHeight();i++)ret.rows[i]=top.rows[i];
545    for(int i=0;i<bottom.getHeight();i++)ret.rows[i+top.getHeight()]=bottom.rows[i];
546
547    return ret;
548  }
549  /**
550     Takes two matrices with the same number of rows and construct
551     a new matrix which has the columns of the matrix left on the left and
552     the columns of the matrix right on the right. The return value is
553     the constructed matrix.
554   */
555  friend Matrix combineLeftRight(Matrix const &left, Matrix const &right)
556  {
557    assert(left.getHeight()==right.getHeight());
558    Matrix ret(left.getHeight(),left.getWidth()+right.getWidth());
559    for(int i=0;i<left.getHeight();i++)
560      {
561        for(int j=0;j<left.getWidth();j++)ret.rows[i][j]=left.rows[i][j];
562        for(int j=0;j<right.getWidth();j++)ret.rows[i][j+left.getWidth()]=right.rows[i][j];
563      }
564    return ret;
565  }
566};
567
568typedef Matrix<Integer> ZMatrix;
569typedef Matrix<Rational> QMatrix;
570typedef Matrix<int> IntMatrix;
571
572inline QMatrix ZToQMatrix(ZMatrix const &m)
573{
574  QMatrix ret(m.getHeight(),m.getWidth());
575  for(int i=0;i<m.getHeight();i++)ret[i]=ZToQVector(m[i]);
576  return ret;
577}
578
579inline ZMatrix QToZMatrixPrimitive(QMatrix const &m)
580{
581  ZMatrix ret(m.getHeight(),m.getWidth());
582  for(int i=0;i<m.getHeight();i++)ret[i]=QToZVectorPrimitive(m[i]);
583  return ret;
584}
585
586
587inline IntMatrix ZToIntMatrix(ZMatrix const &m)
588{
589  IntMatrix ret(m.getHeight(),m.getWidth());
590  for(int i=0;i<m.getHeight();i++)ret[i]=ZToIntVector(m[i]);
591  return ret;
592}
593
594
595inline ZMatrix IntToZMatrix(IntMatrix const &m)
596{
597  ZMatrix ret(m.getHeight(),m.getWidth());
598  for(int i=0;i<m.getHeight();i++)ret[i]=IntToZVector(m[i]);
599  return ret;
600}
601
602inline QMatrix canonicalizeSubspace(QMatrix const &m)
603{
604  QMatrix temp=m;
605  temp.reduce();
606  temp.REformToRREform();
607  temp.removeZeroRows();
608  return temp;
609}
610
611inline ZMatrix canonicalizeSubspace(ZMatrix const &m)
612{
613  return QToZMatrixPrimitive(canonicalizeSubspace(ZToQMatrix(m)));
614}
615
616
617inline QMatrix kernel(QMatrix const &m)
618{
619  QMatrix temp=m;
620  return temp.reduceAndComputeKernel();
621}
622
623inline ZMatrix kernel(ZMatrix const &m)
624{
625  return QToZMatrixPrimitive(kernel(ZToQMatrix(m)));
626}
627
628}
629
630
631#endif /* LIB_ZMATRIX_H_ */
Note: See TracBrowser for help on using the repository browser.