Loading...
Searching...
No Matches
symmetric_band_matrix.c
Go to the documentation of this file.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/gmath.h>
#include <grass/glocale.h>
/**
* \brief Convert a symmetrix matrix into a symmetric band matrix
*
* \verbatim
Symmetric matrix with bandwidth of 3
5 2 1 0
2 5 2 1
1 2 5 2
0 1 2 5
will be converted into the symmetric band matrix
5 2 1
5 2 1
5 2 0
5 0 0
\endverbatim
* \param A (double**) the symmetric matrix
* \param rows (int)
* \param bandwidth (int)
* \return B (double**) new created symmetric band matrix
* */
double **G_math_matrix_to_sband_matrix(double **A, int rows, int bandwidth)
{
int i, j;
double **B = G_alloc_matrix(rows, bandwidth);
double tmp;
for (i = 0; i < rows; i++) {
for (j = 0; j < bandwidth; j++) {
if (i + j < rows) {
tmp = A[i][i + j];
B[i][j] = tmp;
}
else {
B[i][j] = 0.0;
}
}
}
return B;
}
/**
* \brief Convert a symmetric band matrix into a symmetric matrix
*
* \verbatim
Such a symmetric band matrix with banwidth 3
5 2 1
5 2 1
5 2 0
5 0 0
Will be converted into this symmetric matrix
5 2 1 0
2 5 2 1
1 2 5 2
0 1 2 5
\endverbatim
* \param A (double**) the symmetric band matrix
* \param rows (int)
* \param bandwidth (int)
* \return B (double**) new created symmetric matrix
* */
double **G_math_sband_matrix_to_matrix(double **A, int rows, int bandwidth)
{
int i, j;
double **B = G_alloc_matrix(rows, rows);
double tmp;
for (i = 0; i < rows; i++) {
for (j = 0; j < bandwidth; j++) {
tmp = A[i][j];
if (i + j < rows) {
B[i][i + j] = tmp;
}
}
}
/*Symmetry */
for (i = 0; i < rows; i++) {
for (j = i; j < rows; j++) {
B[j][i] = B[i][j];
}
}
return B;
}
/*!
* \brief Compute the matrix - vector product
* of symmetric band 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 bandwidth (int)
* \return (void)
*
* */
void G_math_Ax_sband(double **A, double *x, double *y, int rows, int bandwidth)
{
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 < bandwidth; j++) {
if ((i + j) < rows)
tmp += A[i][j] * x[i + j];
}
y[i] = tmp;
}
#pragma omp single
{
for (i = 0; i < rows; i++) {
tmp = 0;
for (j = 1; j < bandwidth; j++) {
if (i + j < rows)
y[i + j] += A[i][j] * x[i];
}
}
}
return;
}
G_alloc_matrix
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:55
G_math_Ax_sband
void G_math_Ax_sband(double **A, double *x, double *y, int rows, int bandwidth)
Compute the matrix - vector product of symmetric band matrix A and vector x.
Definition symmetric_band_matrix.c:123
G_math_matrix_to_sband_matrix
double ** G_math_matrix_to_sband_matrix(double **A, int rows, int bandwidth)
Convert a symmetrix matrix into a symmetric band matrix.
Definition symmetric_band_matrix.c:34
G_math_sband_matrix_to_matrix
double ** G_math_sband_matrix_to_matrix(double **A, int rows, int bandwidth)
Convert a symmetric band matrix into a symmetric matrix.
Definition symmetric_band_matrix.c:80
x
#define x