1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2
3 /*@
4 MatGetColumnVector - Gets the values from a given column of a matrix.
5
6 Not Collective
7
8 Input Parameters:
9 + A - the matrix
10 . yy - the vector
11 - col - the column requested (in global numbering)
12
13 Level: advanced
14
15 Notes:
16 If a `MatType` does not implement the operation, each processor for which this is called
17 gets the values for its rows using `MatGetRow()`.
18
19 The vector must have the same parallel row layout as the matrix.
20
21 Contributed by: Denis Vanderstraeten
22
23 .seealso: [](ch_matrices), `Mat`, `MatGetRow()`, `MatGetDiagonal()`, `MatMult()`
24 @*/
MatGetColumnVector(Mat A,Vec yy,PetscInt col)25 PetscErrorCode MatGetColumnVector(Mat A, Vec yy, PetscInt col)
26 {
27 PetscScalar *y;
28 const PetscScalar *v;
29 PetscInt i, j, nz, N, Rs, Re, rs, re;
30 const PetscInt *idx;
31
32 PetscFunctionBegin;
33 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
34 PetscValidHeaderSpecific(yy, VEC_CLASSID, 2);
35 PetscValidLogicalCollectiveInt(A, col, 3);
36 PetscCheck(col >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Requested negative column: %" PetscInt_FMT, col);
37 PetscCall(MatGetSize(A, NULL, &N));
38 PetscCheck(col < N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Requested column %" PetscInt_FMT " larger than number columns in matrix %" PetscInt_FMT, col, N);
39 PetscCall(MatGetOwnershipRange(A, &Rs, &Re));
40 PetscCall(VecGetOwnershipRange(yy, &rs, &re));
41 PetscCheck(Rs == rs && Re == re, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Matrix %" PetscInt_FMT " %" PetscInt_FMT " does not have same ownership range (size) as vector %" PetscInt_FMT " %" PetscInt_FMT, Rs, Re, rs, re);
42
43 if (A->ops->getcolumnvector) PetscUseTypeMethod(A, getcolumnvector, yy, col);
44 else {
45 PetscCall(VecSet(yy, 0.0));
46 PetscCall(VecGetArray(yy, &y));
47 /* TODO for general matrices */
48 for (i = Rs; i < Re; i++) {
49 PetscCall(MatGetRow(A, i, &nz, &idx, &v));
50 if (nz && idx[0] <= col) {
51 /*
52 Should use faster search here
53 */
54 for (j = 0; j < nz; j++) {
55 if (idx[j] >= col) {
56 if (idx[j] == col) y[i - rs] = v[j];
57 break;
58 }
59 }
60 }
61 PetscCall(MatRestoreRow(A, i, &nz, &idx, &v));
62 }
63 PetscCall(VecRestoreArray(yy, &y));
64 }
65 PetscFunctionReturn(PETSC_SUCCESS);
66 }
67
68 /*@
69 MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix.
70
71 Input Parameters:
72 + A - the matrix
73 - type - `NORM_2`, `NORM_1` or `NORM_INFINITY`
74
75 Output Parameter:
76 . norms - an array as large as the TOTAL number of columns in the matrix
77
78 Level: intermediate
79
80 Note:
81 Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values,
82 if each process wants only some of the values it should extract the ones it wants from the array.
83
84 .seealso: [](ch_matrices), `Mat`, `NormType`, `MatNorm()`
85 @*/
MatGetColumnNorms(Mat A,NormType type,PetscReal norms[])86 PetscErrorCode MatGetColumnNorms(Mat A, NormType type, PetscReal norms[])
87 {
88 /* NOTE: MatGetColumnNorms() could simply be a macro that calls MatGetColumnReductions().
89 * I've kept this as a function because it allows slightly more in the way of error checking,
90 * erroring out if MatGetColumnNorms() is not called with a valid NormType. */
91
92 PetscFunctionBegin;
93 PetscCheck(type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown NormType");
94 PetscCall(MatGetColumnReductions(A, type, norms));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
98 /*@
99 MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix.
100
101 Input Parameter:
102 . A - the matrix
103
104 Output Parameter:
105 . sums - an array as large as the TOTAL number of columns in the matrix
106
107 Level: intermediate
108
109 Note:
110 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
111 if each process wants only some of the values it should extract the ones it wants from the array.
112
113 .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
114 @*/
MatGetColumnSumsRealPart(Mat A,PetscReal sums[])115 PetscErrorCode MatGetColumnSumsRealPart(Mat A, PetscReal sums[])
116 {
117 PetscFunctionBegin;
118 PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_REALPART, sums));
119 PetscFunctionReturn(PETSC_SUCCESS);
120 }
121
122 /*@
123 MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix.
124
125 Input Parameter:
126 . A - the matrix
127
128 Output Parameter:
129 . sums - an array as large as the TOTAL number of columns in the matrix
130
131 Level: intermediate
132
133 Note:
134 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
135 if each process wants only some of the values it should extract the ones it wants from the array.
136
137 .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
138 @*/
MatGetColumnSumsImaginaryPart(Mat A,PetscReal sums[])139 PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A, PetscReal sums[])
140 {
141 PetscFunctionBegin;
142 PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_IMAGINARYPART, sums));
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
146 /*@
147 MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix.
148
149 Input Parameter:
150 . A - the matrix
151
152 Output Parameter:
153 . sums - an array as large as the TOTAL number of columns in the matrix
154
155 Level: intermediate
156
157 Note:
158 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
159 if each process wants only some of the values it should extract the ones it wants from the array.
160
161 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
162 @*/
MatGetColumnSums(Mat A,PetscScalar sums[])163 PetscErrorCode MatGetColumnSums(Mat A, PetscScalar sums[])
164 {
165 #if defined(PETSC_USE_COMPLEX)
166 PetscInt i, n;
167 PetscReal *work;
168 #endif
169
170 PetscFunctionBegin;
171 #if !defined(PETSC_USE_COMPLEX)
172 PetscCall(MatGetColumnSumsRealPart(A, sums));
173 #else
174 PetscCall(MatGetSize(A, NULL, &n));
175 PetscCall(PetscArrayzero(sums, n));
176 PetscCall(PetscCalloc1(n, &work));
177 PetscCall(MatGetColumnSumsRealPart(A, work));
178 for (i = 0; i < n; i++) sums[i] = work[i];
179 PetscCall(MatGetColumnSumsImaginaryPart(A, work));
180 for (i = 0; i < n; i++) sums[i] += work[i] * PETSC_i;
181 PetscCall(PetscFree(work));
182 #endif
183 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185
186 /*@
187 MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix.
188
189 Input Parameter:
190 . A - the matrix
191
192 Output Parameter:
193 . means - an array as large as the TOTAL number of columns in the matrix
194
195 Level: intermediate
196
197 Note:
198 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
199 if each process wants only some of the values it should extract the ones it wants from the array.
200
201 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
202 @*/
MatGetColumnMeansRealPart(Mat A,PetscReal means[])203 PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[])
204 {
205 PetscFunctionBegin;
206 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means));
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
210 /*@
211 MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix.
212
213 Input Parameter:
214 . A - the matrix
215
216 Output Parameter:
217 . means - an array as large as the TOTAL number of columns in the matrix
218
219 Level: intermediate
220
221 Note:
222 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
223 if each process wants only some of the values it should extract the ones it wants from the array.
224
225 .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
226 @*/
MatGetColumnMeansImaginaryPart(Mat A,PetscReal means[])227 PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[])
228 {
229 PetscFunctionBegin;
230 PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means));
231 PetscFunctionReturn(PETSC_SUCCESS);
232 }
233
234 /*@
235 MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix.
236
237 Input Parameter:
238 . A - the matrix
239
240 Output Parameter:
241 . means - an array as large as the TOTAL number of columns in the matrix
242
243 Level: intermediate
244
245 Note:
246 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
247 if each process wants only some of the values it should extract the ones it wants from the array.
248
249 .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
250 @*/
MatGetColumnMeans(Mat A,PetscScalar means[])251 PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[])
252 {
253 #if defined(PETSC_USE_COMPLEX)
254 PetscInt i, n;
255 PetscReal *work;
256 #endif
257
258 PetscFunctionBegin;
259 #if !defined(PETSC_USE_COMPLEX)
260 PetscCall(MatGetColumnMeansRealPart(A, means));
261 #else
262 PetscCall(MatGetSize(A, NULL, &n));
263 PetscCall(PetscArrayzero(means, n));
264 PetscCall(PetscCalloc1(n, &work));
265 PetscCall(MatGetColumnMeansRealPart(A, work));
266 for (i = 0; i < n; i++) means[i] = work[i];
267 PetscCall(MatGetColumnMeansImaginaryPart(A, work));
268 for (i = 0; i < n; i++) means[i] += work[i] * PETSC_i;
269 PetscCall(PetscFree(work));
270 #endif
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
274 /*@
275 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix.
276
277 Input Parameters:
278 + A - the matrix
279 - type - A constant defined in `NormType` or `ReductionType`: `NORM_2`, `NORM_1`, `NORM_INFINITY`, `REDUCTION_SUM_REALPART`,
280 `REDUCTION_SUM_IMAGINARYPART`, `REDUCTION_MEAN_REALPART`, `REDUCTION_MEAN_IMAGINARYPART`
281
282 Output Parameter:
283 . reductions - an array as large as the TOTAL number of columns in the matrix
284
285 Level: developer
286
287 Note:
288 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values,
289 if each process wants only some of the values it should extract the ones it wants from the array.
290
291 Developer Notes:
292 This routine is primarily intended as a back-end.
293 `MatGetColumnNorms()`, `MatGetColumnSums()`, and `MatGetColumnMeans()` are implemented using this routine.
294
295 .seealso: [](ch_matrices), `Mat`, `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()`
296 @*/
MatGetColumnReductions(Mat A,PetscInt type,PetscReal reductions[])297 PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[])
298 {
299 PetscFunctionBegin;
300 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
301 PetscUseTypeMethod(A, getcolumnreductions, type, reductions);
302 PetscFunctionReturn(PETSC_SUCCESS);
303 }
304