blas_level_3.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>


void G_math_d_aA_B(double **A, double **B, double a, double **C, int rows,
                   int cols)
{
    int i, j;


    /*If B is null, scale the matrix A with th scalar a */
    if (B == NULL) {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = a * A[i][j];

        return;
    }

    /*select special cases */
    if (a == 1.0) {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = A[i][j] + B[i][j];
    }
    else if (a == -1.0) {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = B[i][j] - A[i][j];
    }
    else {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = a * A[i][j] + B[i][j];
    }

    return;
}

void G_math_f_aA_B(float **A, float **B, float a, float **C, int rows,
                   int cols)
{
    int i, j;

    /*If B is null, scale the matrix A with th scalar a */
    if (B == NULL) {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = a * A[i][j];
        return;
    }

    /*select special cases */
    if (a == 1.0) {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = A[i][j] + B[i][j];
    }
    else if (a == -1.0) {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = B[i][j] - A[i][j];
    }
    else {
#pragma omp for schedule (static) private(i, j)
        for (i = rows - 1; i >= 0; i--)
            for (j = cols - 1; j >= 0; j--)
                C[i][j] = a * A[i][j] + B[i][j];
    }

    return;
}


void G_math_d_AB(double **A, double **B, double **C, int rows_A,
                 int cols_A, int rows_B)
{
    int i, j, k;

#pragma omp for schedule (static) private(i, j, k)
    for (i = 0; i < rows_A; i++) {
        for (j = 0; j < rows_B; j++) {
            C[i][j] = 0.0;
            for (k = cols_A - 1; k >= 0; k--) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return;
}

void G_math_f_AB(float **A, float **B, float **C, int rows_A,
                 int cols_A, int rows_B)
{
    int i, j, k;

#pragma omp for schedule (static) private(i, j, k)
    for (i = 0; i < rows_A; i++) {
        for (j = 0; j < rows_B; j++) {
            C[i][j] = 0.0;
            for (k = cols_A - 1; k >= 0; k--) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return;
}