Commit c5acdce0 authored by Nabil Chouika's avatar Nabil Chouika
Browse files

refs#16

In NumA++:
- Added -pedantic flag in CMakeLists for warnings.
- Moved LinAlgUtils to parent folder.
- Added LSMR iterative algorithm for sparse least-squares problems, taken from https://github.com/tvercaut/LSQR-cpp (APACHE licence) and adapted to fit into NumA.

In NumA++_Test:
- Added -pedantic flag in CMakeLists for warnings.
- Added test for LSMR also taken from github.

In PARTONS:
- Added -pedantic flag in CMakeLists for warnings.
- Converted one observable scenario to the new cross configuration.
- Various modifications and additions to the Radon Inverse files. Now, the DoubleDistribution modules are not needed anymore. It was too complicated to deal with... But it's still possible to use them I guess.
- Renamed the Radon Inverse daughter class because the brute-force approach was dropped in favor of a fixed number of rows for the matrix.

In PARTONS_exe:
- Added -pedantic flag in CMakeLists for warnings.
- Updated partons.properties to the new scheme.
- Added a test for the Radon inversion.

In PARTONS_release:
- Added -pedantic flag in CMakeLists for warnings.
- Added a scenario for GPD.
- Close PARTONS properly when exception.
parent 84888091
......@@ -5,7 +5,7 @@ cmake_minimum_required(VERSION 2.6)
project(NumA++ CXX)
# Force C++98
add_definitions(-std=c++98)
add_definitions(-std=c++98 -pedantic)
# placez vos exécutables dans des dossiers portant le nom du type de compilation (Release, Debug...). Ce dernier est défini par la variable CMAKE_BUILD_TYPE (qui n'a d'effet seulement si le générateur choisi est basé sur Make)
......
......@@ -8,8 +8,8 @@
#ifndef LINALGUTILS_H_
#define LINALGUTILS_H_
#include "../matrix/MatrixD.h"
#include "../vector/VectorD.h"
#include "matrix/MatrixD.h"
#include "vector/VectorD.h"
namespace NumA {
class MatrixD;
......@@ -35,7 +35,8 @@ public:
* @param method Solver method.
* @return X Solution: VectorD
*/
static VectorD solve(const MatrixD & A, const VectorD & Y, Method method = ColQR);
static VectorD solve(const MatrixD & A, const VectorD & Y, Method method =
ColQR);
/**
* @brief Computes the rank of the matrix with the given method.
......
This diff is collapsed.
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef LSQR_lsmrDense_h
#define LSQR_lsmrDense_h
#include "../matrix/MatrixD.h"
#include "lsmrBase.h"
namespace NumA {
/** \class lsmrDense
*
* Specific implementation of the solver for a type of dense Matrix.
*
*/
class lsmrDense: public lsmrBase {
public:
lsmrDense();
lsmrDense(const MatrixD& A);
virtual ~lsmrDense();
/**
* computes y = y + A*x without altering x,
* where A is a matrix of dimensions A[m][n].
* The size of the vector x is n.
* The size of the vector y is m.
*/
void Aprod1(unsigned int m, unsigned int n, const double * x,
double * y) const;
/**
* computes x = x + A'*y without altering y,
* where A is a matrix of dimensions A[m][n].
* The size of the vector x is n.
* The size of the vector y is m.
*/
void Aprod2(unsigned int m, unsigned int n, double * x,
const double * y) const;
/** Set the matrix A of the equation to be solved A*x = b. */
void setMatrix(const MatrixD& A);
private:
MatrixD m_A;
};
} /* namespace NumA */
#endif
......@@ -10,7 +10,7 @@
#include <stddef.h>
#include "../linear_algebra/eigen/LinAlgUtils.h"
#include "../linear_algebra/LinAlgUtils.h"
#include "../linear_algebra/matrix/MatrixD.h"
#include "../linear_algebra/vector/VectorD.h"
......
#include "../../../../include/NumA/linear_algebra/eigen/LinAlgUtils.h"
#include "../../../include/NumA/linear_algebra/LinAlgUtils.h"
#include <Eigen/Core>
#include <Eigen/LU>
......@@ -6,13 +6,12 @@
#include <sstream>
#include <stdexcept>
#include "../../../../include/NumA/linear_algebra/eigen/EigenUtils.h"
//#include "../../../../include/NumA/linear_algebra/matrix/MatrixD.h"
#include "../../../include/NumA/linear_algebra/eigen/EigenUtils.h"
namespace NumA {
VectorD LinAlgUtils::solve(const MatrixD & A, const VectorD & Y, Method method) {
VectorD LinAlgUtils::solve(const MatrixD & A, const VectorD & Y,
Method method) {
if (A.cols() != A.rows()) {
std::stringstream formatter;
formatter << "[Solve::solve] Matrix column size = " << A.cols()
......
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "../../../../include/NumA/linear_algebra/least_squares/lsmrBase.h"
#include <algorithm>
#include <cmath>
#include <iterator>
#include <numeric>
#include <vector>
namespace NumA {
inline void daxpy(unsigned int n, double alpha, const double * x, double * y) {
const double * xend = x + n;
while (x != xend) {
*y++ += alpha * *x++;
}
}
#define Sqr(x) ((x)*(x))
lsmrBase::lsmrBase() {
this->eps = 1e-16;
this->atol = 1e-6;
this->btol = 1e-6;
this->conlim = 1.0 / (10 * sqrt(this->eps));
this->itnlim = 10;
this->nout = NULL;
this->istop = 0;
this->itn = 0;
this->normA = 0.0;
this->condA = 0.0;
this->normr = 0.0;
this->normAr = 0.0;
this->normx = 0.0;
this->normb = 0.0;
this->dxmax = 0.0;
this->maxdx = 0;
this->damp = 0.0;
this->damped = false;
this->localSize = 0;
}
lsmrBase::~lsmrBase() {
}
unsigned int lsmrBase::getStoppingReason() const {
return this->istop;
}
std::string lsmrBase::getStoppingReasonMessage() const {
std::string msg;
switch (this->istop) {
case 0:
msg = "The exact solution is x = 0";
break;
case 1:
msg = "Ax - b is small enough, given atol, btol";
break;
case 2:
msg = "The least-squares solution is good enough, given atol";
break;
case 3:
msg = "The estimate of cond(Abar) has exceeded conlim";
break;
case 4:
msg = "Ax - b is small enough for this machine";
break;
case 5:
msg = "The LS solution is good enough for this machine";
break;
case 6:
msg = "Cond(Abar) seems to be too large for this machine";
break;
case 7:
msg = "The iteration limit has been reached";
break;
default:
msg = "Error. Unknown stopping reason";
break;
}
return msg;
}
unsigned int lsmrBase::getNumberOfIterationsPerformed() const {
return this->itn;
}
double lsmrBase::getFrobeniusNormEstimateOfAbar() const {
return this->normA;
}
double lsmrBase::getConditionNumberEstimateOfAbar() const {
return this->condA;
}
double lsmrBase::getFinalEstimateOfNormRbar() const {
return this->normr;
}
double lsmrBase::getFinalEstimateOfNormOfResiduals() const {
return this->normAr;
}
double lsmrBase::getFinalEstimateOfNormOfX() const {
return this->normx;
}
void lsmrBase::setLocalSize(unsigned int n) {
this->localSize = n;
}
void lsmrBase::setEpsilon(double value) {
this->eps = value;
}
void lsmrBase::setDamp(double value) {
this->damp = value;
}
void lsmrBase::setToleranceA(double value) {
this->atol = value;
}
void lsmrBase::setToleranceB(double value) {
this->btol = value;
}
void lsmrBase::setMaximumNumberOfIterations(unsigned int value) {
this->itnlim = value;
}
void lsmrBase::setUpperLimitOnConditional(double value) {
this->conlim = value;
}
void lsmrBase::setOutputStream(std::ostream & os) {
this->nout = &os;
}
/**
* returns sqrt( a**2 + b**2 )
* with precautions to avoid overflow.
*/
double lsmrBase::D2Norm(double a, double b) const {
const double scale = std::abs(a) + std::abs(b);
const double zero = 0.0;
if (scale == zero) {
return zero;
}
const double sa = a / scale;
const double sb = b / scale;
return scale * sqrt(sa * sa + sb * sb);
}
/** Simplified for this use from the BLAS version. */
void lsmrBase::scale(unsigned int n, double factor, double *x) const {
double * xend = x + n;
while (x != xend) {
*x++ *= factor;
}
}
double lsmrBase::Dnrm2(unsigned int n, const double *x) const {
double magnitudeOfLargestElement = 0.0;
double sumOfSquaresScaled = 1.0;
for (unsigned int i = 0; i < n; i++) {
if (x[i] != 0.0) {
double dx = x[i];
const double absxi = std::abs(dx);
if (magnitudeOfLargestElement < absxi) {
// rescale the sum to the range of the new element
dx = magnitudeOfLargestElement / absxi;
sumOfSquaresScaled = sumOfSquaresScaled * (dx * dx) + 1.0;
magnitudeOfLargestElement = absxi;
} else {
// rescale the new element to the range of the sum
dx = absxi / magnitudeOfLargestElement;
sumOfSquaresScaled += dx * dx;
}
}
}
const double norm = magnitudeOfLargestElement * sqrt(sumOfSquaresScaled);
return norm;
}
/**
*
* The array b must have size m
*
*/
void lsmrBase::solve(unsigned int m, unsigned int n, const VectorD& b,
VectorD& x) {
const double zero = 0.0;
const double one = 1.0;
double test1;
double test2;
// Initialize.
unsigned int localVecs = std::min(localSize, std::min(m, n));
if (this->nout) {
(*this->nout)
<< " Enter LSMR. Least-squares solution of Ax = b\n"
<< std::endl;
(*this->nout) << " The matrix A has " << m << " rows and " << n
<< " columns" << std::endl;
(*this->nout) << " damp = " << this->damp << std::endl;
(*this->nout) << " atol = " << this->atol << ", conlim = "
<< this->conlim << std::endl;
(*this->nout) << " btol = " << this->btol << ", itnlim = "
<< this->itnlim << std::endl;
(*this->nout)
<< " localSize (no. of vectors for local reorthogonalization) = "
<< this->localSize << std::endl;
}
int pfreq = 20;
int pcount = 0;
this->damped = (this->damp > zero);
std::vector<double> workBuffer(m + 5 * n + n * localVecs);
double * u = &workBuffer[0];
double * v = u + m;
double * w = v + n;
double * h = w + n;
double * hbar = h + n;
double * localV = hbar + n;
//-------------------------------------------------------------------
// Set up the first vectors u and v for the bidiagonalization.
// These satisfy beta*u = b, alpha*v = A(transpose)*u.
//-------------------------------------------------------------------
std::copy(b.toStdVector().begin(), b.toStdVector().begin() + m, u);
std::fill(v, v + n, zero);
std::fill(w, v + n, zero);
// std::fill(x, x + n, zero);
this->scale(m, (-1.0), u);
this->Aprod1(m, n, &(*(x.toStdVector().begin())), u);
this->scale(m, (-1.0), u);
double alpha = zero;
double beta = this->Dnrm2(m, u);
if (beta > zero) {
this->scale(m, (one / beta), u);
this->Aprod2(m, n, v, u); // v = A'*u
alpha = this->Dnrm2(n, v);
}
if (alpha > zero) {
this->scale(n, (one / alpha), v);
std::copy(v, v + n, w);
}
this->normAr = alpha * beta;
if (this->normAr == zero) {
this->TerminationPrintOut();
return;
}
// Initialization for local reorthogonalization.
bool localOrtho = false;
bool localVQueueFull = false;
unsigned int localPointer = 0;
if (localVecs > 0) {
localOrtho = true;
std::copy(v, v + n, localV);
}
// Initialize variables for 1st iteration.
this->itn = 0;
double zetabar = alpha * beta;
double alphabar = alpha;
double rho = one;
double rhobar = one;
double cbar = one;
double sbar = zero;
std::copy(v, v + n, h);
std::fill(hbar, hbar + n, zero);
// Initialize variables for estimation of ||r||.
double betadd = beta;
double betad = zero;
double rhodold = one;
double tautildeold = zero;
double thetatilde = zero;
double zeta = zero;
double d = zero;
// Initialize variables for estimation of ||A|| and cond(A).
double normA2 = alpha * alpha;
double maxrbar = zero;
double minrbar = 1e+100;
// Items for use in stopping rules.
this->normb = beta;
this->istop = 0;
double ctol = zero;
if (this->conlim > zero) {
ctol = one / this->conlim;
}
this->normr = beta;
if (this->nout) {
if (damped) {
(*this->nout)
<< " Itn x(1) norm rbar Abar'rbar"
" Compatible LS norm Abar cond Abar\n";
} else {
(*this->nout)
<< " Itn x(1) norm r A'r "
" Compatible LS norm A cond A\n";
}
test1 = one;
test2 = alpha / beta;
(*this->nout) << this->itn << ", " << x[0] << ", " << this->normr
<< ", " << this->normA << ", " << test1 << ", " << test2
<< std::endl;
}
// Main iteration loop
do {
this->itn++;
//----------------------------------------------------------------
// Perform the next step of the bidiagonalization to obtain the
// next beta, u, alpha, v. These satisfy
// beta*u = A*v - alpha*u,
// alpha*v = A'*u - beta*v.
//----------------------------------------------------------------
this->scale(m, (-alpha), u);
this->Aprod1(m, n, v, u); // u = u + A * v
beta = this->Dnrm2(m, u);
if (beta > zero) {
this->scale(m, (one / beta), u);
if (localOrtho) {
if (localPointer + 1 < localVecs) {
localPointer = localPointer + 1;
} else {
localPointer = 0;
localVQueueFull = true;
}
std::copy(v, v + n, localV + localPointer * n);
}
this->scale(n, (-beta), v);
this->Aprod2(m, n, v, u); // v = A'*u
if (localOrtho) {
unsigned int localOrthoLimit =
localVQueueFull ? localVecs : localPointer + 1;
for (unsigned int localOrthoCount = 0;
localOrthoCount < localOrthoLimit; ++localOrthoCount) {
double d = std::inner_product(v, v + n,
localV + n * localOrthoCount, 0.0);
daxpy(n, -d, localV + localOrthoCount * n, v);
}
}
alpha = this->Dnrm2(n, v);
if (alpha > zero) {
this->scale(n, (one / alpha), v);
}
}
// At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
//----------------------------------------------------------------
// Construct rotation Qhat_{k,2k+1}.
double alphahat = this->D2Norm(alphabar, damp);
double chat = alphabar / alphahat;
double shat = damp / alphahat;
// Use a plane rotation (Q_i) to turn B_i to R_i.
double rhoold = rho;
rho = D2Norm(alphahat, beta);
double c = alphahat / rho;
double s = beta / rho;
double thetanew = s * alpha;
alphabar = c * alpha;
// Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar.
double rhobarold = rhobar;
double zetaold = zeta;
double thetabar = sbar * rho;
double rhotemp = cbar * rho;
rhobar = this->D2Norm(cbar * rho, thetanew);
cbar = cbar * rho / rhobar;
sbar = thetanew / rhobar;
zeta = cbar * zetabar;
zetabar = -sbar * zetabar;
// Update h, h_hat, x.
for (unsigned int i = 0; i < n; ++i) {
hbar[i] = h[i] - (thetabar * rho / (rhoold * rhobarold)) * hbar[i];
x[i] = x[i] + (zeta / (rho * rhobar)) * hbar[i];
h[i] = v[i] - (thetanew / rho) * h[i];
}
// Estimate ||r||.
// Apply rotation Qhat_{k,2k+1}.
double betaacute = chat * betadd;
double betacheck = -shat * betadd;
// Apply rotation Q_{k,k+1}.
double betahat = c * betaacute;
betadd = -s * betaacute;
// Apply rotation Qtilde_{k-1}.
// betad = betad_{k-1} here.
double thetatildeold = thetatilde;
double rhotildeold = this->D2Norm(rhodold, thetabar);
double ctildeold = rhodold / rhotildeold;
double stildeold = thetabar / rhotildeold;
thetatilde = stildeold * rhobar;
rhodold = ctildeold * rhobar;
betad = -stildeold * betad + ctildeold * betahat;
// betad = betad_k here.
// rhodold = rhod_k here.
tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold;