source: git/gfanlib/gfanlib_matrix.h @ 26b713

jengelh-datetimespielwiese
Last change on this file since 26b713 was 26b713, checked in by Yue Ren <ren@…>, 10 years ago
chg: updated gfanlib package to newest version in master
  • Property mode set to 100644
File size: 17.7 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):width(a.getWidth()),height(a.getHeight()),rows(a.rows){
24  }
25  Matrix(int height_, int width_):width(width_),height(height_),rows(height_){
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, bool makePivotsOne=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            if(makePivotsOne)
326              {//THE PIVOT SHOULD BE SET TO ONE IF INTEGRAL IS FALSE
327                if(!rows[currentRow][i].sign()<0)retSwaps++;
328                typ inverse=typ(1)/rows[currentRow][i];
329                //                if(!rows[currentRow][i].isOne())
330                  {
331                    for(int k=0;k<width;k++)
332                      if(!rows[currentRow][k].isZero())
333                        rows[currentRow][k]*=inverse;
334                  }
335              }
336            for(int j=currentRow+1;j<height;j++)
337              if(integral)
338                {
339                  if(!rows[j][i].isZero())
340                    {
341                      typ s;
342                      typ t;
343
344                      typ g=typ::gcd(rows[currentRow][i],rows[j][i],s,t);
345                      typ u=-rows[j][i]/g;
346                      typ v=rows[currentRow][i]/g;
347                        /* We want the (s,t) vector to be as small as possible.
348                         * We are allowed to adjust by multiples of (u,v).
349                         * The following computes the correct multiplier (in most cases).
350                         */
351/*                        {
352                          FieldElement multiplier=(s*u+t*v)*((u*u+v*v).inverse());
353                          double d=mpq_get_d(*(multiplier.getGmpRationalTemporaryPointer()));
354                          multiplier=multiplier.getField()->zHomomorphism(-(((int)(d+10000.5))-10000));
355                          s.madd(multiplier,u);
356                          t.madd(multiplier,v);
357                        }*/
358                        for(int k=0;k<width;k++)
359                          {
360                            typ A=rows[currentRow][k];
361                            typ B=rows[j][k];
362
363                            rows[currentRow][k]=s*A+t*B;
364                            rows[j][k]=u*A+v*B;
365                          }
366                      }
367                  }
368                else
369                  {
370                    if(!rows[j][i].isZero())
371                      madd(currentRow,-rows[j][i]/rows[currentRow][i],j);
372                  }
373              currentRow++;
374            }
375          else
376            if(returnIfZeroDeterminant)return -1;
377        }
378
379      return retSwaps;
380    }
381  /**
382     Computes a reduced row echelon form from a row echelon form. In
383     Groebner basis terms this is the same as tranforming a minimal
384     Groebner basis to a reduced one except that we do not force
385     pivots to be 1 unless the scalePivotsToOne parameter is set.
386   */
387  int REformToRREform(bool scalePivotsToOne=false)
388  {
389    int ret=0;
390    int pivotI=-1;
391    int pivotJ=-1;
392    while(nextPivot(pivotI,pivotJ))
393      {
394        if(scalePivotsToOne)
395          rows[pivotI]=rows[pivotI]/rows[pivotI][pivotJ];
396        for(int i=0;i<pivotI;i++)
397          if(!rows[i][pivotJ].isZero())
398            madd(pivotI,-rows[i][pivotJ]/rows[pivotI][pivotJ],i);
399      }
400    return ret;
401  }
402  /**
403     This function may be called if the FieldMatrix is in Row Echelon
404     Form. The input is a FieldVector which is rewritten modulo the
405     rows of the matrix. The result is unique and is the same as the
406     normal form of the input vector modulo the Groebner basis of the
407     linear ideal generated by the rows of the matrix.
408  */
409  Vector<typ> canonicalize(Vector<typ> v)const
410  {
411    assert(typ::isField());
412    assert(v.size()==getWidth());
413
414    int pivotI=-1;
415    int pivotJ=-1;
416
417    while(nextPivot(pivotI,pivotJ))
418      if(!v[pivotJ].isZero())
419      {
420        typ s=-v[pivotJ]/rows[pivotI][pivotJ];
421
422        for(int k=0;k<width;k++)
423          if(!rows[pivotI][k].isZero())
424            v[k].madd(rows[pivotI][k],s);
425      }
426    return v;
427  }
428  /**
429     Calls reduce() and constructs matrix whose rows forms a basis of
430     the kernel of the linear map defined by the original matrix. The
431     return value is the new matrix.
432   */
433  Matrix reduceAndComputeKernel()
434  {
435    Matrix ret(width-reduceAndComputeRank(),width);
436
437    REformToRREform();
438
439    int k=0;
440    int pivotI=-1;
441    int pivotJ=-1;
442    bool pivotExists=nextPivot(pivotI,pivotJ);
443    for(int j=0;j<width;j++)
444      {
445        if(pivotExists && (pivotJ==j))
446          {
447            pivotExists=nextPivot(pivotI,pivotJ);
448            continue;
449          }
450        int pivot2I=-1;
451        int pivot2J=-1;
452        while(nextPivot(pivot2I,pivot2J))
453          {
454            ret[k][pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J];
455          }
456        ret[k][j]=typ(-1);
457        k++;
458      }
459    return ret;
460  }
461  /**
462     Assumes that the matrix has a kernel of dimension 1.
463     Calls reduce() and returns a non-zero vector in the kernel.
464     If the matrix is an (n-1)x(n) matrix then the returned vector has
465     the property that if it was appended as a row to the original matrix
466     then the determinant of that matrix would be positive. Of course
467     this property, as described here, only makes sense for ordered fields.
468     Only allowed for fields at the moment.
469   */
470  Vector<typ> reduceAndComputeVectorInKernel()
471  {
472    assert(typ::isField());
473    // TODO: (optimization) if the field is ordered, then it is better to just keep track of signs rather than
474    // multiplying by sign*diagonalProduct*lastEntry at the end.
475    int nswaps=this->reduce();
476    typ sign=typ(1-2*(nswaps&1));
477    int rank=reduceAndComputeRank();
478    assert(rank+1==width);
479
480    REformToRREform();
481
482    Vector<typ> ret(width);
483
484    typ diagonalProduct(1);
485    {
486      int pivot2I=-1;
487      int pivot2J=-1;
488      while(nextPivot(pivot2I,pivot2J))
489        {
490          diagonalProduct*=rows[pivot2I][pivot2J];
491        }
492    }
493    {
494      int j=nonPivotColumns().front();
495      int pivot2I=-1;
496      int pivot2J=-1;
497      ret[j]=typ(-1);
498      // Pretend that we are appending ret to the matrix, and reducing this
499      // new row by the previous ones. The last entry of the resulting matrix
500      // is lastEntry.
501      typ lastEntry=ret[j];
502      while(nextPivot(pivot2I,pivot2J))
503        {
504          ret[pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J];
505          lastEntry-=ret[pivot2J]*ret[pivot2J];
506        }
507      ret=(sign*(diagonalProduct*lastEntry))*ret;
508    }
509
510    return ret;
511  }
512
513  /**
514     Calls reduce() and returns the rank of the matrix.
515   */
516  int reduceAndComputeRank()
517  {
518    if (typ::isField())
519      reduce();
520    else
521      reduce(false,true,false);
522    int ret=0;
523    int pivotI=-1;
524    int pivotJ=-1;
525    while(nextPivot(pivotI,pivotJ))ret++;
526    return ret;
527  }
528  /**
529   * Sort the rows of the matrix.
530   */
531  void sortRows()
532  {
533    std::sort(rows.begin(),rows.end());
534  }
535  /**
536   * Sort the rows of the matrix and remove duplicate rows.
537   */
538  void sortAndRemoveDuplicateRows()
539  {
540    sortRows();
541    if(getHeight()==0)return;
542    Matrix B(0,getWidth());
543    B.appendRow((*this)[0]);
544    for(int i=1;i<getHeight();i++)
545      if(rows[i]!=rows[i-1])B.appendRow((*this)[i]);
546    *this=B;
547  }
548  /**
549     Takes two matrices with the same number of columns and construct
550     a new matrix which has the rows of the matrix top on the top and
551     the rows of the matrix bottom at the bottom. The return value is
552     the constructed matrix.
553   */
554  friend Matrix combineOnTop(Matrix const &top, Matrix const &bottom)
555  {
556    assert(bottom.getWidth()==top.getWidth());
557    Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth());
558    for(int i=0;i<top.getHeight();i++)ret.rows[i]=top.rows[i];
559    for(int i=0;i<bottom.getHeight();i++)ret.rows[i+top.getHeight()]=bottom.rows[i];
560
561    return ret;
562  }
563  /**
564     Takes two matrices with the same number of rows and construct
565     a new matrix which has the columns of the matrix left on the left and
566     the columns of the matrix right on the right. The return value is
567     the constructed matrix.
568   */
569  friend Matrix combineLeftRight(Matrix const &left, Matrix const &right)
570  {
571    assert(left.getHeight()==right.getHeight());
572    Matrix ret(left.getHeight(),left.getWidth()+right.getWidth());
573    for(int i=0;i<left.getHeight();i++)
574      {
575        for(int j=0;j<left.getWidth();j++)ret.rows[i][j]=left.rows[i][j];
576        for(int j=0;j<right.getWidth();j++)ret.rows[i][j+left.getWidth()]=right.rows[i][j];
577      }
578    return ret;
579  }
580};
581
582typedef Matrix<Integer> ZMatrix;
583typedef Matrix<Rational> QMatrix;
584typedef Matrix<int> IntMatrix;
585
586inline QMatrix ZToQMatrix(ZMatrix const &m)
587{
588  QMatrix ret(m.getHeight(),m.getWidth());
589  for(int i=0;i<m.getHeight();i++)ret[i]=ZToQVector(m[i]);
590  return ret;
591}
592
593inline ZMatrix QToZMatrixPrimitive(QMatrix const &m)
594{
595  ZMatrix ret(m.getHeight(),m.getWidth());
596  for(int i=0;i<m.getHeight();i++)ret[i]=QToZVectorPrimitive(m[i]);
597  return ret;
598}
599
600
601inline IntMatrix ZToIntMatrix(ZMatrix const &m)
602{
603  IntMatrix ret(m.getHeight(),m.getWidth());
604  for(int i=0;i<m.getHeight();i++)ret[i]=ZToIntVector(m[i]);
605  return ret;
606}
607
608
609inline ZMatrix IntToZMatrix(IntMatrix const &m)
610{
611  ZMatrix ret(m.getHeight(),m.getWidth());
612  for(int i=0;i<m.getHeight();i++)ret[i]=IntToZVector(m[i]);
613  return ret;
614}
615
616inline QMatrix canonicalizeSubspace(QMatrix const &m)
617{
618  QMatrix temp=m;
619  temp.reduce();
620  temp.REformToRREform();
621  temp.removeZeroRows();
622  return temp;
623}
624
625inline ZMatrix canonicalizeSubspace(ZMatrix const &m)
626{
627  return QToZMatrixPrimitive(canonicalizeSubspace(ZToQMatrix(m)));
628}
629
630
631inline QMatrix kernel(QMatrix const &m)
632{
633  QMatrix temp=m;
634  return temp.reduceAndComputeKernel();
635}
636
637inline ZMatrix kernel(ZMatrix const &m)
638{
639  return QToZMatrixPrimitive(kernel(ZToQMatrix(m)));
640}
641
642}
643
644
645#endif /* LIB_ZMATRIX_H_ */
Note: See TracBrowser for help on using the repository browser.