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()`
 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;

 37:   MatGetSize(A, NULL, &N);
 39:   MatGetOwnershipRange(A, &Rs, &Re);
 40:   VecGetOwnershipRange(yy, &rs, &re);

 43:   if (A->ops->getcolumnvector) PetscUseTypeMethod(A, getcolumnvector, yy, col);
 44:   else {
 45:     VecSet(yy, 0.0);
 46:     VecGetArray(yy, &y);
 47:     /* TODO for general matrices */
 48:     for (i = Rs; i < Re; i++) {
 49:       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:       MatRestoreRow(A, i, &nz, &idx, &v);
 62:     }
 63:     VecRestoreArray(yy, &y);
 64:   }
 65:   return 0;
 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: `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:   if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) {
 93:     MatGetColumnReductions(A, type, norms);
 94:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown NormType");
 95:   return 0;
 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: `MatGetColumnSumsImaginaryPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
114: @*/
115: PetscErrorCode MatGetColumnSumsRealPart(Mat A, PetscReal sums[])
116: {
117:   MatGetColumnReductions(A, REDUCTION_SUM_REALPART, sums);
118:   return 0;
119: }

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

124:    Input Parameter:
125: .  A - the matrix

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

130:    Level: intermediate

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

136: .seealso: `MatGetColumnSumsRealPart()`, `VecSum()`, `MatGetColumnMeans()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
137: @*/
138: PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A, PetscReal sums[])
139: {
140:   MatGetColumnReductions(A, REDUCTION_SUM_IMAGINARYPART, sums);
141:   return 0;
142: }

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

147:    Input Parameter:
148: .  A - the matrix

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

153:    Level: intermediate

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

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


169: #if !defined(PETSC_USE_COMPLEX)
170:   MatGetColumnSumsRealPart(A, sums);
171: #else
172:   MatGetSize(A, NULL, &n);
173:   PetscArrayzero(sums, n);
174:   PetscCalloc1(n, &work);
175:   MatGetColumnSumsRealPart(A, work);
176:   for (i = 0; i < n; i++) sums[i] = work[i];
177:   MatGetColumnSumsImaginaryPart(A, work);
178:   for (i = 0; i < n; i++) sums[i] += work[i] * PETSC_i;
179:   PetscFree(work);
180: #endif
181:   return 0;
182: }

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

187:    Input Parameter:
188: .  A - the matrix

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

193:    Level: intermediate

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

199: .seealso: `MatGetColumnMeansImaginaryPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
200: @*/
201: PetscErrorCode MatGetColumnMeansRealPart(Mat A, PetscReal means[])
202: {
203:   MatGetColumnReductions(A, REDUCTION_MEAN_REALPART, means);
204:   return 0;
205: }

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

210:    Input Parameter:
211: .  A - the matrix

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

216:    Level: intermediate

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

222: .seealso: `MatGetColumnMeansRealPart()`, `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
223: @*/
224: PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A, PetscReal means[])
225: {
226:   MatGetColumnReductions(A, REDUCTION_MEAN_IMAGINARYPART, means);
227:   return 0;
228: }

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

233:    Input Parameter:
234: .  A - the matrix

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

239:    Level: intermediate

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

245: .seealso: `VecSum()`, `MatGetColumnSums()`, `MatGetColumnNorms()`, `MatGetColumnReductions()`
246: @*/
247: PetscErrorCode MatGetColumnMeans(Mat A, PetscScalar means[])
248: {
249: #if defined(PETSC_USE_COMPLEX)
250:   PetscInt   i, n;
251:   PetscReal *work;
252: #endif


255: #if !defined(PETSC_USE_COMPLEX)
256:   MatGetColumnMeansRealPart(A, means);
257: #else
258:   MatGetSize(A, NULL, &n);
259:   PetscArrayzero(means, n);
260:   PetscCalloc1(n, &work);
261:   MatGetColumnMeansRealPart(A, work);
262:   for (i = 0; i < n; i++) means[i] = work[i];
263:   MatGetColumnMeansImaginaryPart(A, work);
264:   for (i = 0; i < n; i++) means[i] += work[i] * PETSC_i;
265:   PetscFree(work);
266: #endif
267:   return 0;
268: }

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

273:   Input Parameters:
274: +  A - the matrix
275: -  type - A constant defined in `NormType` or `ReductionType`: `NORM_2`, `NORM_1`, `NORM_INFINITY`, `REDUCTION_SUM_REALPART`,
276:           `REDUCTION_SUM_IMAGINARYPART`, `REDUCTION_MEAN_REALPART`, `REDUCTION_MEAN_IMAGINARYPART`

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

281:    Level: developer

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

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

291: .seealso: `ReductionType`, `NormType`, `MatGetColumnNorms()`, `MatGetColumnSums()`, `MatGetColumnMeans()`
292: @*/
293: PetscErrorCode MatGetColumnReductions(Mat A, PetscInt type, PetscReal reductions[])
294: {
296:   PetscUseTypeMethod(A, getcolumnreductions, type, reductions);
297:   return 0;
298: }