Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Path: gap4r8 / pkg / NormalizInterface-1.0.2 / Normaliz.git / DST / include / libQnormaliz / Qmatrix.h
Views: 418426/*1* Normaliz2* Copyright (C) 2007-2014 Winfried Bruns, Bogdan Ichim, Christof Soeger3* This program is free software: you can redistribute it and/or modify4* it under the terms of the GNU General Public License as published by5* the Free Software Foundation, either version 3 of the License, or6* (at your option) any later version.7*8* This program is distributed in the hope that it will be useful,9* but WITHOUT ANY WARRANTY; without even the implied warranty of10* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the11* GNU General Public License for more details.12*13* You should have received a copy of the GNU General Public License14* along with this program. If not, see <http://www.gnu.org/licenses/>.15*16* As an exception, when this program is distributed through (i) the App Store17* by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play18* by Google Inc., then that store may impose any digital rights management,19* device limits and/or redistribution restrictions that are required by its20* terms of service.21*/2223//---------------------------------------------------------------------------24#ifndef MATRIX_HPP25#define MATRIX_HPP26//---------------------------------------------------------------------------272829#include <vector>30#include <list>31#include <iostream>32#include <string>3334#include <libQnormaliz/libQnormaliz.h>35#include <libQnormaliz/Qinteger.h>36#include <libQnormaliz/Qconvert.h>3738//---------------------------------------------------------------------------3940namespace libQnormaliz {41using std::list;42using std::vector;43using std::string;4445template<typename Number> class Matrix {4647template<typename> friend class Matrix;48template<typename> friend class Lineare_Transformation;49template<typename> friend class Sublattice_Representation;5051// public:5253size_t nr;54size_t nc;55vector< vector<Number> > elem;5657//---------------------------------------------------------------------------58// Private routines, used in the public routines59//---------------------------------------------------------------------------6061//---------------------------------------------------------------------------62// Rows and columns exchange63//---------------------------------------------------------------------------6465void exchange_rows(const size_t& row1, const size_t& row2); //row1 is exchanged with row266void exchange_columns(const size_t& col1, const size_t& col2); // col1 is exchanged with col26768//---------------------------------------------------------------------------69// Row and column reduction70//---------------------------------------------------------------------------71// return value false undicates failure because of overflow72// for all the routines below7374// reduction via integer division and elemntary transformations75bool reduce_row(size_t corner); //reduction by the corner-th row76bool reduce_row (size_t row, size_t col); // corner at position (row,col)7778// replaces two columns by linear combinations of them79bool linear_comb_columns(const size_t& col,const size_t& j,80const Number& u,const Number& w,const Number& v,const Number& z);8182// column reduction with gcd method83bool gcd_reduce_column (size_t corner, Matrix<Number>& Right);8485//---------------------------------------------------------------------------86// Work horses87//---------------------------------------------------------------------------8889// takes |product of the diagonal elem|90Number compute_vol(bool& success);9192// Solve system with coefficients and right hand side in one matrix, using elementary transformations93// solution replaces right hand side94bool solve_destructive_inner(bool ZZinvertible, Number& denom);9596// asembles the matrix of the system (left side the submatrix of mother given by key97// right side from column vectors pointed to by RS98// both in a single matrix99void solve_system_submatrix_outer(const Matrix<Number>& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,100Number& denom, bool ZZ_invertible, bool transpose, size_t red_col, size_t sign_col,101bool compute_denom=true, bool make_sol_prime=false);102103size_t row_echelon_inner_elem(bool& success); // does the work and checks for overflows104// size_t row_echelon_inner_bareiss(bool& success, Number& det);105// NOTE: Bareiss cannot be used if z-invertible transformations are needed106107size_t row_echelon(bool& success); // transforms this into row echelon form and returns rank108size_t row_echelon(bool& success, Number& det); // computes also |det|109size_t row_echelon(bool& success, bool do_compute_vol, Number& det); // chooses elem (or bareiss in former time)110111// reduces the rows a matrix in row echelon form upwards, from left to right112bool reduce_rows_upwards();113size_t row_echelon_reduce(bool& success); // combines row_echelon and reduce_rows_upwards114115// computes rank and index simultaneously, returns rank116Number full_rank_index(bool& success);117118vector<key_t> max_rank_submatrix_lex_inner(bool& success) const;119120// A version of invert that circumvents protection and leaves it to the calling routine121Matrix invert_unprotected(Number& denom, bool& sucess) const;122123bool SmithNormalForm_inner(size_t& rk, Matrix<Number>& Right);124125126//---------------------------------------------------------------------------127// Pivots for rows/columns operations128//---------------------------------------------------------------------------129130vector<long> pivot(size_t corner); //Find the position of an element x with131//0<abs(x)<=abs(y) for all y!=0 in the right-lower submatrix of this132//described by an int corner133134long pivot_column(size_t col); //Find the position of an element x with135//0<abs(x)<=abs(y) for all y!=0 in the lower half of the column of this136//described by an int col137138long pivot_column(size_t row,size_t col); //in column col starting from row139140//---------------------------------------------------------------------------141// Helpers for linear systems142//---------------------------------------------------------------------------143144Matrix bundle_matrices(const Matrix<Number>& Right_side)const;145Matrix extract_solution() const;146vector<vector<Number>* > row_pointers();147void customize_solution(size_t dim, Number& denom, size_t red_col,148size_t sign_col, bool make_sol_prime);149150public:151152size_t row_echelon_inner_bareiss(bool& success, Number& det);153154vector<vector<Number>* > submatrix_pointers(const vector<key_t>& key);155156//---------------------------------------------------------------------------157158//---------------------------------------------------------------------------159// Construction and destruction160//---------------------------------------------------------------------------161162Matrix();163Matrix(size_t dim); //constructor of unit matrix164Matrix(size_t row, size_t col); //main constructor, all entries 0165Matrix(size_t row, size_t col, Number value); //constructor, all entries set to value166Matrix(const vector< vector<Number> >& elem); //constuctor, elem=elem167Matrix(const list< vector<Number> >& elems);168169//---------------------------------------------------------------------------170// Data access171//---------------------------------------------------------------------------172173void write_column(size_t col, const vector<Number>& data); //write a column174void print(const string& name, const string& suffix) const; // writes matrix into name.suffix175void print_append(const string& name,const string& suffix) const; // the same, but appends matrix176void print(std::ostream& out) const; // writes matrix to the stream177void pretty_print(std::ostream& out, bool with_row_nr=false) const; // writes matrix in a nice format to the stream // read a row178size_t nr_of_rows() const; // returns nr179size_t nr_of_columns() const; // returns nc180void set_nr_of_columns(size_t c);181/* generates a pseudo random matrix for tests, entries form 0 to mod-1 */182void random(int mod=3);183184void set_zero(); // sets all entries to 0185186/* returns a submatrix with rows corresponding to indices given by187* the entries of rows, Numbering from 0 to n-1 ! */188Matrix submatrix(const vector<key_t>& rows) const;189Matrix submatrix(const vector<int>& rows) const;190Matrix submatrix(const vector<bool>& rows) const;191192void swap (Matrix<Number>& x);193194// returns the permutation created by sorting the rows with a grading function195// or by 1-norm if computed is false196vector<key_t> perm_sort_by_degree(const vector<key_t>& key, const vector<Number>& grading, bool computed) const;197vector<key_t> perm_by_weights(const Matrix<Number>& Weights, vector<bool> absolute);198199void select_submatrix(const Matrix<Number>& mother, const vector<key_t>& rows);200void select_submatrix_trans(const Matrix<Number>& mother, const vector<key_t>& rows);201202Matrix& remove_zero_rows(); // remove zero rows, modifies this203204// resizes the matrix to the given number of rows/columns205// if the size shrinks it will keep all its allocated memory206// useful when the size varies207void resize(size_t nr_rows);208void resize(size_t nr_rows, size_t nr_cols);209void resize_columns(size_t nr_cols);210void Shrink_nr_rows(size_t new_nr_rows);211212vector<Number> diagonal() const; //returns the diagonale of this213//this should be a quadratic matrix214size_t maximal_decimal_length() const; //return the maximal number of decimals215//needed to write an entry216217vector<size_t> maximal_decimal_length_columnwise() const; // the same per column218219void append(const Matrix& M); // appends the rows of M to this220void append(const vector<vector<Number> >& M); // the same, but for another type of matrix221void append(const vector<Number>& v); // append the row v to this222void append_column(const vector<Number>& v); // append the column v to this223void remove_row(const vector<Number>& row); // removes all appearances of this row, not very efficient!224void remove_duplicate_and_zero_rows();225226inline const Number& get_elem(size_t row, size_t col) const {227return elem[row][col];228}229inline const vector< vector<Number> >& get_elements() const {230return elem;231}232inline vector<Number> const& operator[] (size_t row) const {233return elem[row];234}235inline vector<Number>& operator[] (size_t row) {236return elem[row];237}238void set_nc(size_t col){239nc=col;240}241void set_nr(size_t rows){242nc=rows;243}244245//---------------------------------------------------------------------------246// Basic matrices operations247//---------------------------------------------------------------------------248249Matrix add(const Matrix& A) const; // returns this+A250Matrix multiplication(const Matrix& A) const; // returns this*A251Matrix multiplication(const Matrix& A, long m) const;// returns this*A (mod m)252Matrix<Number> multiplication_cut(const Matrix<Number>& A, const size_t& c) const; // returns253// this*(first c columns of A)254bool equal(const Matrix& A) const; // returns this==A255bool equal(const Matrix& A, long m) const; // returns this==A (mod m)256Matrix transpose() const; // returns the transpose of this257258bool is_diagonal() const;259260//---------------------------------------------------------------------------261// Scalar operations262//---------------------------------------------------------------------------263264void scalar_multiplication(const Number& scalar); //this=this*scalar265void scalar_division(const Number& scalar);266//this=this div scalar, all the elem of this must be divisible with the scalar267void reduction_modulo(const Number& modulo); //this=this mod scalar268Number matrix_gcd() const; //returns the gcd of all elem269void simplify_rows(); //each row of this is reduced by its gcd,270// vector of gcds returned271void make_cols_prime(size_t from_col, size_t to_col);272// the columns of this in the specified range are reduced by their gcd273274Matrix multiply_rows(const vector<Number>& m) const; //returns matrix were row i is multiplied by m[i]275276//---------------------------------------------------------------------------277// Vector operations278//---------------------------------------------------------------------------279280void MxV(vector<Number>& result, const vector<Number>& v) const;//result = this*V281vector<Number> MxV(const vector<Number>& v) const;//returns this*V282vector<Number> VxM(const vector<Number>& v) const;//returns V*this283vector<Number> VxM_div(const vector<Number>& v, const Number& divisor,bool& success) const; // additionally divides by divisor284285//---------------------------------------------------------------------------286// Matrix operations287// --- these are more complicated algorithms ---288//---------------------------------------------------------------------------289290// Normal forms291292// converts this to row echelon form over ZZ and returns rank, GMP protected, uses only elementary transformations over ZZ293size_t row_echelon();294295// public version of row_echelon_reduce, GMP protected, uses only elementary transformations over ZZ296size_t row_echelon_reduce();297298// transforms matrix into lower triangular form via column transformations299// assumes that rk is the rank and that the matrix is zero after the first rk rows300// Right = Right*(column transformation of this call)301bool column_trigonalize(size_t rk, Matrix<Number>& Right);302303// combines row_echelon_reduce and column_trigonalize304// returns column transformation matrix305Matrix<Number> row_column_trigonalize(size_t& rk, bool& success);306307// rank and determinant308309size_t rank() const; //returns rank310Number full_rank_index() const; // returns index of full rank sublattice311size_t rank_submatrix(const vector<key_t>& key) const; //returns rank of submarix defined by key312313// returns rank of submatrix of mother. "this" is used as work space314size_t rank_submatrix(const Matrix<Number>& mother, const vector<key_t>& key);315316// vol stands for |det|317Number vol() const;318Number vol_submatrix(const vector<key_t>& key) const;319Number vol_submatrix(const Matrix<Number>& mother, const vector<key_t>& key);320321// find linearly indepenpendent submatrix of maximal rank322323vector<key_t> max_rank_submatrix_lex() const; //returns a vector with entries324//the indices of the first rows in lexicographic order of this forming325//a submatrix of maximal rank.326327// Solution of linear systems with square matrix328329// In the following routines, denom is the absolute value of the determinant of the330// left side matrix.331// If the diagonal is asked for, ZZ-invertible transformations are used.332// Otherwise ther is no restriction on the used algorithm333334//The diagonal of left hand side after transformation into an upper triangular matrix335//is saved in diagonal, denom is |determinant|.336337// System with "this" as left side338Matrix solve(const Matrix& Right_side, Number& denom) const;339Matrix solve(const Matrix& Right_side, vector< Number >& diagonal, Number& denom) const;340// solve the system this*Solution=denom*Right_side.341342// system is defined by submatrix of mother given by key (left side) and column vectors pointed to by RS (right side)343// NOTE: this is used as the matrix for the work344void solve_system_submatrix(const Matrix& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,345vector< Number >& diagonal, Number& denom, size_t red_col, size_t sign_col);346void solve_system_submatrix(const Matrix& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,347Number& denom, size_t red_col, size_t sign_col, bool compute_denom=true, bool make_sol_prime=false);348// the left side gets transposed349void solve_system_submatrix_trans(const Matrix& mother, const vector<key_t>& key, const vector<vector<Number>* >& RS,350Number& denom, size_t red_col, size_t sign_col);351352353// For non-square matrices354355// The next two solve routines do not require the matrix to be square.356// However, we want rank = number of columns, ensuring unique solvability357358vector<Number> solve_rectangular(const vector<Number>& v, Number& denom) const;359// computes solution vector for right side v, solution over the rationals360// matrix needs not be quadratic, but must have rank = number of columns361// with denominator denom.362// gcd of denom and solution is extracted !!!!!363364vector<Number> solve_ZZ(const vector<Number>& v) const;365// computes solution vector for right side v366// insists on integrality of the solution367368// homogenous linear systems369370Matrix<Number> kernel () const;371// computes a ZZ-basis of the solutions of (*this)x=0372// the basis is formed by the ROWS of the returned matrix373374// inverse matrix375376//this*Solution=denom*I. "this" should be a quadratic matrix with nonzero determinant.377Matrix invert(Number& denom) const;378379void invert_submatrix(const vector<key_t>& key, Number& denom, Matrix<Number>& Inv,380bool compute_denom=true, bool make_sol_prime=false) const;381382// find linear form that is constant on the rows383384vector<Number> find_linear_form () const;385// Tries to find a linear form which gives the same value an all rows of this386// this should be a m x n matrix (m>=n) of maxinal rank387// returns an empty vector if there does not exist such a linear form388389vector<Number> find_linear_form_low_dim () const;390//same as find_linear_form but also works with not maximal rank391//uses a linear transformation to get a full rank matrix392393// normal forms394395Matrix AlmostHermite(size_t& rk);396// Converts "this" into lower trigonal column Hermite normal form, returns column397// transformation matrix398// Almost: elements left of diagonal are not reduced mod diagonal399400// Computes Smith normal form and returns column transformation matrix401Matrix SmithNormalForm(size_t& rk);402403//for simplicial subcones404405// computes support hyperplanes and volume406void simplex_data(const vector<key_t>& key, Matrix<Number>& Supp, Number& vol, bool compute_vol) const;407408// Sorting of rows409410Matrix& sort_by_weights(const Matrix<Number>& Weights, vector<bool> absolute);411Matrix& sort_lex();412void order_rows_by_perm(const vector<key_t>& perm);413414// solve homogeneous congruences415416Matrix<Number> solve_congruences(bool& zero_modulus) const;417418// saturate sublattice419420void saturate();421422// find the indices of the rows in which the linear form L takes its max and min values423424vector<key_t> max_and_min(const vector<Number>& L, const vector<Number>& norm) const;425426// try to sort the rows in such a way that the extreme points of the polytope spanned by the rows come first427428size_t extreme_points_first(const vector<Number> norm=vector<Number>(0));429430// find an inner point in the cone spanned by the rows of the matrix431432vector<Number> find_inner_point();433};434//class end *****************************************************************435436//---------------------------------------------------------------------------437// Utilities438//---------------------------------------------------------------------------439440template<typename Number> class order_helper {441442public:443444vector<Number> weight;445key_t index;446vector<Number>* v;447};448449template<typename Number>450vector<vector<Number> > to_matrix(const vector<Number>& v){451452vector<vector<Number> > mat(1);453mat[0]=v;454return mat;455}456457template<typename Number>458Matrix<Number> readMatrix(const string project);459460//---------------------------------------------------------------------------461// Conversion between integer types462//---------------------------------------------------------------------------463464template<typename ToType, typename FromType>465void convert(Matrix<ToType>& to_mat, const Matrix<FromType>& from_mat);466467template<typename Number>468void mat_to_mpz(const Matrix<Number>& mat, Matrix<mpz_class>& mpz_mat);469470template<typename Number>471void mat_to_Int(const Matrix<mpz_class>& mpz_mat, Matrix<Number>& mat);472473template<typename Number>474void mpz_submatrix(Matrix<mpz_class>& sub, const Matrix<Number>& mother, const vector<key_t>& selection);475476template<typename Number>477void mpz_submatrix_trans(Matrix<mpz_class>& sub, const Matrix<Number>& mother, const vector<key_t>& selection);478479} // namespace480481//---------------------------------------------------------------------------482#endif483//---------------------------------------------------------------------------484485486