Commit 1d695006 authored by Pierre Kraemer's avatar Pierre Kraemer

remove old OpenNL + update algos for new version

parent d8eb0172
......@@ -73,6 +73,7 @@ SET (EXT_INCLUDES
# define libs for external libs
SET (EXT_LIBS
PythonQt
nl
${OPENGL_LIBRARY}
${GLEW_LIBRARIES}
${ZLIB_LIBRARIES}
......
......@@ -8,14 +8,15 @@
#include "Container/fakeAttribute.h"
#include "OpenNL/linear_solver.h"
#include "NL/nl.h"
#include "Algo/LinearSolving/basic.h"
#include "Eigen/Dense"
using namespace CGoGN;
using namespace SCHNApps;
// namespace CGoGN { namespace Utils { class Drawer; } }
namespace CGoGN { namespace Utils { class Drawer; } }
enum SelectionMode
......@@ -51,7 +52,8 @@ struct PerMapParameterSet
VertexAttribute<unsigned int> vIndex;
unsigned int nb_vertices;
LinearSolver<PFP2::REAL>* solver;
// LinearSolver<PFP2::REAL>* solver;
NLContext nlContext;
};
struct ParameterSet
......@@ -143,7 +145,7 @@ private:
SurfaceDeformationDockTab* m_dockTab;
QHash<View*, ParameterSet*> h_viewParams;
// Utils::Drawer* m_drawer;
Utils::Drawer* m_drawer;
bool b_refreshingUI;
......
......@@ -6,7 +6,7 @@
#include "Algo/Geometry/normal.h"
#include "Algo/Geometry/laplacian.h"
//#include "Utils/drawer.h"
#include "Utils/drawer.h"
#include <QKeyEvent>
#include <QMouseEvent>
......@@ -48,16 +48,38 @@ PerMapParameterSet::PerMapParameterSet(MapHandlerGen* mh) :
lockingMarker = new CellMarker<VERTEX>(*map);
handleMarker = new CellMarker<VERTEX>(*map);
positionInit = mh->getAttribute<PFP2::VEC3, VERTEX>("positionInit") ;
if(!positionInit.isValid())
positionInit = mh->addAttribute<PFP2::VEC3, VERTEX>("positionInit") ;
map->copyAttribute(positionInit, positionAttribute) ;
vertexNormal = mh->getAttribute<PFP2::VEC3, VERTEX>("vertexNormal");
if(!vertexNormal.isValid())
vertexNormal = mh->addAttribute<PFP2::VEC3, VERTEX>("vertexNormal");
edgeAngle = mh->getAttribute<PFP2::REAL, EDGE>("edgeAngle");
if(!edgeAngle.isValid())
edgeAngle = mh->addAttribute<PFP2::REAL, EDGE>("edgeAngle");
edgeWeight = mh->getAttribute<PFP2::REAL, EDGE>("edgeWeight");
if(!edgeWeight.isValid())
edgeWeight = mh->addAttribute<PFP2::REAL, EDGE>("edgeWeight");
vertexArea = mh->getAttribute<PFP2::REAL, VERTEX>("vertexArea");
if(!vertexArea.isValid())
vertexArea = mh->addAttribute<PFP2::REAL, VERTEX>("vertexArea");
diffCoord = mh->getAttribute<PFP2::VEC3, VERTEX>("diffCoord");
if(!diffCoord.isValid())
diffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("diffCoord");
vertexRotationMatrix = mh->addAttribute<Eigen_Matrix3f, VERTEX>("vertexRotationMatrix") ;
rotatedDiffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("rotatedDiffCoord") ;
vertexRotationMatrix = mh->getAttribute<Eigen_Matrix3f, VERTEX>("vertexRotationMatrix") ;
if(!vertexRotationMatrix.isValid())
vertexRotationMatrix = mh->addAttribute<Eigen_Matrix3f, VERTEX>("vertexRotationMatrix");
rotatedDiffCoord = mh->getAttribute<PFP2::VEC3, VERTEX>("rotatedDiffCoord") ;
if(!rotatedDiffCoord.isValid())
rotatedDiffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("rotatedDiffCoord");
Algo::Surface::Geometry::computeNormalVertices<PFP2>(*map, positionAttribute, vertexNormal) ;
Algo::Surface::Geometry::computeAnglesBetweenNormalsOnEdges<PFP2>(*map, positionAttribute, edgeAngle) ;
......@@ -68,16 +90,23 @@ PerMapParameterSet::PerMapParameterSet(MapHandlerGen* mh) :
for(unsigned int i = vertexRotationMatrix.begin(); i != vertexRotationMatrix.end() ; vertexRotationMatrix.next(i))
vertexRotationMatrix[i] = Eigen::Matrix3f::Identity();
vIndex = mh->getAttribute<unsigned int, VERTEX>("vIndex");
if(!vIndex.isValid())
vIndex = mh->addAttribute<unsigned int, VERTEX>("vIndex");
nb_vertices = static_cast<AttribMap*>(map)->computeIndexCells<VERTEX>(vIndex);
solver = new LinearSolver<PFP2::REAL>(nb_vertices);
nlContext = nlNewContext();
nlMakeCurrent(nlContext) ;
nlSolverParameteri(NL_NB_VARIABLES, nb_vertices) ;
nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE) ;
nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT) ;
}
PerMapParameterSet::~PerMapParameterSet()
{
delete lockingMarker;
delete handleMarker;
nlDeleteContext(nlContext);
}
......@@ -95,22 +124,22 @@ bool SurfaceDeformationPlugin::enable()
connect(m_window, SIGNAL(viewAndPluginUnlinked(View*, Plugin*)), this, SLOT(viewUnlinked(View*, Plugin*)));
connect(m_window, SIGNAL(currentViewChanged(View*)), this, SLOT(currentViewChanged(View*)));
// m_drawer = new Utils::Drawer();
// registerShader(m_drawer->getShader());
m_drawer = new Utils::Drawer();
registerShader(m_drawer->getShader());
return true;
}
void SurfaceDeformationPlugin::disable()
{
// delete m_drawer;
delete m_drawer;
}
void SurfaceDeformationPlugin::redraw(View* view)
{
if(selecting)
{
/* glDisable(GL_LIGHTING) ;
glDisable(GL_LIGHTING) ;
m_drawer->newList(GL_COMPILE_AND_EXECUTE) ;
m_drawer->lineWidth(2.0f) ;
m_drawer->begin(GL_LINES) ;
......@@ -129,7 +158,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
m_drawer->vertex(selectionCenter + selectionRadius * PFP2::VEC3(0,0,-1)) ;
m_drawer->end() ;
m_drawer->endList() ;
*/ }
}
ParameterSet* params = h_viewParams[view];
MapHandlerGen* mh = params->selectedMap;
......@@ -139,7 +168,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
if(!perMap->locked_vertices.empty() || !perMap->handle_vertices.empty())
{
/* glDisable(GL_LIGHTING) ;
glDisable(GL_LIGHTING) ;
m_drawer->newList(GL_COMPILE_AND_EXECUTE) ;
m_drawer->pointSize(4.0f) ;
m_drawer->begin(GL_POINTS) ;
......@@ -154,7 +183,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
m_drawer->vertex(perMap->positionAttribute[perMap->handle_vertices[i]]) ;
m_drawer->end() ;
m_drawer->endList() ;
*/ }
}
}
}
......@@ -166,6 +195,12 @@ void SurfaceDeformationPlugin::keyPress(View* view, QKeyEvent* event)
selecting = true;
view->updateGL();
}
else if(event->key() == Qt::Key_M)
{
ParameterSet* params = h_viewParams[view];
MapHandlerGen* map = params->selectedMap;
matchDiffCoord(view, map);
}
}
void SurfaceDeformationPlugin::keyRelease(View* view, QKeyEvent* event)
......@@ -221,7 +256,8 @@ void SurfaceDeformationPlugin::mousePress(View* view, QMouseEvent* event)
}
}
LinearSolving::resetSolver(perMap->solver, false) ;
nlMakeCurrent(perMap->nlContext) ;
nlReset(NL_FALSE) ;
view->updateGL() ;
}
......@@ -277,7 +313,7 @@ void SurfaceDeformationPlugin::mouseMove(View* view, QMouseEvent* event)
dragPrevious = q ;
// matchDiffCoord(view, map);
// for(unsigned int i = 0; i < 2; ++i)
for(unsigned int i = 0; i < 2; ++i)
asRigidAsPossible(view, map);
params->selectedMap->updateVBO(perMap->positionAttribute);
......@@ -414,11 +450,12 @@ void SurfaceDeformationPlugin::changeSelectedMap(View* view, MapHandlerGen* map)
if(view->isCurrentView())
{
if(prev)
disconnect(map, SIGNAL(attributeAdded()), this, SLOT(attributeAdded()));
disconnect(prev, SIGNAL(attributeAdded()), this, SLOT(attributeAdded()));
if(map)
{
connect(map, SIGNAL(attributeAdded()), this, SLOT(attributeAdded()));
selectionRadius = map->getBBdiagSize() / 50.0;
}
m_dockTab->refreshUI(params);
view->updateGL();
......@@ -431,10 +468,7 @@ void SurfaceDeformationPlugin::changePositionAttribute(View* view, MapHandlerGen
params->perMap[map->getName()]->positionAttribute = attribute;
if(view->isCurrentView())
{
m_dockTab->refreshUI(params);
// view->updateGL();
}
}
void SurfaceDeformationPlugin::changeVerticesSelectionMode(View* view, MapHandlerGen* map, SelectionMode m)
......@@ -443,10 +477,7 @@ void SurfaceDeformationPlugin::changeVerticesSelectionMode(View* view, MapHandle
params->perMap[map->getName()]->verticesSelectionMode = m;
if(view->isCurrentView())
{
m_dockTab->refreshUI(params);
// view->updateGL();
}
}
void SurfaceDeformationPlugin::cb_selectedMapChanged()
......@@ -494,17 +525,19 @@ void SurfaceDeformationPlugin::matchDiffCoord(View* view, MapHandlerGen* mh)
PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
PerMapParameterSet* perMap = h_viewParams[view]->perMap[mh->getName()];
LinearSolving::initSolver(perMap->solver, perMap->nb_vertices, true, true) ;
nlMakeCurrent(perMap->nlContext);
if(nlGetCurrentState() == NL_STATE_INITIAL)
nlBegin(NL_SYSTEM) ;
for(int coord = 0; coord < 3; ++coord)
{
LinearSolving::setupVariables<PFP2>(*map, perMap->solver, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
LinearSolving::startMatrix(perMap->solver) ;
LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
// LinearSolving::addRowsRHS_Laplacian_Topo<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->diffCoord, coord);
LinearSolving::endMatrix(perMap->solver) ;
LinearSolving::solve(perMap->solver) ;
LinearSolving::getResult<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->positionAttribute, coord) ;
LinearSolving::resetSolver(perMap->solver, true) ;
LinearSolving::setupVariables<PFP2>(*map, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
nlBegin(NL_MATRIX) ;
LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
nlEnd(NL_MATRIX) ;
nlEnd(NL_SYSTEM) ;
nlSolve() ;
LinearSolving::getResult<PFP2>(*map, perMap->vIndex, perMap->positionAttribute, coord) ;
nlReset(NL_TRUE) ;
}
}
......@@ -620,16 +653,19 @@ void SurfaceDeformationPlugin::asRigidAsPossible(View* view, MapHandlerGen* mh)
}
}
LinearSolving::initSolver(perMap->solver, perMap->nb_vertices, true, true) ;
nlMakeCurrent(perMap->nlContext);
if(nlGetCurrentState() == NL_STATE_INITIAL)
nlBegin(NL_SYSTEM) ;
for(int coord = 0; coord < 3; ++coord)
{
LinearSolving::setupVariables<PFP2>(*map, perMap->solver, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
LinearSolving::startMatrix(perMap->solver) ;
LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->rotatedDiffCoord, coord) ;
LinearSolving::endMatrix(perMap->solver) ;
LinearSolving::solve(perMap->solver) ;
LinearSolving::getResult<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->positionAttribute, coord) ;
LinearSolving::resetSolver(perMap->solver, true) ;
LinearSolving::setupVariables<PFP2>(*map, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
nlBegin(NL_MATRIX) ;
LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->rotatedDiffCoord, coord) ;
nlEnd(NL_MATRIX) ;
nlEnd(NL_SYSTEM) ;
nlSolve() ;
LinearSolving::getResult<PFP2>(*map, perMap->vIndex, perMap->positionAttribute, coord) ;
nlReset(NL_TRUE) ;
}
}
......
importPlugin = schnapps.loadPlugin("ImportSurface");
renderPlugin = schnapps.loadPlugin("Render");
renderVectorPlugin = schnapps.loadPlugin("RenderVector");
differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
#renderVectorPlugin = schnapps.loadPlugin("RenderVector");
#differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
#subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
surfaceDeformationPlugin = schnapps.loadPlugin("SurfaceDeformation");
#obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_good_9k.off");
obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_good_9k.off");
#v = schnapps.getView("view_0");
v = schnapps.getView("view_0");
#schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
#schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
#schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
#schnapps.linkViewAndMap(v.getName(), obj.getName());
schnapps.linkViewAndMap(v.getName(), obj.getName());
#differentialPropertiesPlugin.computeNormal(obj.getName());
#differentialPropertiesPlugin.computeCurvature(obj.getName());
......@@ -15,7 +15,7 @@ ELSE (${CMAKE_CURRENT_BINARY_DIR} MATCHES "(.*)Debug")
ENDIF (${CMAKE_CURRENT_BINARY_DIR} MATCHES "(.*)Debug")
IF(WIN32)
SET(LIBRARY_OUTPUT_PATH ${CGoGN_ROOT_DIR}/lib)#release added by visual
SET(LIBRARY_OUTPUT_PATH ${CGoGN_ROOT_DIR}/lib) #release added by visual
INCLUDE_DIRECTORIES(${CGoGN_ROOT_DIR}/windows_dependencies/include/)
ELSE(WIN32)
SET(LIBRARY_OUTPUT_PATH ${CGoGN_ROOT_DIR}/lib/${CMAKE_BUILD_TYPE})
......@@ -40,6 +40,9 @@ add_subdirectory(Tools Tools/build)
add_subdirectory(PythonQt PythonQt/build)
INSTALL (DIRECTORY PythonQt/src/ DESTINATION ${CGoGN_ROOT_DIR}/ThirdParty/include/PythonQt FILES_MATCHING PATTERN "*.h")
add_subdirectory(OpenNL OpenNL/build)
INSTALL (DIRECTORY OpenNL/src/NL/ DESTINATION ${CGoGN_ROOT_DIR}/ThirdParty/include/NL FILES_MATCHING PATTERN "*.h")
IF (WITH_ZINRI)
add_subdirectory(Zinri Zinri/build)
INSTALL (DIRECTORY Zinri/ DESTINATION ${CGoGN_ROOT_DIR}/ThirdParty/include/Zinri FILES_MATCHING PATTERN "*.h")
......
This diff is collapsed.
--------------------- Open Numerical Library <T> -------------------------
General information
===================
This is OpenNL<T>, a library to easily construct and solve sparse linear systems.
* OpenNL<T> is supplied with a set of built-in iterative solvers
(Conjugate gradient and BICGSTAB)
* OpenNL can also use other solvers (by writing a SolverTraits class)
* ../lscm_with_generic_api.cpp implements the LSCM parameterization method with OpenNL<T>
Sorry, there is no more doc than that ...
/*
* OpenNL<T>: Generic Numerical Library
* Copyright (C) 2004 Bruno Levy
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* In addition, as a special exception, the INRIA gives permission to link the
* code of this program with the CGAL library, and distribute linked combinations
* including the two. You must obey the GNU General Public License in all respects
* for all of the code used other than CGAL.
*
* If you modify this software, you should include a notice giving the
* name of the person performing the modification, the date of modification,
* and the reason for such modification.
*
* Contact: Bruno Levy
*
* levy@loria.fr
*
* ISA-ALICE Project
* LORIA, INRIA Lorraine,
* Campus Scientifique, BP 239
* 54506 VANDOEUVRE LES NANCY CEDEX
* FRANCE
*
* Note that the GNU General Public License does not permit incorporating
* the Software into proprietary programs.
*/
#ifndef __BICGSTAB__
#define __BICGSTAB__
#include <cmath>
#include <cassert>
#include "blas.h"
/**
* 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 matrix, a function
* mult(const MATRIX& M, 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 VECTOR>
class Solver_BICGSTAB {
public:
typedef MATRIX Matrix ;
typedef VECTOR Vector ;
typedef typename Vector::CoeffType CoeffType ;
Solver_BICGSTAB() {
epsilon_ = 1e-5 ;
max_iter_ = 0 ;
}
void set_epsilon(CoeffType eps) { epsilon_ = eps ; }
void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; }
void solve(const MATRIX &A, const VECTOR& b, VECTOR& x) {
assert(A.n() == A.m()) ;
assert(A.has_symmetric_storage()) ;
assert(A.n() == b.dimension()) ;
assert(A.n() == x.dimension()) ;
unsigned int n = A.n() ;
unsigned int max_iter = max_iter_ ;
if(max_iter == 0) {
max_iter = 5 * n ;
}
Vector rT(n) ;
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 ;
CoeffType err = epsilon_ * epsilon_ * BLAS<Vector>::dot(b,b) ;
Vector r(n) ;
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) ;
assert( BLAS<Vector>::dot(rT,rT)>1e-40 ) ;
rTh=BLAS<Vector>::dot(rT,h) ;
rTr=BLAS<Vector>::dot(r,r) ;
while ( rTr > err && its < max_iter) {
mult(A,d,Ad) ;
rTAd = BLAS<Vector>::dot(rT,Ad) ;
assert( fabs(rTAd)>1e-40 ) ;
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 ( fabs(st) < 1e-40 || fabs(tt) < 1e-40 )
omega = 0 ;
else
omega = st/tt ;
BLAS<Vector>::axpy(-omega,t,r) ;
BLAS<Vector>::axpy(-alpha,d,x) ;
BLAS<Vector>::axpy(-omega,s,x) ;
BLAS<Vector>::copy(s,h) ;
BLAS<Vector>::axpy(-omega,t,h) ;
beta = (alpha/omega) / rTh ;
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) ;
rTr = BLAS<Vector>::dot(r,r) ;
its++ ;
}
}
private:
CoeffType epsilon_ ;
unsigned int max_iter_ ;
} ;
#endif
/*
* OpenNL<T>: Generic Numerical Library
* Copyright (C) 2004 Bruno Levy
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* In addition, as a special exception, the INRIA gives permission to link the
* code of this program with the CGAL library, and distribute linked combinations
* including the two. You must obey the GNU General Public License in all respects
* for all of the code used other than CGAL.
*
* If you modify this software, you should include a notice giving the
* name of the person performing the modification, the date of modification,
* and the reason for such modification.
*
* Contact: Bruno Levy
*
* levy@loria.fr
*
* ISA-ALICE Project
* LORIA, INRIA Lorraine,
* Campus Scientifique, BP 239
* 54506 VANDOEUVRE LES NANCY CEDEX
* FRANCE
*
* Note that the GNU General Public License does not permit incorporating
* the Software into proprietary programs.
*/
#ifndef __BLAS__
#define __BLAS__
#include <cassert>
/*
* Basic Linear Algebra Subroutines
*
* VectorType should have the following members :
* CoeffType -> type of the coefficients used in the vector
* operator[](unsigned int i)
* dimension()
*
*/
template <class VectorType>
class BLAS {
public:
typedef typename VectorType::CoeffType CoeffType ;
/** y <- y + a*x */
static void axpy(CoeffType a, const VectorType& x, VectorType& y) {
assert(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) {
assert(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 ;
assert(x.dimension() == y.dimension()) ;
for(unsigned int i=0; i<x.dimension(); i++) {
result += y[i] * x[i] ;
}
return result ;
}
} ;
#endif
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps *
* version 0.1 *
* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg *
* *
* This library is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by the *
* Free Software Foundation; either version 2.1 of the License, or (at your *
* option) any later version. *
* *
* This library is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this library; if not, write to the Free Software Foundation, *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. *
* *
* Web site: http://cgogn.unistra.fr/ *
* Contact information: cgogn@unistra.fr *
* *
*******************************************************************************/
#ifndef __GPU_CONJUGATE_GRADIENT__
#define __GPU_CONJUGATE_GRADIENT__
#include "CNC/cnc_gpu_solver.h"
template< class MATRIX, class VECTOR >
class Solver_CG_GPU
{
public:
typedef MATRIX Matrix ;
typedef VECTOR Vector ;
typedef typename Vector::CoeffType CoeffType ;
Solver_CG_GPU() {
epsilon_ = 1e-5 ;
max_iter_ = 2000 ;
}
void set_epsilon(CoeffType eps) { epsilon_ = eps ; }
void set_max_iter(unsigned int max_iter) { max_iter_ = max_iter ; }