source: git/gfanlib/gfanlib_matrix.h

spielwiese
Last change on this file was 15813d, checked in by Hans Schoenemann <hannes@…>, 8 years ago
format
  • Property mode set to 100644
File size: 23.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;
20  std::vector<typ> data;
21public:
22  // rowIterator;
23 // std::vector<Vector<typ> >::iterator rowsBegin(){return rows.begin();}
24//  std::vector<Vector<typ> >::iterator rowsEnd(){return rows.end();}
25  inline int getHeight()const{return height;};
26  inline int getWidth()const{return width;};
27  Matrix(const Matrix &a):width(a.getWidth()),height(a.getHeight()),data(a.data){
28  }
29  Matrix(int height_, int width_):width(width_),height(height_),data(width_*height_){
30    assert(height>=0);
31    assert(width>=0);
32  };
33  ~Matrix(){
34  };
35  Matrix():width(0),height(0){
36  };
37  static Matrix rowVectorMatrix(Vector<typ> const &v)
38  {
39    Matrix ret(1,v.size());
40    for(int i=0;i<v.size();i++)ret[0][i]=v[i];
41    return ret;
42  }
43  Vector<typ> column(int i)const
44    {
45      assert(i>=0);
46      assert(i<getWidth());
47      Vector<typ> ret(getHeight());
48      for(int j=0;j<getHeight();j++)ret[j]=(*this)[j][i];
49      return ret;
50    }
51  Matrix transposed()const
52    {
53      Matrix ret(getWidth(),getHeight());
54      for(int i=0;i<getWidth();i++)
55          for(int j=0;j<getHeight();j++)
56                  ret[i][j]=(*this)[j][i];
57      return ret;
58    }
59  static Matrix identity(int n)
60    {
61      Matrix m(n,n);
62      for(int i=0;i<n;i++)m[i][i]=typ(1);
63      return m;
64    }
65  void append(Matrix const &m)
66    {
67      assert(m.getWidth()==width);
68          data.resize((height+m.height)*width);
69          int oldHeight=height;
70      height+=m.height;
71      for(int i=0;i<m.height;i++)
72        {
73          for(int j=0;j<m.width;j++)
74                  (*this)[i+oldHeight][j]=m[i][j];
75        }
76    }
77  void appendRow(Vector<typ> const &v)
78    {
79          assert(v.size()==width);
80          data.resize((height+1)*width);
81          height++;
82          for(int j=0;j<width;j++)
83                  (*this)[height-1][j]=v[j];
84    }
85  void eraseLastRow()
86  {
87    assert(height>0);
88    data.resize((height-1)*width);
89    height--;
90  }
91  /*IntegerVector vectormultiply(IntegerVector const &v)const
92    {
93      assert(v.size()==width);
94      IntegerVector ret(height);
95      for(int i=0;i<height;i++)
96        ret[i]=dot(rows[i],v);
97      return ret;
98    }*/
99  /**
100   * Decides if v is in the kernel of the matrix.
101   */
102/*  bool inKernel(IntegerVector const &v)const
103    {
104      assert(v.size()==width);
105      for(int i=0;i<height;i++)
106          if(dotLong(rows[i],v)!=0)return false;
107      return true;
108    }
109*/
110  friend Matrix operator*(const typ &s, const Matrix& q)
111    {
112      Matrix p=q;
113      for(int i=0;i<q.height;i++)
114          for(int j=0;j<q.width;j++)
115                  p[i][j]=s*(q[i][j]);
116      return p;
117    }
118  friend Matrix operator*(const Matrix& a, const Matrix& b)
119    {
120      assert(a.width==b.height);
121      Matrix ret(b.width,a.height);
122      for(int i=0;i<b.width;i++)
123        ret[i]=a.vectormultiply(b.column(i));
124      return ret.transposed();
125    }
126  /*  template<class T>
127    Matrix<T>(const Matrix<T>& c):v(c.size()){
128    for(int i=0;i<size();i++)v[i]=typ(c[i]);}
129  */
130  friend Matrix operator-(const Matrix &b)
131  {
132    Matrix ret(b.height,b.width);
133    for(int i=0;i<b.height;i++)ret[i]=-b[i];
134    return ret;
135  }
136
137  /**
138     Returns the specified submatrix. The endRow and endColumn are not included.
139   */
140  Matrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const
141  {
142    assert(startRow>=0);
143    assert(startColumn>=0);
144    assert(endRow>=startRow);
145    assert(endColumn>=startColumn);
146    assert(endRow<=height);
147    assert(endColumn<=width);
148    Matrix ret(endRow-startRow,endColumn-startColumn);
149    for(int i=startRow;i<endRow;i++)
150      for(int j=startColumn;j<endColumn;j++)
151        ret[i-startRow][j-startColumn]=(*this)[i][j];
152    return ret;
153  }
154
155  class RowRef;
156  class const_RowRef{
157    int rowNumM;
158    Matrix const &matrix;
159    friend class RowRef;
160  public:
161  inline const_RowRef(const Matrix  &matrix_, int rowNum_)__attribute__((always_inline)):
162    rowNumM(rowNum_*matrix_.width),
163      matrix(matrix_)
164      {
165      }
166  inline typ const &operator[](int j)const __attribute__((always_inline))
167    {
168        assert(j>=0);
169        assert(j<matrix.width);
170        return matrix.data[rowNumM+j];
171    }
172  inline typ const &UNCHECKEDACCESS(int j)const __attribute__((always_inline))
173    {
174        return matrix.data[rowNumM+j];
175    }
176    const Vector<typ> toVector()const
177    {
178      Vector<typ> ret(matrix.width);
179      for(int j=0;j<matrix.width;j++)
180              ret[j]=matrix.data[rowNumM+j];
181      return ret;
182    }
183    operator Vector<typ>()const
184                {
185                        return toVector();
186                }
187    bool operator==(Vector<typ> const &b)const
188                {
189                        return toVector()==b;
190                }
191/*    typ dot(Vector<typ> const &b)const
192                {
193                        return dot(toVector(),b);
194                }*/
195    Vector<typ> operator-()const
196    {
197            return -toVector();
198    }
199  };
200  class RowRef{
201    int rowNumM;
202    Matrix &matrix;
203  public:
204  inline RowRef(Matrix &matrix_, int rowNum_):
205    rowNumM(rowNum_*matrix_.width),
206      matrix(matrix_)
207      {
208      }
209    inline typ &operator[](int j) __attribute__((always_inline))
210      {
211            assert(j>=0);
212            assert(j<matrix.width);
213            return matrix.data[rowNumM+j];
214      }
215    inline typ &UNCHECKEDACCESS(int j)
216      {
217            return matrix.data[rowNumM+j];
218      }
219    RowRef &operator=(Vector<typ> const &v)
220    {
221        assert(v.size()==matrix.width);
222        for(int j=0;j<matrix.width;j++)
223                matrix.data[rowNumM+j]=v[j];
224
225            return *this;
226    }
227    RowRef &operator=(RowRef const &v)
228    {
229        assert(v.matrix.width==matrix.width);
230        for(int j=0;j<matrix.width;j++)
231                matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j];
232
233            return *this;
234    }
235/*    RowRef &operator+=(Vector<typ> const &v)
236    {
237        assert(v.size()==matrix.width);
238        for(int j=0;j<matrix.width;j++)
239                matrix.data[rowNumM+j]+=v.v[j];
240
241            return *this;
242    }*/
243    RowRef &operator+=(RowRef const &v)
244    {
245        assert(v.matrix.width==matrix.width);
246        for(int j=0;j<matrix.width;j++)
247                matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j];
248
249            return *this;
250    }
251    RowRef &operator+=(const_RowRef const &v)
252    {
253        assert(v.matrix.width==matrix.width);
254        for(int j=0;j<matrix.width;j++)
255                matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j];
256
257            return *this;
258    }
259    RowRef &operator=(const_RowRef const &v)
260    {
261        assert(v.matrix.width==matrix.width);
262        for(int j=0;j<matrix.width;j++)
263                matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j];
264
265            return *this;
266    }
267    const Vector<typ> toVector()const
268    {
269      Vector<typ> ret(matrix.width);
270      for(int j=0;j<matrix.width;j++)
271              ret[j]=matrix.data[rowNumM+j];
272      return ret;
273    }
274    operator Vector<typ>()const
275                {
276                        return toVector();
277                }
278/*    typ dot(Vector<typ> const &b)const
279                {
280                        return dot(toVector(),b);
281                }*/
282    bool isZero()const
283        {
284          for(int j=0;j<matrix.width;j++)if(!(matrix.data[rowNumM+j].isZero()))return false;
285          return true;
286        }
287  };
288
289
290  inline RowRef operator[](int i) __attribute__((always_inline))
291  {
292    assert(i>=0);
293    assert(i<height);
294    return RowRef(*this,i);
295  }
296  inline const_RowRef operator[](int i)const __attribute__((always_inline))
297  {
298    assert(i>=0);
299    assert(i<height);
300    return const_RowRef(*this,i);
301  }
302
303
304  const typ& UNCHECKEDACCESS(int i,int j)const __attribute__((always_inline)){
305/*            assert(i>=0);
306            assert(i<height);
307            assert(j>=0);
308            assert(j<width);*/
309          return data[i*width+j];}
310  typ& UNCHECKEDACCESS(int i,int j) __attribute__((always_inline)){
311/*            assert(i>=0);
312            assert(i<height);
313            assert(j>=0);
314            assert(j<width);*/
315            return data[i*width+j];}
316
317
318
319  bool operator<(const Matrix & b)const
320  {
321    if(getWidth()<b.getWidth())return true;
322    if(b.getWidth()<getWidth())return false;
323    if(getHeight()<b.getHeight())return true;
324    if(b.getHeight()<getHeight())return false;
325
326    for(int i=0;i<getHeight();i++)
327      {
328        if((*this)[i].toVector()<b[i].toVector())return true;
329        if(b[i].toVector()<(*this)[i].toVector())return false;
330      }
331    return false;
332  }
333  /**
334     Adds a times the i th row to the j th row.
335  */
336  void madd(int i, typ a, int j)
337  {
338    assert(i!=j);
339    assert(i>=0 && i<height);
340    assert(j>=0 && j<height);
341
342    if(!a.isZero())
343    for(int k=0;k<width;k++)
344      if(!(*this)[i][k].isZero())
345              (*this)[j][k].madd((*this)[i][k],a);
346  }
347
348  friend std::ostream &operator<<(std::ostream &f, Matrix const &a){
349    f<<"{";
350    for(int i=0;i<a.getHeight();i++)
351      {
352        if(i)f<<","<<std::endl;
353        f<<a[i].toVector();
354      }
355    f<<"}"<<std::endl;
356    return f;
357  }
358
359  std::string toString()const
360  {
361          std::stringstream f;
362          f<<*this;
363          return f.str();
364  }
365
366  /**
367     Swaps the i th and the j th row.
368   */
369  void swapRows(int i, int j)
370  {
371    for(int a=0;a<width;a++)std::swap((*this)[i][a],(*this)[j][a]);
372  }
373  /**
374     This method is used for iterating through the pivots in a matrix
375     in row echelon form. To find the first pivot put i=-1 and
376     j=-1 and call this routine. The (i,j) th entry of the matrix is a
377     pivot. Call the routine again to get the next pivot. When no more
378     pivots are found the routine returns false.
379  */
380  bool nextPivot(int &i, int &j)const
381  {
382    i++;
383    if(i>=height)return false;
384    while(++j<width)
385      {
386        if(!(*this)[i][j].isZero()) return true;
387      }
388    return false;
389  }
390  /**
391     Returns the indices of the columns containing a pivot.
392     The returned list is sorted.
393     The matrix must be in row echelon form.
394   */
395  std::vector<int> pivotColumns()const
396  {
397    std::vector<int> ret;
398    int pivotI=-1;
399    int pivotJ=-1;
400    while(nextPivot(pivotI,pivotJ))ret.push_back(pivotJ);
401    return ret;
402  }
403  /**
404     Returns the indices of the columns not containing a pivot.
405     The returned list is sorted.
406     The matrix must be in row echelon form.
407   */
408  std::vector<int> nonPivotColumns()const
409  {
410    std::vector<int> ret;
411    int pivotI=-1;
412    int pivotJ=-1;
413    int firstPossiblePivot=0;
414    while(nextPivot(pivotI,pivotJ))
415      {
416        for(int j=firstPossiblePivot;j<pivotJ;j++)
417          ret.push_back(j);
418        firstPossiblePivot=pivotJ+1;
419      }
420    for(int j=firstPossiblePivot;j<getWidth();j++)
421      ret.push_back(j);
422
423    return ret;
424  }
425  /**
426     This routine removes the zero rows of the matrix.
427   */
428  void removeZeroRows()
429  {
430    int nonZeros=0;
431    for(int i=0;i<height;i++)if(!(*this)[i].isZero())nonZeros++;
432    if(nonZeros==height)return;
433
434    Matrix b(nonZeros,width);
435
436    int j=0;
437    for(int i=0;i<height;i++)
438      {
439        if(!(*this)[i].isZero())
440          {
441            b[j]=(*this)[i];
442            j++;
443          }
444      }
445    *this=b;
446  }
447  /**
448     Returns the index of a row whose index is at least currentRow and
449     which has a non-zero element on the column th entry. If no such
450     row exists then -1 is returned. This routine is used in the Gauss
451     reduction. To make the reduction more efficient the routine
452     chooses its row with as few non-zero entries as possible.
453   */
454  int findRowIndex(int column, int currentRow)const
455  {
456    int best=-1;
457    int bestNumberOfNonZero=0;
458    for(int i=currentRow;i<height;i++)
459      if(!(*this)[i][column].isZero())
460        {
461          int nz=0;
462          for(int k=column+1;k<width;k++)
463            if(!(*this)[i][k].isZero())nz++;
464          if(best==-1)
465            {
466              best=i;
467              bestNumberOfNonZero=nz;
468            }
469          else if(nz<bestNumberOfNonZero)
470            {
471              best=i;
472              bestNumberOfNonZero=nz;
473            }
474        }
475    return best;
476  }
477  /**
478     Performs a Gauss reduction and returns the number of row swaps (and negative scalings)
479     done. The result is a matrix in row echelon form. The pivots may
480     not be all 1.  In terms of Groebner bases, what is computed is a
481     minimal (not necessarily reduced) Groebner basis of the linear
482     ideal generated by the rows.  The number of swaps is need if one
483     wants to compute the determinant afterwards. In this case it is
484     also a good idea to set the flag returnIfZeroDeterminant which
485     make the routine terminate before completion if it discovers that
486     the determinant is zero.
487  */
488  int reduce(bool returnIfZeroDeterminant=false, bool integral=false, bool makePivotsOne=false)
489  {
490    assert(integral || typ::isField());
491    assert(!makePivotsOne || !integral);
492
493    int retSwaps=0;
494    int currentRow=0;
495
496    for(int i=0;i<width;i++)
497      {
498        int s=findRowIndex(i,currentRow);
499
500        if(s!=-1)
501          {
502            if(s!=currentRow)
503              {
504                swapRows(currentRow,s);
505                retSwaps++;
506              }
507            if(makePivotsOne)
508              {//THE PIVOT SHOULD BE SET TO ONE IF INTEGRAL IS FALSE
509                if((*this)[currentRow][i].sign()>=0)retSwaps++;
510                typ inverse=typ(1)/(*this)[currentRow][i];
511                //                if(!rows[currentRow][i].isOne())
512                  {
513                    for(int k=0;k<width;k++)
514                      if(!(*this)[currentRow][k].isZero())
515                        (*this)[currentRow][k]*=inverse;
516                  }
517              }
518            for(int j=currentRow+1;j<height;j++)
519              if(integral)
520                {
521                  if(!(*this)[j][i].isZero())
522                    {
523                      typ s;
524                      typ t;
525
526                      typ g=typ::gcd((*this)[currentRow][i],(*this)[j][i],s,t);
527                      typ u=-(*this)[j][i]/g;
528                      typ v=(*this)[currentRow][i]/g;
529                        /* We want the (s,t) vector to be as small as possible.
530                         * We are allowed to adjust by multiples of (u,v).
531                         * The following computes the correct multiplier (in most cases).
532                         */
533/*                        {
534                          FieldElement multiplier=(s*u+t*v)*((u*u+v*v).inverse());
535                          double d=mpq_get_d(*(multiplier.getGmpRationalTemporaryPointer()));
536                          multiplier=multiplier.getField()->zHomomorphism(-(((int)(d+10000.5))-10000));
537                          s.madd(multiplier,u);
538                          t.madd(multiplier,v);
539                        }*/
540                        for(int k=0;k<width;k++)
541                          {
542                            typ A=(*this)[currentRow][k];
543                            typ B=(*this)[j][k];
544
545                            (*this)[currentRow][k]=s*A+t*B;
546                            (*this)[j][k]=u*A+v*B;
547                          }
548                      }
549                  }
550                else
551                  {
552                    if(!(*this)[j][i].isZero())
553                      madd(currentRow,-(*this)[j][i]/(*this)[currentRow][i],j);
554                  }
555              currentRow++;
556            }
557          else
558            if(returnIfZeroDeterminant)return -1;
559        }
560
561      return retSwaps;
562    }
563  /**
564     Computes a reduced row echelon form from a row echelon form. In
565     Groebner basis terms this is the same as tranforming a minimal
566     Groebner basis to a reduced one except that we do not force
567     pivots to be 1 unless the scalePivotsToOne parameter is set.
568   */
569  int REformToRREform(bool scalePivotsToOne=false)
570  {
571    int ret=0;
572    int pivotI=-1;
573    int pivotJ=-1;
574    while(nextPivot(pivotI,pivotJ))
575      {
576            if(scalePivotsToOne)
577          (*this)[pivotI]=(*this)[pivotI].toVector()/(*this)[pivotI][pivotJ];
578        for(int i=0;i<pivotI;i++)
579          if(!(*this)[i][pivotJ].isZero())
580            madd(pivotI,-(*this)[i][pivotJ]/(*this)[pivotI][pivotJ],i);
581        }
582    return ret;
583  }
584  /**
585     This function may be called if the FieldMatrix is in Row Echelon
586     Form. The input is a FieldVector which is rewritten modulo the
587     rows of the matrix. The result is unique and is the same as the
588     normal form of the input vector modulo the Groebner basis of the
589     linear ideal generated by the rows of the matrix.
590  */
591  Vector<typ> canonicalize(Vector<typ> v)const
592  {
593    assert(typ::isField());
594    assert((int)v.size()==getWidth());
595
596    int pivotI=-1;
597    int pivotJ=-1;
598
599    while(nextPivot(pivotI,pivotJ))
600      if(!v[pivotJ].isZero())
601      {
602        typ s=-v[pivotJ]/(*this)[pivotI][pivotJ];
603
604        for(int k=0;k<width;k++)
605          if(!(*this)[pivotI][k].isZero())
606            v[k].madd((*this)[pivotI][k],s);
607      }
608    return v;
609  }
610  /**
611     Calls reduce() and constructs matrix whose rows forms a basis of
612     the kernel of the linear map defined by the original matrix. The
613     return value is the new matrix.
614   */
615  Matrix reduceAndComputeKernel()
616  {
617    Matrix ret(width-reduceAndComputeRank(),width);
618
619    REformToRREform();
620
621    int k=0;
622    int pivotI=-1;
623    int pivotJ=-1;
624    bool pivotExists=nextPivot(pivotI,pivotJ);
625    for(int j=0;j<width;j++)
626      {
627        if(pivotExists && (pivotJ==j))
628          {
629            pivotExists=nextPivot(pivotI,pivotJ);
630            continue;
631          }
632        int pivot2I=-1;
633        int pivot2J=-1;
634        while(nextPivot(pivot2I,pivot2J))
635          {
636            ret[k][pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J];
637          }
638        ret[k][j]=typ(-1);
639        k++;
640      }
641    return ret;
642  }
643  /**
644     Assumes that the matrix has a kernel of dimension 1.
645     Calls reduce() and returns a non-zero vector in the kernel.
646     If the matrix is an (n-1)x(n) matrix then the returned vector has
647     the property that if it was appended as a row to the original matrix
648     then the determinant of that matrix would be positive. Of course
649     this property, as described here, only makes sense for ordered fields.
650     Only allowed for fields at the moment.
651   */
652  Vector<typ> reduceAndComputeVectorInKernel()
653  {
654    assert(typ::isField());
655    // TODO: (optimization) if the field is ordered, then it is better to just keep track of signs rather than
656    // multiplying by sign*diagonalProduct*lastEntry at the end.
657    int nswaps=this->reduce();
658    typ sign=typ(1-2*(nswaps&1));
659    int rank=reduceAndComputeRank();
660    assert(rank+1==width);
661
662    REformToRREform();
663
664    Vector<typ> ret(width);
665
666    typ diagonalProduct(1);
667    {
668      int pivot2I=-1;
669      int pivot2J=-1;
670      while(nextPivot(pivot2I,pivot2J))
671        {
672          diagonalProduct*=(*this)[pivot2I][pivot2J];
673        }
674    }
675    {
676      int j=nonPivotColumns().front();
677      int pivot2I=-1;
678      int pivot2J=-1;
679      ret[j]=typ(-1);
680      // Pretend that we are appending ret to the matrix, and reducing this
681      // new row by the previous ones. The last entry of the resulting matrix
682      // is lastEntry.
683      typ lastEntry=ret[j];
684      while(nextPivot(pivot2I,pivot2J))
685        {
686          ret[pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J];
687          lastEntry-=ret[pivot2J]*ret[pivot2J];
688        }
689      ret=(sign*(diagonalProduct*lastEntry))*ret;
690    }
691
692    return ret;
693  }
694
695  /**
696     Calls reduce() and returns the rank of the matrix.
697   */
698  int reduceAndComputeRank()
699  {
700    reduce(false,!typ::isField(),false);
701    int ret=0;
702    int pivotI=-1;
703    int pivotJ=-1;
704    while(nextPivot(pivotI,pivotJ))ret++;
705    return ret;
706  }
707  /**
708   * Sort the rows of the matrix.
709   */
710  struct rowComparer{
711    bool operator()(std::pair<Matrix*,int> i, std::pair<Matrix*,int> j) {return ((*i.first)[i.second].toVector()<(*j.first)[j.second].toVector());}
712  } theRowComparer;
713  void sortRows()
714  {
715          std::vector<std::pair<Matrix*,int> > v;
716          for(int i=0;i<height;i++)v.push_back(std::pair<Matrix*,int>(this,i));
717          std::sort(v.begin(),v.end(),theRowComparer);
718          Matrix result(height,width);
719          for(int i=0;i<height;i++)
720                  result[i]=(*this)[v[i].second].toVector();
721          data=result.data;
722  }
723  /**
724   * Sort the rows of the matrix and remove duplicate rows.
725   */
726  void sortAndRemoveDuplicateRows()
727  {
728    sortRows();
729    if(getHeight()==0)return;
730    Matrix B(0,getWidth());
731    B.appendRow((*this)[0]);
732    for(int i=1;i<getHeight();i++)
733      if((*this)[i].toVector()!=(*this)[i-1].toVector())B.appendRow((*this)[i].toVector());
734    *this=B;
735  }
736  /**
737     Takes two matrices with the same number of columns and construct
738     a new matrix which has the rows of the matrix top on the top and
739     the rows of the matrix bottom at the bottom. The return value is
740     the constructed matrix.
741   */
742  friend Matrix combineOnTop(Matrix const &top, Matrix const &bottom)
743  {
744    assert(bottom.getWidth()==top.getWidth());
745    Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth());
746    for(int i=0;i<top.getHeight();i++)ret[i]=top[i];
747    for(int i=0;i<bottom.getHeight();i++)ret[i+top.getHeight()]=bottom[i];
748
749    return ret;
750  }
751  /**
752     Takes two matrices with the same number of rows and construct
753     a new matrix which has the columns of the matrix left on the left and
754     the columns of the matrix right on the right. The return value is
755     the constructed matrix.
756   */
757  friend Matrix combineLeftRight(Matrix const &left, Matrix const &right)
758  {
759    assert(left.getHeight()==right.getHeight());
760    Matrix ret(left.getHeight(),left.getWidth()+right.getWidth());
761    for(int i=0;i<left.getHeight();i++)
762      {
763        for(int j=0;j<left.getWidth();j++)ret[i][j]=left[i][j];
764        for(int j=0;j<right.getWidth();j++)ret[i][j+left.getWidth()]=right[i][j];
765      }
766    return ret;
767  }
768};
769
770typedef Matrix<Integer> ZMatrix;
771typedef Matrix<Rational> QMatrix;
772typedef Matrix<int> IntMatrix;
773
774inline QMatrix ZToQMatrix(ZMatrix const &m)
775{
776  QMatrix ret(m.getHeight(),m.getWidth());
777  for(int i=0;i<m.getHeight();i++)ret[i]=ZToQVector(m[i]);
778  return ret;
779}
780
781inline ZMatrix QToZMatrixPrimitive(QMatrix const &m)
782{
783  ZMatrix ret(m.getHeight(),m.getWidth());
784  for(int i=0;i<m.getHeight();i++)ret[i]=QToZVectorPrimitive(m[i]);
785  return ret;
786}
787
788
789inline IntMatrix ZToIntMatrix(ZMatrix const &m)
790{
791  IntMatrix ret(m.getHeight(),m.getWidth());
792  for(int i=0;i<m.getHeight();i++)ret[i]=ZToIntVector(m[i]);
793  return ret;
794}
795
796
797inline ZMatrix IntToZMatrix(IntMatrix const &m)
798{
799  ZMatrix ret(m.getHeight(),m.getWidth());
800  for(int i=0;i<m.getHeight();i++)ret[i]=IntToZVector(m[i]);
801  return ret;
802}
803
804inline QMatrix canonicalizeSubspace(QMatrix const &m)
805{
806  QMatrix temp=m;
807  temp.reduce();
808  temp.REformToRREform();
809  temp.removeZeroRows();
810  return temp;
811}
812
813inline ZMatrix canonicalizeSubspace(ZMatrix const &m)
814{
815  return QToZMatrixPrimitive(canonicalizeSubspace(ZToQMatrix(m)));
816}
817
818
819inline QMatrix kernel(QMatrix const &m)
820{
821  QMatrix temp=m;
822  return temp.reduceAndComputeKernel();
823}
824
825inline ZMatrix kernel(ZMatrix const &m)
826{
827  return QToZMatrixPrimitive(kernel(ZToQMatrix(m)));
828}
829
830}
831
832
833#endif /* LIB_ZMATRIX_H_ */
Note: See TracBrowser for help on using the repository browser.