Commit 03989c58 authored by Nabil Chouika's avatar Nabil Chouika
Browse files

refs#16

In NumA++:
- Initialize Eigen vectors to zero in LSMRSolver, to avoid random results (and maybe nan).

In PARTONS_Fits:
- Propagate changes in parameters names for DVCSConstantCFFModel.
parent 7604c6f0
......@@ -87,7 +87,6 @@ private:
std::string m_hdg1;
std::string m_hdg2;
unsigned int m_pfreq; // print frequency (for repeating the heading)
size_t m_pcount; // print counter
/**
* Routine of rotation from scipy.
......
......@@ -4,10 +4,7 @@
#include <cmath>
#include <ctime>
#include <limits>
#include <string>
#include <vector>
#include <ElementaryUtils/string_utils/Formatter.h>
#include <ctime>
#include "../../../../include/NumA/linear_algebra/eigen/EigenUtils.h"
......@@ -17,8 +14,7 @@ LSMRSolver::LSMRSolver(double damp, double atol, double btol, double conlim,
size_t maxiter, ElemUtils::Formatter* output) :
m_damp(damp), m_atol(atol), m_btol(btol), m_conlim(conlim), m_maxiter(
maxiter), m_output(output), m_istop(0), m_itn(0), m_normr(0.), m_normar(
0.), m_normA(0.), m_condA(0.), m_normx(0.), m_pfreq(20), m_pcount(
0) {
0.), m_normA(0.), m_condA(0.), m_normx(0.), m_pfreq(20) {
m_msg.assign(8, "");
m_msg[0] = "The exact solution is x = 0.";
m_msg[1] = "Ax - b is small enough, given atol, btol.";
......@@ -68,6 +64,7 @@ VectorD LSMRSolver::solve(const MatrixD& A, const VectorD& B) {
// VectorD v(n);
Eigen::VectorXd v(n);
v.setZero();
double alpha = 0.;
if (beta > 0) {
......@@ -83,6 +80,7 @@ VectorD LSMRSolver::solve(const MatrixD& A, const VectorD& B) {
// Initialize variables for 1st iteration.
unsigned int pcount = 0; // print counter
m_itn = 0;
double zetabar = alpha * beta;
double alphabar = alpha;
......@@ -96,7 +94,9 @@ VectorD LSMRSolver::solve(const MatrixD& A, const VectorD& B) {
// VectorD x(n);
Eigen::VectorXd h(v);
Eigen::VectorXd hbar(n);
hbar.setZero();
Eigen::VectorXd x(n);
x.setZero();
// Initialize variables for estimation of ||r||.
......@@ -141,7 +141,7 @@ VectorD LSMRSolver::solve(const MatrixD& A, const VectorD& B) {
*m_output << m_hdg1 << " " << m_hdg2 << "\n";
test1 = 1.;
test2 = alpha / beta;
*m_output << m_itn << " " << x[0] << " " << m_normr << " " << m_normar
*m_output << m_itn << " " << x(0) << " " << m_normr << " " << m_normar
<< " " << test1 << " " << test2 << "\n";
}
......@@ -249,7 +249,7 @@ VectorD LSMRSolver::solve(const MatrixD& A, const VectorD& B) {
// Test for convergence.
// Compute norms for convergence testing.
m_normar = std::abs(zetabar);
m_normar = fabs(zetabar);
m_normx = x.norm();
// Now use these norms to estimate certain other quantities,
......@@ -298,16 +298,31 @@ VectorD LSMRSolver::solve(const MatrixD& A, const VectorD& B) {
or (m_itn % 10 == 0) or (test3 <= 1.1 * ctol)
or (test2 <= 1.1 * m_atol) or (test1 <= 1.1 * rtol)
or (m_istop != 0)) {
if (m_pcount >= m_pfreq) {
m_pcount = 0;
if (pcount >= m_pfreq) {
pcount = 0;
*m_output << "\n";
*m_output << m_hdg1 << m_hdg2 << "\n";
}
m_pcount++;
*m_output << m_itn << " " << x[0] << " " << m_normr << " "
<< m_normar << " " << test1 << " " << test2 << " "
<< m_normA << " " << m_condA << " " << duration << "s"
<< "\n";
pcount++;
*m_output << m_itn;
*m_output << " ";
*m_output << x(0);
*m_output << " ";
*m_output << m_normr;
*m_output << " ";
*m_output << m_normar;
*m_output << " ";
*m_output << test1;
*m_output << " ";
*m_output << test2;
*m_output << " ";
*m_output << m_normA;
*m_output << " ";
*m_output << m_condA;
*m_output << " ";
*m_output << duration;
*m_output << "s";
*m_output << "\n";
}
}
......@@ -420,12 +435,12 @@ void LSMRSolver::_sym_ortho(double a, double b, double& c, double& s,
if (b == 0) {
c = (a > 0) - (a < 0);
s = 0;
r = std::abs(a);
r = fabs(a);
} else if (a == 0) {
c = 0;
s = (b > 0) - (b < 0);
r = std::abs(b);
} else if (std::abs(b) > std::abs(a)) {
r = fabs(b);
} else if (fabs(b) > fabs(a)) {
double tau = a / b;
s = ((b > 0) - (b < 0)) / std::sqrt(1. + tau * tau);
c = s * tau;
......
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