From 902dd43db1073708b493a23649cde1b1ec3a299f Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 12 Jul 2016 09:09:23 +0200
Subject: [PATCH] removed unneeded source code
---
MathLib/LinAlg/MatrixTools.h | 118 ---------
MathLib/Nonlinear/NewtonRaphson-impl.h | 139 -----------
MathLib/Nonlinear/NewtonRaphson.h | 123 ----------
MathLib/Nonlinear/Picard-impl.h | 111 ---------
MathLib/Nonlinear/Picard.cpp | 30 ---
MathLib/Nonlinear/Picard.h | 107 --------
MathLib/vector_io.h | 108 --------
Tests/MathLib/TestLocalMatrixFunctions.cpp | 33 ---
Tests/MathLib/TestNonlinearNewton.cpp | 271 ---------------------
Tests/MathLib/TestNonlinearPicard.cpp | 174 -------------
10 files changed, 1214 deletions(-)
delete mode 100644 MathLib/LinAlg/MatrixTools.h
delete mode 100644 MathLib/Nonlinear/NewtonRaphson-impl.h
delete mode 100644 MathLib/Nonlinear/NewtonRaphson.h
delete mode 100644 MathLib/Nonlinear/Picard-impl.h
delete mode 100644 MathLib/Nonlinear/Picard.cpp
delete mode 100644 MathLib/Nonlinear/Picard.h
delete mode 100644 MathLib/vector_io.h
delete mode 100644 Tests/MathLib/TestLocalMatrixFunctions.cpp
delete mode 100644 Tests/MathLib/TestNonlinearNewton.cpp
delete mode 100644 Tests/MathLib/TestNonlinearPicard.cpp
diff --git a/MathLib/LinAlg/MatrixTools.h b/MathLib/LinAlg/MatrixTools.h
deleted file mode 100644
index 0c049b9c3d2..00000000000
--- a/MathLib/LinAlg/MatrixTools.h
+++ /dev/null
@@ -1,118 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef MATHLIB_MATRIXTOOLS_H_
-#define MATHLIB_MATRIXTOOLS_H_
-
-#include <type_traits>
-
-#include <Eigen/Eigen>
-
-
-namespace MathLib
-{
-
-namespace details
-{
-
-template <class T, int N, typename = void>
-struct Determinant
-{
- static_assert(N != Eigen::Dynamic, "MathLib::Determinant() for a fixed-size matrix was called for a dynamic one");
- static double eval(T const& mat)
- {
- return mat.determinant();
- }
-};
-
-template <class T, int N>
-struct Determinant<T, N, typename std::enable_if<N == Eigen::Dynamic>::type>
-{
- static_assert(N == Eigen::Dynamic, "MathLib::Determinant() for a dynamic-size matrix was called for a static one");
- static double eval(T const& mat)
- {
- if (mat.rows()==1)
- return mat(0,0);
- else if (mat.rows()==2)
- return mat(0,0) * mat(1,1) - mat(0,1) * mat(1,0);
- else if (mat.rows()==3)
- return mat(0,0) * (mat(1,1) * mat(2,2) - mat(2,1) * mat(1,2))
- + mat(2,0) * (mat(0,1) * mat(1,2) - mat(1,1) * mat(0,2))
- + mat(1,0) * (mat(0,2) * mat(2,1) - mat(2,2) * mat(0,1));
- else
- return mat.determinant();
- }
-};
-
-
-template <class T, int N, typename = void>
-struct Inverse
-{
- static_assert(N != Eigen::Dynamic, "MathLib::Inverse() for a fixed-size matrix was called for a dynamic one");
- static void eval(T const& mat, double /*detJ*/, T &result)
- {
- result.noalias() = mat.inverse();
- }
-};
-
-template <class T, int N>
-struct Inverse<T, N, typename std::enable_if<N == Eigen::Dynamic>::type>
-{
- static_assert(N == Eigen::Dynamic, "MathLib::Inverse() for a dynamic-size matrix was called for a static one");
- static void eval(T const& mat, double detJ, T &result)
- {
- if (mat.rows()==1) {
- result(0,0) = 1./detJ;
- } else if (mat.rows()==2) {
- result(0,0) = mat(1,1);
- result(0,1) = -mat(0,1);
- result(1,0) = -mat(1,0);
- result(1,1) = mat(0,0);
- result *= 1./detJ;
- } else if (mat.rows()==3) {
- result(0,0) = mat(1,1) * mat(2,2) - mat(2,1) * mat(1,2);
- result(0,1) = mat(0,2) * mat(2,1) - mat(0,1) * mat(2,2);
- result(0,2) = mat(0,1) * mat(1,2) - mat(0,2) * mat(1,1);
- result(1,0) = mat(1,2) * mat(2,0) - mat(2,2) * mat(1,0);
- result(1,1) = mat(0,0) * mat(2,2) - mat(2,0) * mat(0,2);
- result(1,2) = mat(0,2) * mat(1,0) - mat(1,2) * mat(0,0);
- result(2,0) = mat(1,0) * mat(2,1) - mat(2,0) * mat(1,1);
- result(2,1) = mat(0,1) * mat(2,0) - mat(2,1) * mat(0,0);
- result(2,2) = mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
- result *= 1./detJ;
- } else {
- result.noalias() = mat.inverse();
- }
- }
-};
-
-} // details
-
-/// returns determinant of a given Eigen dense matrix
-template <typename Derived>
-double determinant(Eigen::MatrixBase<Derived> const& mat)
-{
- using MatrixType = Eigen::MatrixBase<Derived>;
- return details::Determinant<MatrixType, MatrixType::SizeAtCompileTime>::eval(mat);
-}
-
-/// computes inverse of a given Eigen dense matrix
-/// \param mat a dense matrix whose inverse is computed
-/// \param det_mat determinant of the given matrix
-/// \param result inverted matrix whose memory should be allocated beforehand
-template <typename Derived>
-void inverse(Eigen::MatrixBase<Derived> const& mat, double det_mat, Eigen::MatrixBase<Derived> & result)
-{
- using MatrixType = Eigen::MatrixBase<Derived>;
- return details::Inverse<MatrixType, MatrixType::SizeAtCompileTime>::eval(mat, det_mat, result);
-}
-
-} // namespace
-
-#endif /* MATHLIB_MATRIXTOOLS_H_ */
diff --git a/MathLib/Nonlinear/NewtonRaphson-impl.h b/MathLib/Nonlinear/NewtonRaphson-impl.h
deleted file mode 100644
index 6fed4038533..00000000000
--- a/MathLib/Nonlinear/NewtonRaphson-impl.h
+++ /dev/null
@@ -1,139 +0,0 @@
-/**
- * \author Norihiro Watanabe
- * \date 2012-06-25
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <limits>
-
-#include <logog/include/logog.hpp>
-
-namespace MathLib
-{
-
-namespace Nonlinear
-{
-
-NewtonRaphson::NewtonRaphson()
-: _normType(VecNormType::INFINITY_N),
- _r_abs_tol(std::numeric_limits<double>::max()), _r_rel_tol(1e-6), _dx_rel_tol(.0),
- _max_itr(25), _printErrors(true), _n_iterations(0), _r_abs_error(.0), _r_rel_error(.0), _dx_rel_error(.0)
-{
-}
-
-template<class F_RESIDUAL, class F_DX, class T_VALUE>
-bool NewtonRaphson::solve(F_RESIDUAL &f_residual, F_DX &f_dx, const T_VALUE &x0, T_VALUE &x_new)
-{
- const bool checkAbsResidualError = (_r_abs_tol < std::numeric_limits<double>::max());
- const bool checkRelResidualError = (_r_rel_tol < std::numeric_limits<double>::max());
- const bool checkRelDxError = (_dx_rel_tol > .0);
- const bool needXNorm = (checkRelResidualError || checkRelDxError);
-
- INFO("------------------------------------------------------------------");
- INFO("*** NEWTON-RAPHSON nonlinear solver");
- INFO("-> iteration started");
-
- // evaluate initial residual
- x_new = x0;
- T_VALUE r(x0);
- f_residual(x_new, r);
-
- // check convergence
- double r_norm = norm(r, _normType);
-
- T_VALUE dx(x0);
- double dx_norm = norm(dx, _normType);
-
- double x_norm = -1.;
- if (needXNorm)
- x_norm = norm(x_new, _normType);
-
- bool converged = ((r_norm < _r_abs_tol && r_norm < _r_rel_tol*x_norm)
- || (checkRelDxError && dx_norm < _dx_rel_tol*x_norm));
- if (_printErrors)
- INFO("-> %d: ||r||=%1.3e, ||dx||=%1.3e, ||x||=%1.3e, ||dx||/||x||=%1.3e", 0, r_norm, dx_norm, x_norm, x_norm==0 ? dx_norm : dx_norm/x_norm);
-
- std::size_t itr_cnt = 0;
- if (!converged) {
- for (itr_cnt=1; itr_cnt<_max_itr; itr_cnt++) {
- // solve dx=-J^-1*r
- f_dx(x_new, r, dx);
- x_new += dx;
- // evaluate residual
- f_residual(x_new, r);
-#ifdef DEBUG_NEWTON_RAPHSON
- printout(std::cout, itr_cnt, x_new, r, dx);
-#endif
- // check convergence
- r_norm = norm(r, _normType);
- dx_norm = norm(dx, _normType);
- if (needXNorm)
- x_norm = norm(x_new, _normType);
- converged = ((r_norm < _r_abs_tol && r_norm < _r_rel_tol*x_norm)
- || (checkRelDxError && dx_norm < _dx_rel_tol*x_norm));
- if (_printErrors)
- INFO("-> %d: ||r||=%1.3e, ||dx||=%1.3e, ||x||=%1.3e, ||dx||/||x||=%1.3e", itr_cnt, r_norm, dx_norm, x_norm, x_norm==0 ? dx_norm : dx_norm/x_norm);
-
- if (converged)
- break;
- }
- }
-
- INFO("-> iteration finished");
- if (_max_itr==1) {
- INFO("status : iteration not required");
- } else {
- INFO("status : %s", (converged ? "CONVERGED" : "***DIVERGED***"));
- }
- INFO("iteration : %d/%d", itr_cnt, _max_itr);
- if (checkAbsResidualError)
- INFO("abs. res. : %1.3e (tolerance=%1.3e)", r_norm, _r_abs_tol);
- if (checkRelResidualError)
- INFO("rel. res. : %1.3e (tolerance=%1.3e)", x_norm==0?r_norm:r_norm/x_norm, _r_rel_tol);
- if (checkRelDxError)
- INFO("dx : %1.3e (tolerance=%1.3e)", x_norm==0?dx_norm:dx_norm/x_norm, _dx_rel_tol);
- INFO("norm type : %s", convertVecNormTypeToString(_normType).c_str());
- INFO("------------------------------------------------------------------");
-
- this->_n_iterations = itr_cnt;
- this->_r_abs_error = r_norm;
- this->_r_rel_error = (x_norm==0 ? r_norm : r_norm / x_norm);
- this->_dx_rel_error = (x_norm==0 ? dx_norm : dx_norm / x_norm);
-
- return converged;
-}
-
-
-#ifdef DEBUG_NEWTON_RAPHSON
-template<class T_VALUE>
-inline void NewtonRaphson::printout(std::ostream& os, std::size_t i, T_VALUE& x_new, T_VALUE& r, T_VALUE& dx)
-{
- os << "-> " << i <<": x=(";
- for (std::size_t i=0; i<x_new.size(); i++)
- os << x_new[i] << " ";
- os << "), r=(";
- for (std::size_t i=0; i<dx.size(); i++)
- os << r[i] << " ";
- os << "), dx=(";
- for (std::size_t i=0; i<dx.size(); i++)
- os << dx[i] << " ";
- os << ")\n";
-}
-
-// in case of double
-template<>
-inline void NewtonRaphson::printout(std::ostream& os, std::size_t i, double& x_new, double& r, double& dx)
-{
- os << "-> " << i <<": x=" << x_new << ", r=" << r << ", dx=" << dx << "\n";
-}
-#endif
-
-}
-
-} //end
diff --git a/MathLib/Nonlinear/NewtonRaphson.h b/MathLib/Nonlinear/NewtonRaphson.h
deleted file mode 100644
index 332db56c947..00000000000
--- a/MathLib/Nonlinear/NewtonRaphson.h
+++ /dev/null
@@ -1,123 +0,0 @@
-/**
- * \author Norihiro Watanabe
- * \date 2012-06-25
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef MATHLIB_NONLINEAR_NEWTONRAPHSON_H_
-#define MATHLIB_NONLINEAR_NEWTONRAPHSON_H_
-
-#include "MathLib/LinAlg/VectorNorms.h"
-
-namespace MathLib
-{
-
-namespace Nonlinear
-{
-
-/**
- * \brief Newton-Raphson method
- */
-class NewtonRaphson
-{
-public:
- /// Default constructor with norm type INFINITY_N, relative tolerance 1e-6
- /// for residual, and the maximum number of iterations 25
- NewtonRaphson();
-
- /// set a vector norm type
- void setNormType(VecNormType normType) {_normType = normType;}
-
- /// set the maximum number of iterations
- void setMaxIterations(std::size_t max_itr) {_max_itr = max_itr;}
-
- /// set the absolute residual tolerance used by the stopping criteria
- void setAbsResidualTolerance(double abs_tol) {_r_abs_tol = abs_tol;}
-
- /// set the relative residual tolerance used by the stopping criteria
- void setRelResidualTolerance(double rel_tol) {_r_rel_tol = rel_tol;}
-
- /// set the relative solution increment tolerance used by the stopping criteria
- void setRelDxTolerance(double rel_tol) {_dx_rel_tol = rel_tol;}
-
- /// print errors during iterations
- void printErrors(bool flag) {_printErrors = flag;}
-
- /**
- * solve a nonlinear problem
- *
- * \tparam F_RESIDUAL Function object returning a residual of an equation.
- * The object should have an operator ()(\f$x_k\f$, \f$r_k}\f$).
- * \f$x_k\f$ is current solution. \f$r_k\f$ is calculated residual and an output
- * of this function.
- * \tparam F_DX Function object returning a solution increment.
- * The object should have an operator ()(\f$x_k\f$, \f$r_k\f$, \f$\Delta x_{k+1}}\f$).
- * \f$\Delta x_{k+1}}\f$ is a calculated solution increment and an output of
- * this function.
- * \tparam T_VALUE Data type of \f$x_k\f$
- * Both scalar and vector types are available as far as the following conditions
- * are satisfied
- * - T_VALUE has a copy constructor (for non-primitive data types)
- * - MathLib::norm(T_VALUE) exists
- * \param f_residual Residual function object \f$r(x)\f$
- * \param f_dx Solution increment function object \f$dx(x)=J^{-1}r\f$
- * \param x0 Initial guess
- * \param x_new Solution
- * \return true if converged
- */
- template<class F_RESIDUAL, class F_DX, class T_VALUE>
- bool solve(F_RESIDUAL &f_residual, F_DX &f_dx, const T_VALUE &x0, T_VALUE &x_new);
-
- /// return the number of iterations
- std::size_t getNumberOfIterations() const {return _n_iterations; }
-
- /// return absolute error in the last iteration
- double getAbsResidualError() const {return _r_abs_error; }
-
- /// return relative error in the last iteration
- double getRelResidualError() const {return _r_rel_error; }
-
- /// return relative error in the last iteration
- double getRelDxError() const {return _dx_rel_error; }
-
-private:
-#ifdef DEBUG_NEWTON_RAPHSON
- /// print out for debugging
- template<class T_VALUE>
- inline void printout(std::ostream& os, std::size_t i, T_VALUE& x_new, T_VALUE& r, T_VALUE& dx);
-#endif
-
- /// vector norm type used to evaluate errors
- VecNormType _normType;
- /// absolute tolerance for residual
- double _r_abs_tol;
- /// relative tolerance for residual
- double _r_rel_tol;
- /// relative tolerance for dx
- double _dx_rel_tol;
- /// the maximum allowed iteration number
- std::size_t _max_itr;
- /// print iteration errors or not
- bool _printErrors;
- /// the number of iterations in the last calculation
- std::size_t _n_iterations;
- /// absolute residual error in the last calculation
- double _r_abs_error;
- /// relative residual error in the last calculation
- double _r_rel_error;
- /// relative dx error in the last calculation
- double _dx_rel_error;
-};
-
-} // Nonlinear
-} // MathLib
-
-#include "NewtonRaphson-impl.h"
-
-#endif // MATHLIB_NONLINEAR_NEWTONRAPHSON_H_
diff --git a/MathLib/Nonlinear/Picard-impl.h b/MathLib/Nonlinear/Picard-impl.h
deleted file mode 100644
index 1a42c45e941..00000000000
--- a/MathLib/Nonlinear/Picard-impl.h
+++ /dev/null
@@ -1,111 +0,0 @@
-/**
- * \author Norihiro Watanabe
- * \date 2012-06-25
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <limits>
-
-#include <logog/include/logog.hpp>
-
-namespace MathLib
-{
-
-namespace Nonlinear
-{
-
-template<class T_FUNCTOR, class T_VALUE>
-bool Picard::solve(T_FUNCTOR &functor, const T_VALUE &x0, T_VALUE &x_new)
-{
- T_VALUE x_old(x0);
- T_VALUE dx(x0);
-
- const bool checkAbsError = (_abs_tol<std::numeric_limits<double>::max());
- const bool checkRelError = (_rel_tol<std::numeric_limits<double>::max());
-
- INFO("------------------------------------------------------------------");
- INFO("*** PICARD nonlinear solver");
- INFO("-> iteration started");
- bool converged = false;
- std::size_t itr_cnt = 0;
- double x_norm = -1.;
- double abs_error = -1.;
- double rel_error = -1.;
- for (itr_cnt=0; itr_cnt<_max_itr; itr_cnt++) {
- functor(x_old, x_new);
- dx = x_new;
- dx -= x_old;
-
- abs_error = norm(dx, _normType);
- if (checkRelError) {
- x_norm = norm(x_new, _normType);
- if (x_norm>.0)
- rel_error = abs_error / x_norm;
- else
- rel_error = abs_error;
- }
- converged = (abs_error < _abs_tol && rel_error < _rel_tol);
- if (_printErrors)
- INFO("-> %d: ||dx||=%1.3e, ||x||=%1.3e, ||dx||/||x||=%1.3e", itr_cnt, abs_error, x_norm, rel_error);
-
-#ifdef DEBUG_PICARD
- printout(std::cout, itr_cnt, x_new, dx);
-#endif
- if (converged) {
- break;
- }
- x_old = x_new;
- }
-
- INFO("-> iteration finished");
- if (_max_itr==1) {
- INFO("status : iteration not required");
- } else {
- INFO("status : %s", (converged ? "CONVERGED" : "***DIVERGED***"));
- }
- INFO("iteration : %d/%d", itr_cnt, _max_itr);
- if (checkAbsError)
- INFO("abs error = %1.3e (tolerance=%1.3e)", abs_error, _abs_tol);
- if (checkRelError)
- INFO("rel error = %1.3e (tolerance=%1.3e)", rel_error, _rel_tol);
- INFO("norm type : %s", convertVecNormTypeToString(_normType).c_str());
- INFO("------------------------------------------------------------------");
-
- this->_n_iterations = itr_cnt;
- this->_abs_error = abs_error;
- this->_rel_error = rel_error;
-
- return converged;
-}
-
-
-#ifdef DEBUG_PICARD
-template<class T_VALUE>
-inline void Picard::printout(std::ostream& os, std::size_t i, T_VALUE& x_new, T_VALUE& dx)
-{
- os << "-> " << i <<": x=(";
- for (std::size_t i=0; i<x_new.size(); i++)
- os << x_new[i] << " ";
- os << "), dx=(";
- for (std::size_t i=0; i<dx.size(); i++)
- os << dx[i] << " ";
- os << ")\n";
-}
-
-// in case of double
-template<>
-inline void Picard::printout(std::ostream& os, std::size_t i, double& x_new, double& dx)
-{
- os << "-> " << i <<": x=" << x_new << ", dx=" << dx << "\n";
-}
-#endif
-
-}
-
-} //end
diff --git a/MathLib/Nonlinear/Picard.cpp b/MathLib/Nonlinear/Picard.cpp
deleted file mode 100644
index 856c35b8fa1..00000000000
--- a/MathLib/Nonlinear/Picard.cpp
+++ /dev/null
@@ -1,30 +0,0 @@
-/*
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "Picard.h"
-
-namespace MathLib
-{
-namespace Nonlinear
-{
-Picard::Picard()
- : _normType(VecNormType::INFINITY_N),
- _abs_tol(std::numeric_limits<double>::max()),
- _rel_tol(1e-6),
- _max_itr(25),
- _printErrors(false),
- _n_iterations(0),
- _abs_error(.0),
- _rel_error(.0)
-{
-}
-
-} // namespace Nonlinear
-
-} // namespace MathLib
diff --git a/MathLib/Nonlinear/Picard.h b/MathLib/Nonlinear/Picard.h
deleted file mode 100644
index 3790ed8fff1..00000000000
--- a/MathLib/Nonlinear/Picard.h
+++ /dev/null
@@ -1,107 +0,0 @@
-/**
- * \author Norihiro Watanabe
- * \date 2012-06-25
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef MATHLIB_NONLINEAR_PICARD_H_
-#define MATHLIB_NONLINEAR_PICARD_H_
-
-#include "MathLib/LinAlg/VectorNorms.h"
-
-namespace MathLib
-{
-
-namespace Nonlinear
-{
-
-/**
- * \brief Picard method (Fixed-point iteration method)
- */
-class Picard
-{
-public:
- /// Default constructor with norm type INFINITY_N, relative tolerance 1e-6,
- /// and the maximum number of iterations 25
- Picard();
-
- /// set a vector norm type
- void setNormType(VecNormType normType) {_normType = normType;}
-
- /// set the maximum number of iterations
- void setMaxIterations(std::size_t max_itr) {_max_itr = max_itr;}
-
- /// set the absolute tolerance used by the stopping criteria
- void setAbsTolerance(double abs_tol) {_abs_tol = abs_tol;}
-
- /// set the relative tolerance used by the stopping criteria
- void setRelTolerance(double rel_tol) {_rel_tol = rel_tol;}
-
- /// print errors during iterations
- void printErrors(bool flag) {_printErrors = flag;}
-
- /**
- * solve a nonlinear problem
- *
- * \tparam T_FUNCTOR Function object type which supports an
- * operator ()(\f$x_k\f$, \f$x_{k+1}\f$). \f$x_k\f$ is a previous iteration
- * step value. \f$x_{k+1}\f$ is a new solution.
- * \tparam T_VALUE Data type of \f$x_k\f$ and \f$x_{k+1}\f$.
- * Both scalar and vector types are available as far as the following conditions
- * are satisfied
- * - T_VALUE has a copy constructor (for non-primitive data types)
- * - MathLib::norm(T_VALUE) exists
- * \param functor Fixed point function object (\f$x_{k+1}=g(x_k)\f$)
- * \param x0 Initial guess
- * \param x_new Solution
- * \return true if converged
- */
- template<class T_FUNCTOR, class T_VALUE>
- bool solve(T_FUNCTOR &functor, const T_VALUE &x0, T_VALUE &x_new);
-
- /// return the number of iterations
- std::size_t getNumberOfIterations() const {return _n_iterations; }
-
- /// return absolute error in the last iteration
- double getAbsError() const {return _abs_error; }
-
- /// return relative error in the last iteration
- double getRelError() const {return _rel_error; }
-
-private:
-#ifdef DEBUG_PICARD
- /// print out for debugging
- template<class T_VALUE>
- inline void printout(std::ostream& os, std::size_t i, T_VALUE& x_new, T_VALUE& dx);
-#endif
-
- /// vector norm type used to evaluate errors
- VecNormType _normType;
- /// absolute tolerance for dx
- double _abs_tol;
- /// relative tolerance for dx
- double _rel_tol;
- /// the maximum allowed iteration number
- std::size_t _max_itr;
- /// print iteration errors or not
- bool _printErrors;
- /// the number of iterations in the last calculation
- std::size_t _n_iterations;
- /// absolute error in the last calculation
- double _abs_error;
- /// relative error in the last calculation
- double _rel_error;
-};
-
-} // Nonlinear
-} // MathLib
-
-#include "Picard-impl.h"
-
-#endif // MATHLIB_NONLINEAR_PICARD_H_
diff --git a/MathLib/vector_io.h b/MathLib/vector_io.h
deleted file mode 100644
index bff8f54d95d..00000000000
--- a/MathLib/vector_io.h
+++ /dev/null
@@ -1,108 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date no date
- * \brief Definition of vector IO functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef VECTOR_IO_H
-#define VECTOR_IO_H
-
-#include <fstream>
-#include <string>
-
-// von http://www.zfx.info/Display/Thread.php?TID=177779
-#include <sstream>
-
-#ifdef UNICODE
-typedef std::wstringstream unistringstream;
-typedef std::wstring unistring;
-#else
-typedef std::stringstream unistringstream;
-typedef std::string unistring;
-#endif
-
-template<typename A, typename T> inline const A lexical_cast(const T& source)
-{
- unistringstream s;
-
- s << source;
-
- A destination;
- s >> destination;
-
- return (destination);
-}
-
-/** reads the number of lines of non-binary stream */
-unsigned readLines ( std::istream &in )
-{
- unsigned k = 0;
- while (!in.eof()) {
- std::string str;
- getline (in, str);
- k++;
- }
-
- return k;
-}
-
-/** reads N values of non-binary stream */
-template <class T> void read ( std::istream &in, unsigned N, T *a )
-{
- unsigned k = 0;
- std::string ws (" \t");
-
- while (!in.eof () and k <= N) {
- std::string t;
- getline (in, t);
- std::size_t i1;
- if (t.length () != 0) {
- i1 = t.find_first_not_of (ws, 0);
- if (i1 != std::string::npos) {
- std::size_t i2 = t.find_first_of (ws, i1);
- if (i2 != std::string::npos) {
- a[k++] = lexical_cast<T> ( t.substr(i1, i2-i1) );
- } else {
- a[k++] = lexical_cast<T> ( t.substr(i1, t.size()-i1));
- }
- }
- }
- }
-}
-
-template <class T> int readArrays ( const std::string &fname, std::size_t &s,
- std::size_t n, T* &arr )
-{
- // open stream
- std::ifstream in (fname.c_str());
- // read number of rows
- if (in) {
- s = readLines(in)-1;
- in.close();
- } else {
- std::cout << "could not open " << fname << std::endl;
- return 1;
- }
- if (s > 0) {
- arr = new T[s * n];
- std::size_t j(0), l(0);
- // read values
- in.open (fname.c_str(), std::ios::in);
- for (j=0; j<s; ++j) {
- for (l=0; l<n; l++) in >> arr[l*s+j];
- }
- in.close ();
- }
- return 0;
-}
-
-#endif
-
diff --git a/Tests/MathLib/TestLocalMatrixFunctions.cpp b/Tests/MathLib/TestLocalMatrixFunctions.cpp
deleted file mode 100644
index e577578e5ad..00000000000
--- a/Tests/MathLib/TestLocalMatrixFunctions.cpp
+++ /dev/null
@@ -1,33 +0,0 @@
-/**
- * @copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/LICENSE.txt
- */
-
-#include <gtest/gtest.h>
-
-#include <Eigen/Eigen>
-
-#include "MathLib/LinAlg/MatrixTools.h"
-
-#include "Tests/TestTools.h"
-
-TEST(MathLib, LocalMatrixDeterminantInverse_Eigen)
-{
- Eigen::Matrix3d fMat, fInv;
- fMat << 1, 2, 3,
- 0, 1, 4,
- 5, 6, 0;
- double fMat_det = MathLib::determinant(fMat);
- MathLib::inverse(fMat, fMat_det, fInv);
-
- Eigen::MatrixXd dMat(3,3), dInv(3,3);
- dMat = fMat;
- double dMat_det = MathLib::determinant(dMat);
- MathLib::inverse(dMat, dMat_det, dInv);
-
- ASSERT_NEAR(fMat_det, dMat_det, std::numeric_limits<double>::epsilon());
- ASSERT_ARRAY_NEAR(fInv.data(), dInv.data(), fInv.size(), std::numeric_limits<double>::epsilon());
-}
diff --git a/Tests/MathLib/TestNonlinearNewton.cpp b/Tests/MathLib/TestNonlinearNewton.cpp
deleted file mode 100644
index 4a58dc722bc..00000000000
--- a/Tests/MathLib/TestNonlinearNewton.cpp
+++ /dev/null
@@ -1,271 +0,0 @@
-/**
- * \author Norihiro Watanabe
- * \date 2012-06-25
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- */
-
-#include <gtest/gtest.h>
-
-#include <valarray>
-
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
-#include "MathLib/Nonlinear/NewtonRaphson.h"
-#include "Tests/TestTools.h"
-
-namespace
-{
-
-typedef MathLib::DenseMatrix<double> MatrixType;
-typedef std::valarray<double> VectorType;
-typedef MathLib::GaussAlgorithm<MatrixType,VectorType> DenseSolverType;
-
-template<class F_JACOBIAN>
-class ScalarDx
-{
-public:
- explicit ScalarDx(F_JACOBIAN &f_J) : _f_Jacobian(f_J) {}
- // dx = - r/J
- void operator()(const double &x, const double &r, double &dx)
- {
- double j;
- _f_Jacobian(x, j);
- dx = - r/j;
- }
-
-private:
- F_JACOBIAN &_f_Jacobian;
-};
-
-template<class F_JACOBIAN>
-class VectorDx
-{
-public:
- VectorDx(F_JACOBIAN &f_J, MatrixType &matJ) : _f_Jacobian(f_J), _matJ(matJ) {}
-
- // dx = - r/J
- void operator()(const VectorType &x, const VectorType &r, VectorType &dx)
- {
- _f_Jacobian(x, _matJ);
- DenseSolverType solver;
- VectorType rhs(r);
- rhs *= -1.;
- solver.solve(_matJ, rhs, dx);
- }
-
-private:
- F_JACOBIAN &_f_Jacobian;
- MatrixType &_matJ;
-};
-
-//##############################################################################
-// Example problem 1 (one variable)
-// f(x) = x*x -4 = 0
-// x = 2,-2
-//##############################################################################
-namespace Example1
-{
-
-class Residual
-{
-public:
- void operator()(const double &x, double &r) { r = x*x-4.; }
-};
-
-class Jacobian
-{
-public:
- void operator()(const double &x, double &j) { j = 2*x; }
-};
-
-} // Example1
-
-
-//##############################################################################
-// Example problem 2 (two variables)
-// 3x-y=-2
-// 2x^2-y=0
-// (x,y) = (-1/2, 1/2) and (2, 8)
-//##############################################################################
-namespace Example2
-{
-
-class Residual
-{
-public:
- void operator()(const VectorType &x, VectorType &r)
- {
- r[0] = 3*x[0]-x[1]+2.;
- r[1] = 2*x[0]*x[0]-x[1];
- }
-};
-
-class Jacobian
-{
-public:
- void operator()(const VectorType &x, MatrixType &j)
- {
- j(0,0) = 3.;
- j(0,1) = -1.0;
- j(1,0) = 4.*x[0];
- j(1,1) = -1.0;
- }
-};
-
-} // Example2
-
-
-//##############################################################################
-// Example problem 3 (10 variables)
-//##############################################################################
-
-namespace Example3
-{
-
-class Residual
-{
-public:
- void operator()(const VectorType &x, VectorType &r)
- {
- double P = 1.;
- double R = 10.;
- double s = sqrt(2.);
- r[1-1]= (9*P*x[1-1])/4 + (9*x[2-1]*x[3-1])/(8*s) + (P*R*x[7-1])/s;
- r[2-1]= (81*P*x[2-1])/4 + (9*x[1-1]*x[3-1])/(8*s) + (P*R*x[8-1])/s;
- r[3-1]= (-9*x[1-1]*x[2-1])/(4*s) + 9*P*x[3-1] + s*P*R*x[9-1];
- r[4-1]= 36*P*x[4-1] + s*P*R*x[10-1];
- r[5-1]= -2*x[5-1] + (x[2-1]*x[7-1])/(2*s) + (x[1-1]*x[8-1])/(2*s) - (x[4-1]*x[9-1])/s + s*x[4-1]*x[9-1] - (x[3-1]*x[10-1])/s + s*x[3-1]*x[10-1];
- r[6-1]= -8*x[6-1] - (x[1-1]*x[7-1])/s - s*x[3-1]*x[9-1];
- r[7-1]= -(x[1-1]/s) - (x[2-1]*x[5-1])/(2*s) + (x[1-1]*x[6-1])/s - (3*x[7-1])/2.0 + (3*x[3-1]*x[8-1])/(4*s) + (3*x[2-1]*x[9-1])/(4*s);
- r[8-1]= -(x[2-1]/s) - (x[1-1]*x[5-1])/(2*s) - (3*x[3-1]*x[7-1])/(4*s) - (9*x[8-1])/2.0 - (3*x[1-1]*x[9-1])/(4*s);
- r[9-1]= -(s*x[3-1]) - (x[4-1]*x[5-1])/s + s*x[3-1]*x[6-1] - (3*x[2-1]*x[7-1])/(4*s) + (3*x[1-1]*x[8-1])/(4*s) - 3*x[9-1];
- r[10-1]= -(s*x[4-1]) - (x[3-1]*x[5-1])/s - 6*x[10-1];
- }
-};
-
-class Jacobian
-{
-public:
- void operator()(const VectorType &x, MatrixType &j)
- {
- double P = 1.;
- double R = 10.;
- double s = sqrt(2.);
- j = .0;
- j(1-1,1-1) = (9*P)/4.0;
- j(7-1,1-1) = -(1/s)+x[6-1]/s;
- j(1-1,2-1) = (9*x[3-1])/(8*s);
- j(7-1,2-1) = -x[5-1]/(2*s) + (3*x[9-1])/(4*s);
- j(1-1,3-1) = (9*x[2-1])/(8*s);
- j(7-1,3-1) = (3*x[8-1])/(4*s);
- j(1-1,7-1) = (P*R)/s;
- j(7-1,5-1) = -x[2-1]/(2*s);
- j(2-1,1-1) = (9*x[3-1])/(8*s);
- j(7-1,6-1) = x[1-1]/s;
- j(2-1,2-1) = (81*P)/4.0;
- j(7-1,7-1) = -1.5;
- j(2-1,3-1) = (9*x[1-1])/(8*s);
- j(7-1,8-1) = (3*x[3-1])/(4*s);
- j(2-1,8-1) = (P*R)/s;
- j(7-1,9-1) = (3*x[2-1])/(4*s);
- j(3-1,1-1) = (-9*x[2-1])/(4*s);
- j(8-1,1-1) = -x[5-1]/(2*s) - (3*x[9-1])/(4*s);
- j(3-1,2-1) = (-9*x[1-1])/(4*s);
- j(8-1,2-1) = -(1/s);
- j(3-1,3-1) = 9*P;
- j(8-1,3-1) = (-3*x[7-1])/(4*s);
- j(3-1,9-1) = s*P*R;
- j(8-1,5-1) = -x[1-1]/(2*s);
- j(4-1,4-1) = 36*P;
- j(8-1,7-1) = (-3*x[3-1])/(4*s);
- j(4-1,10-1)= s*P*R;
- j(8-1,8-1) = -4.5;
- j(5-1,1-1) = x[8-1]/(2*s);
- j(8-1,9-1) = (-3*x[1-1])/(4*s);
- j(5-1,2-1) = x[7-1]/(2*s);
- j(9-1,1-1) = (3*x[8-1])/(4*s);
- j(5-1,3-1) = -(x[10-1]/s) + s*x[10-1];
- j(9-1,2-1) = (-3*x[7-1])/(4*s);
- j(5-1,4-1) = -(x[9-1]/s) + s*x[9-1];
- j(9-1,3-1) = -s + s*x[6-1];
- j(5-1,5-1) = -2.0;
- j(9-1,4-1) = -(x[5-1]/s);
- j(5-1,7-1) = x[2-1]/(2*s);
- j(9-1,5-1) = -(x[4-1]/s);
- j(5-1,8-1) = x[1-1]/(2*s);
- j(9-1,6-1) = s*x[3-1];
- j(5-1,9-1) = -(x[4-1]/s) + s*x[4-1];
- j(9-1,7-1) = (-3*x[2-1])/(4*s);
- j(5-1,10-1)= -(x[3-1]/s) + s*x[3-1];
- j(9-1,8-1) = (3*x[1-1])/(4*s);
- j(6-1,1-1) = -(x[7-1]/s);
- j(9-1,9-1) = -3.0;
- j(6-1,3-1) = -(s*x[9-1]);
- j(10-1,3-1) = -(x[5-1]/s);
- j(6-1,6-1) = -8.0;
- j(10-1,4-1) = -s;
- j(6-1,7-1) = -(x[1-1]/s);
- j(10-1,5-1) = -(x[3-1]/s);
- j(6-1,9-1) = -(s*x[3-1]);
- j(10-1,10-1)= -6.0;
- }
-};
-
-} // Example 3
-
-} //namespace
-
-//##############################################################################
-// Tests
-//##############################################################################
-TEST(MathLib, NonlinearNR_double)
-{
- Example1::Residual f_r;
- Example1::Jacobian f_j;
- ScalarDx<Example1::Jacobian> f_dx(f_j);
- double x0 = 6.0;
- double x = .0;
- MathLib::Nonlinear::NewtonRaphson nr;
- nr.solve(f_r, f_dx, x0, x);
-
- ASSERT_NEAR(2.0, x, 1e-5);
-}
-
-TEST(MathLib, NonlinearNR_dense)
-{
- Example2::Residual f_r;
- Example2::Jacobian f_j;
- MatrixType matJ(2, 2);
- VectorDx<Example2::Jacobian> f_dx(f_j, matJ);
- VectorType x0(2), x(2);
- x0 = 6.0;
- x = .0;
- MathLib::Nonlinear::NewtonRaphson nr;
- nr.solve(f_r, f_dx, x0, x);
-
- double my_expect[] = {2., 8.};
- ASSERT_ARRAY_NEAR(my_expect, x, 2, 1e-5);
-}
-
-TEST(MathLib, NonlinearNR_dense2)
-{
- Example3::Residual f_r;
- Example3::Jacobian f_j;
- const std::size_t n = 10;
- MatrixType matJ(n, n, .0);
- VectorDx<Example3::Jacobian> f_dx(f_j, matJ);
- VectorType x0(n), x(n);
- x0 = 1.;
- x = 0.;
- MathLib::Nonlinear::NewtonRaphson nr;
- nr.solve(f_r, f_dx, x0, x);
-
- double my_expect[] = {3.39935, 3.70074e-018, -1.42576e-017, 1.4903e-021, 4.35602e-018, 0.325, -1.08167, -5.61495e-018, 7.58394e-018, -3.79368e-021};
- ASSERT_ARRAY_NEAR(my_expect, x, n, 1e-5);
-}
-
diff --git a/Tests/MathLib/TestNonlinearPicard.cpp b/Tests/MathLib/TestNonlinearPicard.cpp
deleted file mode 100644
index dcfde74741a..00000000000
--- a/Tests/MathLib/TestNonlinearPicard.cpp
+++ /dev/null
@@ -1,174 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- *
- * \file TestNonlinearPicard.cpp
- *
- * Created on 2012-08-03 by Norihiro Watanabe
- */
-
-#include <gtest/gtest.h>
-
-#include <valarray>
-
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
-#include "MathLib/Nonlinear/Picard.h"
-#include "Tests/TestTools.h"
-
-namespace
-{
-
-typedef MathLib::DenseMatrix<double> MatrixType;
-typedef std::valarray<double> VectorType;
-typedef MathLib::GaussAlgorithm<MatrixType,VectorType> DenseSolverType;
-
-
-//##############################################################################
-// Example problem 1 (one variable)
-// f(x) = x*x -4 = 0
-// x = 2,-2
-// the above function can be transformed from Newton-Rahpson equation as follows
-// x = g(x) = x - 1/f'(x) * f(x) = (x*x+4)/(2x)
-//##############################################################################
-class Example1
-{
-public:
- void operator()(const double &x, double &x_new) { x_new = (x*x+4.)/(2.*x); }
-};
-
-
-//##############################################################################
-// Example problem 2 (two variables)
-// 3x-y=-2
-// 2x^2-y=0
-// (x,y) = (-1/2, 1/2) and (2, 8)
-//##############################################################################
-class Example2
-{
-public:
- explicit Example2(std::size_t n) : A(n, n), b(n) {}
-
- void operator()(VectorType &x, VectorType &x_new)
- {
- A(0,0) = 3;
- A(0,1) = -1.0;
- A(1,0) = 2*x[0];
- A(1,1) = -1;
- b[0] = -2;
- b[1] = 0.;
-
- DenseSolverType solver;
- solver.solve(A, b, x_new);
- }
-private:
- MatrixType A;
- VectorType b;
-};
-
-} //namespace
-
-//##############################################################################
-// Tests
-//##############################################################################
-TEST(MathLib, NonlinearPicard_double)
-{
- Example1 f1;
- double x0 = 6.0;
- double x = .0;
- MathLib::Nonlinear::Picard picard;
- //picard.printErrors(true);
-
- // abs tol, reach max. iterations
- picard.setMaxIterations(3);
- picard.setAbsTolerance(1e-5);
- picard.setRelTolerance(std::numeric_limits<double>::max());
- ASSERT_FALSE(picard.solve(f1, x0, x));
- ASSERT_EQ(3u, picard.getNumberOfIterations());
-
- // abs tol, converge
- picard.setMaxIterations(10);
- ASSERT_TRUE(picard.solve(f1, x0, x));
- ASSERT_NEAR(2.0, x, 1e-5);
- ASSERT_EQ(5u, picard.getNumberOfIterations());
-
- // rel tol, converge
- picard.setMaxIterations(100);
- picard.setAbsTolerance(std::numeric_limits<double>::max());
- picard.setRelTolerance(1e-5);
- ASSERT_TRUE(picard.solve(f1, x0, x));
- ASSERT_NEAR(2.0, x, 1e-5);
- ASSERT_EQ(5u, picard.getNumberOfIterations());
-}
-
-TEST(MathLib, NonlinearPicard_vector_x0)
-{
- Example2 f2(2);
- VectorType x0(2), x(2);
-
- // initial guess1
- x0[0] = 2.;
- x0[1] = 9.;
- x = 0.0;
- MathLib::Nonlinear::Picard picard;
- ASSERT_TRUE(picard.solve(f2, x0, x));
-
- double my_expect1[] = {2., 8.};
- ASSERT_ARRAY_NEAR(my_expect1, x, 2, 1e-5);
-
- // initial guess2
- x0 = 6.0;
- x = 0.0;
- picard.solve(f2, x0, x);
-
- double my_expect2[] = {-0.5, 0.5};
- ASSERT_ARRAY_NEAR(my_expect2, x, 2, 1e-5);
-}
-
-TEST(MathLib, NonlinearPicard_vector_norms)
-{
- Example2 f2(2);
- VectorType x0(2), x(2);
- MathLib::Nonlinear::Picard picard;
- //picard.printErrors(true);
-
- x0[0] = 2.;
- x0[1] = 9.;
- x = 0.0;
- double my_expect1[] = {2., 8.};
-
- // check relative errors with different norm types
- picard.setMaxIterations(1);
- picard.setNormType(MathLib::VecNormType::NORM1);
- ASSERT_FALSE(picard.solve(f2, x0, x));
- ASSERT_NEAR(0.1, picard.getRelError(), 1e-3);
-
- picard.setNormType(MathLib::VecNormType::NORM2);
- picard.setMaxIterations(1);
- ASSERT_FALSE(picard.solve(f2, x0, x));
- ASSERT_NEAR(0.1213, picard.getRelError(), 1e-3);
-
- picard.setNormType(MathLib::VecNormType::INFINITY_N);
- picard.setMaxIterations(1);
- ASSERT_FALSE(picard.solve(f2, x0, x));
- ASSERT_NEAR(0.125, picard.getRelError(), 1e-3);
-
- // solution should be converged with any norm types
- picard.setMaxIterations(5);
- picard.setNormType(MathLib::VecNormType::NORM1);
- ASSERT_TRUE(picard.solve(f2, x0, x));
- ASSERT_ARRAY_NEAR(my_expect1, x, 2, 1e-5);
-
- picard.setNormType(MathLib::VecNormType::NORM2);
- ASSERT_TRUE(picard.solve(f2, x0, x));
- ASSERT_ARRAY_NEAR(my_expect1, x, 2, 1e-5);
-
- picard.setNormType(MathLib::VecNormType::INFINITY_N);
- ASSERT_TRUE(picard.solve(f2, x0, x));
- ASSERT_ARRAY_NEAR(my_expect1, x, 2, 1e-5);
-
-}
--
GitLab