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 `MatType` 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: [](ch_matrices), `Mat`, `MatGetRow()`, `MatGetDiagonal()`, `MatMult()`
 25: @*/
 26: PetscErrorCode MatGetColumnVector(Mat A, Vec yy, PetscInt col)
 27: {
 28:   PetscScalar       *y;
 29:   const PetscScalar *v;
 30:   PetscInt           i, j, nz, N, Rs, Re, rs, re;
 31:   const PetscInt    *idx;

 33:   PetscFunctionBegin;
 37:   PetscCheck(col >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Requested negative column: %" PetscInt_FMT, col);
 38:   PetscCall(MatGetSize(A, NULL, &N));
 39:   PetscCheck(col < N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Requested column %" PetscInt_FMT " larger than number columns in matrix %" PetscInt_FMT, col, N);
 40:   PetscCall(MatGetOwnershipRange(A, &Rs, &Re));
 41:   PetscCall(VecGetOwnershipRange(yy, &rs, &re));
 42:   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);

 44:   if (A->ops->getcolumnvector) PetscUseTypeMethod(A, getcolumnvector, yy, col);
 45:   else {
 46:     PetscCall(VecSet(yy, 0.0));
 47:     PetscCall(VecGetArray(yy, &y));
 48:     /* TODO for general matrices */
 49:     for (i = Rs; i < Re; i++) {
 50:       PetscCall(MatGetRow(A, i, &nz, &idx, &v));
 51:       if (nz && idx[0] <= col) {
 52:         /*
 53:           Should use faster search here
 54:         */
 55:         for (j = 0; j < nz; j++) {
 56:           if (idx[j] >= col) {
 57:             if (idx[j] == col) y[i - rs] = v[j];
 58:             break;
 59:           }
 60:         }
 61:       }
 62:       PetscCall(MatRestoreRow(A, i, &nz, &idx, &v));
 63:     }
 64:     PetscCall(VecRestoreArray(yy, &y));
 65:   }
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: /*@
 70:    MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix.

 72:    Input Parameters:
 73: +  A - the matrix
 74: -  type - `NORM_2`, `NORM_1` or `NORM_INFINITY`

 76:    Output Parameter:
 77: .  norms - an array as large as the TOTAL number of columns in the matrix

 79:    Level: intermediate

 81:    Note:
 82:     Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values,
 83:     if each process wants only some of the values it should extract the ones it wants from the array.

 85: .seealso: [](ch_matrices), `Mat`, `NormType`, `MatNorm()`
 86: @*/
 87: PetscErrorCode MatGetColumnNorms(Mat A, NormType type, PetscReal norms[])
 88: {
 89:   /* NOTE: MatGetColumnNorms() could simply be a macro that calls MatGetColumnReductions().
 90:    * I've kept this as a function because it allows slightly more in the way of error checking,
 91:    * erroring out if MatGetColumnNorms() is not called with a valid NormType. */

 93:   PetscFunctionBegin;
 94:   if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) {
 95:     PetscCall(MatGetColumnReductions(A, type, norms));
 96:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown NormType");
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: /*@
101:    MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix.

103:    Input Parameter:
104: .  A - the matrix

106:    Output Parameter:
107: .  sums - an array as large as the TOTAL number of columns in the matrix

109:    Level: intermediate

111:    Note:
112:     Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
113:     if each process wants only some of the values it should extract the ones it wants from the array.

115: .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
116: @*/
117: PetscErrorCode MatGetColumnSumsRealPart(Mat A, PetscReal sums[])
118: {
119:   PetscFunctionBegin;
120:   PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_REALPART, sums));
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

124: /*@
125:    MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix.

127:    Input Parameter:
128: .  A - the matrix

130:    Output Parameter:
131: .  sums - an array as large as the TOTAL number of columns in the matrix

133:    Level: intermediate

135:    Note:
136:     Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
137:     if each process wants only some of the values it should extract the ones it wants from the array.

139: .seealso: [](ch_matrices), `Mat`, `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
140: @*/
141: PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A, PetscReal sums[])
142: {
143:   PetscFunctionBegin;
144:   PetscCall(MatGetColumnReductions(A, REDUCTION_SUM_IMAGINARYPART, sums));
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

148: /*@
149:    MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix.

151:    Input Parameter:
152: .  A - the matrix

154:    Output Parameter:
155: .  sums - an array as large as the TOTAL number of columns in the matrix

157:    Level: intermediate

159:    Note:
160:     Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values,
161:     if each process wants only some of the values it should extract the ones it wants from the array.

163: .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
164: @*/
165: PetscErrorCode MatGetColumnSums(Mat A, PetscScalar sums[])
166: {
167: #if defined(PETSC_USE_COMPLEX)
168:   PetscInt   i, n;
169:   PetscReal *work;
170: #endif

172:   PetscFunctionBegin;

174: #if !defined(PETSC_USE_COMPLEX)
175:   PetscCall(MatGetColumnSumsRealPart(A, sums));
176: #else
177:   PetscCall(MatGetSize(A, NULL, &n));
178:   PetscCall(PetscArrayzero(sums, n));
179:   PetscCall(PetscCalloc1(n, &work));
180:   PetscCall(MatGetColumnSumsRealPart(A, work));
181:   for (i = 0; i < n; i++) sums[i] = work[i];
182:   PetscCall(MatGetColumnSumsImaginaryPart(A, work));
183:   for (i = 0; i < n; i++) sums[i] += work[i] * PETSC_i;
184:   PetscCall(PetscFree(work));
185: #endif
186:   PetscFunctionReturn(PETSC_SUCCESS);
187: }

189: /*@
190:    MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix.

192:    Input Parameter:
193: .  A - the matrix

195:    Output Parameter:
196: .  sums - an array as large as the TOTAL number of columns in the matrix

198:    Level: intermediate

200:    Note:
201:     Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
202:     if each process wants only some of the values it should extract the ones it wants from the array.

204: .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
205: @*/
206: PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[])
207: {
208:   PetscFunctionBegin;
209:   PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means));
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: /*@
214:    MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix.

216:    Input Parameter:
217: .  A - the matrix

219:    Output Parameter:
220: .  sums - an array as large as the TOTAL number of columns in the matrix

222:    Level: intermediate

224:    Note:
225:     Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
226:     if each process wants only some of the values it should extract the ones it wants from the array.

228: .seealso: [](ch_matrices), `Mat`, `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
229: @*/
230: PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[])
231: {
232:   PetscFunctionBegin;
233:   PetscCall(MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means));
234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }

237: /*@
238:    MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix.

240:    Input Parameter:
241: .  A - the matrix

243:    Output Parameter:
244: .  means - an array as large as the TOTAL number of columns in the matrix

246:    Level: intermediate

248:    Note:
249:     Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values,
250:     if each process wants only some of the values it should extract the ones it wants from the array.

252: .seealso: [](ch_matrices), `Mat`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
253: @*/
254: PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[])
255: {
256: #if defined(PETSC_USE_COMPLEX)
257:   PetscInt   i, n;
258:   PetscReal *work;
259: #endif

261:   PetscFunctionBegin;

263: #if !defined(PETSC_USE_COMPLEX)
264:   PetscCall(MatGetColumnMeansRealPart(A, means));
265: #else
266:   PetscCall(MatGetSize(A, NULL, &n));
267:   PetscCall(PetscArrayzero(means, n));
268:   PetscCall(PetscCalloc1(n, &work));
269:   PetscCall(MatGetColumnMeansRealPart(A, work));
270:   for (i = 0; i < n; i++) means[i] = work[i];
271:   PetscCall(MatGetColumnMeansImaginaryPart(A, work));
272:   for (i = 0; i < n; i++) means[i] += work[i] * PETSC_i;
273:   PetscCall(PetscFree(work));
274: #endif
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: /*@
279:     MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix.

281:   Input Parameters:
282: +  A - the matrix
283: -  type - A constant defined in `NormType` or `ReductionType`: `NORM_2`, `NORM_1`, `NORM_INFINITY`, `REDUCTION_SUM_REALPART`,
284:           `REDUCTION_SUM_IMAGINARYPART`, `REDUCTION_MEAN_REALPART`, `REDUCTION_MEAN_IMAGINARYPART`

286:   Output Parameter:
287: .  reductions - an array as large as the TOTAL number of columns in the matrix

289:    Level: developer

291:    Note:
292:     Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values,
293:     if each process wants only some of the values it should extract the ones it wants from the array.

295:   Developer Note:
296:     This routine is primarily intended as a back-end.
297:     `MatGetColumnNorms()`, `MatGetColumnSums()`, and `MatGetColumnMeans()` are implemented using this routine.

299: .seealso: [](ch_matrices), `Mat`, `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()`
300: @*/
301: PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[])
302: {
303:   PetscFunctionBegin;
305:   PetscUseTypeMethod(A, getcolumnreductions, type, reductions);
306:   PetscFunctionReturn(PETSC_SUCCESS);
307: }