Actual source code: getcolv.c
2: #include <petsc/private/matimpl.h>
4: /*@
5: MatGetColumnVector - Gets the values from a given column of a matrix.
7: Not Collective
9: Input Parameters:
10: + A - the matrix
11: . yy - the vector
12: - col - the column requested (in global numbering)
14: Level: advanced
16: Notes:
17: If a Mat type does not implement the operation, each processor for which this is called
18: gets the values for its rows using MatGetRow().
20: The vector must have the same parallel row layout as the matrix.
22: Contributed by: Denis Vanderstraeten
24: .seealso: MatGetRow(), MatGetDiagonal(), MatMult()
26: @*/
27: PetscErrorCode MatGetColumnVector(Mat A,Vec yy,PetscInt col)
28: {
29: PetscScalar *y;
30: const PetscScalar *v;
31: PetscErrorCode ierr;
32: PetscInt i,j,nz,N,Rs,Re,rs,re;
33: const PetscInt *idx;
39: if (col < 0) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col);
40: MatGetSize(A,NULL,&N);
41: if (col >= N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N);
42: MatGetOwnershipRange(A,&Rs,&Re);
43: VecGetOwnershipRange(yy,&rs,&re);
44: if (Rs != rs || Re != re) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Matrix %D %D does not have same ownership range (size) as vector %D %D",Rs,Re,rs,re);
46: if (A->ops->getcolumnvector) {
47: (*A->ops->getcolumnvector)(A,yy,col);
48: } else {
49: VecSet(yy,0.0);
50: VecGetArray(yy,&y);
51: /* TODO for general matrices */
52: for (i=Rs; i<Re; i++) {
53: MatGetRow(A,i,&nz,&idx,&v);
54: if (nz && idx[0] <= col) {
55: /*
56: Should use faster search here
57: */
58: for (j=0; j<nz; j++) {
59: if (idx[j] >= col) {
60: if (idx[j] == col) y[i-rs] = v[j];
61: break;
62: }
63: }
64: }
65: MatRestoreRow(A,i,&nz,&idx,&v);
66: }
67: VecRestoreArray(yy,&y);
68: }
69: return(0);
70: }
72: /*@
73: MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix.
75: Input Parameters:
76: + A - the matrix
77: - type - NORM_2, NORM_1 or NORM_INFINITY
79: Output Parameter:
80: . norms - an array as large as the TOTAL number of columns in the matrix
82: Level: intermediate
84: Notes:
85: Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values,
86: if each process wants only some of the values it should extract the ones it wants from the array.
88: .seealso: NormType, MatNorm()
90: @*/
91: PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[])
92: {
93: /* NOTE: MatGetColumnNorms() could simply be a macro that calls MatGetColumnReductions().
94: * I've kept this as a function because it allows slightly more in the way of error checking,
95: * erroring out if MatGetColumnNorms() is not called with a valid NormType. */
99: if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) {
100: MatGetColumnReductions(A,type,norms);
101: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
102: return(0);
103: }
105: /*@
106: MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix.
108: Input Parameter:
109: . A - the matrix
111: Output Parameter:
112: . sums - an array as large as the TOTAL number of columns in the matrix
114: Level: intermediate
116: Notes:
117: Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
118: if each process wants only some of the values it should extract the ones it wants from the array.
120: .seealso: MatGetColumnSumsImaginaryPart(), VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions()
122: @*/
123: PetscErrorCode MatGetColumnSumsRealPart(Mat A,PetscReal sums[])
124: {
128: MatGetColumnReductions(A,REDUCTION_SUM_REALPART,sums);
129: return(0);
130: }
132: /*@
133: MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix.
135: Input Parameter:
136: . A - the matrix
138: Output Parameter:
139: . sums - an array as large as the TOTAL number of columns in the matrix
141: Level: intermediate
143: Notes:
144: Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
145: if each process wants only some of the values it should extract the ones it wants from the array.
147: .seealso: MatGetColumnSumsRealPart(), VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions()
149: @*/
150: PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A,PetscReal sums[])
151: {
155: MatGetColumnReductions(A,REDUCTION_SUM_IMAGINARYPART,sums);
156: return(0);
157: }
159: /*@
160: MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix.
162: Input Parameter:
163: . A - the matrix
165: Output Parameter:
166: . sums - an array as large as the TOTAL number of columns in the matrix
168: Level: intermediate
170: Notes:
171: Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
172: if each process wants only some of the values it should extract the ones it wants from the array.
174: .seealso: VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions()
176: @*/
177: PetscErrorCode MatGetColumnSums(Mat A,PetscScalar sums[])
178: {
180: #if defined(PETSC_USE_COMPLEX)
181: PetscInt i,n;
182: PetscReal *work;
183: #endif
187: #if !defined(PETSC_USE_COMPLEX)
188: MatGetColumnSumsRealPart(A,sums);
189: #else
190: MatGetSize(A,NULL,&n);
191: PetscArrayzero(sums,n);
192: PetscCalloc1(n,&work);
193: MatGetColumnSumsRealPart(A,work);
194: for (i=0; i<n; i++) sums[i] = work[i];
195: MatGetColumnSumsImaginaryPart(A,work);
196: for (i=0; i<n; i++) sums[i] += work[i]*PETSC_i;
197: PetscFree(work);
198: #endif
199: return(0);
200: }
202: /*@
203: MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix.
205: Input Parameter:
206: . A - the matrix
208: Output Parameter:
209: . sums - an array as large as the TOTAL number of columns in the matrix
211: Level: intermediate
213: Notes:
214: Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
215: if each process wants only some of the values it should extract the ones it wants from the array.
217: .seealso: MatGetColumnMeansImaginaryPart(), VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions()
219: @*/
220: PetscErrorCode MatGetColumnMeansRealPart(Mat A,PetscReal means[])
221: {
225: MatGetColumnReductions(A,REDUCTION_MEAN_REALPART,means);
226: return(0);
227: }
229: /*@
230: MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix.
232: Input Parameter:
233: . A - the matrix
235: Output Parameter:
236: . sums - an array as large as the TOTAL number of columns in the matrix
238: Level: intermediate
240: Notes:
241: Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
242: if each process wants only some of the values it should extract the ones it wants from the array.
244: .seealso: MatGetColumnMeansRealPart(), VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions()
246: @*/
247: PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A,PetscReal means[])
248: {
252: MatGetColumnReductions(A,REDUCTION_MEAN_IMAGINARYPART,means);
253: return(0);
254: }
256: /*@
257: MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix.
259: Input Parameter:
260: . A - the matrix
262: Output Parameter:
263: . means - an array as large as the TOTAL number of columns in the matrix
265: Level: intermediate
267: Notes:
268: Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
269: if each process wants only some of the values it should extract the ones it wants from the array.
271: .seealso: VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions()
273: @*/
274: PetscErrorCode MatGetColumnMeans(Mat A,PetscScalar means[])
275: {
277: #if defined(PETSC_USE_COMPLEX)
278: PetscInt i,n;
279: PetscReal *work;
280: #endif
284: #if !defined(PETSC_USE_COMPLEX)
285: MatGetColumnMeansRealPart(A,means);
286: #else
287: MatGetSize(A,NULL,&n);
288: PetscArrayzero(means,n);
289: PetscCalloc1(n,&work);
290: MatGetColumnMeansRealPart(A,work);
291: for (i=0; i<n; i++) means[i] = work[i];
292: MatGetColumnMeansImaginaryPart(A,work);
293: for (i=0; i<n; i++) means[i] += work[i]*PETSC_i;
294: PetscFree(work);
295: #endif
296: return(0);
297: }
299: /*@
300: MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix.
302: Input Parameters:
303: + A - the matrix
304: - type - A constant defined in NormType or ReductionType: NORM_2, NORM_1, NORM_INFINITY, REDUCTION_SUM_REALPART,
305: REDUCTION_SUM_IMAGINARYPART, REDUCTION_MEAN_REALPART, REDUCTION_MEAN_IMAGINARYPART
307: Output Parameter:
308: . reductions - an array as large as the TOTAL number of columns in the matrix
310: Level: developer
312: Notes:
313: Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values,
314: if each process wants only some of the values it should extract the ones it wants from the array.
316: Developer Note:
317: This routine is primarily intended as a back-end.
318: MatGetColumnNorms(), MatGetColumnSums(), and MatGetColumnMeans() are implemented using this routine.
320: .seealso: ReductionType, NormType, MatGetColumnNorms(), MatGetColumnSums(), MatGetColumnMeans()
322: @*/
323: PetscErrorCode MatGetColumnReductions(Mat A,PetscInt type,PetscReal reductions[])
324: {
329: if (A->ops->getcolumnreductions) {
330: (*A->ops->getcolumnreductions)(A,type,reductions);
331: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type");
332: return(0);
333: }