Commit 73308d32 authored by Nabil Chouika's avatar Nabil Chouika
Browse files

refs#16

In ElementaryUtils (trunk):
- Change GenericType and Formatter precision to 16 instead of 17, to avoid messy numbers in formatting.

In NumA++ (trunk):
- Change some std::runtime_error to CustonException. A lot more to change...
- Time monitoring in LSMR.

In PARTONS:
- Inherit IncompleteGPDModule from GPDModule to be able to use GPDService. Lot of useless code commented.
- Some changes in Radon inversion modules...
parent c5acdce0
#include "../../../../include/NumA/integration/one_dimension/ChebyshevAIntegrator1D.h" #include "../../../../include/NumA/integration/one_dimension/ChebyshevAIntegrator1D.h"
#include <stdexcept>
#include <vector> #include <vector>
#include <ElementaryUtils/logger/CustomException.h>
#include <math.h> #include <math.h>
namespace NumA { namespace NumA {
ChebyshevAIntegrator1D::ChebyshevAIntegrator1D(unsigned int N) : ChebyshevAIntegrator1D::ChebyshevAIntegrator1D(unsigned int N) :
...@@ -41,16 +40,16 @@ void ChebyshevAIntegrator1D::setNodesAndWeights() { ...@@ -41,16 +40,16 @@ void ChebyshevAIntegrator1D::setNodesAndWeights() {
// Check number of steps // Check number of steps
if (m_N < 1) { if (m_N < 1) {
throw std::runtime_error( throw ElemUtils::CustomException("ChebyshevAIntegrator1D", __func__,
"[ChebyshevAIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 1"); "Invalid number of nodes. Must be > 0");
} else { } else {
m_nodes.assign(m_N, 0.); m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.); m_weights.assign(m_N, 0.);
for (unsigned int i = 0; i <m_N; i++) { for (unsigned int i = 0; i < m_N; i++) {
phi = M_PI*double(m_N-i)/double(m_N+1); phi = M_PI * double(m_N - i) / double(m_N + 1);
m_nodes.at(i) = cos(phi); m_nodes.at(i) = cos(phi);
m_weights.at(i) = M_PI/double(m_N+1)*(sin(phi)*sin(phi)); m_weights.at(i) = M_PI / double(m_N + 1) * (sin(phi) * sin(phi));
} }
} }
} }
......
#include "../../../../include/NumA/integration/one_dimension/ChebyshevBIntegrator1D.h" #include "../../../../include/NumA/integration/one_dimension/ChebyshevBIntegrator1D.h"
#include <stdexcept>
#include <vector> #include <vector>
#include <ElementaryUtils/logger/CustomException.h>
#include <math.h> #include <math.h>
namespace NumA { namespace NumA {
ChebyshevBIntegrator1D::ChebyshevBIntegrator1D(unsigned int N) : ChebyshevBIntegrator1D::ChebyshevBIntegrator1D(unsigned int N) :
...@@ -41,16 +40,16 @@ void ChebyshevBIntegrator1D::setNodesAndWeights() { ...@@ -41,16 +40,16 @@ void ChebyshevBIntegrator1D::setNodesAndWeights() {
// Check number of steps // Check number of steps
if (m_N < 1) { if (m_N < 1) {
throw std::runtime_error( throw ElemUtils::CustomException("ChebyshevBIntegrator1D", __func__,
"[ChebyshevBIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 1"); "Invalid number of nodes. Must be > 0");
} else { } else {
m_nodes.assign(m_N, 0.); m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.); m_weights.assign(m_N, 0.);
for (unsigned int i = 0; i < m_N; i++) { for (unsigned int i = 0; i < m_N; i++) {
phi = M_PI*double(2*m_N-1-2*i)/double(2*m_N); phi = M_PI * double(2 * m_N - 1 - 2 * i) / double(2 * m_N);
m_nodes.at(i) = cos(phi); m_nodes.at(i) = cos(phi);
m_weights.at(i) = M_PI/double(m_N); m_weights.at(i) = M_PI / double(m_N);
} }
} }
} }
......
#include "../../../../include/NumA/integration/one_dimension/GaussLegendreIntegrator1D.h" #include "../../../../include/NumA/integration/one_dimension/GaussLegendreIntegrator1D.h"
#include <stdexcept>
#include <vector> #include <vector>
#include <ElementaryUtils/logger/CustomException.h>
#include <math.h> #include <math.h>
//#include "../../../../include/NumA/utils/GenericType.h" //#include "../../../../include/NumA/utils/GenericType.h"
...@@ -41,8 +41,8 @@ void GaussLegendreIntegrator1D::setNodesAndWeights() { ...@@ -41,8 +41,8 @@ void GaussLegendreIntegrator1D::setNodesAndWeights() {
// Check number of steps // Check number of steps
if (m_N <= 0) { if (m_N <= 0) {
throw std::runtime_error( throw ElemUtils::CustomException("GaussLegendreIntegrator1D", __func__,
"[GaussLegendreIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 0"); "Invalid number of nodes. Must be > 0");
} else { } else {
m_nodes.assign(m_N, 0.); m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.); m_weights.assign(m_N, 0.);
......
...@@ -2,13 +2,10 @@ ...@@ -2,13 +2,10 @@
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
//#include <functional>
//#include <iostream>
#include <limits> #include <limits>
//#include <numeric> //#include <stdexcept>
//#include <sstream>
#include <stdexcept>
#include <vector> #include <vector>
#include <ElementaryUtils/logger/CustomException.h>
//#include "../../../../include/NumA/utils/GenericType.h" //#include "../../../../include/NumA/utils/GenericType.h"
//#include "../../../../include/NumA/utils/Parameters.h" //#include "../../../../include/NumA/utils/Parameters.h"
...@@ -47,46 +44,50 @@ void GaussLegendreSeStIntegrator1D::setNodesAndWeights() { ...@@ -47,46 +44,50 @@ void GaussLegendreSeStIntegrator1D::setNodesAndWeights() {
// Check number of steps // Check number of steps
if (m_N <= 0) { if (m_N <= 0) {
throw std::runtime_error( throw ElemUtils::CustomException("GaussLegendreSeStIntegrator1D",
"[GaussLegendreSeStIntegrator1D::setNodesAndWeights()] Invalid number of nodes. Must be > 0"); __func__, "Invalid number of nodes. Must be > 0");
} else { } else {
m_nodes.assign(m_N, 0.); m_nodes.assign(m_N, 0.);
m_weights.assign(m_N, 0.); m_weights.assign(m_N, 0.);
std::vector< double > c(m_N,0.); std::vector<double> c(m_N, 0.);
double cc = 1.; double cc = 1.;
double r = 0.; double r = 0.;
double xtemp = 0.; double xtemp = 0.;
unsigned int const maxit = 10; unsigned int const maxit = 10;
double p0, dp0, p1, dp1, p2, dp2, d; double p0, dp0, p1, dp1, p2, dp2, d;
for (unsigned int i = 0; i <m_N; i++) { for (unsigned int i = 0; i < m_N; i++) {
c[i] = double(i*i)/double((2*i+1)*(2*i-1)); c[i] = double(i * i) / double((2 * i + 1) * (2 * i - 1));
} }
for (unsigned int i = 1; i <m_N; i++) { for (unsigned int i = 1; i < m_N; i++) {
cc *= c[i]; cc *= c[i];
} }
cc = 2*cc; cc = 2 * cc;
for (unsigned int i = 0; i <m_N; i++) { for (unsigned int i = 0; i < m_N; i++) {
if (i == 0) { if (i == 0) {
r = 2.78/(4.0+double(m_N*m_N)); r = 2.78 / (4.0 + double(m_N * m_N));
xtemp = 1.0-r; xtemp = 1.0 - r;
} else if (i == 1) { } else if (i == 1) {
r = 1.0+0.06*((double) m_N-8)/((double) m_N); r = 1.0 + 0.06 * ((double) m_N - 8) / ((double) m_N);
xtemp = xtemp-4.1*r*(1.0-xtemp); xtemp = xtemp - 4.1 * r * (1.0 - xtemp);
} else if (i == 2) { } else if (i == 2) {
r = 1.0+0.22*((double) m_N-8)/((double) m_N); r = 1.0 + 0.22 * ((double) m_N - 8) / ((double) m_N);
xtemp = xtemp-1.67*r*(m_nodes[0]-xtemp); xtemp = xtemp - 1.67 * r * (m_nodes[0] - xtemp);
} else if (i < m_N-2) { } else if (i < m_N - 2) {
xtemp = 3.0*m_nodes[i-1]-3.0*m_nodes[i-2]+m_nodes[i-3]; xtemp = 3.0 * m_nodes[i - 1] - 3.0 * m_nodes[i - 2]
} else if (i == m_N-2) { + m_nodes[i - 3];
r = 1.0/(1.0+0.639*((double) m_N-4)/(1.0+0.71*((double) m_N-4))); } else if (i == m_N - 2) {
xtemp = xtemp+r*(xtemp-m_nodes[i-2])/0.766; r = 1.0
} else if (i == m_N-1) { / (1.0
r = 1.0/(1.0+0.22*((double) m_N-8)/((double) m_N)); + 0.639 * ((double) m_N - 4)
xtemp = xtemp+r*(xtemp-m_nodes[i-2])/1.67; / (1.0 + 0.71 * ((double) m_N - 4)));
xtemp = xtemp + r * (xtemp - m_nodes[i - 2]) / 0.766;
} else if (i == m_N - 1) {
r = 1.0 / (1.0 + 0.22 * ((double) m_N - 8) / ((double) m_N));
xtemp = xtemp + r * (xtemp - m_nodes[i - 2]) / 1.67;
} }
for (unsigned int j = 1; j <= maxit; j++) { for (unsigned int j = 1; j <= maxit; j++) {
...@@ -99,19 +100,22 @@ void GaussLegendreSeStIntegrator1D::setNodesAndWeights() { ...@@ -99,19 +100,22 @@ void GaussLegendreSeStIntegrator1D::setNodesAndWeights() {
dp0 = dp1; dp0 = dp1;
p1 = p2; p1 = p2;
dp1 = dp2; dp1 = dp2;
p2 = xtemp*p1-c[k]*p0; p2 = xtemp * p1 - c[k] * p0;
dp2 = xtemp*dp1+p1-c[k]*dp0; dp2 = xtemp * dp1 + p1 - c[k] * dp0;
} }
d = p2/dp2; d = p2 / dp2;
xtemp = xtemp-d; xtemp = xtemp - d;
if (std::abs(d) < std::numeric_limits<double>::epsilon()*(std::abs(d)+1.)) break; if (std::abs(d)
< std::numeric_limits<double>::epsilon()
* (std::abs(d) + 1.))
break;
} }
m_nodes[i] = xtemp; m_nodes[i] = xtemp;
m_weights[i] = cc/dp2/p1; m_weights[i] = cc / dp2 / p1;
} }
std::reverse(m_nodes.begin(), m_nodes.end()); std::reverse(m_nodes.begin(), m_nodes.end());
std::reverse(m_weights.begin(), m_weights.end()); std::reverse(m_weights.begin(), m_weights.end());
......
...@@ -23,6 +23,7 @@ ...@@ -23,6 +23,7 @@
#include <iterator> #include <iterator>
#include <numeric> #include <numeric>
#include <vector> #include <vector>
#include <ctime>
namespace NumA { namespace NumA {
...@@ -353,10 +354,17 @@ void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b, ...@@ -353,10 +354,17 @@ void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b,
<< std::endl; << std::endl;
} }
// Test time of an iteration
std::clock_t start;
double duration;
// Main iteration loop // Main iteration loop
do { do {
this->itn++; this->itn++;
// Test time of an iteration
start = std::clock();
//---------------------------------------------------------------- //----------------------------------------------------------------
// Perform the next step of the bidiagonalization to obtain the // Perform the next step of the bidiagonalization to obtain the
// next beta, u, alpha, v. These satisfy // next beta, u, alpha, v. These satisfy
...@@ -522,6 +530,9 @@ void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b, ...@@ -522,6 +530,9 @@ void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b,
if (test1 <= rtol) if (test1 <= rtol)
this->istop = 1; this->istop = 1;
// Test time of an iteration
duration = (std::clock() - start) / (double) CLOCKS_PER_SEC;
//---------------------------------------------------------------- //----------------------------------------------------------------
// See if it is time to print something. // See if it is time to print something.
//---------------------------------------------------------------- //----------------------------------------------------------------
...@@ -561,7 +572,7 @@ void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b, ...@@ -561,7 +572,7 @@ void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b,
(*this->nout) << this->itn << ", " << x[0] << ", " (*this->nout) << this->itn << ", " << x[0] << ", "
<< this->normr << ", " << this->normAr << ", " << test1 << this->normr << ", " << this->normAr << ", " << test1
<< ", " << test2 << ", " << this->normA << ", " << ", " << test2 << ", " << this->normA << ", "
<< this->condA << std::endl; << this->condA << ", " << duration << "s" << std::endl;
} }
} }
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment