// Copyright (c) 1997-2000 // Utrecht University (The Netherlands), // ETH Zurich (Switzerland), // INRIA Sophia-Antipolis (France), // Max-Planck-Institute Saarbruecken (Germany), // and Tel-Aviv University (Israel). All rights reserved. // // This file is part of CGAL (www.cgal.org) // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Kernel_d/include/CGAL/Kernel_d/Linear_algebraHd_impl.h $ // $Id: Linear_algebraHd_impl.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial // // // Author(s) : Michael Seel //--------------------------------------------------------------------- // file generated by notangle from Linear_algebra.lw // please debug or modify noweb file // based on LEDA architecture by S. Naeher, C. Uhrig // coding: K. Mehlhorn, M. Seel // debugging and templatization: M. Seel //--------------------------------------------------------------------- #ifndef CGAL_LINEAR_ALGEBRAHD_C #define CGAL_LINEAR_ALGEBRAHD_C namespace CGAL { template bool Linear_algebraHd:: linear_solver(const Matrix& A, const Vector& b, Vector& x, RT& D, Matrix& spanning_vectors, Vector& c) { bool solvable = true; int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j= 0; i--) { RT_ h = C(i,cols) * D; for (j = i + 1; j < rank; j++) { h -= C(i,j)*x[var[j]]; } x[var[i]]= h / C(i,i); } #ifdef CGAL_LA_SELFTEST CGAL_assertion( (M*x).is_zero() ); #endif int defect = cols - rank; // dimension of kernel spanning_vectors = Matrix(cols,defect); if (defect > 0) { for(int l=0; l < defect; l++) { spanning_vectors(var[rank + l],l)=D; for(i = rank - 1; i >= 0 ; i--) { RT_ h = - C(i,rank + l)*D; for ( j= i + 1; j RT_ Linear_algebraHd:: determinant(const Matrix& A) { if (A.row_dimension() != A.column_dimension()) CGAL_error_msg( "determinant(): only square matrices are legal inputs."); Vector b(A.row_dimension()); // zero - vector int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j RT_ Linear_algebraHd:: determinant(const Matrix& A, Matrix& Ld, Matrix& Ud, std::vector& q, Vector& c) { if (A.row_dimension() != A.column_dimension()) CGAL_error_msg( "determinant(): only square matrices are legal inputs."); Vector b(A.row_dimension()); // zero - vector int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j int Linear_algebraHd:: sign_of_determinant(const Matrix& M) { return CGAL_NTS sign(determinant(M)); } template bool Linear_algebraHd:: verify_determinant(const Matrix& A, RT_ D, Matrix& L, Matrix& U, const std::vector& q, Vector& c) { if ((int)q.size() != A.column_dimension()) CGAL_error_msg( "verify_determinant: q should be a permutation array \ with index range [0,A.column_dimension() - 1]."); int n = A.row_dimension(); int i,j; if (D == 0) { /* we have $c^T \cdot A = 0$ */ Vector zero(n); return (transpose(A) * c == zero); } else { /* we check the conditions on |L| and |U| */ if (L(0,0) != 1) return false; for (i = 0; i 0 && L(i,i) != U(i - 1,i - 1)) return false; for (j = i + 1; j < n; j++) if (L(i,j) != 0) return false; } /* check whether $L \cdot A \cdot Q = U$ */ Matrix LA = L * A; for (j = 0; j < n; j++) if (LA.column(q[j]) != U.column(j)) return false; /* compute sign |s| of |Q| */ int sign = 1; /* we chase the cycles of |q|. An even length cycle contributes - 1 and vice versa */ std::vector already_considered(n); for (i = 0; i < n; i++) already_considered[i] = false; for (i = 0; i < n; i++) already_considered[q[i]] = true; for (i = 0; i < n; i++) if (! already_considered[i]) CGAL_error_msg("verify_determinant:q is not a permutation."); else already_considered[i] = false; for (i = 0; i < n; i++) { if (already_considered[i]) continue; /* we have found a new cycle with minimal element $i$. */ int k = q[i]; already_considered[i] =true; while (k != i) { sign = - sign; already_considered[k]= true; k = q[k]; } } return (D == RT(sign) * U(n - 1,n - 1)); } } template int Linear_algebraHd:: independent_columns(const Matrix& A, std::vector& columns) { Vector b(A.row_dimension()); // zero - vector int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j(rank); for(i = 0; i < rank; i++) columns[i] = var[i]; return rank; } template int Linear_algebraHd:: rank(const Matrix& A) { Vector b(A.row_dimension()); // zero - vector int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j bool Linear_algebraHd:: inverse(const Matrix& A, Matrix& inverse, RT_& D, Vector& c) { if (A.row_dimension() != A.column_dimension()) CGAL_error_msg("inverse: only square matrices are legal inputs."); Vector b(A.row_dimension()); // zero - vector int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j= 0; j--) { h = L (j,i) * D; for (int l = j + 1; l int Linear_algebraHd:: homogeneous_linear_solver(const Matrix &A, Matrix& spanning_vectors) /* returns the dimension of the solution space of the homogeneous system $Ax = 0$. The columns of spanning\_vectors span the solution space. */ { Vector b(A.row_dimension()); // zero - vector RT_ D; int i,j,k; // indices to step through the matrix int rows = A.row_dimension(); int cols = A.column_dimension(); /* at this point one might want to check whether the computation can be carried out with doubles, see section \ref{ optimization }. */ CGAL_assertion_msg(b.dimension() == rows, "linear_solver: b has wrong dimension"); Matrix C(rows,cols + 1); // the matrix in which we will calculate ($C = (A \vert b)$) /* copy |A| and |b| into |C| and L becomes the identity matrix */ Matrix L(rows); // zero initialized for(i=0; i var(cols); // column $j$ of |C| represents the |var[j]| - th variable // the array is indexed between |0| and |cols - 1| for(j=0; j= 0; i--) { RT_ h = C(i,cols) * D; for (j = i + 1; j < rank; j++) { h -= C(i,j)*x[var[j]]; } x[var[i]]= h / C(i,i); } #ifdef CGAL_LA_SELFTEST CGAL_assertion( (M*x).is_zero() ); #endif int defect = cols - rank; // dimension of kernel spanning_vectors = Matrix(cols,defect); if (defect > 0) { for(int l=0; l < defect; l++) { spanning_vectors(var[rank + l],l)=D; for(i = rank - 1; i >= 0 ; i--) { RT_ h = - C(i,rank + l)*D; for ( j= i + 1; j bool Linear_algebraHd:: homogeneous_linear_solver(const Matrix& A, Vector& x) /* returns true if the homogeneous system $Ax = 0$ has a non - trivial solution and false otherwise. */ { x = Vector(A.row_dimension()); Matrix spanning_vectors; int dimension = homogeneous_linear_solver(A,spanning_vectors); if (dimension == 0) return false; /* return first column of |spanning_vectors| */ for (int i = 0; i < spanning_vectors.row_dimension(); i++) x[i] = spanning_vectors(i,0); return true; } template typename Linear_algebraHd::Matrix Linear_algebraHd::transpose(const Matrix& M) { int d1 = M.row_dimension(); int d2 = M.column_dimension(); Matrix result(d2,d1); for(int i = 0; i < d2; i++) for(int j = 0; j < d1; j++) result(i,j) = M(j,i); return result; } } //namespace CGAL #endif //CGAL_LINEAR_ALGEBRAHD_C