Loading...
Searching...
No Matches
blas_level_2.c
Go to the documentation of this file.
/*****************************************************************************
*
* MODULE: Grass numerical math interface
* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
* soerengebbert <at> googlemail <dot> com
*
* PURPOSE: grass blas implementation
* part of the gmath library
*
* COPYRIGHT: (C) 2010 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>
#define EPSILON 0.00000000000000001
/*!
* \brief Compute the matrix - vector product
* of matrix A and vector x.
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
* y = A * x
*
*
* \param A (double ** )
* \param x (double *)
* \param y (double *)
* \param rows (int)
* \param cols (int)
* \return (void)
*
* */
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;
}
/*!
* \brief Compute the matrix - vector product
* of matrix A and vector x.
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
* y = A * x
*
*
* \param A (float ** )
* \param x (float *)
* \param y (float *)
* \param rows (int)
* \param cols (int)
* \return (void)
*
* */
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;
}
/*!
* \brief Compute the dyadic product of two vectors.
* The result is stored in the matrix A.
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
* A = x * y^T
*
*
* \param x (double *)
* \param y (double *)
* \param A (float **) -- matrix of size rows*cols
* \param rows (int) -- length of vector x
* \param cols (int) -- length of vector y
* \return (void)
*
* */
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;
}
/*!
* \brief Compute the dyadic product of two vectors.
* The result is stored in the matrix A.
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
* A = x * y^T
*
*
* \param x (float *)
* \param y (float *)
* \param A (float **= -- matrix of size rows*cols
* \param rows (int) -- length of vector x
* \param cols (int) -- length of vector y
* \return (void)
*
* */
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;
}
/*!
* \brief Compute the scaled matrix - vector product
* of matrix double **A and vector x and y.
*
* z = a * A * x + b * y
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
*
* \param A (double **)
* \param x (double *)
* \param y (double *)
* \param a (double)
* \param b (double)
* \param z (double *)
* \param rows (int)
* \param cols (int)
* \return (void)
*
* */
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;
}
/*!
* \brief Compute the scaled matrix - vector product
* of matrix A and vectors x and y.
*
* z = a * A * x + b * y
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
*
* \param A (float **)
* \param x (float *)
* \param y (float *)
* \param a (float)
* \param b (float)
* \param z (float *)
* \param rows (int)
* \param cols (int)
* \return (void)
*
* */
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;
}
/*!
* \fn int G_math_d_A_T(double **A, int rows)
*
* \brief Compute the transposition of matrix A.
* Matrix A will be overwritten.
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
* Returns 0.
*
* \param A (double **)
* \param rows (int)
* \return int
*/
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;
}
/*!
* \fn int G_math_f_A_T(float **A, int rows)
*
* \brief Compute the transposition of matrix A.
* Matrix A will be overwritten.
*
* This function is multi-threaded with OpenMP and can be called within a
* parallel OpenMP region.
*
* Returns 0.
*
* \param A (float **)
* \param rows (int)
* \return int
*/
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;
}
G_math_d_A_T
int G_math_d_A_T(double **A, int rows)
Compute the transposition of matrix A. Matrix A will be overwritten.
Definition blas_level_2.c:341
G_math_d_aAx_by
void G_math_d_aAx_by(double **A, double *x, double *y, double a, double b, double *z, int rows, int cols)
Compute the scaled matrix - vector product of matrix double **A and vector x and y.
Definition blas_level_2.c:182
G_math_d_Ax
void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
Definition blas_level_2.c:46
G_math_f_aAx_by
void G_math_f_aAx_by(float **A, float *x, float *y, float a, float b, float *z, int rows, int cols)
Compute the scaled matrix - vector product of matrix A and vectors x and y.
Definition blas_level_2.c:265
G_math_f_x_dyad_y
void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
Compute the dyadic product of two vectors. The result is stored in the matrix A.
Definition blas_level_2.c:147
G_math_f_Ax
void G_math_f_Ax(float **A, float *x, float *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
Definition blas_level_2.c:81
G_math_d_x_dyad_y
void G_math_d_x_dyad_y(double *x, double *y, double **A, int rows, int cols)
Compute the dyadic product of two vectors. The result is stored in the matrix A.
Definition blas_level_2.c:116
G_math_f_A_T
int G_math_f_A_T(float **A, int rows)
Compute the transposition of matrix A. Matrix A will be overwritten.
Definition blas_level_2.c:374
b
double b
Definition driver/set_window.c:5
x
#define x