File: matrix.h

package info (click to toggle)
regina-normal 5.1-6
  • links: PTS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 54,488 kB
  • sloc: cpp: 142,029; ansic: 19,218; xml: 9,844; objc: 7,729; perl: 1,190; python: 623; sh: 614; makefile: 34
file content (455 lines) | stat: -rw-r--r-- 20,849 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
/*
 * Normaliz
 * Copyright (C) 2007-2014  Winfried Bruns, Bogdan Ichim, Christof Soeger
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * As an exception, when this program is distributed through (i) the App Store
 * by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play
 * by Google Inc., then that store may impose any digital rights management,
 * device limits and/or redistribution restrictions that are required by its
 * terms of service.
 */

//---------------------------------------------------------------------------
#ifndef MATRIX_HPP
#define MATRIX_HPP
//---------------------------------------------------------------------------


#include <vector>
#include <list>
#include <iostream>
#include <string>

#include <libnormaliz/libnormaliz.h>
#include <libnormaliz/integer.h>
#include <libnormaliz/convert.h>

//---------------------------------------------------------------------------

namespace libnormaliz {
using std::list;
using std::vector;
using std::string;

template<typename Integer> class Matrix {

    template<typename> friend class Matrix;
    template<typename> friend class Lineare_Transformation;
    template<typename> friend class Sublattice_Representation;
    
    // public:

    size_t nr;
    size_t nc;
    vector< vector<Integer> > elem;

//---------------------------------------------------------------------------
//              Private routines, used in the public routines
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
//                      Rows and columns exchange
//---------------------------------------------------------------------------

    void exchange_rows(const size_t& row1, const size_t& row2);      //row1 is exchanged with row2
    void exchange_columns(const size_t& col1, const size_t& col2); // col1 is exchanged with col2

//---------------------------------------------------------------------------
//              Row and column reduction
//---------------------------------------------------------------------------
    // return value false undicates failure because of overflow
    // for all the routines below
    
    // reduction via integer division and elemntary transformations
    bool reduce_row(size_t corner);      //reduction by the corner-th row
    bool reduce_row (size_t row, size_t col); // corner at position (row,col)
            
    // replaces two columns by linear combinations of them
    bool linear_comb_columns(const size_t& col,const size_t& j,
            const Integer& u,const Integer& w,const Integer& v,const Integer& z);
    
    // column reduction with gcd method
    bool gcd_reduce_column (size_t corner, Matrix<Integer>& Right);
    
//---------------------------------------------------------------------------
//                      Work horses
//---------------------------------------------------------------------------

    // takes |product of the diagonal elem|
    Integer compute_vol(bool& success);  
        
    // Solve system with coefficients and right hand side in one matrix, using elementary transformations
    // solution replaces right hand side
    bool solve_destructive_inner(bool ZZinvertible, Integer& denom);

    // asembles the matrix of the system (left side the submatrix of mother given by key
    // right side from column vectors pointed to by RS
    // both in a single matrix    
    void solve_system_submatrix_outer(const Matrix<Integer>& mother, const vector<key_t>& key, const vector<vector<Integer>* >& RS,
         Integer& denom, bool ZZ_invertible, bool transpose, size_t red_col, size_t sign_col);
                    
    size_t row_echelon_inner_elem(bool& success); // does the work and checks for overflows
    // size_t row_echelon_inner_bareiss(bool& success, Integer& det);
    // NOTE: Bareiss cannot be used if z-invertible transformations are needed
    
    size_t row_echelon(bool& success); // transforms this into row echelon form and returns rank
    size_t row_echelon(bool& success, Integer& det); // computes also |det|
    size_t row_echelon(bool& success, bool do_compute_vol, Integer& det); // chooses elem or bareiss
    
    // reduces the rows a matrix in row echelon form upwards, from left to right
    bool reduce_rows_upwards();
    size_t row_echelon_reduce(bool& success); // combines row_echelon and reduce_rows_upwards
    
    // computes rank and index simultaneously, returns rank
    Integer full_rank_index(bool& success);
    
    vector<key_t> max_rank_submatrix_lex_inner(bool& success) const;
    
    // A version of invert that circumvents protection and leaves it to the calling routine
    Matrix invert_unprotected(Integer& denom, bool& sucess) const;
    
    bool SmithNormalForm_inner(size_t& rk, Matrix<Integer>& Right);
    

//---------------------------------------------------------------------------
//                      Pivots for rows/columns operations
//---------------------------------------------------------------------------

    vector<long> pivot(size_t corner); //Find the position of an element x with
    //0<abs(x)<=abs(y) for all y!=0 in the right-lower submatrix of this
    //described by an int corner

    long pivot_column(size_t col);  //Find the position of an element x with
    //0<abs(x)<=abs(y) for all y!=0 in the lower half of the column of this
    //described by an int col
    
    long pivot_column(size_t row,size_t col); //in column col starting from row
    
//---------------------------------------------------------------------------
//                     Helpers for linear systems
//---------------------------------------------------------------------------

    Matrix bundle_matrices(const Matrix<Integer>& Right_side)const;
    Matrix extract_solution() const;
    vector<vector<Integer>* > row_pointers();
                    
public:

size_t row_echelon_inner_bareiss(bool& success, Integer& det);

    vector<vector<Integer>* > submatrix_pointers(const vector<key_t>& key);     
  
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
//                      Construction and destruction
//---------------------------------------------------------------------------

    Matrix();
    Matrix(size_t dim);                           //constructor of unit matrix
    Matrix(size_t row, size_t col);                 //main constructor, all entries 0
    Matrix(size_t row, size_t col, Integer value); //constructor, all entries set to value
    Matrix(const vector< vector<Integer> >& elem); //constuctor, elem=elem
    Matrix(const list< vector<Integer> >& elems);

//---------------------------------------------------------------------------
//                             Data access
//---------------------------------------------------------------------------

    void write(std::istream& in = std::cin);                // to be modified, just for tests
    void write(size_t row, const vector<Integer>& data); //write a row
    void write(size_t row, const vector<int>& data); //write a row
    void write_column(size_t col, const vector<Integer>& data); //write a column
    void write(size_t row, size_t col, Integer data);  // write data at (row,col)
    void print(const string& name, const string& suffix) const;         //  writes matrix into name.suffix
    void print_append(const string& name,const string& suffix) const;  // the same, but appends matrix
    void print(std::ostream& out) const;          // writes matrix to the stream
    void pretty_print(std::ostream& out, bool with_row_nr=false) const;  // writes matrix in a nice format to the stream
    void read() const;                 // to be modified, just for tests
    vector<Integer> read(size_t row) const;                   // read a row
    Integer read (size_t row, size_t col) const;         // read data at (row,col)
    size_t nr_of_rows() const;                       // returns nr
    size_t nr_of_columns() const;                   // returns nc
    /* generates a pseudo random matrix for tests, entries form 0 to mod-1 */
    void random(int mod=3);

    void set_zero(); // sets all entries to 0

    /* returns a submatrix with rows corresponding to indices given by
     * the entries of rows, Numbering from 0 to n-1 ! */
    Matrix submatrix(const vector<key_t>& rows) const;
    Matrix submatrix(const vector<int>& rows) const;
    Matrix submatrix(const vector<bool>& rows) const;

	// returns the permutation created by sorting the rows with a grading function
    // or by 1-norm if computed is false
    vector<key_t> perm_sort_by_degree(const vector<key_t>& key, const vector<Integer>& grading, bool computed) const;
    vector<key_t> perm_by_weights(const Matrix<Integer>& Weights, vector<bool> absolute);
    
    void select_submatrix(const Matrix<Integer>& mother, const vector<key_t>& rows);
    void select_submatrix_trans(const Matrix<Integer>& mother, const vector<key_t>& rows);

    Matrix& remove_zero_rows(); // remove zero rows, modifies this

    // resizes the matrix to the given number of rows/columns
    // if the size shrinks it will keep all its allocated memory
    // useful when the size varies
    void resize(size_t nr_rows);
    void resize(size_t nr_rows, size_t nr_cols);
    void resize_columns(size_t nr_cols);

    vector<Integer> diagonal() const;     //returns the diagonale of this
                                  //this should be a quadratic matrix
    size_t maximal_decimal_length() const;    //return the maximal number of decimals
                                      //needed to write an entry
                                      
    vector<size_t> maximal_decimal_length_columnwise() const; // the same per column

    void append(const Matrix& M); // appends the rows of M to this
    void append(const vector<vector<Integer> >& M); // the same, but for another type of matrix
    void append(const vector<Integer>& v); // append the row v to this
    void append_column(const vector<Integer>& v); // append the column v to this
    void remove_row(const vector<Integer>& row); // removes all appearances of this row, not very efficient!
    void remove_duplicate_and_zero_rows();

    inline const Integer& get_elem(size_t row, size_t col) const {
        return elem[row][col];
    }
    inline const vector< vector<Integer> >& get_elements() const {
        return elem;
    }
    inline vector<Integer> const& operator[] (size_t row) const {
        return elem[row];
    }
    inline vector<Integer>& operator[] (size_t row) { 
        return elem[row];
    }
    void set_nc(size_t col){
        nc=col;
    }
    void set_nr(size_t rows){
        nc=rows;
    }

//---------------------------------------------------------------------------
//                  Basic matrices operations
//---------------------------------------------------------------------------

    Matrix add(const Matrix& A) const;                       // returns this+A
    Matrix multiplication(const Matrix& A) const;          // returns this*A
    Matrix multiplication(const Matrix& A, long m) const;// returns this*A (mod m)
    Matrix<Integer> multiplication_cut(const Matrix<Integer>& A, const size_t& c) const; // returns 
    // this*(first c columns of A)
    bool equal(const Matrix& A) const;             // returns this==A
    bool equal(const Matrix& A, long m) const;     // returns this==A (mod m)
    Matrix transpose() const;                     // returns the transpose of this
    
    bool is_diagonal() const;

//---------------------------------------------------------------------------
//                          Scalar operations
//---------------------------------------------------------------------------

    void scalar_multiplication(const Integer& scalar);  //this=this*scalar
    void scalar_division(const Integer& scalar);
    //this=this div scalar, all the elem of this must be divisible with the scalar
    void reduction_modulo(const Integer& modulo);     //this=this mod scalar
    Integer matrix_gcd() const; //returns the gcd of all elem
    vector<Integer> make_prime();         //each row of this is reduced by its gcd
    //return a vector containing the gcd of the rows

    Matrix multiply_rows(const vector<Integer>& m) const;  //returns matrix were row i is multiplied by m[i]

//---------------------------------------------------------------------------
//                          Vector operations
//---------------------------------------------------------------------------

   void MxV(vector<Integer>& result, const vector<Integer>& v) const;//result = this*V
   vector<Integer> MxV(const vector<Integer>& v) const;//returns this*V
   vector<Integer> VxM(const vector<Integer>& v) const;//returns V*this
   vector<Integer> VxM_div(const vector<Integer>& v, const Integer& divisor,bool& success) const; // additionally divides by divisor

//---------------------------------------------------------------------------
//                          Matrix operations
//           --- these are more complicated algorithms ---
//---------------------------------------------------------------------------

// Normal forms

    // converts this to row echelon form over QQ and returns rank, GMP protected
    size_t row_echelon();

    // public version of row_echelon_reduce (), GMP protected, uses only elementary transformations over ZZ
    void row_echelon_reduce();

    // transforms matrix into lower triangular form via column transformations
    // assumes that rk is the rank and that the matrix is zero after the first rk rows
    // Right = Right*(column transformation of this call)
    bool column_trigonalize(size_t rk, Matrix<Integer>& Right);
    
    // combines row_echelon_reduce and column_trigonalize
    // returns column transformation matrix
    Matrix<Integer> row_column_trigonalize(size_t& rk, bool& success);
    
// rank and determinant

    size_t rank() const; //returns rank
    Integer full_rank_index() const; // returns index of full rank sublattice
    size_t rank_submatrix(const vector<key_t>& key) const; //returns rank of submarix defined by key
    
    // returns rank of submatrix of mother. "this" is used as work space    
    size_t rank_submatrix(const Matrix<Integer>& mother, const vector<key_t>& key);
 
    // vol stands for |det|
    Integer vol() const;
    Integer vol_submatrix(const vector<key_t>& key) const;
    Integer vol_submatrix(const Matrix<Integer>& mother, const vector<key_t>& key);
    
// find linearly indepenpendent submatrix of maximal rank

    vector<key_t>  max_rank_submatrix_lex() const; //returns a vector with entries
    //the indices of the first rows in lexicographic order of this forming
    //a submatrix of maximal rank.
    
// Solution of linear systems with square matrix
  
    // In the following routines, denom is the absolute value of the determinant of the
    // left side matrix.
    // If the diagonal is asked for, ZZ-invertible transformations are used.
    // Otherwise ther is no restriction on the used algorithm
    
    //The diagonal of left hand side after transformation into an upper triangular matrix
    //is saved in diagonal, denom is |determinant|.
    
    // System with "this" as left side
    Matrix solve(const Matrix& Right_side, Integer& denom) const;
    Matrix solve(const Matrix& Right_side, vector< Integer >& diagonal, Integer& denom) const;
    // solve the system this*Solution=denom*Right_side. 

    // system is defined by submatrix of mother given by key (left side) and column vectors pointed to by RS (right side)
    // NOTE: this is used as the matrix for the woek     
    void solve_system_submatrix(const Matrix& mother, const vector<key_t>& key, const vector<vector<Integer>* >& RS,
         vector< Integer >& diagonal, Integer& denom, size_t red_col, size_t sign_col);
    void solve_system_submatrix(const Matrix& mother, const vector<key_t>& key, const vector<vector<Integer>* >& RS,
         Integer& denom, size_t red_col, size_t sign_col);
    // the left side gets transposed
    void solve_system_submatrix_trans(const Matrix& mother, const vector<key_t>& key, const vector<vector<Integer>* >& RS,
         Integer& denom, size_t red_col, size_t sign_col);
        
                    
// For non-square matrices
                    
    // The next two solve routines do not require the matrix to be square.
    // However, we want rank = number of columns, ensuring unique solvability
    
    vector<Integer> solve_rectangular(const vector<Integer>& v, Integer& denom) const;
    // computes solution vector for right side v, solution over the rationals
    // matrix needs not be quadratic, but must have rank = number of columns
    // with denominator denom. 
    // gcd of denom and solution is extracted !!!!!
    
    vector<Integer> solve_ZZ(const vector<Integer>& v) const;
    // computes solution vector for right side v
    // insists on integrality of the solution
                    
// homogenous linear systems

    Matrix<Integer> kernel () const;
    // computes a ZZ-basis of the solutions of (*this)x=0
    // the basis is formed by the ROWS of the returned matrix
                    
// inverse matrix
                    
    //this*Solution=denom*I. "this" should be a quadratic matrix with nonzero determinant. 
    Matrix invert(Integer& denom) const;
    
    void invert_submatrix(const vector<key_t>& key, Integer& denom, Matrix<Integer>& Inv) const;
                    
// find linear form that is constant on the rows 

    vector<Integer> find_linear_form () const;
    // Tries to find a linear form which gives the same value an all rows of this
    // this should be a m x n matrix (m>=n) of maxinal rank
    // returns an empty vector if there does not exist such a linear form
  
    vector<Integer> find_linear_form_low_dim () const;
    //same as find_linear_form but also works with not maximal rank
    //uses a linear transformation to get a full rank matrix
    
// normal forms
        
    Matrix AlmostHermite(size_t& rk);
    // Converts "this" into lower trigonal column Hermite normal form, returns column 
    // transformation matrix
    // Almost: elements left of diagonal are not reduced mod diagonal 
    
    // Computes Smith normal form and returns column transformation matrix
    Matrix SmithNormalForm(size_t& rk);
    
//for simplicial subcones

    // computes support hyperplanes and volume
    void simplex_data(const vector<key_t>& key, Integer& vol, Matrix& Supp) const; 
    
// Sorting of rows
    
    Matrix& sort_by_weights(const Matrix<Integer>& Weights, vector<bool> absolute);
    Matrix& sort_lex();
    void order_rows_by_perm(const vector<key_t>& perm);
    
// solve homogeneous congruences
    
    Matrix<Integer> solve_congruences(bool& zero_modulus) const;
    
// saturate sublattice
    
    void saturate();

};
//class end *****************************************************************

template<typename Integer> class order_helper {
    
public:
    
    vector<Integer> weight;
    key_t index;
    vector<Integer>* v;
};

//---------------------------------------------------------------------------
//                  Conversion between integer types
//---------------------------------------------------------------------------

template<typename ToType, typename FromType>
void convert(Matrix<ToType>& to_mat, const Matrix<FromType>& from_mat);

template<typename Integer>
void mat_to_mpz(const Matrix<Integer>& mat, Matrix<mpz_class>& mpz_mat);

template<typename Integer>
void mat_to_Int(const Matrix<mpz_class>& mpz_mat, Matrix<Integer>& mat);

template<typename Integer>
void mpz_submatrix(Matrix<mpz_class>& sub, const Matrix<Integer>& mother, const vector<key_t>& selection);

template<typename Integer>
void mpz_submatrix_trans(Matrix<mpz_class>& sub, const Matrix<Integer>& mother, const vector<key_t>& selection);

} // namespace

//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------