blas_level_2.c

Go to the documentation of this file.

/*****************************************************************************
 *
 * MODULE:       Grass PDE Numerical Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2007
 *              soerengebbert <at> gmx <dot> de
 *               
 * PURPOSE:      linear equation system solvers
 *              part of the gpde library
 *               
 * COPYRIGHT:    (C) 2007 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 <math.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <grass/gmath.h>
#include <grass/gis.h>
#include <grass/gisdefs.h>

#define EPSILON 0.00000000000000001

void G_math_Ax_sparse(G_math_spvector ** Asp, double *x, double *y, int rows)
{
    int i, j;

    double tmp;

#pragma omp for schedule (static) private(i, j, tmp)
    for (i = 0; i < rows; i++) {
        tmp = 0;
        for (j = 0; j < Asp[i]->cols; j++) {
            tmp += Asp[i]->values[j] * x[Asp[i]->index[j]];
        }
        y[i] = tmp;
    }
    return;
}

void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
{
    int i, j;

    double tmp;

#pragma omp for schedule (static) private(i, j, tmp)
    for (i = 0; i < rows; i++) {
        tmp = 0;
        for (j = cols - 1; j >= 0; j--) {
            tmp += A[i][j] * x[j];
        }
        y[i] = tmp;
    }
    return;
}

void G_math_f_Ax(float **A, float *x, float *y, int rows, int cols)
{
    int i, j;

    float tmp;

#pragma omp for schedule (static) private(i, j, tmp)
    for (i = 0; i < rows; i++) {
        tmp = 0;
        for (j = cols - 1; j >= 0; j--) {
            tmp += A[i][j] * x[j];
        }
        y[i] = tmp;
    }
    return;
}

void G_math_d_x_dyad_y(double *x, double *y, double **A, int rows, int cols)
{
    int i, j;

#pragma omp for schedule (static) private(i, j)
    for (i = 0; i < rows; i++) {
        for (j = cols - 1; j >= 0; j--) {
            A[i][j] = x[i] * y[j];
        }
    }
    return;
}

void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
{
    int i, j;

#pragma omp for schedule (static) private(i, j)
    for (i = 0; i < rows; i++) {
        for (j = cols - 1; j >= 0; j--) {
            A[i][j] = x[i] * y[j];
        }
    }
    return;
}

void G_math_d_aAx_by(double **A, double *x, double *y, double a, double b,
                     double *z, int rows, int cols)
{
    int i, j;

    double tmp;

    /*catch specific cases */
    if (a == b) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += A[i][j] * x[j] + y[j];
            }
            z[i] = a * tmp;
        }
    }
    else if (b == -1.0) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += a * A[i][j] * x[j] - y[j];
            }
            z[i] = tmp;
        }
    }
    else if (b == 0.0) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += A[i][j] * x[j];
            }
            z[i] = a * tmp;
        }
    }
    else if (a == -1.0) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += b * y[j] - A[i][j] * x[j];
            }
            z[i] = tmp;
        }
    }
    else {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += a * A[i][j] * x[j] + b * y[j];
            }
            z[i] = tmp;
        }
    }
    return;
}

void G_math_f_aAx_by(float **A, float *x, float *y, float a, float b,
                     float *z, int rows, int cols)
{
    int i, j;

    float tmp;

    /*catch specific cases */
    if (a == b) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += A[i][j] * x[j] + y[j];
            }
            z[i] = a * tmp;
        }
    }
    else if (b == -1.0) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += a * A[i][j] * x[j] - y[j];
            }
            z[i] = tmp;
        }
    }
    else if (b == 0.0) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += A[i][j] * x[j];
            }
            z[i] = a * tmp;
        }
    }
    else if (a == -1.0) {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += b * y[j] - A[i][j] * x[j];
            }
            z[i] = tmp;
        }
    }
    else {
#pragma omp for schedule (static) private(i, j, tmp)
        for (i = 0; i < rows; i++) {
            tmp = 0;
            for (j = cols - 1; j >= 0; j--) {
                tmp += a * A[i][j] * x[j] + b * y[j];
            }
            z[i] = tmp;
        }
    }
    return;
}



int G_math_d_A_T(double **A, int rows)
{
    int i, j;

    double tmp;

#pragma omp for schedule (static) private(i, j, tmp)
    for (i = 0; i < rows; i++)
        for (j = 0; j < i; j++) {
            tmp = A[i][j];

            A[i][j] = A[j][i];
            A[j][i] = tmp;
        }

    return 0;
}

int G_math_f_A_T(float **A, int rows)
{
    int i, j;

    float tmp;

#pragma omp for schedule (static) private(i, j, tmp)
    for (i = 0; i < rows; i++)
        for (j = 0; j < i; j++) {
            tmp = A[i][j];

            A[i][j] = A[j][i];
            A[j][i] = tmp;
        }

    return 0;
}