test_blas1.c

/*****************************************************************************
 *
 * MODULE:       Grass PDE Numerical Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 *                 soerengebbert <at> gmx <dot> de
 *
 * PURPOSE:      Unit tests for les creation
 *
 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 *****************************************************************************/

#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/gmath.h>
#include <math.h>
#include "test_gmath_lib.h"

#define EPSILON 0.000001

/* prototypes */
static int test_blas_level_1_double(void);
static int test_blas_level_1_float(void);
static int test_blas_level_1_int(void);

/* *************************************************************** */
/* Perfrome the blas level 1 unit tests ************************** */
/* *************************************************************** */
int unit_test_blas_level_1(void)
{
    int sum = 0;

    G_message(_("\n++ Running blas level 1 unit tests ++"));

    sum += test_blas_level_1_double();
    sum += test_blas_level_1_float();
    sum += test_blas_level_1_int();

    if (sum > 0)
        G_warning(_("\n-- blas level 1 unit tests failure --"));
    else
        G_message(_("\n-- blas level 1 unit tests finished successfully --"));

    return sum;
}

/* *************************************************************** */
/* ************** D O U B L E ************************************ */
/* *************************************************************** */
int test_blas_level_1_double(void)
{

    int sum = 0;
    int rows = 10000;
    double *x, *y, *z, a = 0.0, b = 0.0, c = 0.0, d = 0.0, e = 0.0;

    x = G_alloc_vector(rows);
    y = G_alloc_vector(rows);
    z = G_alloc_vector(rows);

    fill_d_vector_scalar(x, 1, rows);
    fill_d_vector_scalar(y, 2, rows);

    /*test the grass implementation */
    G_math_d_x_dot_y(x, y, &a, rows);
    G_math_d_asum_norm(x, &b, rows);
    G_math_d_euclid_norm(x, &c, rows);

    if (a != 2.0 * rows) {
        G_message("Error in G_math_d_x_dot_y %f != %f", 2.0 * rows, a);
        sum++;
    }

    if (b != rows) {
        G_message("Error in G_math_d_asum_norm");
        sum++;
    }

    if (c != sqrt(rows)) {
        G_message("Error in G_math_d_euclid_norm");
        sum++;
    }

    /*test the ATLAS implementation */
    a = G_math_dnrm2(x, rows);
    b = G_math_dasum(x, rows);
    c = G_math_ddot(x, y, rows);

    if (a != sqrt(rows)) {
        G_message("Error in G_math_dnrm2 %f != %f", sqrt(rows), a);
        sum++;
    }

    if (b != rows) {
        G_message("Error in G_math_dasum %f != %f", 2.0 * rows, b);
        sum++;
    }

    if (c != 2.0 * rows) {
        G_message("Error in G_math_ddot %f != %f", 2.0 * rows, c);
        sum++;
    }

    fill_d_vector_range_1(x, 1.0, rows);
    fill_d_vector_range_2(y, 1.0, rows);

    /*grass function */
    G_math_d_max_norm(x, &a, rows);
    /*atlas function */
    b = G_math_idamax(x, rows);

    if (a != 1.0 * (rows - 1)) {
        G_message("Error in G_math_d_max_norm: %f != %f",
                  (double)1.0 * (rows - 1), a);
        sum++;
    }

    if (b != 1.0 * (rows - 1)) {
        G_message("Error in G_math_idamax: %f != %f", (double)1.0 * (rows - 1),
                  b);
        sum++;
    }

#pragma omp parallel default(shared)
    {
        G_math_d_ax_by(x, y, z, 1.0, 1.0, rows);
    }
    G_math_d_asum_norm(z, &a, rows);

#pragma omp parallel default(shared)
    {
        G_math_d_ax_by(x, y, z, 1.0, -1.0, rows);
    }
    G_math_d_asum_norm(z, &b, rows);

#pragma omp parallel default(shared)
    {
        G_math_d_ax_by(x, y, z, 2.0, 1.0, rows);
    }
    G_math_d_asum_norm(z, &c, rows);

    if (a != 1.0 * (rows - 1) * rows) {
        G_message("Error in G_math_d_ax_by: %f != %f",
                  (double)1.0 * (rows - 1) * rows, a);
        sum++;
    }

    if (b != 5.0 * (rows) * (rows / 10)) {
        G_message("Error in G_math_d_ax_by: %f != %f",
                  (double)5.0 * (rows) * (rows / 10), b);
        sum++;
    }

    if (c != 149985000) {
        G_message("Error in G_math_d_ax_by: 149985000 != %f", c);
        sum++;
    }

#pragma omp parallel default(shared)
    {
        /*scale x with 1 */
        G_math_d_ax_by(x, z, z, 1.0, 0.0, rows);
    }
    G_math_d_asum_norm(x, &a, rows);
    G_math_d_asum_norm(z, &b, rows);
    /*scale x with -1 */
#pragma omp parallel default(shared)
    {
        G_math_d_ax_by(x, z, z, -1.0, 0.0, rows);
    }
    G_math_d_asum_norm(z, &c, rows);

    /*ATLAS implementation */
    G_math_dscal(x, 1.0, rows);
    G_math_d_asum_norm(x, &d, rows);

    /*ATLAS implementation */
    fill_d_vector_range_1(x, 1.0, rows);
    fill_d_vector_scalar(z, 0.0, rows);
    G_math_daxpy(x, z, 1.0, rows);
    G_math_d_asum_norm(z, &e, rows);

    if (a != 49995000 || a != b || b != c) {
        G_message("Error in G_math_d_ax: 49995000 != %f", a);
        sum++;
    }

    if (49995000 != d) {
        G_message("Error in G_math_dscal: 49995000 != %f", d);
        sum++;
    }

    if (49995000 != e) {
        G_message("Error in G_math_daxpy: 49995000 != %f", e);
        sum++;
    }

    fill_d_vector_scalar(z, 0, rows);

    G_math_d_copy(x, z, rows);
    G_math_d_asum_norm(x, &a, rows);
    G_math_dcopy(x, z, rows);
    G_math_d_asum_norm(x, &b, rows);

    if (a != 49995000) {
        G_message("Error in G_math_d_copy: 49995000 != %f", a);
        sum++;
    }

    if (b != 49995000) {
        G_message("Error in G_math_dcopy: 49995000 != %f", a);
        sum++;
    }

    G_free_vector(x);
    G_free_vector(y);
    G_free_vector(z);

    return sum;
}

/* *************************************************************** */
/* ************** F L O A T ************************************** */
/* *************************************************************** */
int test_blas_level_1_float(void)
{

    int sum = 0;
    int rows = 1000;
    float *x, *y, *z, a = 0.0, b = 0.0, c = 0.0, d = 0.0, e = 0.0;

    x = G_alloc_fvector(rows);
    y = G_alloc_fvector(rows);
    z = G_alloc_fvector(rows);

    fill_f_vector_scalar(x, 1, rows);
    fill_f_vector_scalar(y, 2, rows);

    /*test the grass implementation */
    G_math_f_x_dot_y(x, y, &a, rows);
    G_math_f_asum_norm(x, &b, rows);
    G_math_f_euclid_norm(x, &c, rows);

    if (a != 2.0 * rows) {
        G_message("Error in G_math_f_x_dot_y %f != %f", 2.0 * rows, a);
        sum++;
    }

    if (b != rows) {
        G_message("Error in G_math_f_asum_norm");
        sum++;
    }

    if (fabs(c - (float)sqrt(rows)) > EPSILON) {
        G_message("Error in G_math_f_euclid_norm");
        sum++;
    }

    /*test the ATLAS implementation */
    a = G_math_snrm2(x, rows);
    b = G_math_sasum(x, rows);
    c = G_math_sdot(x, y, rows);

    if (fabs(a - sqrt(rows)) > EPSILON) {
        G_message("Error in G_math_snrm2 %f != %f", sqrt(rows), a);
        sum++;
    }

    if (b != rows) {
        G_message("Error in G_math_sasum %f != %f", 2.0 * rows, b);
        sum++;
    }

    if (c != 2.0 * rows) {
        G_message("Error in G_math_sdot %f != %f", 2.0 * rows, c);
        sum++;
    }

    fill_f_vector_range_1(x, 1.0, rows);
    fill_f_vector_range_2(y, 1.0, rows);

    /*grass function */
    G_math_f_max_norm(x, &a, rows);
    /*atlas function */
    b = G_math_isamax(x, rows);

    if (a != 1.0 * (rows - 1)) {
        G_message("Error in G_math_f_max_norm: %f != %f",
                  (float)1.0 * (rows - 1), a);
        sum++;
    }

    if (b != 1.0 * (rows - 1)) {
        G_message("Error in G_math_isamax: %f != %f", (float)1.0 * (rows - 1),
                  b);
        sum++;
    }

#pragma omp parallel default(shared)
    {
        G_math_f_ax_by(x, y, z, 1.0, 1.0, rows);
    }
    G_math_f_asum_norm(z, &a, rows);
#pragma omp parallel default(shared)
    {
        G_math_f_ax_by(x, y, z, 1.0, -1.0, rows);
    }
    G_math_f_asum_norm(z, &b, rows);
#pragma omp parallel default(shared)
    {
        G_math_f_ax_by(x, y, z, 2.0, 1.0, rows);
    }
    G_math_f_asum_norm(z, &c, rows);

    if (fabs(a - 1.0 * (rows - 1) * rows) > EPSILON) {
        G_message("Error in G_math_f_ax_by 1: %f != %f",
                  (float)1.0 * (rows - 1) * rows, a);
        sum++;
    }

    if (fabs(b - 5.0 * (rows) * (rows / 10)) > EPSILON) {
        G_message("Error in G_math_f_ax_by 2: %f != %f",
                  (float)5.0 * (rows) * (rows / 10), b);
        sum++;
    }

    if (fabs(c - 1498500) > EPSILON) {
        G_message("Error in G_math_f_ax_by 3: 14998500 != %f", c);
        sum++;
    }

#pragma omp parallel default(shared)
    {
        /*scale x with 1 */
        G_math_f_ax_by(x, z, z, 1.0, 0.0, rows);
    }
    G_math_f_asum_norm(x, &a, rows);
    G_math_f_asum_norm(z, &b, rows);
    /*scale x with -1 */
#pragma omp parallel default(shared)
    {
        G_math_f_ax_by(x, z, z, -1.0, 0.0, rows);
    }
    G_math_f_asum_norm(z, &c, rows);

    /*ATLAS implementation */
    G_math_sscal(x, 1.0, rows);
    G_math_f_asum_norm(x, &d, rows);

    /*ATLAS implementation */
    fill_f_vector_range_1(x, 1.0, rows);
    fill_f_vector_scalar(z, 0.0, rows);
    G_math_saxpy(x, z, 1.0, rows);
    G_math_f_asum_norm(z, &e, rows);

    if (fabs(a - 499500) > EPSILON) {
        G_message("Error in G_math_f_ax_by 4: 4999500 != %f", a);
        sum++;
    }
    if (fabs(b - 499500) > EPSILON) {
        G_message("Error in G_math_f_ax_by 4: 4999500 != %f", b);
        sum++;
    }
    if (fabs(c - 499500) > EPSILON) {
        G_message("Error in G_math_f_ax_by 4: 4999500 != %f", c);
        sum++;
    }
    if (fabs(d - 499500) > EPSILON) {
        G_message("Error in G_math_sscal: 4999500 != %f", d);
        sum++;
    }

    if (fabs(e - 499500) > EPSILON) {
        G_message("Error in G_math_saxpy: 4999500 != %f", e);
        sum++;
    }

    fill_f_vector_range_1(x, 1.0, rows);
    fill_f_vector_scalar(z, 0, rows);

    G_math_f_copy(x, z, rows);
    G_math_f_asum_norm(x, &a, rows);
    G_math_scopy(x, z, rows);
    G_math_f_asum_norm(x, &b, rows);

    if (fabs(a - 499500) > EPSILON) {
        G_message("Error in G_math_f_copy: 4999500 != %f", a);
        sum++;
    }

    if (fabs(b - 499500) > EPSILON) {
        G_message("Error in G_math_scopy: 4999500 != %f", b);
        sum++;
    }

    G_free_fvector(x);
    G_free_fvector(y);
    G_free_fvector(z);

    return sum;
}

/* *************************************************************** */
/* ************** I N T E G E R ********************************** */
/* *************************************************************** */
int test_blas_level_1_int(void)
{

    int sum = 0;
    int rows = 10000;
    int *x, *y, *z, max;
    double a, b, c;

    x = G_alloc_ivector(rows);
    y = G_alloc_ivector(rows);
    z = G_alloc_ivector(rows);

    fill_i_vector_scalar(x, 1, rows);
    fill_i_vector_scalar(y, 2, rows);

    G_math_i_x_dot_y(x, y, &a, rows);
    G_math_i_asum_norm(x, &b, rows);
    G_math_i_euclid_norm(x, &c, rows);

    if (a != 2 * rows) {
        G_message("Error in G_math_i_x_dot_y");
        sum++;
    }

    if (b != rows) {
        G_message("Error in G_math_i_asum_norm");
        sum++;
    }

    if (c != sqrt((double)rows)) {
        G_message("Error in G_math_i_euclid_norm");
        sum++;
    }

    fill_i_vector_range_1(x, 1, rows);
    fill_i_vector_range_2(y, 1, rows);

    G_math_i_max_norm(x, &max, rows);

    if (max != (rows - 1)) {
        G_message("Error in G_math_i_max_norm: %i != %i", (rows - 1), max);
        sum++;
    }

#pragma omp parallel default(shared)
    {
        G_math_i_ax_by(x, y, z, 1, 1, rows);
    }
    G_math_i_asum_norm(z, &a, rows);
#pragma omp parallel default(shared)
    {
        G_math_i_ax_by(x, y, z, 1, -1, rows);
    }
    G_math_i_asum_norm(z, &b, rows);
#pragma omp parallel default(shared)
    {
        G_math_i_ax_by(x, y, z, 2, 1, rows);
    }
    G_math_i_asum_norm(z, &c, rows);

    if (a != 1.0 * (rows - 1) * rows) {
        G_message("Error in G_math_i_ax_by: %f != %f", 1.0 * (rows - 1) * rows,
                  a);
        sum++;
    }

    if (b != 5.0 * (rows) * (rows / 10)) {
        G_message("Error in G_math_i_ax_by: %f != %f",
                  5.0 * (rows) * (rows / 10), b);
        sum++;
    }

    if (c != 149985000) {
        G_message("Error in G_math_i_ax_by: 149985000 != %f", c);
        sum++;
    }

#pragma omp parallel default(shared)
    {
        /*scale x with 1 */
        G_math_i_ax_by(x, z, z, 1, 0, rows);
    }
    G_math_i_asum_norm(x, &a, rows);
    G_math_i_asum_norm(z, &b, rows);

    /*scale a with -1 */
#pragma omp parallel default(shared)
    {
        G_math_i_ax_by(x, z, z, -1, 0, rows);
    }
    G_math_i_asum_norm(z, &c, rows);

    if (a != 49995000 || a != b || b != c) {
        G_message("Error in G_math_i_ax_by: 49995000 != %f", a);
        sum++;
    }

    return sum;
}