21#include <grass/gmath.h>
40 "Add sparse vector %p to the sparse linear equation system at "
42 (
void *)spvector, row);
61 G_math_spvector **spmatrix;
63 G_debug(4,
"Allocate memory for a sparse matrix with %i rows\n", rows);
65 spmatrix = (G_math_spvector **)G_calloc(rows,
sizeof(G_math_spvector *));
79 G_math_spvector *spvector;
81 G_debug(4,
"Allocate memory for a sparse vector with %i cols\n", cols);
83 spvector = (G_math_spvector *)G_calloc(1,
sizeof(G_math_spvector));
85 spvector->cols = cols;
86 spvector->index = (
unsigned int *)G_calloc(cols,
sizeof(
unsigned int));
87 spvector->values = (
double *)G_calloc(cols,
sizeof(
double));
102 if (spvector->values)
127 for (i = 0; i < rows; i++)
152 for (i = 0; i < rows; i++) {
153 for (j = 0; j < rows; j++) {
155 for (k = 0; k < Asp[i]->cols; k++) {
156 if (Asp[i]->index[k] == (
unsigned int)j) {
157 fprintf(stdout,
"%4.5f ", Asp[i]->values[k]);
162 fprintf(stdout,
"%4.5f ", 0.0);
164 fprintf(stdout,
"\n");
190#pragma omp parallel for schedule(static) private(i, j)
191 for (i = 0; i < rows; i++) {
192 for (j = 0; j < Asp[i]->cols; j++) {
193 A[i][Asp[i]->index[j]] = Asp[i]->values[j];
231 assert(rows >= 0 && bandwidth >= 0);
235 for (i = 0; i < (
unsigned int)rows; i++) {
236 for (j = 0; j < Asp[i]->cols; j++) {
237 if (Asp[i]->index[j] == i) {
238 A[i][0] = Asp[i]->values[j];
240 else if (Asp[i]->index[j] > i) {
241 A[i][Asp[i]->index[j] - i] = Asp[i]->values[j];
264 int nonull,
count = 0;
266 G_math_spvector **Asp =
NULL;
270#pragma omp parallel for schedule(static) private(i, j, nonull, count)
271 for (i = 0; i < rows; i++) {
274 for (j = 0; j < rows; j++) {
275 if (A[i][j] > epsilon)
282 for (j = 0; j < rows; j++) {
283 if (A[i][j] > epsilon) {
285 v->values[
count] = A[i][j];
310 int bandwidth,
double epsilon)
314 int nonull,
count = 0;
316 G_math_spvector **Asp =
NULL;
320 for (i = 0; i < rows; i++) {
323 for (j = 0; j < bandwidth; j++) {
324 if (A[i][j] > epsilon)
333 if (A[i][0] > epsilon) {
335 v->values[
count] = A[i][0];
339 for (j = 1; j < bandwidth; j++) {
340 if (A[i][j] > epsilon && i + j < rows) {
341 v->index[
count] = i + j;
342 v->values[
count] = A[i][j];
376#pragma omp for schedule(static) private(i, j, tmp)
377 for (i = 0; i < rows; i++) {
379 for (j = 0; j < Asp[i]->cols; j++) {
380 tmp += Asp[i]->values[j] *
x[Asp[i]->index[j]];
void G_free(void *buf)
Free allocated memory.
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
int G_debug(int level, const char *msg,...)
Print debugging message.
#define assert(condition)
void G_math_Ax_sparse(G_math_spvector **Asp, double *x, double *y, int rows)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
double ** G_math_Asp_to_sband_matrix(G_math_spvector **Asp, int rows, int bandwidth)
Convert a symmetric sparse matrix into a symmetric band matrix.
void G_math_print_spmatrix(G_math_spvector **Asp, int rows)
print the sparse matrix Asp to stdout
double ** G_math_Asp_to_A(G_math_spvector **Asp, int rows)
Convert a sparse matrix into a quadratic matrix.
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
void G_math_free_spvector(G_math_spvector *spvector)
Release the memory of the sparse vector.
void G_math_free_spmatrix(G_math_spvector **Asp, int rows)
Release the memory of the sparse matrix.
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse matrix at position row.
G_math_spvector ** G_math_sband_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)
Convert a symmetric band matrix into a sparse matrix.
G_math_spvector ** G_math_A_to_Asp(double **A, int rows, double epsilon)
Convert a quadratic matrix into a sparse matrix.
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.