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:   PetscInt          i,j,nz,N,Rs,Re,rs,re;
 32:   const PetscInt    *idx;

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

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

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

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

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

 80:    Level: intermediate

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

 86: .seealso: NormType, MatNorm()

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

 95:   if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) {
 96:     MatGetColumnReductions(A,type,norms);
 97:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
 98:   return 0;
 99: }

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

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

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

110:    Level: intermediate

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

116: .seealso: MatGetColumnSumsImaginaryPart(), VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions()

118: @*/
119: PetscErrorCode MatGetColumnSumsRealPart(Mat A,PetscReal sums[])
120: {
121:   MatGetColumnReductions(A,REDUCTION_SUM_REALPART,sums);
122:   return 0;
123: }

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

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

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

134:    Level: intermediate

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

140: .seealso: MatGetColumnSumsRealPart(), VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions()

142: @*/
143: PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A,PetscReal sums[])
144: {
145:   MatGetColumnReductions(A,REDUCTION_SUM_IMAGINARYPART,sums);
146:   return 0;
147: }

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

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

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

158:    Level: intermediate

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

164: .seealso: VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions()

166: @*/
167: PetscErrorCode MatGetColumnSums(Mat A,PetscScalar sums[])
168: {
169: #if defined(PETSC_USE_COMPLEX)
170:   PetscInt       i,n;
171:   PetscReal      *work;
172: #endif


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

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

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

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

199:    Level: intermediate

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

205: .seealso: MatGetColumnMeansImaginaryPart(), VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions()

207: @*/
208: PetscErrorCode MatGetColumnMeansRealPart(Mat A,PetscReal means[])
209: {
210:   MatGetColumnReductions(A,REDUCTION_MEAN_REALPART,means);
211:   return 0;
212: }

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

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

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

223:    Level: intermediate

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

229: .seealso: MatGetColumnMeansRealPart(), VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions()

231: @*/
232: PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A,PetscReal means[])
233: {
234:   MatGetColumnReductions(A,REDUCTION_MEAN_IMAGINARYPART,means);
235:   return 0;
236: }

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

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

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

247:    Level: intermediate

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

253: .seealso: VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions()

255: @*/
256: PetscErrorCode MatGetColumnMeans(Mat A,PetscScalar means[])
257: {
258: #if defined(PETSC_USE_COMPLEX)
259:   PetscInt       i,n;
260:   PetscReal      *work;
261: #endif


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

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

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

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

290:    Level: developer

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

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

300: .seealso: ReductionType, NormType, MatGetColumnNorms(), MatGetColumnSums(), MatGetColumnMeans()

302: @*/
303: PetscErrorCode MatGetColumnReductions(Mat A,PetscInt type,PetscReal reductions[])
304: {
306:   if (A->ops->getcolumnreductions) {
307:     (*A->ops->getcolumnreductions)(A,type,reductions);
308:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type");
309:   return 0;
310: }