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: }