From 86502db00a85ae05d4656d7c43bb56bbedd5d992 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 20 Oct 2023 14:44:17 +0100 Subject: [PATCH] Parameterization: cleanup --- OpenNL/include/CGAL/OpenNL/bicgstab.h | 332 ------------------ OpenNL/include/CGAL/OpenNL/blas.h | 67 ---- .../include/CGAL/OpenNL/conjugate_gradient.h | 240 ------------- OpenNL/include/CGAL/OpenNL/full_vector.h | 148 -------- OpenNL/include/CGAL/OpenNL/linear_solver.h | 179 ---------- OpenNL/include/CGAL/OpenNL/preconditioner.h | 214 ----------- OpenNL/include/CGAL/OpenNL/sparse_matrix.h | 256 -------------- .../LSCM_parameterizer_3.h | 7 +- 8 files changed, 1 insertion(+), 1442 deletions(-) delete mode 100644 OpenNL/include/CGAL/OpenNL/bicgstab.h delete mode 100644 OpenNL/include/CGAL/OpenNL/blas.h delete mode 100644 OpenNL/include/CGAL/OpenNL/conjugate_gradient.h delete mode 100644 OpenNL/include/CGAL/OpenNL/full_vector.h delete mode 100644 OpenNL/include/CGAL/OpenNL/preconditioner.h delete mode 100644 OpenNL/include/CGAL/OpenNL/sparse_matrix.h diff --git a/OpenNL/include/CGAL/OpenNL/bicgstab.h b/OpenNL/include/CGAL/OpenNL/bicgstab.h deleted file mode 100644 index a9cef2c92e8..00000000000 --- a/OpenNL/include/CGAL/OpenNL/bicgstab.h +++ /dev/null @@ -1,332 +0,0 @@ -// Copyright (c) 2005-2008 Inria Loria (France). -/* - * author: Bruno Levy, INRIA, project ALICE - * website: https://www.loria.fr/~levy/software - * - * This file is part of CGAL (www.cgal.org) - * - * Scientific work that use this software can reference the website and - * the following publication: - * - * @INPROCEEDINGS {levy:NMDGP:05, - * AUTHOR = Bruno Levy, - * TITLE = Numerical Methods for Digital Geometry Processing, - * BOOKTITLE =Israel Korea Bi-National Conference, - * YEAR=November 2005, - * URL=https://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics - * } - * - * Laurent Saboret 2005-2008: Changes for CGAL: - * - Added OpenNL namespace - * - solve() returns true on success - * - check divisions by zero - * - added comments and traces - * - copied BICGSTAB algorithm WITH preconditioner from Graphite 1.9 code - * - * $URL$ - * $Id$ - * SPDX-License-Identifier: LGPL-3.0-or-later - */ - -#ifndef __OPENNL_BICGSTAB__ -#define __OPENNL_BICGSTAB__ - -#include -#include - -#include -#include -#include -#include - -namespace OpenNL { - - -/** - * The BICGSTAB algorithm without preconditioner: - * Ashby, Manteuffel, Saylor - * A taxononmy for conjugate gradient methods - * SIAM J Numer Anal 27, 1542-1568 (1990) - * - * This implementation is inspired by the lsolver library, - * by Christian Badura, available from: - * http://www.mathematik.uni-freiburg.de/IAM/Research/projectskr/lin_solver/ - * - * @param A generic square matrix; a function - * mult(const MATRIX& M, const double* x, double* y) - * and a member function - * int dimension() const - * must to be defined. - * @param b right hand side of the system. - * @param x initial value. - * @param eps threshold for the residual. - * @param max_iter maximum number of iterations. - */ - -template class Solver_BICGSTAB { -public: - typedef MATRIX Matrix ; - typedef VECTOR Vector ; - typedef typename Vector::CoeffType CoeffType ; - -public: - Solver_BICGSTAB() { - epsilon_ = 1e-6 ; - max_iter_ = 0 ; - } - - // Default copy constructor, operator =() and destructor are fine - - void set_epsilon(CoeffType eps) { epsilon_ = eps ; } - void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; } - - // Solve the sparse linear system "A*x = b". Return true on success. - bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) - { -#ifdef DEBUG_TRACE - std::cerr << " Call BICGSTAB" << std::endl; -#endif - CGAL_assertion(A.dimension() > 0); - unsigned int n = A.dimension() ; // (Square) matrix dimension - - unsigned int max_iter = max_iter_ ; // Max number of iterations - if(max_iter == 0) { - max_iter = 5 * n ; - } - - Vector rT(n) ; // Initial residue rT=Ax-b - Vector d(n) ; - Vector h(n) ; - Vector u(n) ; - Vector Ad(n) ; - Vector t(n) ; - Vector& s = h ; - CoeffType rTh, rTAd, rTr, alpha, beta, omega, st, tt; - unsigned int its=0; // Loop counter - CoeffType err=epsilon_*epsilon_*BLAS::dot(b,b); // Error to reach - Vector r(n) ; // residue r=A*x-b - mult(A,x,r); - BLAS::axpy(-1,b,r); - BLAS::copy(r,d); - BLAS::copy(d,h); - BLAS::copy(h,rT); - //CGAL_assertion(BLAS::dot(rT,rT) == 0.0); // may happen for small systems - rTh=BLAS::dot(rT,h); // rTh = (rT|h) - rTr=BLAS::dot(r,r); // error rTr = (r|r) - - while ( rTr>err && its < max_iter) { - mult(A,d,Ad); - rTAd=BLAS::dot(rT,Ad); - CGAL_assertion(rTAd != 0.0); - alpha=rTh/rTAd; - BLAS::axpy(-alpha,Ad,r); - BLAS::copy(h,s); - BLAS::axpy(-alpha,Ad,s); - mult(A,s,t); - BLAS::axpy(1,t,u); - BLAS::scal(alpha,u); - st=BLAS::dot(s,t); - tt=BLAS::dot(t,t); - if (st == 0.0 || tt == 0.0) - omega = 0 ; - else - omega = st/tt; - BLAS::axpy(-omega,t,r); - BLAS::axpy(-alpha,d,x); - BLAS::axpy(-omega,s,x); - rTr=BLAS::dot(r,r); - BLAS::copy(s,h); - BLAS::axpy(-omega,t,h); - if (omega == 0.0) // stop if omega==0 (success) - { -#ifdef DEBUG_TRACE - std::cerr << " BICGSTAB: omega=0" << std::endl; -#endif - break; - } - if (rTh == 0.0) // stop if rTh==0 (failure) - { -#ifdef DEBUG_TRACE - std::cerr << " BICGSTAB: rTh=0" << std::endl; -#endif - break; - } - beta=(alpha/omega)/rTh; // beta = (rTh/"old rTh") * (alpha/omega) - rTh=BLAS::dot(rT,h); - beta*=rTh; - BLAS::scal(beta,d); - BLAS::axpy(1,h,d); - BLAS::axpy(-beta*omega,Ad,d); - ++its; - } - - bool success = (rTr <= err); - return success; - } - -private: - // Test if a floating point number is (close to) 0.0 - static bool IsZero(CoeffType a) - { - return (std::fabs(a) < 10.0 * (std::numeric_limits::min)()); - } - -private: - CoeffType epsilon_ ; - unsigned int max_iter_ ; -} ; - - -/** - * The BICGSTAB algorithm WITH preconditioner: - * Ashby, Manteuffel, Saylor - * A taxononmy for conjugate gradient methods - * SIAM J Numer Anal 27, 1542-1568 (1990) - * - * This implementation is inspired by the lsolver library, - * by Christian Badura, available from: - * http://www.mathematik.uni-freiburg.de/IAM/Research/projectskr/lin_solver/ - * - * @param A generic square matrix; a function - * mult(const MATRIX& M, const double* x, double* y) - * and a member function - * int dimension() const - * must to be defined. - * @param C preconditioner; a function - * mult(const PC_MATRIX& C, const double* x, double* y) - * needs to be defined. - * @param b right hand side of the system. - * @param x initial value. - * @param eps threshold for the residual. - * @param max_iter maximum number of iterations. - */ - -template< class MATRIX, class PC_MATRIX, class VECTOR > -class Solver_preconditioned_BICGSTAB -{ -public: - typedef MATRIX Matrix ; - typedef PC_MATRIX Preconditioner ; - typedef VECTOR Vector ; - typedef typename Vector::CoeffType CoeffType ; - -public: - Solver_preconditioned_BICGSTAB() { - epsilon_ = 1e-6 ; - max_iter_ = 0 ; - } - - // Default copy constructor, operator =() and destructor are fine - - void set_epsilon(CoeffType eps) { epsilon_ = eps ; } - void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; } - - // Solve the sparse linear system "A*x = b". Return true on success. - bool solve(const MATRIX &A, const PC_MATRIX &C, const VECTOR& b, VECTOR& x) - { -#ifdef DEBUG_TRACE - std::cerr << " Call BICGSTAB with preconditioner" << std::endl; -#endif - CGAL_assertion(A.dimension() > 0); - unsigned int n = A.dimension() ; // (Square) matrix dimension - - unsigned int max_iter = max_iter_ ; // Max number of iterations - if(max_iter == 0) { - max_iter = 5 * n ; - } - - Vector rT(n) ; // Initial residue rT=Ax-b - Vector d(n) ; - Vector h(n) ; - Vector u(n) ; - Vector Sd(n) ; - Vector t(n) ; - Vector aux(n) ; - Vector& s = h ; - CoeffType rTh, rTSd, rTr, alpha, beta, omega, st, tt; - unsigned int its=0; // Loop counter - CoeffType err=epsilon_*epsilon_*BLAS::dot(b,b); // Error to reach - Vector r(n) ; // residue r=A*x-b - mult(A,x,r); - BLAS::axpy(-1,b,r); - mult(C,r,d); - BLAS::copy(d,h); - BLAS::copy(h,rT); - //CGAL_assertion(BLAS::dot(rT,rT) == 0.0); // may happen for small systems - rTh=BLAS::dot(rT,h); // rTh = (rT|h) - rTr=BLAS::dot(r,r); // error rTr = (r|r) - - while (rTr>err && its < max_iter) { - mult(A,d,aux); - mult(C,aux,Sd); - rTSd=BLAS::dot(rT,Sd); - if (rTSd == 0.0) // stop if rTSd==0 (failure) - { -#ifdef DEBUG_TRACE - std::cerr << " BICGSTAB with preconditioner: rTSd=0" << std::endl; -#endif - break; - } - alpha=rTh/rTSd; - BLAS::axpy(-alpha,aux,r); - BLAS::copy(h,s); - BLAS::axpy(-alpha,Sd,s); - mult(A,s,aux); - mult(C,aux,t); - BLAS::axpy(1,t,u); - BLAS::scal(alpha,u); - st=BLAS::dot(s,t); - tt=BLAS::dot(t,t); - if (st == 0.0 || tt == 0.0) - omega = 0 ; - else - omega = st/tt; - BLAS::axpy(-omega,aux,r); - BLAS::axpy(-alpha,d,x); - BLAS::axpy(-omega,s,x); - rTr=BLAS::dot(r,r); - BLAS::copy(s,h); - BLAS::axpy(-omega,t,h); - if (omega == 0.0) // stop if omega==0 (success) - { -#ifdef DEBUG_TRACE - std::cerr << " BICGSTAB with preconditioner: omega=0" << std::endl; -#endif - break; - } - if (rTh == 0.0) // stop if rTh==0 (failure) - { -#ifdef DEBUG_TRACE - std::cerr << " BICGSTAB with preconditioner: rTh=0" << std::endl; -#endif - break; - } - beta=(alpha/omega)/rTh; // beta = (rTh/"old rTh") * (alpha/omega) - rTh=BLAS::dot(rT,h); - beta*=rTh; - BLAS::scal(beta,d); // d = h + beta * (d - omega * Sd); - BLAS::axpy(1,h,d); - BLAS::axpy(-beta*omega,Sd,d); - ++its; - } - - bool success = (rTr <= err); - return success; - } - -private: - // Test if a floating point number is (close to) 0.0 - static bool IsZero(CoeffType a) - { - return (std::fabs(a) < 10.0 * (std::numeric_limits::min)()); - } - -private: - CoeffType epsilon_ ; - unsigned int max_iter_ ; -} ; - - -} // namespace OpenNL - -#endif // __OPENNL_BICGSTAB__ diff --git a/OpenNL/include/CGAL/OpenNL/blas.h b/OpenNL/include/CGAL/OpenNL/blas.h deleted file mode 100644 index 8cea9da67ad..00000000000 --- a/OpenNL/include/CGAL/OpenNL/blas.h +++ /dev/null @@ -1,67 +0,0 @@ -// Copyright (c) 2005-2008 Inria Loria (France). -/* - * author: Bruno Levy, INRIA, project ALICE - * website: https://www.loria.fr/~levy/software - * - * This file is part of CGAL (www.cgal.org) - * - * Scientific work that use this software can reference the website and - * the following publication: - * - * @INPROCEEDINGS {levy:NMDGP:05, - * AUTHOR = Bruno Levy, - * TITLE = Numerical Methods for Digital Geometry Processing, - * BOOKTITLE =Israel Korea Bi-National Conference, - * YEAR=November 2005, - * URL=https://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics - * } - * - * Laurent Saboret 01/2005: Change for CGAL: - * - Added OpenNL namespace - * Andreas Meyer 2007 changes for CGAL: - * - replaced assert with CGAL_assertion/CGAL_error etc. - * - * $URL$ - * $Id$ - * SPDX-License-Identifier: LGPL-3.0-or-later -*/ - -#ifndef __OPENNL_BLAS__ -#define __OPENNL_BLAS__ - -#include - -namespace OpenNL { - - -/** Basic Linear Algebra Subroutines */ -template class BLAS { -public: - typedef VECTOR VectorType ; - typedef typename VECTOR::CoeffType CoeffType ; - - /** y <- y + a*x */ - static void axpy(CoeffType /*a*/, const VectorType& /*x*/, VectorType& /*y*/) { - CGAL_error(); - } - - /** x <- a*x */ - static void scal(CoeffType /*a*/, VectorType& /*x*/) { - CGAL_error(); - } - - /** y <- x */ - static void copy(const VectorType& /*x*/, VectorType& /*y*/) { - CGAL_error(); - } - - /** returns x^t * y */ - static CoeffType dot(const VectorType& /*x*/, const VectorType& /*y*/) { - CGAL_error(); - } -} ; - - -} // namespace OpenNL - -#endif diff --git a/OpenNL/include/CGAL/OpenNL/conjugate_gradient.h b/OpenNL/include/CGAL/OpenNL/conjugate_gradient.h deleted file mode 100644 index 6f2e6f5b2e8..00000000000 --- a/OpenNL/include/CGAL/OpenNL/conjugate_gradient.h +++ /dev/null @@ -1,240 +0,0 @@ -// Copyright (c) 2005-2008 Inria Loria (France). -/* - * author: Bruno Levy, INRIA, project ALICE - * website: https://www.loria.fr/~levy/software - * - * This file is part of CGAL (www.cgal.org) - * - * Scientific work that use this software can reference the website and - * the following publication: - * - * @INPROCEEDINGS {levy:NMDGP:05, - * AUTHOR = Bruno Levy, - * TITLE = Numerical Methods for Digital Geometry Processing, - * BOOKTITLE =Israel Korea Bi-National Conference, - * YEAR=November 2005, - * URL=https://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics - * } - * - * Laurent Saboret 2005-2006: Changes for CGAL: - * - Added OpenNL namespace - * - solve() returns true on success - * - check divisions by zero - * - added comments - * - copied Conjugate Gradient algorithm WITH preconditioner from Graphite 1.9 code - * - * $URL$ - * $Id$ - * SPDX-License-Identifier: LGPL-3.0-or-later - */ - -#ifndef __OPENNL_CONJUGATE_GRADIENT__ -#define __OPENNL_CONJUGATE_GRADIENT__ - -#include -#include - -#include -#include -#include - -namespace OpenNL { - - -/** - * The Conjugate Gradient algorithm WITHOUT preconditioner: - * Ashby, Manteuffel, Saylor - * A taxononmy for conjugate gradient methods - * SIAM J Numer Anal 27, 1542-1568 (1990) - * - * This implementation is inspired by the lsolver library, - * by Christian Badura, available from: - * http://www.mathematik.uni-freiburg.de/IAM/Research/projectskr/lin_solver/ - * - * @param A generic square matrix; a function - * mult(const MATRIX& M, const double* x, double* y) - * and a member function - * int dimension() const - * must to be defined. - * @param b right hand side of the system. - * @param x initial value. - * @param eps threshold for the residual. - * @param max_iter maximum number of iterations. - */ - -template class Solver_CG { -public: - typedef MATRIX Matrix ; - typedef VECTOR Vector ; - typedef typename Vector::CoeffType CoeffType ; - -public: - Solver_CG() { - epsilon_ = 1e-6 ; - max_iter_ = 0 ; - } - - // Default copy constructor, operator =() and destructor are fine - - void set_epsilon(CoeffType eps) { epsilon_ = eps ; } - void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; } - - // Solve the sparse linear system "A*x = b" for A symmetric positive definite - // Return true on success - bool solve(const MATRIX &A, const VECTOR& b, VECTOR& x) - { -#ifdef DEBUG_TRACE - std::cerr << " Call Conjugate Gradient" << std::endl; -#endif - CGAL_assertion(A.dimension() > 0); - unsigned int n = A.dimension() ; // (Square) matrix dimension - - unsigned int max_iter = max_iter_ ; // Max number of iterations - if(max_iter == 0) { - max_iter = 5 * n ; - } - - Vector g(n) ; - Vector r(n) ; - Vector p(n) ; - unsigned int its=0; // Loop counter - CoeffType t, tau, sig, rho, gam; - CoeffType bnorm2 = BLAS::dot(b,b) ; - CoeffType err=epsilon_*epsilon_*bnorm2 ; // Error to reach - // residue g=b-A*x - mult(A,x,g); - BLAS::axpy(-1,b,g); - BLAS::scal(-1,g); - // Initially, r=g=b-A*x - BLAS::copy(g,r); // r = g - - CoeffType gg=BLAS::dot(g,g); // error gg = (g|g) - while ( gg>err && its < max_iter) { - mult(A,r,p); - rho=BLAS::dot(p,p); - sig=BLAS::dot(r,p); - tau=BLAS::dot(g,r); - CGAL_assertion(sig != 0.0); - t=tau/sig; - BLAS::axpy(t,r,x); - BLAS::axpy(-t,p,g); - CGAL_assertion(tau != 0.0); - gam=(t*t*rho-tau)/tau; - BLAS::scal(gam,r); - BLAS::axpy(1,g,r); - gg=BLAS::dot(g,g); // Update error gg = (g|g) - ++its; - } - - bool success = (gg <= err); - return success; - } - -private: - CoeffType epsilon_ ; - unsigned int max_iter_ ; -} ; - - -/** - * The Conjugate Gradient algorithm WITH preconditioner: - * Ashby, Manteuffel, Saylor - * A taxononmy for conjugate gradient methods - * SIAM J Numer Anal 27, 1542-1568 (1990) - * - * This implementation is inspired by the lsolver library, - * by Christian Badura, available from: - * http://www.mathematik.uni-freiburg.de/IAM/Research/projectskr/lin_solver/ - * - * @param A generic square matrix; a function - * mult(const MATRIX& M, const double* x, double* y) - * and a member function - * int dimension() const - * must to be defined. - * @param C preconditioner; a function - * mult(const PC_MATRIX& C, const double* x, double* y) - * needs to be defined. - * @param b right hand side of the system. - * @param x initial value. - * @param eps threshold for the residual. - * @param max_iter maximum number of iterations. - */ - -template< class MATRIX, class PC_MATRIX, class VECTOR > -class Solver_preconditioned_CG -{ -public: - typedef MATRIX Matrix ; - typedef PC_MATRIX Preconditioner ; - typedef VECTOR Vector ; - typedef typename Vector::CoeffType CoeffType ; - -public: - Solver_preconditioned_CG() { - epsilon_ = 1e-6 ; - max_iter_ = 0 ; - } - - // Default copy constructor, operator =() and destructor are fine - - void set_epsilon(CoeffType eps) { epsilon_ = eps ; } - void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; } - - // Solve the sparse linear system "A*x = b" for A symmetric positive definite - // Return true on success - bool solve(const MATRIX &A, const PC_MATRIX &C, const VECTOR& b, VECTOR& x) - { -#ifdef DEBUG_TRACE - std::cerr << " Call Conjugate Gradient with preconditioner" << std::endl; -#endif - CGAL_assertion(A.dimension() > 0); - unsigned int n = A.dimension() ; // (Square) matrix dimension - - unsigned int max_iter = max_iter_ ; // Max number of iterations - if(max_iter == 0) { - max_iter = 5 * n ; - } - - Vector r(n) ; // residue r=Ax-b - Vector d(n) ; - Vector h(n) ; - Vector Ad = h ; - unsigned int its=0; // Loop counter - CoeffType rh, alpha, beta; - CoeffType bnorm2 = BLAS::dot(b,b) ; - CoeffType err=epsilon_*epsilon_*bnorm2 ; // Error to reach - mult(A,x,r); - BLAS::axpy(-1,b,r); - mult(C,r,d); - BLAS::copy(d,h); - rh=BLAS::dot(r,h); - - CoeffType rr=BLAS::dot(r,r); // error rr = (r|r) - while ( rr>err && its < max_iter) { - mult(A,d,Ad); - CGAL_assertion(BLAS::dot(d,Ad) != 0.0); - alpha=rh/BLAS::dot(d,Ad); - BLAS::axpy(-alpha,d,x); - BLAS::axpy(-alpha,Ad,r); - mult(C,r,h); - CGAL_assertion(rh != 0.0); - beta=1/rh; rh=BLAS::dot(r,h); beta*=rh; - BLAS::scal(beta,d); - BLAS::axpy(1,h,d); - rr=BLAS::dot(r,r); // Update error rr = (r|r) - ++its; - } - - bool success = (rr <= err); - return success; - } - -private: - CoeffType epsilon_ ; - unsigned int max_iter_ ; -} ; - - -} // namespace OpenNL - -#endif // __OPENNL_CONJUGATE_GRADIENT__ diff --git a/OpenNL/include/CGAL/OpenNL/full_vector.h b/OpenNL/include/CGAL/OpenNL/full_vector.h deleted file mode 100644 index 0459808b638..00000000000 --- a/OpenNL/include/CGAL/OpenNL/full_vector.h +++ /dev/null @@ -1,148 +0,0 @@ -// Copyright (c) 2005-2008 Inria Loria (France). -/* - * author: Bruno Levy, INRIA, project ALICE - * website: https://www.loria.fr/~levy/software - * - * This file is part of CGAL (www.cgal.org) - * - * Scientific work that use this software can reference the website and - * the following publication: - * - * @INPROCEEDINGS {levy:NMDGP:05, - * AUTHOR = Bruno Levy, - * TITLE = Numerical Methods for Digital Geometry Processing, - * BOOKTITLE =Israel Korea Bi-National Conference, - * YEAR=November 2005, - * URL=https://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics - * } - * - * Laurent Saboret 01/2005: Change for CGAL: - * - Added OpenNL namespace - * - FullVector is now a model of the SparseLinearAlgebraTraits_d::Vector concept - * - Coefficients are initialized with zeros - * - * $URL$ - * $Id$ - * SPDX-License-Identifier: LGPL-3.0-or-later - */ - - -#ifndef __OPENNL_FULL_VECTOR__ -#define __OPENNL_FULL_VECTOR__ - -#include -#include - -#include - -namespace OpenNL { - - -// Class FullVector -// Model of the SparseLinearAlgebraTraits_d::Vector concept -template class FullVector -{ -// Public types -public: - typedef T CoeffType ; - - // for SparseLinearAlgebraTraits_d::Vector concept - typedef T NT; - - // Create a vector initialized with zeros - FullVector(unsigned int dim) { - dimension_ = dim ; coeff_ = new T[dim] ; - - // init with zeros (for SparseLinearAlgebraTraits_d::Vector concept) - for (unsigned int i=0; i < dimension_; i++) - coeff_[i] = 0; - } - // Copy constructor - FullVector(const FullVector& toCopy) { - dimension_ = toCopy.dimension_ ; - coeff_ = new T[dimension_] ; - for (unsigned int i=0; i < dimension_; i++) - coeff_[i] = toCopy.coeff_[i]; - } - // operator =() - FullVector& operator =(const FullVector& toCopy) { - delete[] coeff_ ; - - dimension_ = toCopy.dimension_ ; - coeff_ = new T[dimension_] ; - for (unsigned int i=0; i < dimension_; i++) - coeff_[i] = toCopy.coeff_[i]; - - return *this; - } - - ~FullVector() { - delete[] coeff_ ; - coeff_ = nullptr ; - } - - // Return the vector's number of coefficients - unsigned int dimension() const { - return dimension_ ; - } - // Read/write access to 1 vector coefficient. - // - // Preconditions: - // * 0 <= i < dimension() - const T& operator[](unsigned int i) const { - CGAL_assertion(i < dimension_) ; - return coeff_[i] ; - } - T& operator[](unsigned int i) { - CGAL_assertion(i < dimension_) ; - return coeff_[i] ; - } - -private: - unsigned int dimension_ ; - T* coeff_ ; -} ; - -template class BLAS< FullVector > { -public: - typedef FullVector VectorType ; - typedef T CoeffType ; - - /** y <- y + a*x */ - static void axpy(CoeffType a, const VectorType& x, VectorType& y) { - CGAL_assertion(x.dimension() == y.dimension()) ; - for(unsigned int i=0; i -#include -#include -#include -#include - #include #include #include @@ -47,128 +41,6 @@ namespace OpenNL { -// Class DefaultLinearSolverTraits -// is a traits class for solving general sparse linear systems. -// It uses BICGSTAB solver with Jacobi preconditioner. -// -// Concept: Model of the SparseLinearAlgebraTraits_d concept. - -template -< - class COEFFTYPE, // type of matrix and vector coefficients - class MATRIX = SparseMatrix, // model of SparseLinearSolverTraits_d::Matrix - class VECTOR = FullVector // model of SparseLinearSolverTraits_d::Vector -> -class DefaultLinearSolverTraits -{ -// Public types -public: - typedef COEFFTYPE CoeffType ; - typedef COEFFTYPE NT; - typedef MATRIX Matrix ; - typedef VECTOR Vector ; - -// Private types -private: - typedef Jacobi_Preconditioner Preconditioner ; - typedef Solver_preconditioned_BICGSTAB - Preconditioned_solver ; - typedef Solver_BICGSTAB Solver ; - -// Public operations -public: - // Default constructor, copy constructor, operator=() and destructor are fine - - // Solve the sparse linear system "A*X = B" - // Return true on success. The solution is then (1/D) * X. - // - // Preconditions: - // - A.row_dimension() == B.dimension() - // - A.column_dimension() == X.dimension() - bool linear_solver (const Matrix& A, const Vector& B, Vector& X, NT& D) - { - D = 1; // OpenNL does not support homogeneous coordinates - - // Solve using BICGSTAB solver with preconditioner - Preconditioned_solver preconditioned_solver ; - NT omega = 1.5; - Preconditioner C(A, omega); - X = B; - if (preconditioned_solver.solve(A, C, B, X)) - return true; - - // On error, solve using BICGSTAB solver without preconditioner -#ifdef DEBUG_TRACE - std::cerr << " Failure of BICGSTAB solver with Jacobi preconditioner. " - << "Trying BICGSTAB." << std::endl; -#endif - Solver solver ; - X = B; - return solver.solve(A, B, X) ; - } -} ; - -// Class SymmetricLinearSolverTraits -// is a traits class for solving symmetric positive definite sparse linear systems. -// It uses Conjugate Gradient solver with Jacobi preconditioner. -// -// Concept: Model of the SparseLinearAlgebraTraits_d concept. - -template -< - class COEFFTYPE, // type of matrix and vector coefficients - class MATRIX = SparseMatrix, // model of SparseLinearSolverTraits_d::Matrix - class VECTOR = FullVector // model of SparseLinearSolverTraits_d::Vector -> -class SymmetricLinearSolverTraits -{ -// Public types -public: - typedef COEFFTYPE CoeffType ; - typedef COEFFTYPE NT; - typedef MATRIX Matrix ; - typedef VECTOR Vector ; - -// Private types -private: - typedef Jacobi_Preconditioner Preconditioner ; - typedef Solver_preconditioned_CG - Preconditioned_solver ; - typedef Solver_CG Solver ; - -// Public operations -public: - // Default constructor, copy constructor, operator=() and destructor are fine - - // Solve the sparse linear system "A*X = B" - // Return true on success. The solution is then (1/D) * X. - // - // Preconditions: - // - A.row_dimension() == B.dimension() - // - A.column_dimension() == X.dimension() - bool linear_solver (const Matrix& A, const Vector& B, Vector& X, NT& D) - { - D = 1; // OpenNL does not support homogeneous coordinates - - // Solve using Conjugate Gradient solver with preconditioner - Preconditioned_solver preconditioned_solver ; - NT omega = 1.5; - Preconditioner C(A, omega); - X = B; - if (preconditioned_solver.solve(A, C, B, X)) - return true; - - // On error, solve using Conjugate Gradient solver without preconditioner -#ifdef DEBUG_TRACE - std::cerr << " Failure of Conjugate Gradient solver with Jacobi preconditioner. " - << "Trying Conjugate Gradient." << std::endl; -#endif - Solver solver ; - X = B; - return solver.solve(A, B, X) ; - } -}; - /* * Solves a linear system or minimizes a quadratic form. @@ -280,11 +152,6 @@ public: bk_ = 0 ; } - void set_right_hand_side(double b) { - check_state(IN_ROW) ; - bk_ = b ; - } - void add_coefficient(unsigned int iv, double a) { check_state(IN_ROW) ; Variable& v = variable(iv) ; @@ -297,34 +164,6 @@ public: } } - void normalize_row(CoeffType weight = 1) { - check_state(IN_ROW) ; - CoeffType norm = 0.0 ; - unsigned int nf = af_.size() ; - for(unsigned int i=0; i1e-40 ); - scale_row(weight / norm) ; - } - - void scale_row(CoeffType s) { - check_state(IN_ROW) ; - unsigned int nf = af_.size() ; - for(unsigned int i=0; i -#include -#include - -#include -#include - -namespace OpenNL { - - -/** - * Base class for some preconditioners. - */ - -template -class Preconditioner { -public: - typedef T CoeffType ; - -public: - /** - * The matrix A should be square. - */ - Preconditioner( - const SparseMatrix& A, CoeffType omega = 1.0 - ) ; - - const SparseMatrix& A() const { return A_ ; } - CoeffType omega() const { return omega_ ; } - - /** - * To use this function, the matrix A should be symmetric. - */ - void mult_upper_inverse(const FullVector& x, FullVector& y) const ; - - /** - * To use this function, the matrix A should be symmetric. - */ - void mult_lower_inverse(const FullVector& x, FullVector& y) const ; - - void mult_diagonal(FullVector& xy) const ; - void mult_diagonal_inverse(FullVector& xy) const ; - -private: - const SparseMatrix& A_ ; - CoeffType omega_ ; -} ; - -template -Preconditioner::Preconditioner( - const SparseMatrix& A, CoeffType omega -) : A_(A), omega_(omega) { - //CGAL_assertion(A.is_square()) ; -} - -template -void Preconditioner::mult_lower_inverse( - const FullVector& x, FullVector& y -) const { - //CGAL_assertion(A_.has_symmetric_storage()) ; - //CGAL_assertion(A_.rows_are_stored()) ; - int n = A_.dimension() ; - for(int i=0; i::Row& Ri = A_.row(i) ; - for(int ij=0; ij < Ri.size(); ij++) { - const typename SparseMatrix::Coeff& c = Ri[ij] ; - if (c.index < i) // traverse only lower half matrix - S += c.a * y[c.index] ; - } - y[i] = (x[i] - S) * omega_ / A_.get_coef(i,i) ; - } -} -template -void Preconditioner::mult_upper_inverse( - const FullVector& x, FullVector& y -) const { - //CGAL_assertion(A_.has_symmetric_storage()) ; - //CGAL_assertion(A_.columns_are_stored()) ; - int n = A_.dimension() ; - for(int i=n-1; i>=0; i--) { - double S = 0 ; - const typename SparseMatrix::Row& Ci = A_.row(i) ; // column i == row i - for(int ij=0; ij < Ci.size(); ij++) { - const typename SparseMatrix::Coeff& c = Ci[ij] ; - if(c.index > i) // traverse only upper half matrix - S += c.a * y[c.index] ; - } - y[i] = (x[i] - S) * omega_ / A_.get_coef(i,i) ; - } -} - -template -void Preconditioner::mult_diagonal(FullVector& xy) const { - int n = A_.dimension() ; - for(int i=0; i -void Preconditioner::mult_diagonal_inverse(FullVector& xy) const { - int n = A_.dimension() ; - for(int i=0; i -class Jacobi_Preconditioner : public Preconditioner { -public: - typedef T CoeffType ; - -public: - Jacobi_Preconditioner( - const SparseMatrix& A, CoeffType omega = 1.0 - ) ; -} ; - -template -Jacobi_Preconditioner::Jacobi_Preconditioner( - const SparseMatrix& A, CoeffType omega -) : Preconditioner(A, omega) { -} - -template -void mult(const Jacobi_Preconditioner& M, const FullVector& x, FullVector& y) { - BLAS< FullVector >::copy(x, y) ; - M.mult_diagonal_inverse(y) ; -} - - -/** - * The SSOR preconditioner, sharing storage with the matrix. - */ - -template -class SSOR_Preconditioner : public Preconditioner { -public: - typedef T CoeffType ; - -public: - /** - * The matrix A should be symmetric. - */ - SSOR_Preconditioner( - const SparseMatrix& A, CoeffType omega = 1.0 - ) ; -} ; - -template -SSOR_Preconditioner::SSOR_Preconditioner( - const SparseMatrix& A, CoeffType omega -) : Preconditioner(A, omega) { -} - -/** y <- M*x */ -template -void mult(const SSOR_Preconditioner& M, const FullVector& x, FullVector& y) { - - CGAL_STATIC_THREAD_LOCAL_VARIABLE(FullVector, work,0) ; - - const SparseMatrix& A = M.A() ; - int n = A.dimension() ; - - if(work.dimension() != n) { - work = FullVector(n) ; - } - - M.mult_lower_inverse(x, work) ; - M.mult_diagonal(work) ; - M.mult_upper_inverse(work, y) ; - - BLAS< FullVector >::scal(2 - M.omega(), y) ; -} - - -} // namespace OpenNL - -#endif // __OPENNL_PRECONDITIONER__ diff --git a/OpenNL/include/CGAL/OpenNL/sparse_matrix.h b/OpenNL/include/CGAL/OpenNL/sparse_matrix.h deleted file mode 100644 index 4e4c01f46d4..00000000000 --- a/OpenNL/include/CGAL/OpenNL/sparse_matrix.h +++ /dev/null @@ -1,256 +0,0 @@ -// Copyright (c) 2005-2008 Inria Loria (France). -/* - * author: Bruno Levy, INRIA, project ALICE - * website: https://www.loria.fr/~levy/software - * - * This file is part of CGAL (www.cgal.org) - * - * Scientific work that use this software can reference the website and - * the following publication: - * - * @INPROCEEDINGS {levy:NMDGP:05, - * AUTHOR = Bruno Levy, - * TITLE = Numerical Methods for Digital Geometry Processing, - * BOOKTITLE =Israel Korea Bi-National Conference, - * YEAR=November 2005, - * URL=https://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics - * } - * - * Laurent Saboret 01/2005: Change for CGAL: - * - Added OpenNL namespace - * - SparseMatrix is now a model of the SparseLinearAlgebraTraits_d::Matrix concept - * - * $URL$ - * $Id$ - * SPDX-License-Identifier: LGPL-3.0-or-later -*/ - -#ifndef __OPENNL_SPARSE_MATRIX__ -#define __OPENNL_SPARSE_MATRIX__ - -#include -#include -#include - -#include -#include - -namespace OpenNL { - - -//________________________________________________________________ -// Class SparseMatrix -// Model of the SparseLinearAlgebraTraits_d::Matrix concept -template class SparseMatrix -{ -// Public types -public: - - typedef T CoeffType ; - - // added for SparseLinearAlgebraTraits_d::Matrix concept - typedef T NT; - - struct Coeff { - Coeff() { } - Coeff(unsigned int i, T val) : index(i), a(val) { } - unsigned int index ; - T a ; - } ; - -//__________________________________________________ - - /** - * A row or a column of a SparseMatrix. The row/column is - * compressed, and stored in the form of a list of - * (value,index) couples. - */ - class Row : public std::vector { - typedef typename std::vector superclass ; - public: - /** a_{index} <- a_{index} + val */ - void add_coef(unsigned int index, T val) - { - // search for coefficient in superclass vector - for(typename superclass::iterator it = superclass::begin() ; - it != superclass::end() ; - it++) - { - if(it->index == index) { - it->a += val ; // += - return ; - } - } - // coefficient doesn't exist yet if we reach this point - superclass::push_back(Coeff(index, val)) ; - } - - // a_{index} <- val - // (added for SparseLinearAlgebraTraits_d::Matrix concept) - // - // Optimization: - // - Caller can optimize this call by setting 'new_coef' to true - // if the coefficient does not already exists in the matrix. - void set_coef(unsigned int index, T val, bool new_coef) - { - if (!new_coef) - { - // search for coefficient in superclass vector - for(typename superclass::iterator it = superclass::begin() ; - it != superclass::end() ; - it++) - { - if(it->index == index) { - it->a = val ; // = - return ; - } - } - } - // coefficient doesn't exist yet if we reach this point - superclass::push_back(Coeff(index, val)) ; - } - - // return a_{index} (0 by default) - // (added for SparseLinearAlgebraTraits_d::Matrix concept) - T get_coef(unsigned int index) const - { - // search for coefficient in superclass vector - for(typename superclass::const_iterator it = superclass::begin() ; - it != superclass::end() ; - it++) - { - if(it->index == index) - return it->a ; // return value - } - // coefficient doesn't exist if we reach this point - return 0 ; - } - } ; - -// Public operations -public: - - //__________ constructors / destructor _____ - - // Create a square matrix initialized with zeros - SparseMatrix(unsigned int dim) { - CGAL_assertion(dim > 0); - dimension_ = dim ; - row_ = new Row[dimension_] ; - } - // Create a rectangular matrix initialized with zeros - // (added for SparseLinearAlgebraTraits_d::Matrix concept) - // WARNING: this class supports square matrices only - SparseMatrix (unsigned int rows, unsigned int columns ) { - CGAL_USE(rows); - CGAL_assertion(rows == columns); - CGAL_assertion(columns > 0); - dimension_ = columns ; - row_ = new Row[dimension_] ; - } - - ~SparseMatrix() { - delete[] row_ ; - row_ = nullptr ; - } - - //___________ access ________________________ - - // Return the matrix dimension - unsigned int dimension() const { return dimension_ ; } - - // Added for SparseLinearAlgebraTraits_d::Matrix concept: - // Return the matrix number of rows - unsigned int row_dimension() const { return dimension(); } - // Return the matrix number of columns - unsigned int column_dimension() const { return dimension(); } - - Row& row(unsigned int i) { - CGAL_assertion(i < dimension_) ; - return row_[i] ; - } - - const Row& row(unsigned int i) const { - CGAL_assertion(i < dimension_) ; - return row_[i] ; - } - - // Read access to 1 matrix coefficient - // (added for SparseLinearAlgebraTraits_d::Matrix concept) - // - // Preconditions: - // * 0 <= i < row_dimension() - // * 0 <= j < column_dimension() - NT get_coef (unsigned int i, unsigned int j) const { - CGAL_assertion(i < dimension_) ; - CGAL_assertion(j < dimension_) ; - return row(i).get_coef(j) ; - } - - // Write access to 1 matrix coefficient: a_ij <- a_ij + val - // - // Preconditions: - // * 0 <= i < row_dimension() - // * 0 <= j < column_dimension() - void add_coef(unsigned int i, unsigned int j, T val) { - CGAL_assertion(i < dimension_) ; - CGAL_assertion(j < dimension_) ; - row(i).add_coef(j, val) ; - } - - // Write access to 1 matrix coefficient: a_ij <- val - //(added for SparseLinearAlgebraTraits_d::Matrix concept) - // - // Optimization: - // - Caller can optimize this call by setting 'new_coef' to true - // if the coefficient does not already exists in the matrix. - // - // Preconditions: - // - 0 <= i < row_dimension(). - // - 0 <= j < column_dimension(). - void set_coef(unsigned int i, unsigned int j, NT val, bool new_coef = false) { - CGAL_assertion(i < dimension_) ; - CGAL_assertion(j < dimension_) ; - row(i).set_coef(j, val, new_coef) ; - } - - /** - * removes all the coefficients and frees the allocated - * space. - */ - void clear() { - for(unsigned int i=0; i -void mult(const SparseMatrix& M, const FullVector& x, FullVector& y) { - unsigned int N = M.dimension() ; - CGAL_assertion(x.dimension() == N) ; - CGAL_assertion(y.dimension() == N) ; - for(unsigned int i=0; i::Row& R = M.row(i) ; - for(unsigned int jj=0; jj > > /// \endcode -/// Otherwise, it uses CGAL's wrapping function to the OpenNL library: -/// \code -/// OpenNL::SymmetricLinearSolverTraits -/// \endcode +/// /// /// \sa `CGAL::Surface_mesh_parameterization::Two_vertices_parameterizer_3` /// @@ -103,8 +100,6 @@ public: // is always used... CGAL::Eigen_solver_traits< Eigen::SimplicialLDLT::EigenType> > - #else - OpenNL::SymmetricLinearSolverTraits #endif >::type Solver_traits; #else