songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Parameterization: cleanup

This commit is contained in:
Andreas Fabri 2023-10-20 14:44:17 +01:00
parent a0ca27e40a
commit 86502db00a
8 changed files with 1 additions and 1442 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

332
OpenNL/include/CGAL/OpenNL/bicgstab.h
View File

@ -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 <CGAL/OpenNL/blas.h>
#include <CGAL/assertions.h>
#include <cmath>
#include <cfloat>
#include <climits>
#include <limits>
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 MATRIX, class VECTOR> 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<Vector>::dot(b,b); // Error to reach
Vector r(n) ; // residue r=A*x-b
mult(A,x,r);
BLAS<Vector>::axpy(-1,b,r);
BLAS<Vector>::copy(r,d);
BLAS<Vector>::copy(d,h);
BLAS<Vector>::copy(h,rT);
//CGAL_assertion(BLAS<Vector>::dot(rT,rT) == 0.0); // may happen for small systems
rTh=BLAS<Vector>::dot(rT,h); // rTh = (rT|h)
rTr=BLAS<Vector>::dot(r,r); // error rTr = (r|r)
while ( rTr>err && its < max_iter) {
mult(A,d,Ad);
rTAd=BLAS<Vector>::dot(rT,Ad);
CGAL_assertion(rTAd != 0.0);
alpha=rTh/rTAd;
BLAS<Vector>::axpy(-alpha,Ad,r);
BLAS<Vector>::copy(h,s);
BLAS<Vector>::axpy(-alpha,Ad,s);
mult(A,s,t);
BLAS<Vector>::axpy(1,t,u);
BLAS<Vector>::scal(alpha,u);
st=BLAS<Vector>::dot(s,t);
tt=BLAS<Vector>::dot(t,t);
if (st == 0.0 || tt == 0.0)
omega = 0 ;
else
omega = st/tt;
BLAS<Vector>::axpy(-omega,t,r);
BLAS<Vector>::axpy(-alpha,d,x);
BLAS<Vector>::axpy(-omega,s,x);
rTr=BLAS<Vector>::dot(r,r);
BLAS<Vector>::copy(s,h);
BLAS<Vector>::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<Vector>::dot(rT,h);
beta*=rTh;
BLAS<Vector>::scal(beta,d);
BLAS<Vector>::axpy(1,h,d);
BLAS<Vector>::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<CoeffType>::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<Vector>::dot(b,b); // Error to reach
Vector r(n) ; // residue r=A*x-b
mult(A,x,r);
BLAS<Vector>::axpy(-1,b,r);
mult(C,r,d);
BLAS<Vector>::copy(d,h);
BLAS<Vector>::copy(h,rT);
//CGAL_assertion(BLAS<Vector>::dot(rT,rT) == 0.0); // may happen for small systems
rTh=BLAS<Vector>::dot(rT,h); // rTh = (rT|h)
rTr=BLAS<Vector>::dot(r,r); // error rTr = (r|r)
while (rTr>err && its < max_iter) {
mult(A,d,aux);
mult(C,aux,Sd);
rTSd=BLAS<Vector>::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<Vector>::axpy(-alpha,aux,r);
BLAS<Vector>::copy(h,s);
BLAS<Vector>::axpy(-alpha,Sd,s);
mult(A,s,aux);
mult(C,aux,t);
BLAS<Vector>::axpy(1,t,u);
BLAS<Vector>::scal(alpha,u);
st=BLAS<Vector>::dot(s,t);
tt=BLAS<Vector>::dot(t,t);
if (st == 0.0 || tt == 0.0)
omega = 0 ;
else
omega = st/tt;
BLAS<Vector>::axpy(-omega,aux,r);
BLAS<Vector>::axpy(-alpha,d,x);
BLAS<Vector>::axpy(-omega,s,x);
rTr=BLAS<Vector>::dot(r,r);
BLAS<Vector>::copy(s,h);
BLAS<Vector>::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<Vector>::dot(rT,h);
beta*=rTh;
BLAS<Vector>::scal(beta,d); // d = h + beta * (d - omega * Sd);
BLAS<Vector>::axpy(1,h,d);
BLAS<Vector>::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<CoeffType>::min)());
}
private:
CoeffType epsilon_ ;
unsigned int max_iter_ ;
} ;
} // namespace OpenNL
#endif // __OPENNL_BICGSTAB__

67
OpenNL/include/CGAL/OpenNL/blas.h
View File

