Commit e9e9c0b5 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

eigen_solution_of_symmetric_matrix(): changed type parameter eig_val to Varray<double> &.

parent 941aa212
Pipeline #4862 passed with stages
in 16 minutes and 43 seconds
......@@ -513,7 +513,7 @@ EOFs(void *process)
/* SOLVE THE EIGEN PROBLEM */
if (Options::Timer) timer_start(timer_eig);
auto eig_val = eofdata[varID][levelID].eig_val.data();
auto &eig_val = eofdata[varID][levelID].eig_val;
if (eigen_mode == JACOBI)
// TODO: use return status (>0 okay, -1 did not converge at all)
cdo::parallel_eigen_solution_of_symmetric_matrix(covar, eig_val, n, __func__);
......
......@@ -314,9 +314,9 @@ EOF3d(void *process)
if (Options::cdoVerbose) cdoPrint("Processed correlation matrix for var %2i | npack: %zu", varID, n);
if (eigen_mode == JACOBI)
cdo::parallel_eigen_solution_of_symmetric_matrix(cov, eigv.data(), n, __func__);
cdo::parallel_eigen_solution_of_symmetric_matrix(cov, eigv, n, __func__);
else
cdo::eigen_solution_of_symmetric_matrix(cov, eigv.data(), n, __func__);
cdo::eigen_solution_of_symmetric_matrix(cov, eigv, n, __func__);
// NOW: cov contains the eigenvectors, eigv the eigenvalues
if (Options::cdoVerbose) cdoPrint("Processed SVD decomposition for var %d from %d x %d matrix", varID, n, n);
......
......@@ -56,7 +56,7 @@ namespace cdo
{
void
heap_sort(double *eig_val, double **a, int n)
heap_sort(Varray<double> &eig_val, double **a, int n)
{
int j_next;
......@@ -94,8 +94,8 @@ heap_sort(double *eig_val, double **a, int n)
}
}
void
make_symmetric_matrix_triangular(double **a, int n, double *d, double *e)
static void
make_symmetric_matrix_triangular(double **a, int n, Varray<double> &d, Varray<double> &e)
{
int i, j, k;
double f, g, h, hh, scale;
......@@ -186,8 +186,8 @@ pythagoras(double a, double b)
#define MAX_ITER 1000
void
eigen_solution_of_triangular_matrix(double *d, double *e, int n, double **a, const char *prompt)
static void
eigen_solution_of_triangular_matrix(Varray<double> &d, Varray<double> &e, int n, double **a, const char *prompt)
{
constexpr double eps = 1.e-6;
int i, k, l, m, iter;
......@@ -257,13 +257,13 @@ eigen_solution_of_triangular_matrix(double *d, double *e, int n, double **a, con
}
void
eigen_solution_of_symmetric_matrix(double **a, double *eig_val, int n, const char *prompt)
eigen_solution_of_symmetric_matrix(double **a, Varray<double> &eig_val, int n, const char *prompt)
// After return the rows (!!!) of a are the eigenvectors
{
{
std::vector<double> e(n);
make_symmetric_matrix_triangular(a, n, eig_val, e.data());
eigen_solution_of_triangular_matrix(eig_val, e.data(), n, a, prompt);
make_symmetric_matrix_triangular(a, n, eig_val, e);
eigen_solution_of_triangular_matrix(eig_val, e, n, a, prompt);
}
for (int i = 0; i < n; i++)
......@@ -424,8 +424,8 @@ fft(double *restrict real, double *restrict imag, int n, int sign)
{
const auto w_r = std::cos(M_PI / step);
const auto w_i = std::sin(M_PI / step) * sign;
double ww_r = 1;
double ww_i = 0;
double ww_r = 1.0;
double ww_i = 0.0;
for (int i = 0; i < step; i++)
{
double temp_r, temp_i;
......@@ -444,7 +444,7 @@ fft(double *restrict real, double *restrict imag, int n, int sign)
}
}
const auto norm = 1. / std::sqrt(n);
const auto norm = 1.0 / std::sqrt(n);
for (int i = 0; i < n; i++) real[i] *= norm;
for (int i = 0; i < n; i++) imag[i] *= norm;
}
......@@ -471,10 +471,10 @@ ft(double *real, double *imag, int n, int sign)
{
const auto w_r = std::cos(2 * M_PI * k / n);
const auto w_i = std::sin(2 * M_PI * k / n) * sign;
double ww_r = 1;
double ww_i = 0;
double sum_r = 0;
double sum_i = 0;
double ww_r = 1.0;
double ww_i = 0.0;
double sum_r = 0.0;
double sum_i = 0.0;
for (int j = 0; j < n; j++)
{
sum_r += real[j] * ww_r - imag[j] * ww_i;
......@@ -503,10 +503,10 @@ ft_r(double *restrict real, double *restrict imag, int n, int sign, double *rest
{
const auto w_r = std::cos(2 * M_PI * k / n);
const auto w_i = std::sin(2 * M_PI * k / n) * sign;
double ww_r = 1;
double ww_i = 0;
double sum_r = 0;
double sum_i = 0;
double ww_r = 1.0;
double ww_i = 0.0;
double sum_r = 0.0;
double sum_i = 0.0;
for (int j = 0; j < n; j++)
{
sum_r += real[j] * ww_r - imag[j] * ww_i;
......@@ -519,7 +519,7 @@ ft_r(double *restrict real, double *restrict imag, int n, int sign, double *rest
work_i[k] = sum_i;
}
const auto norm = 1. / std::sqrt(n);
const auto norm = 1.0 / std::sqrt(n);
for (int k = 0; k < n; k++) real[k] = work_r[k] * norm;
for (int k = 0; k < n; k++) imag[k] = work_i[k] * norm;
}
......@@ -545,7 +545,7 @@ gamma_help_1(double a, double x, const char *prompt)
double gln = lngamma(a);
double ap = a;
double sum, del = sum = 1 / a;
double sum, del = sum = 1.0 / a;
for (int i = 1; i <= 100; i++)
{
......@@ -577,7 +577,7 @@ gamma_help_2(double a, double x, const char *prompt)
for (int i = 1; i <= 100; i++)
{
an = -i * (i - a);
b += 2;
b += 2.0;
d = an * d + b;
if (std::fabs(d) < very_small) d = very_small;
c = b + an / c;
......@@ -597,21 +597,21 @@ gamma_help_2(double a, double x, const char *prompt)
double
incomplete_gamma(double a, double x, const char *prompt)
{
if (x < 0 || a <= 0)
if (x < 0.0 || a <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"incomplete_gamma\")\n", prompt);
exit(4);
}
if (x < (a + 1))
if (x < (a + 1.0))
return gamma_help_1(a, x, prompt);
else
return 1 - gamma_help_2(a, x, prompt);
return 1.0 - gamma_help_2(a, x, prompt);
}
double
beta(double a, double b, const char *prompt)
{
if (a <= 0 || b <= 0)
if (a <= 0.0 || b <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"beta\")\n", prompt);
exit(4);
......@@ -623,15 +623,15 @@ beta(double a, double b, const char *prompt)
static double
beta_help(double a, double b, double x, const char *prompt)
{
constexpr double very_small = 1000 * DBL_MIN;
constexpr double very_small = 1000.0 * DBL_MIN;
constexpr double eps = 3.e-07;
double aa, del;
const auto qab = a + b;
const auto qap = a + 1;
const auto qam = a - 1;
double c = 1;
double d = 1 - qab * x / qap;
double c = 1.0;
double d = 1.0 - qab * x / qap;
if (std::fabs(d) < very_small) d = very_small;
d = 1 / d;
double h = d;
......@@ -641,16 +641,16 @@ beta_help(double a, double b, double x, const char *prompt)
aa = m * (b - m) * x / ((qam + m2) * (a + m2));
d = 1 + aa * d;
if (std::fabs(d) < very_small) d = very_small;
c = 1 + aa / c;
c = 1.0 + aa / c;
if (std::fabs(c) < very_small) c = very_small;
d = 1 / d;
d = 1.0 / d;
h *= d * c;
aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
d = 1 + aa * d;
d = 1.0 + aa * d;
if (std::fabs(d) < very_small) d = very_small;
c = 1 + aa / c;
c = 1.0 + aa / c;
if (std::fabs(c) < very_small) c = very_small;
d = 1 / d;
d = 1.0 / d;
del = d * c;
h *= del;
if (std::fabs(del - 1.0) < eps) return h;
......@@ -665,49 +665,49 @@ beta_help(double a, double b, double x, const char *prompt)
double
incomplete_beta(double a, double b, double x, const char *prompt)
{
if (a <= 0 || b <= 0)
if (a <= 0.0 || b <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"incomplete_beta\")\n", prompt);
exit(4);
}
if (x < 0 || x > 1)
if (x < 0.0 || x > 1.0)
{
fprintf(stderr, "%s: Value out of range (0-1)!\n", prompt);
exit(4);
}
double c = (DBL_IS_EQUAL(x, 0.) || DBL_IS_EQUAL(x, 1.))
? 0
? 0.0
: std::exp(lngamma(a + b) - lngamma(a) - lngamma(b) + a * std::log(x) + b * std::log(1 - x));
if (x < (a + 1) / (a + b + 2))
if (x < (a + 1) / (a + b + 2.0))
return c * beta_help(a, b, x, prompt) / a;
else
return 1 - c * beta_help(b, a, 1 - x, prompt) / b;
return 1.0 - c * beta_help(b, a, 1.0 - x, prompt) / b;
}
double
normal_density(double x)
{
return M_2_SQRTPI / 2 / M_SQRT2 * std::exp(-x * x / 2);
return M_2_SQRTPI / 2.0 / M_SQRT2 * std::exp(-x * x / 2.0);
}
double
normal(double x, const char *prompt)
{
return x > 0 ? 0.5 * (1 + incomplete_gamma(0.5, x * x / 2, prompt))
: x < 0 ? 0.5 * (1 - incomplete_gamma(0.5, x * x / 2, prompt)) : 0.5;
return x > 0.0 ? 0.5 * (1.0 + incomplete_gamma(0.5, x * x / 2.0, prompt))
: x < 0.0 ? 0.5 * (1 - incomplete_gamma(0.5, x * x / 2.0, prompt)) : 0.5;
}
double
normal_inv(double p, const char *prompt)
{
constexpr double eps = 1.e-10;
static double last_p = 0.5, last_x = 0;
static double last_p = 0.5, last_x = 0.0;
double x, xx;
if (p <= 0 || p >= 1)
if (p <= 0.0 || p >= 1.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"normal_inv\")\n", prompt);
exit(4);
......@@ -720,7 +720,7 @@ normal_inv(double p, const char *prompt)
return -normal_inv(1 - p, prompt);
else
{
x = 0;
x = 0.0;
while (true)
{
xx = x - (normal(x, prompt) - p) / normal_density(x);
......@@ -736,7 +736,7 @@ normal_inv(double p, const char *prompt)
double
student_t_density(double n, double x, const char *prompt)
{
if (n <= 0)
if (n <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"student_t_density\")\n", prompt);
exit(4);
......@@ -748,16 +748,16 @@ student_t_density(double n, double x, const char *prompt)
double
student_t(double n, double x, const char *prompt)
{
if (n <= 0)
if (n <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"student_t\")\n", prompt);
exit(4);
}
if (x > 0)
return 1 - 0.5 * incomplete_beta(n / 2, 0.5, n / (n + x * x), prompt);
return 1.0 - 0.5 * incomplete_beta(n / 2.0, 0.5, n / (n + x * x), prompt);
else if (x < 0)
return 0.5 * incomplete_beta(n / 2, 0.5, n / (n + x * x), prompt);
return 0.5 * incomplete_beta(n / 2.0, 0.5, n / (n + x * x), prompt);
else
return 0.5;
}
......@@ -766,10 +766,10 @@ double
student_t_inv(double n, double p, const char *prompt)
{
constexpr double eps = 1.e-10;
static double last_n = 1, last_p = 0.5, last_x = 0;
static double last_n = 1, last_p = 0.5, last_x = 0.0;
double x, xx;
if (n <= 0 || p <= 0 || p >= 1)
if (n <= 0.0 || p <= 0.0 || p >= 1.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"student_t_inv\")\n", prompt);
exit(4);
......@@ -777,12 +777,12 @@ student_t_inv(double n, double p, const char *prompt)
if (DBL_IS_EQUAL(n, last_n) && DBL_IS_EQUAL(p, last_p)) return last_x;
if (DBL_IS_EQUAL(p, 0.5))
return 0;
return 0.0;
else if (p < 0.5)
return -student_t_inv(n, 1 - p, prompt);
return -student_t_inv(n, 1.0 - p, prompt);
else
{
x = 0;
x = 0.0;
while (true)
{
xx = x - (student_t(n, x, prompt) - p) / student_t_density(n, x, prompt);
......@@ -799,7 +799,7 @@ student_t_inv(double n, double p, const char *prompt)
double
chi_square_density(double n, double x, const char *prompt)
{
if (n <= 0)
if (n <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"chi_square_density\")\n", prompt);
exit(4);
......@@ -811,24 +811,24 @@ chi_square_density(double n, double x, const char *prompt)
double
chi_square(double n, double x, const char *prompt)
{
if (n <= 0)
if (n <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"chi_square\")\n", prompt);
exit(4);
}
return x <= 0 ? 0 : incomplete_gamma(n / 2, x / 2, prompt);
return x <= 0.0 ? 0.0 : incomplete_gamma(n / 2.0, x / 2.0, prompt);
}
double
chi_square_inv(double n, double p, const char *prompt)
{
constexpr double eps = 1.e-10;
static double last_n = -1, last_p = -1, last_x = -1;
static double last_last_n = -1, last_last_p = -1, last_last_x = -1;
static double last_n = -1.0, last_p = -1.0, last_x = -1.0;
static double last_last_n = -1.0, last_last_p = -1.0, last_last_x = -1.0;
double x, xx;
if (n <= 0 || p <= 0 || p >= 1)
if (n <= 0.0 || p <= 0.0 || p >= 1.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"chi_square_inv\")\n", prompt);
exit(4);
......@@ -844,7 +844,7 @@ chi_square_inv(double n, double p, const char *prompt)
xx = x - (chi_square(n, x, prompt) - p) / chi_square_density(n, x, prompt);
if (std::fabs(xx - x) < x * eps) break;
if (xx < 0)
x /= 2;
x /= 2.0;
else
x = xx;
}
......@@ -865,7 +865,7 @@ chi_square_constants(double n, double p, double *c1, double *c2, const char *pro
constexpr double eps = 1.e-10;
static double last_n, last_p, last_c1, last_c2;
if (n <= 0 || p <= 0 || p >= 1)
if (n <= 0.0 || p <= 0.0 || p >= 1.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"chi_square_constants\")\n", prompt);
exit(4);
......@@ -895,13 +895,13 @@ chi_square_constants(double n, double p, double *c1, double *c2, const char *pro
const auto delta_c2 = (b2 * a11 - b1 * a21) / det;
if (std::fabs(delta_c1) < *c1 * eps && std::fabs(delta_c2) < *c2 * eps) break;
if (*c1 + delta_c1 >= n)
*c1 = (n + *c1) / 2;
*c1 = (n + *c1) / 2.0;
else if (*c1 + delta_c1 <= 0)
*c1 /= 2;
*c1 /= 2.0;
else
*c1 += delta_c1;
if (*c2 + delta_c2 <= n)
*c2 = (n + *c2) / 2;
*c2 = (n + *c2) / 2.0;
else
*c2 += delta_c2;
}
......@@ -915,7 +915,7 @@ chi_square_constants(double n, double p, double *c1, double *c2, const char *pro
double
beta_distr_density(double a, double b, double x, const char *prompt)
{
if (a <= 0 || b <= 0)
if (a <= 0.0 || b <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"beta_distr_density\")\n", prompt);
exit(4);
......@@ -934,11 +934,11 @@ double
beta_distr_inv(double a, double b, double p, const char *prompt)
{
constexpr double eps = 1.e-10;
static double last_a = -1, last_b, last_p = -1, last_x = -1;
static double last_last_a = -1, last_last_b = -1, last_last_p = -1, last_last_x = -1;
static double last_a = -1.0, last_b, last_p = -1.0, last_x = -1.0;
static double last_last_a = -1.0, last_last_b = -1.0, last_last_p = -1.0, last_last_x = -1.0;
double xx, x;
if (a <= 0 || b <= 0 || p <= 0 || p >= 1)
if (a <= 0.0 || b <= 0.0 || p <= 0.0 || p >= 1.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"beta_distr_inv\")\n", prompt);
exit(4);
......@@ -953,16 +953,16 @@ beta_distr_inv(double a, double b, double p, const char *prompt)
{
xx = x - (beta_distr(a, b, x, prompt) - p) / beta_distr_density(a, b, x, prompt);
if (std::fabs(xx - x) < x * eps) break;
if (xx <= 0)
if (xx <= 0.0)
x /= 2;
else if (xx >= 1)
x = (1 + x) / 2;
else if (xx >= 1.0)
x = (1.0 + x) / 2.0;
else
x = xx;
}
#if 0
for (x_l = 0, x_r = 1; fabs (x_l - x_r) > eps;
x = (x_l+x_r) / 2, beta_distr (a, b, x, prompt) < p ? (x_l=x):(x_r=x));
x = (x_l+x_r) / 2.0, beta_distr (a, b, x, prompt) < p ? (x_l=x):(x_r=x));
#endif
last_last_a = last_a;
......@@ -983,7 +983,7 @@ beta_distr_constants(double a, double b, double p, double *c1, double *c2, const
constexpr double eps = 1.e-10;
static double last_a, last_b, last_p, last_c1, last_c2;
if (a <= 0 || b <= 0 || p <= 0 || p >= 1)
if (a <= 0.0 || b <= 0.0 || p <= 0.0 || p >= 1.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"beta_distr_constants\")\n", prompt);
exit(4);
......@@ -1000,8 +1000,8 @@ beta_distr_constants(double a, double b, double p, double *c1, double *c2, const
*c1 = a / (a + b);
*c2 = a / (a + b);
#endif
*c1 = beta_distr_inv(a, b, p / 2, prompt);
*c2 = beta_distr_inv(a, b, 1 - p / 2, prompt);
*c1 = beta_distr_inv(a, b, p / 2.0, prompt);
*c2 = beta_distr_inv(a, b, 1.0 - p / 2.0, prompt);
while (true)
{
......@@ -1017,15 +1017,15 @@ beta_distr_constants(double a, double b, double p, double *c1, double *c2, const
const auto delta_c2 = (b2 * a11 - b1 * a21) / det;
if (std::fabs(delta_c1) < *c1 * eps && std::fabs(delta_c2) < *c2 * eps) break;
if (*c1 + delta_c1 >= a / (a + b))
*c1 = (a / (a + b) + *c1) / 2;
*c1 = (a / (a + b) + *c1) / 2.0;
else if (*c1 + delta_c1 <= 0)
*c1 /= 2;
*c1 /= 2.0;
else
*c1 += delta_c1;
if (*c2 + delta_c2 >= 1)
*c2 = (1 + *c2) / 2;
if (*c2 + delta_c2 >= 1.0)
*c2 = (1.0 + *c2) / 2.0;
else if (*c2 + delta_c2 <= a / (a + b))
*c2 = (a / (a + b) + *c2) / 2;
*c2 = (a / (a + b) + *c2) / 2.0;
else
*c2 += delta_c2;
}
......@@ -1040,13 +1040,13 @@ beta_distr_constants(double a, double b, double p, double *c1, double *c2, const
double
fisher(double m, double n, double x, const char *prompt)
{
if (m <= 0 || n <= 0)
if (m <= 0.0 || n <= 0.0)
{
fprintf(stderr, "%s: IMPLEMENTATION ERROR! (Invalid argument in function \"fisher\")\n", prompt);
exit(4);
}
return incomplete_beta(m / 2, n / 2, n / (n + m * x), prompt);
return incomplete_beta(m / 2.0, n / 2.0, n / (n + m * x), prompt);
}
/* ******************************************************************************** */
......@@ -1093,10 +1093,10 @@ annihilate_1side(double **M, long i, long j, long n)
return;
}
const auto zeta = (beta - alpha) / (2. * gamma); // tan(2*theta)
auto tk = 1. / (std::fabs(zeta) + std::sqrt(1. + zeta * zeta));
const auto zeta = (beta - alpha) / (2.0 * gamma); // tan(2*theta)
auto tk = 1.0 / (std::fabs(zeta) + std::sqrt(1. + zeta * zeta));
tk = zeta > 0 ? tk : -tk; // = cot(2*theta)
const auto ck = 1. / std::sqrt(1. + tk * tk); // = cos(theta)
const auto ck = 1.0 / std::sqrt(1. + tk * tk); // = cos(theta)
const auto sk = ck * tk; // = sin(theta)
// calculate a_i,j - tilde
......@@ -1110,7 +1110,7 @@ annihilate_1side(double **M, long i, long j, long n)
}
static int
jacobi_1side(double **M, double *A, long n)
jacobi_1side(double **M, Varray<double> &A, long n)
{
Varray2D<int> annihilations(2);
annihilations[0].resize(n * n);
......@@ -1222,7 +1222,7 @@ jacobi_1side(double **M, double *A, long n)
for (long i = 0; i < n; i++)
{
memset(M[i], 0, n * sizeof(double));
memset(A, 0, n * sizeof(double));
memset(A.data(), 0, n * sizeof(double));
}
return -1;
}
......@@ -1249,7 +1249,7 @@ jacobi_1side(double **M, double *A, long n)
/* ******************************************************************************** */
void
parallel_eigen_solution_of_symmetric_matrix(double **M, double *A, int n, const char func[])
parallel_eigen_solution_of_symmetric_matrix(double **M, Varray<double> &A, int n, const char func[])
{
(void) (func); // CDO_UNUSED
......
......@@ -7,7 +7,7 @@
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; version 2 of the License.
the Free Software Fou<ndation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
......@@ -17,10 +17,12 @@
#ifndef STATISTIC_H
#define STATISTIC_H
#include "array.h"
namespace cdo
{
void eigen_solution_of_symmetric_matrix(double **a, double *eig_val, int n, const char *prompt);
void eigen_solution_of_symmetric_matrix(double **a, Varray<double> &eig_val, int n, const char *prompt);
void fft(double *real, double *imag, int n, int sign);
void ft_r(double *real, double *imag, int n, int sign, double *work_r, double *work_i);
double lngamma(double x);
......@@ -44,7 +46,7 @@ void beta_distr_constants(double a, double b, double p, double *c1, double *c2,
double fisher(double m, double n, double x, const char *prompt);
// make parallel eigen solution accessible for eigen value computation in EOF3d.c
void parallel_eigen_solution_of_symmetric_matrix(double **M, double *A, int n, const char func[]);
void parallel_eigen_solution_of_symmetric_matrix(double **M, Varray<double> &A, int n, const char func[]);
}
......
Markdown is supported
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