Commit 34aa6f15 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

eigen_solution_of_symmetric_matrix(): changed type of parameter n to size_t.

parent 191bfe69
Pipeline #4864 failed with stages
in 16 minutes and 30 seconds
......@@ -117,7 +117,6 @@ EOFs(void *process)
size_t nmiss;
int varID, levelID;
int nts = 1;
size_t n = 0;
int grid_space = 0, time_space = 0;
int timer_cov = 0, timer_eig = 0;
......@@ -213,6 +212,7 @@ EOFs(void *process)
}
// reset the requested number of eigen-function to the maximum if neccessary
size_t n = 0;
if (time_space)
{
if (n_eig > nts)
......
......@@ -61,7 +61,6 @@ EOF3d(void *process)
bool missval_warning = false;
size_t nmiss;
int ngrids;
int n = 0;
int timer_cov = 0, timer_eig = 0;
int calendar = CALENDAR_STANDARD;
......@@ -129,7 +128,7 @@ EOF3d(void *process)
n_eig = nts;
}
n = nts;
size_t n = nts;
if (Options::cdoVerbose) cdoPrint("counted %i timesteps", n);
......
......@@ -51,7 +51,7 @@ fourier_fftw(int sign, int varID, int levelID, int nts, size_t gridsize, double
std::vector<FourierMemory> &fourierMemory)
{
#ifdef HAVE_LIBFFTW3
const auto norm = 1. / std::sqrt(nts);
const auto norm = 1.0 / std::sqrt(nts);
#ifdef _OPENMP
#pragma omp parallel for default(shared)
......
......@@ -50,7 +50,7 @@ constexpr int MAX_JACOBI_ITER = 12;
static double fnorm_precision;
static int max_jacobi_iter;
static int n_finished;
static size_t n_finished;
namespace cdo
{
......@@ -187,10 +187,10 @@ pythagoras(double a, double b)
#define MAX_ITER 1000
static void
eigen_solution_of_triangular_matrix(Varray<double> &d, Varray<double> &e, long n, Varray2D<double> &a, const char *prompt)
eigen_solution_of_triangular_matrix(Varray<double> &d, Varray<double> &e, size_t n, Varray2D<double> &a, const char *prompt)
{
constexpr double eps = 1.e-6;
long i, k, l, m, iter;
size_t i, k, l, m, iter;
double b, c, f, g, p, r, s;
for (i = 1; i < n; i++) e[i - 1] = e[i];
......@@ -257,7 +257,7 @@ eigen_solution_of_triangular_matrix(Varray<double> &d, Varray<double> &e, long n
}
void
eigen_solution_of_symmetric_matrix(Varray2D<double> &a, Varray<double> &eig_val, long n, const char *prompt)
eigen_solution_of_symmetric_matrix(Varray2D<double> &a, Varray<double> &eig_val, size_t n, const char *prompt)
// After return the rows (!!!) of a are the eigenvectors
{
{
......@@ -266,8 +266,8 @@ eigen_solution_of_symmetric_matrix(Varray2D<double> &a, Varray<double> &eig_val,
eigen_solution_of_triangular_matrix(eig_val, e, n, a, prompt);
}
for (long i = 0; i < n; i++)
for (long j = i + 1; j < n; j++) std::swap(a[i][j], a[j][i]);
for (size_t i = 0; i < n; i++)
for (size_t j = i + 1; j < n; j++) std::swap(a[i][j], a[j][i]);
heap_sort(eig_val, a, n);
}
......@@ -402,15 +402,15 @@ inverse_of_matrix(double **a, double **b, int n)
}
*/
void
fft(double *restrict real, double *restrict imag, long n, int sign)
fft(double *restrict real, double *restrict imag, int n, int sign)
{
// n must be a power of 2
// sign should be 1 (FT) or -1 (reverse FT)
long j, j1, j2;
long bit;
int j, j1, j2;
int bit;
// Bit reversal part
for (long i = j = 0; i < n; i++) // The bit pattern of i and j are reverse
for (int i = j = 0; i < n; i++) // The bit pattern of i and j are reverse
{
if (i > j) std::swap(real[i], real[j]);
if (i > j) std::swap(imag[i], imag[j]);
......@@ -420,13 +420,13 @@ fft(double *restrict real, double *restrict imag, long n, int sign)
}
// Danielson-Lanczos Part
for (long step = 1; step < n; step <<= 1)
for (int step = 1; step < n; step <<= 1)
{
const auto w_r = std::cos(M_PI / step);
const auto w_i = std::sin(M_PI / step) * sign;
double ww_r = 1.0;
double ww_i = 0.0;
for (long i = 0; i < step; i++)
for (int i = 0; i < step; i++)
{
double temp_r, temp_i;
for (j1 = i, j2 = i + step; j2 < n; j1 += 2 * step, j2 += 2 * step)
......@@ -445,12 +445,12 @@ fft(double *restrict real, double *restrict imag, long n, int sign)
}
const auto norm = 1.0 / std::sqrt(n);
for (long i = 0; i < n; i++) real[i] *= norm;
for (long i = 0; i < n; i++) imag[i] *= norm;
for (int i = 0; i < n; i++) real[i] *= norm;
for (int i = 0; i < n; i++) imag[i] *= norm;
}
void
ft(double *real, double *imag, long n, int sign)
ft(double *real, double *imag, int n, int sign)
{
// sign should be 1 (FT) or -1 (reverse FT)
static double *work_r = 0, *work_i = 0;
......@@ -467,7 +467,7 @@ ft(double *real, double *imag, long n, int sign)
/* free_at_exit (work_i); */
}
for (long k = 0; k < n; k++)
for (int k = 0; k < n; k++)
{
const auto w_r = std::cos(2 * M_PI * k / n);
const auto w_i = std::sin(2 * M_PI * k / n) * sign;
......@@ -475,7 +475,7 @@ ft(double *real, double *imag, long n, int sign)
double ww_i = 0.0;
double sum_r = 0.0;
double sum_i = 0.0;
for (long j = 0; j < n; j++)
for (int j = 0; j < n; j++)
{
sum_r += real[j] * ww_r - imag[j] * ww_i;
sum_i += real[j] * ww_i + imag[j] * ww_r;
......@@ -488,18 +488,18 @@ ft(double *real, double *imag, long n, int sign)
}
const auto norm = 1. / std::sqrt(n);
for (long k = 0; k < n; k++) real[k] = work_r[k] * norm;
for (long k = 0; k < n; k++) imag[k] = work_i[k] * norm;
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;
}
// reentrant version of ft
void
ft_r(double *restrict real, double *restrict imag, long n, int sign, double *restrict work_r, double *restrict work_i)
ft_r(double *restrict real, double *restrict imag, int n, int sign, double *restrict work_r, double *restrict work_i)
{
// sign should be 1 (FT) or -1 (reverse FT)
double temp_r;
for (long k = 0; k < n; k++)
for (int k = 0; k < n; k++)
{
const auto w_r = std::cos(2 * M_PI * k / n);
const auto w_i = std::sin(2 * M_PI * k / n) * sign;
......@@ -507,7 +507,7 @@ ft_r(double *restrict real, double *restrict imag, long n, int sign, double *res
double ww_i = 0.0;
double sum_r = 0.0;
double sum_i = 0.0;
for (long j = 0; j < n; j++)
for (int j = 0; j < n; j++)
{
sum_r += real[j] * ww_r - imag[j] * ww_i;
sum_i += real[j] * ww_i + imag[j] * ww_r;
......@@ -520,8 +520,8 @@ ft_r(double *restrict real, double *restrict imag, long n, int sign, double *res
}
const auto norm = 1.0 / std::sqrt(n);
for (long k = 0; k < n; k++) real[k] = work_r[k] * norm;
for (long k = 0; k < n; k++) imag[k] = work_i[k] * norm;
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;
}
double
......@@ -1058,7 +1058,7 @@ fisher(double m, double n, double x, const char *prompt)
/* the same time. */
/* ******************************************************************************** */
static void
annihilate_1side(Varray2D<double> &M, long i, long j, long n)
annihilate_1side(Varray2D<double> &M, size_t i, size_t j, size_t n)
{
i--;
j--;
......@@ -1072,7 +1072,7 @@ annihilate_1side(Varray2D<double> &M, long i, long j, long n)
#ifdef HAVE_OPENMP4
#pragma omp simd reduction(+ : alpha) reduction(+ : beta) reduction(+ : gamma)
#endif
for (long r = 0; r < n; r++)
for (size_t r = 0; r < n; r++)
{
alpha += Mj[r] * Mj[r];
beta += Mi[r] * Mi[r];
......@@ -1100,7 +1100,7 @@ annihilate_1side(Varray2D<double> &M, long i, long j, long n)
const auto sk = ck * tk; // = sin(theta)
// calculate a_i,j - tilde
for (long r = 0; r < n; r++)
for (size_t r = 0; r < n; r++)
{
const auto mi = Mi[r];
const auto mj = Mj[r];
......@@ -1110,20 +1110,20 @@ annihilate_1side(Varray2D<double> &M, long i, long j, long n)
}
static int
jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
jacobi_1side(Varray2D<double> &M, Varray<double> &A, size_t n)
{
Varray2D<int> annihilations(2);
annihilations[0].resize(n * n);
annihilations[1].resize(n * n);
long count = 0;
for (long k = 1; k < n + 1; k++)
size_t count = 0;
for (size_t k = 1; k < n + 1; k++)
{
if (k < n)
{
{
const auto nmax = (long) std::ceil(1. / 2. * (n - k));
for (long i = 1; i <= nmax; i++)
const auto nmax = (size_t) std::ceil(1. / 2. * (n - k));
for (size_t i = 1; i <= nmax; i++)
{
const auto j = n - k + 2 - i;
annihilations[0][count] = i;
......@@ -1133,8 +1133,8 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
}
if (k > 2)
{
const auto nmax = n - (long) std::floor(1. / 2. * k);
for (long i = n - k + 2; i <= nmax; i++)
const auto nmax = n - (size_t) std::floor(1. / 2. * k);
for (size_t i = n - k + 2; i <= nmax; i++)
{
const auto j = 2 * n - k + 2 - i;
annihilations[0][count] = i;
......@@ -1145,8 +1145,8 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
}
else if (k == n)
{
const auto nmax = (long) std::ceil(1. / 2. * n);
for (long i = 2; i <= nmax; i++)
const auto nmax = (size_t) std::ceil(1. / 2. * n);
for (size_t i = 2; i <= nmax; i++)
{
const auto j = n + 2 - i;
annihilations[0][count] = i;
......@@ -1169,12 +1169,12 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
n_finished = 0;
if (n % 2 == 1)
{
for (long m = 0; m < n; m++)
for (size_t m = 0; m < n; m++)
{
#ifdef _OPENMP
#pragma omp parallel for shared(M, annihilations, n)
#endif
for (long i = 0; i < n / 2; i++)
for (size_t i = 0; i < n / 2; i++)
{
const auto idx = m * (n / 2) + i;
const auto i_ann = annihilations[0][idx];
......@@ -1185,12 +1185,12 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
}
else
{ // n%2 == 0
for (long m = 0; m < n; m++)
for (size_t m = 0; m < n; m++)
{
#ifdef _OPENMP
#pragma omp parallel for shared(M, annihilations, n)
#endif
for (long i = 0; i < n / 2 - (m % 2); i++)
for (size_t i = 0; i < n / 2 - (m % 2); i++)
{
auto idx = m / 2 * (n / 2 + n / 2 - 1);
if (m % 2) idx += n / 2;
......@@ -1219,7 +1219,7 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
{
// Do not overwrite results in case of insufficient convergence
cdoWarning("Setting Matrix and Eigenvalues to 0 before return");
for (long i = 0; i < n; i++)
for (size_t i = 0; i < n; i++)
{
memset(M[i].data(), 0, n * sizeof(double));
memset(A.data(), 0, n * sizeof(double));
......@@ -1228,12 +1228,12 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
}
}
// calculate eigen values as std::sqrt(||m_i||)
for (long i = 0; i < n; i++)
for (size_t i = 0; i < n; i++)
{
A[i] = 0.0;
for (long r = 0; r < n; r++) A[i] += M[i][r] * M[i][r];
for (size_t r = 0; r < n; r++) A[i] += M[i][r] * M[i][r];
A[i] = std::sqrt(A[i]);
for (long r = 0; r < n; r++) M[i][r] /= A[i];
for (size_t r = 0; r < n; r++) M[i][r] /= A[i];
}
heap_sort(A, M, n);
......@@ -1249,7 +1249,7 @@ jacobi_1side(Varray2D<double> &M, Varray<double> &A, long n)
/* ******************************************************************************** */
void
parallel_eigen_solution_of_symmetric_matrix(Varray2D<double> &M, Varray<double> &A, long n, const char func[])
parallel_eigen_solution_of_symmetric_matrix(Varray2D<double> &M, Varray<double> &A, size_t n, const char func[])
{
(void) (func); // CDO_UNUSED
......
......@@ -22,9 +22,9 @@
namespace cdo
{
void eigen_solution_of_symmetric_matrix(Varray2D<double> &a, Varray<double> &eig_val, long n, const char *prompt);
void fft(double *real, double *imag, long n, int sign);
void ft_r(double *real, double *imag, long n, int sign, double *work_r, double *work_i);
void eigen_solution_of_symmetric_matrix(Varray2D<double> &a, Varray<double> &eig_val, size_t 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);
double beta(double a, double b, const char *prompt);
double incomplete_gamma(double a, double x, const char *prompt);
......@@ -46,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(Varray2D<double> &M, Varray<double> &A, long n, const char func[]);
void parallel_eigen_solution_of_symmetric_matrix(Varray2D<double> &M, Varray<double> &A, size_t 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