@ -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 <CGAL/assertions.h>
namespace OpenNL {
/** Basic Linear Algebra Subroutines */
template <class VECTOR> 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

240
OpenNL/include/CGAL/OpenNL/conjugate_gradient.h
View File

@ -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 <CGAL/OpenNL/blas.h>
#include <CGAL/assertions.h>
#include <cmath>
#include <cfloat>
#include <climits>
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 MATRIX, class VECTOR> 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<Vector>::dot(b,b) ;
CoeffType err=epsilon_*epsilon_*bnorm2 ; // Error to reach
// residue g=b-A*x
mult(A,x,g);
BLAS<Vector>::axpy(-1,b,g);
BLAS<Vector>::scal(-1,g);
// Initially, r=g=b-A*x
BLAS<Vector>::copy(g,r); // r = g
CoeffType gg=BLAS<Vector>::dot(g,g); // error gg = (g|g)
while ( gg>err && its < max_iter) {
mult(A,r,p);
rho=BLAS<Vector>::dot(p,p);
sig=BLAS<Vector>::dot(r,p);
tau=BLAS<Vector>::dot(g,r);
CGAL_assertion(sig != 0.0);
t=tau/sig;
BLAS<Vector>::axpy(t,r,x);
BLAS<Vector>::axpy(-t,p,g);
CGAL_assertion(tau != 0.0);
gam=(t*t*rho-tau)/tau;
BLAS<Vector>::scal(gam,r);
BLAS<Vector>::axpy(1,g,r);
gg=BLAS<Vector>::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<Vector>::dot(b,b) ;
CoeffType err=epsilon_*epsilon_*bnorm2 ; // Error to reach
mult(A,x,r);
BLAS<Vector>::axpy(-1,b,r);
mult(C,r,d);
BLAS<Vector>::copy(d,h);
rh=BLAS<Vector>::dot(r,h);
CoeffType rr=BLAS<Vector>::dot(r,r); // error rr = (r|r)
while ( rr>err && its < max_iter) {
mult(A,d,Ad);
CGAL_assertion(BLAS<Vector>::dot(d,Ad) != 0.0);
alpha=rh/BLAS<Vector>::dot(d,Ad);
BLAS<Vector>::axpy(-alpha,d,x);
BLAS<Vector>::axpy(-alpha,Ad,r);
mult(C,r,h);
CGAL_assertion(rh != 0.0);
beta=1/rh; rh=BLAS<Vector>::dot(r,h); beta*=rh;
BLAS<Vector>::scal(beta,d);
BLAS<Vector>::axpy(1,h,d);
rr=BLAS<Vector>::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__

148
OpenNL/include/CGAL/OpenNL/full_vector.h
View File

@ -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 <CGAL/OpenNL/blas.h>
#include <CGAL/assertions.h>
#include <cstdlib>
namespace OpenNL {
// Class FullVector
// Model of the SparseLinearAlgebraTraits_d::Vector concept
template <class T> 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 T> class BLAS< FullVector<T> > {
public:
typedef FullVector<T> 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<x.dimension(); i++) {
y[i] += a * x[i] ;
}
}
/** x <- a*x */
static void scal(CoeffType a, VectorType& x) {
for(unsigned int i=0; i<x.dimension(); i++) {
x[i] *= a ;
}
}
/** y <- x */
static void copy(const VectorType& x, VectorType& y) {
CGAL_assertion(x.dimension() == y.dimension()) ;
for(unsigned int i=0; i<x.dimension(); i++) {
y[i] = x[i] ;
}
}
/** returns x^t * y */
static CoeffType dot(const VectorType& x, const VectorType& y) {
CoeffType result = 0 ;
CGAL_assertion(x.dimension() == y.dimension()) ;
for(unsigned int i=0; i<x.dimension(); i++) {
result += y[i] * x[i] ;
}
return result ;
}
} ;
} // namespace OpenNL
#endif // __OPENNL_FULL_VECTOR__

179
OpenNL/include/CGAL/OpenNL/linear_solver.h
View File

@ -31,12 +31,6 @@
#ifndef __OPENNL_LINEAR_SOLVER__
#define __OPENNL_LINEAR_SOLVER__
#include <CGAL/OpenNL/conjugate_gradient.h>
#include <CGAL/OpenNL/bicgstab.h>
#include <CGAL/OpenNL/preconditioner.h>
#include <CGAL/OpenNL/sparse_matrix.h>
#include <CGAL/OpenNL/full_vector.h>
#include <vector>
#include <iostream>
#include <cstdlib>
@ -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<COEFFTYPE>, // model of SparseLinearSolverTraits_d::Matrix
class VECTOR = FullVector<COEFFTYPE> // 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<NT> Preconditioner ;
typedef Solver_preconditioned_BICGSTAB<Matrix, Preconditioner, Vector>
Preconditioned_solver ;
typedef Solver_BICGSTAB<Matrix, Vector> 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<COEFFTYPE>, // model of SparseLinearSolverTraits_d::Matrix
class VECTOR = FullVector<COEFFTYPE> // 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<NT> Preconditioner ;
typedef Solver_preconditioned_CG<Matrix, Preconditioner, Vector>
Preconditioned_solver ;
typedef Solver_CG<Matrix, Vector> 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; i<nf; i++) {
norm += af_[i] * af_[i] ;
}
unsigned int nl = al_.size() ;
for(unsigned int i=0; i<nl; i++) {
norm += al_[i] * al_[i] ;
}
norm = sqrt(norm) ;
CGAL_assertion( fabs(norm)>1e-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<nf; i++) {
af_[i] *= s ;
}
unsigned int nl = al_.size() ;
for(unsigned int i=0; i<nl; i++) {
al_[i] *= s ;
}
bk_ *= s ;
}
void end_row() {
if(least_squares_) {
@ -411,24 +250,6 @@ protected:
// ----------- Finite state automaton (checks that calling sequence is respected) ---------
std::string state_to_string(State s) {
switch(s) {
case INITIAL:
return "initial" ;
case IN_SYSTEM:
return "in system" ;
case IN_ROW:
return "in row" ;
case CONSTRUCTED:
return "constructed" ;
case SOLVED:
return "solved" ;
}
// Should not go there.
CGAL_error();
return "undefined" ;
}
void check_state(State s) {
CGAL_USE(s);
CGAL_assertion(state_ == s) ;

214
OpenNL/include/CGAL/OpenNL/preconditioner.h
View File

@ -1,214 +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 2006: Changes for CGAL:
* - copied Jacobi preconditioner from Graphite 1.9 code
* - Added OpenNL namespace
*
* $URL$
* $Id$
* SPDX-License-Identifier: LGPL-3.0-or-later
*/
#ifndef __OPENNL_PRECONDITIONER__
#define __OPENNL_PRECONDITIONER__
#include <CGAL/OpenNL/blas.h>
#include <CGAL/OpenNL/sparse_matrix.h>
#include <CGAL/OpenNL/full_vector.h>
#include <CGAL/assertions.h>
#include <CGAL/tss.h>
namespace OpenNL {
/**
* Base class for some preconditioners.
*/
template <class T>
class Preconditioner {
public:
typedef T CoeffType ;
public:
/**
* The matrix A should be square.
*/
Preconditioner(
const SparseMatrix<T>& A, CoeffType omega = 1.0
) ;
const SparseMatrix<T>& 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<T>& x, FullVector<T>& y) const ;
/**
* To use this function, the matrix A should be symmetric.
*/
void mult_lower_inverse(const FullVector<T>& x, FullVector<T>& y) const ;
void mult_diagonal(FullVector<T>& xy) const ;
void mult_diagonal_inverse(FullVector<T>& xy) const ;
private:
const SparseMatrix<T>& A_ ;
CoeffType omega_ ;
} ;
template <class T>
Preconditioner<T>::Preconditioner(
const SparseMatrix<T>& A, CoeffType omega
) : A_(A), omega_(omega) {
//CGAL_assertion(A.is_square()) ;
}
template <class T>
void Preconditioner<T>::mult_lower_inverse(
const FullVector<T>& x, FullVector<T>& y
) const {
//CGAL_assertion(A_.has_symmetric_storage()) ;
//CGAL_assertion(A_.rows_are_stored()) ;
int n = A_.dimension() ;
for(int i=0; i<n; i++) {
double S = 0 ;
const typename SparseMatrix<T>::Row& Ri = A_.row(i) ;
for(int ij=0; ij < Ri.size(); ij++) {
const typename SparseMatrix<T>::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 <class T>
void Preconditioner<T>::mult_upper_inverse(
const FullVector<T>& x, FullVector<T>& 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<T>::Row& Ci = A_.row(i) ; // column i == row i
for(int ij=0; ij < Ci.size(); ij++) {
const typename SparseMatrix<T>::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 <class T>
void Preconditioner<T>::mult_diagonal(FullVector<T>& xy) const {
int n = A_.dimension() ;
for(int i=0; i<n; i++) {
xy[i] *= ( A_.get_coef(i,i) / omega_ ) ;
}
}
template <class T>
void Preconditioner<T>::mult_diagonal_inverse(FullVector<T>& xy) const {
int n = A_.dimension() ;
for(int i=0; i<n; i++) {
xy[i] *= ( omega_ / A_.get_coef(i,i) ) ;
}
}
/**
* Jacobi preconditioner
*/
template <class T>
class Jacobi_Preconditioner : public Preconditioner<T> {
public:
typedef T CoeffType ;
public:
Jacobi_Preconditioner(
const SparseMatrix<T>& A, CoeffType omega = 1.0
) ;
} ;
template <class T>
Jacobi_Preconditioner<T>::Jacobi_Preconditioner(
const SparseMatrix<T>& A, CoeffType omega
) : Preconditioner<T>(A, omega) {
}
template <class T>
void mult(const Jacobi_Preconditioner<T>& M, const FullVector<T>& x, FullVector<T>& y) {
BLAS< FullVector<T> >::copy(x, y) ;
M.mult_diagonal_inverse(y) ;
}
/**
* The SSOR preconditioner, sharing storage with the matrix.
*/
template <class T>
class SSOR_Preconditioner : public Preconditioner<T> {
public:
typedef T CoeffType ;
public:
/**
* The matrix A should be symmetric.
*/
SSOR_Preconditioner(
const SparseMatrix<T>& A, CoeffType omega = 1.0
) ;
} ;
template <class T>
SSOR_Preconditioner<T>::SSOR_Preconditioner(
const SparseMatrix<T>& A, CoeffType omega
) : Preconditioner<T>(A, omega) {
}
/** y <- M*x */
template <class T>
void mult(const SSOR_Preconditioner<T>& M, const FullVector<T>& x, FullVector<T>& y) {
CGAL_STATIC_THREAD_LOCAL_VARIABLE(FullVector<T>, work,0) ;
const SparseMatrix<T>& A = M.A() ;
int n = A.dimension() ;
if(work.dimension() != n) {
work = FullVector<T>(n) ;
}
M.mult_lower_inverse(x, work) ;
M.mult_diagonal(work) ;
M.mult_upper_inverse(work, y) ;
BLAS< FullVector<T> >::scal(2 - M.omega(), y) ;
}
} // namespace OpenNL
#endif // __OPENNL_PRECONDITIONER__

256
OpenNL/include/CGAL/OpenNL/sparse_matrix.h
View File

@ -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 <CGAL/OpenNL/full_vector.h>
#include <CGAL/assertions.h>
#include <CGAL/use.h>
#include <vector>
#include <cstdlib>
namespace OpenNL {
//________________________________________________________________
// Class SparseMatrix
// Model of the SparseLinearAlgebraTraits_d::Matrix concept
template <class T> 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<Coeff> {
typedef typename std::vector<Coeff> 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<dimension_; i++) {
row(i).clear() ;
}
}
private:
unsigned int dimension_ ;
Row* row_ ;
// SparseMatrix cannot be copied
// (for the moment, could be implemented if needed).
SparseMatrix(const SparseMatrix& rhs) ;
SparseMatrix& operator=(const SparseMatrix& rhs) ;
} ;
/** y <- M*x */
template <class T>
void mult(const SparseMatrix<T>& M, const FullVector<T>& x, FullVector<T>& y) {
unsigned int N = M.dimension() ;
CGAL_assertion(x.dimension() == N) ;
CGAL_assertion(y.dimension() == N) ;
for(unsigned int i=0; i<N; i++) {
y[i] = 0 ;
const typename SparseMatrix<T>::Row& R = M.row(i) ;
for(unsigned int jj=0; jj<R.size(); jj++) {
unsigned int j = R[jj].index ;
y[i] += R[jj].a * x[j] ;
}
}
}
} // namespace OpenNL
#endif // __OPENNL_SPARSE_MATRIX__

7
Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/LSCM_parameterizer_3.h
View File

@ -77,10 +77,7 @@ namespace Surface_mesh_parameterization {
/// \code
/// CGAL::Eigen_solver_traits<Eigen::SimplicialLDLT< Eigen::SparseMatrix<double> > >
/// \endcode
/// Otherwise, it uses CGAL's wrapping function to the OpenNL library:
/// \code
/// OpenNL::SymmetricLinearSolverTraits<typename TriangleMesh::NT>
/// \endcode
///
///
/// \sa `CGAL::Surface_mesh_parameterization::Two_vertices_parameterizer_3<TriangleMesh, BorderParameterizer, SolverTraits>`
///
@ -103,8 +100,6 @@ public:
// is always used...
CGAL::Eigen_solver_traits<
Eigen::SimplicialLDLT<Eigen_sparse_symmetric_matrix<double>::EigenType> >
#else
OpenNL::SymmetricLinearSolverTraits<typename TriangleMesh::NT>
#endif
>::type Solver_traits;
#else