GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72100
symmetric_band_matrix.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <grass/gis.h>
5 #include <grass/gmath.h>
6 #include <grass/glocale.h>
7
8 /**
9  * \brief Convert a symmetrix matrix into a symmetric band matrix
10  *
11  * \verbatim
12  Symmetric matrix with bandwidth of 3
13
14  5 2 1 0
15  2 5 2 1
16  1 2 5 2
17  0 1 2 5
18
19  will be converted into the symmetric band matrix
20
21  5 2 1
22  5 2 1
23  5 2 0
24  5 0 0
25
26  \endverbatim
27
28  * \param A (double**) the symmetric matrix
29  * \param rows (int)
30  * \param bandwidth (int)
31  * \return B (double**) new created symmetric band matrix
32  * */
33
34 double **G_math_matrix_to_sband_matrix(double **A, int rows, int bandwidth)
35 {
36  int i, j;
37  double **B = G_alloc_matrix(rows, bandwidth);
38  double tmp;
39
40  for(i = 0; i < rows; i++) {
41  for(j = 0; j < bandwidth; j++) {
42  if(i + j < rows) {
43  tmp = A[i][i + j];
44  B[i][j] = tmp;
45  } else {
46  B[i][j] = 0.0;
47  }
48  }
49  }
50
51  return B;
52 }
53
54
55 /**
56  * \brief Convert a symmetric band matrix into a symmetric matrix
57  *
58  * \verbatim
59  Such a symmetric band matrix with banwidth 3
60
61  5 2 1
62  5 2 1
63  5 2 0
64  5 0 0
65
66  Will be converted into this symmetric matrix
67
68  5 2 1 0
69  2 5 2 1
70  1 2 5 2
71  0 1 2 5
72
73  \endverbatim
74  * \param A (double**) the symmetric band matrix
75  * \param rows (int)
76  * \param bandwidth (int)
77  * \return B (double**) new created symmetric matrix
78  * */
79
80 double **G_math_sband_matrix_to_matrix(double **A, int rows, int bandwidth)
81 {
82  int i, j;
83  double **B = G_alloc_matrix(rows, rows);
84  double tmp;
85
86  for(i = 0; i < rows; i++) {
87  for(j = 0; j < bandwidth; j++) {
88  tmp = A[i][j];
89  if(i + j < rows) {
90  B[i][i + j] = tmp;
91  }
92  }
93  }
94
95  /*Symmetry*/
96  for(i = 0; i < rows; i++) {
97  for(j = i; j < rows; j++) {
98  B[j][i] = B[i][j];
99  }
100  }
101
102  return B;
103 }
104
105
106 /*!
107  * \brief Compute the matrix - vector product
108  * of symmetric band matrix A and vector x.
109  *
110  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
111  *
112  * y = A * x
113  *
114  *
115  * \param A (double **)
116  * \param x (double) *)
117  * \param y (double * )
118  * \param rows (int)
119  * \param bandwidth (int)
120  * \return (void)
121  *
122  * */
123 void G_math_Ax_sband(double ** A, double *x, double *y, int rows, int bandwidth)
124 {
125  int i, j;
126  double tmp;
127
128
129 #pragma omp for schedule (static) private(i, j, tmp)
130  for (i = 0; i < rows; i++) {
131  tmp = 0;
132  for (j = 0; j < bandwidth; j++) {
133  if((i + j) < rows)
134  tmp += A[i][j]*x[i + j];
135  }
136  y[i] = tmp;
137  }
138
139 #pragma omp single
140  {
141  for (i = 0; i < rows; i++) {
142  tmp = 0;
143  for (j = 1; j < bandwidth; j++) {
144  if(i + j < rows)
145  y[i + j] += A[i][j]*x[i];
146  }
147  }
148  }
149  return;
150 }
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
double ** G_math_matrix_to_sband_matrix(double **A, int rows, int bandwidth)
Convert a symmetrix matrix into a symmetric band matrix.
#define x
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.
double ** G_math_sband_matrix_to_matrix(double **A, int rows, int bandwidth)
Convert a symmetric band matrix into a symmetric matrix.