Actual source code: getcolv.c
1: #include <petsc/private/matimpl.h>
3: /*@
4: MatGetColumnVector - Gets the values from a given column of a matrix.
6: Not Collective
8: Input Parameters:
9: + A - the matrix
10: . yy - the vector
11: - col - the column requested (in global numbering)
13: Level: advanced
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()`.
19: The vector must have the same parallel row layout as the matrix.
21: Contributed by: Denis Vanderstraeten
23: .seealso: [](ch_matrices), `Mat`, `MatGetRow()`, `MatGetDiagonal()`, `MatMult()`
24: @*/
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;
32: PetscFunctionBegin;
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);
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: }
68: /*@
69: MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix.
71: Input Parameters:
72: + A - the matrix
73: - type - `NORM_2`, `NORM_1` or `NORM_INFINITY`
75: Output Parameter:
76: . norms - an array as large as the TOTAL number of columns in the matrix
78: Level: intermediate
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.
84: .seealso: [](ch_matrices), `Mat`, `NormType`, `MatNorm()`
85: @*/
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. */
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: }
98: /*@
99: MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix.
101: Input Parameter:
102: . A - the matrix
104: Output Parameter:
105: . sums - an array as large as the TOTAL number of columns in the matrix
107: Level: intermediate
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.
113: .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
114: @*/
115: PetscErrorCode MatGetColumnSumsRealPart(Mat A, PetscReal sums[])
116: {
117: PetscFunctionBegin;
118: PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_REALPART, sums));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /*@
123: MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix.
125: Input Parameter:
126: . A - the matrix
128: Output Parameter:
129: . sums - an array as large as the TOTAL number of columns in the matrix
131: Level: intermediate
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.
137: .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
138: @*/
139: PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A, PetscReal sums[])
140: {
141: PetscFunctionBegin;
142: PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_IMAGINARYPART, sums));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*@
147: MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix.
149: Input Parameter:
150: . A - the matrix
152: Output Parameter:
153: . sums - an array as large as the TOTAL number of columns in the matrix
155: Level: intermediate
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.
161: .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
162: @*/
163: PetscErrorCode MatGetColumnSums(Mat A, PetscScalar sums[])
164: {
165: #if defined(PETSC_USE_COMPLEX)
166: PetscInt i, n;
167: PetscReal *work;
168: #endif
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: }
186: /*@
187: MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix.
189: Input Parameter:
190: . A - the matrix
192: Output Parameter:
193: . means - an array as large as the TOTAL number of columns in the matrix
195: Level: intermediate
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.
201: .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
202: @*/
203: PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[])
204: {
205: PetscFunctionBegin;
206: PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*@
211: MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix.
213: Input Parameter:
214: . A - the matrix
216: Output Parameter:
217: . means - an array as large as the TOTAL number of columns in the matrix
219: Level: intermediate
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.
225: .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
226: @*/
227: PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[])
228: {
229: PetscFunctionBegin;
230: PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix.
237: Input Parameter:
238: . A - the matrix
240: Output Parameter:
241: . means - an array as large as the TOTAL number of columns in the matrix
243: Level: intermediate
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.
249: .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
250: @*/
251: PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[])
252: {
253: #if defined(PETSC_USE_COMPLEX)
254: PetscInt i, n;
255: PetscReal *work;
256: #endif
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: }
274: /*@
275: MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix.
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`
282: Output Parameter:
283: . reductions - an array as large as the TOTAL number of columns in the matrix
285: Level: developer
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.
291: Developer Notes:
292: This routine is primarily intended as a back-end.
293: `MatGetColumnNorms()`, `MatGetColumnSums()`, and `MatGetColumnMeans()` are implemented using this routine.
295: .seealso: [](ch_matrices), `Mat`, `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()`
296: @*/
297: PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[])
298: {
299: PetscFunctionBegin;
301: PetscUseTypeMethod(A, getcolumnreductions, type, reductions);
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }