Actual source code: axpy.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/matimpl.h>  /*I   "petscmat.h"  I*/

  6: /*@
  7:    MatAXPY - Computes Y = a*X + Y.

  9:    Logically  Collective on Mat

 11:    Input Parameters:
 12: +  a - the scalar multiplier
 13: .  X - the first matrix
 14: .  Y - the second matrix
 15: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN
 16:          or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's)

 18:    Level: intermediate

 20: .keywords: matrix, add

 22: .seealso: MatAYPX()
 23:  @*/
 24: PetscErrorCode  MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)
 25: {
 27:   PetscInt       m1,m2,n1,n2;

 33:   MatGetSize(X,&m1,&n1);
 34:   MatGetSize(Y,&m2,&n2);
 35:   if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %D %D %D %D",m1,m2,n1,n2);

 37:   PetscLogEventBegin(MAT_AXPY,Y,0,0,0);
 38:   if (Y->ops->axpy) {
 39:     (*Y->ops->axpy)(Y,a,X,str);
 40:   } else {
 41:     MatAXPY_Basic(Y,a,X,str);
 42:   }
 43:   PetscLogEventEnd(MAT_AXPY,Y,0,0,0);
 44: #if defined(PETSC_HAVE_CUSP)
 45:   if (Y->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
 46:     Y->valid_GPU_matrix = PETSC_CUSP_CPU;
 47:   }
 48: #endif
 49:   return(0);
 50: }

 54: PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
 55: {
 56:   PetscInt          i,start,end,j,ncols,m,n;
 57:   PetscErrorCode    ierr;
 58:   const PetscInt    *row;
 59:   PetscScalar       *val;
 60:   const PetscScalar *vals;

 63:   MatGetSize(X,&m,&n);
 64:   MatGetOwnershipRange(X,&start,&end);
 65:   if (a == 1.0) {
 66:     for (i = start; i < end; i++) {
 67:       MatGetRow(X,i,&ncols,&row,&vals);
 68:       MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);
 69:       MatRestoreRow(X,i,&ncols,&row,&vals);
 70:     }
 71:   } else {
 72:     PetscMalloc((n+1)*sizeof(PetscScalar),&val);
 73:     for (i=start; i<end; i++) {
 74:       MatGetRow(X,i,&ncols,&row,&vals);
 75:       for (j=0; j<ncols; j++) {
 76:         val[j] = a*vals[j];
 77:       }
 78:       MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);
 79:       MatRestoreRow(X,i,&ncols,&row,&vals);
 80:     }
 81:     PetscFree(val);
 82:   }
 83:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
 84:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
 85:   return(0);
 86: }

 90: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str)
 91: {
 92:   PetscInt          i,start,end,j,ncols,m,n;
 93:   PetscErrorCode    ierr;
 94:   const PetscInt    *row;
 95:   PetscScalar       *val;
 96:   const PetscScalar *vals;

 99:   MatGetSize(X,&m,&n);
100:   MatGetOwnershipRange(X,&start,&end);
101:   if (a == 1.0) {
102:     for (i = start; i < end; i++) {
103:       MatGetRow(Y,i,&ncols,&row,&vals);
104:       MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
105:       MatRestoreRow(Y,i,&ncols,&row,&vals);

107:       MatGetRow(X,i,&ncols,&row,&vals);
108:       MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
109:       MatRestoreRow(X,i,&ncols,&row,&vals);
110:     }
111:   } else {
112:     PetscMalloc((n+1)*sizeof(PetscScalar),&val);
113:     for (i=start; i<end; i++) {
114:       MatGetRow(Y,i,&ncols,&row,&vals);
115:       MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
116:       MatRestoreRow(Y,i,&ncols,&row,&vals);

118:       MatGetRow(X,i,&ncols,&row,&vals);
119:       for (j=0; j<ncols; j++) {
120:         val[j] = a*vals[j];
121:       }
122:       MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);
123:       MatRestoreRow(X,i,&ncols,&row,&vals);
124:     }
125:     PetscFree(val);
126:   }
127:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
128:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
129:   return(0);
130: }

134: /*@
135:    MatShift - Computes Y =  Y + a I, where a is a PetscScalar and I is the identity matrix.

137:    Neighbor-wise Collective on Mat

139:    Input Parameters:
140: +  Y - the matrices
141: -  a - the PetscScalar

143:    Level: intermediate

145: .keywords: matrix, add, shift

147: .seealso: MatDiagonalSet()
148:  @*/
149: PetscErrorCode  MatShift(Mat Y,PetscScalar a)
150: {
152:   PetscInt       i,start,end;

156:   if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
157:   if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
158:   MatCheckPreallocated(Y,1);

160:   if (Y->ops->shift) {
161:     (*Y->ops->shift)(Y,a);
162:   } else {
163:     PetscScalar alpha = a;
164:     MatGetOwnershipRange(Y,&start,&end);
165:     for (i=start; i<end; i++) {
166:       MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
167:     }
168:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
169:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
170:   }
171: #if defined(PETSC_HAVE_CUSP)
172:   if (Y->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
173:     Y->valid_GPU_matrix = PETSC_CUSP_CPU;
174:   }
175: #endif
176:   return(0);
177: }

181: PetscErrorCode  MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is)
182: {
184:   PetscInt       i,start,end,vstart,vend;
185:   PetscScalar    *v;

188:   VecGetOwnershipRange(D,&vstart,&vend);
189:   MatGetOwnershipRange(Y,&start,&end);
190:   if (vstart != start || vend != end) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %D %D vec %D %D mat",vstart,vend,start,end);
191:   VecGetArray(D,&v);
192:   for (i=start; i<end; i++) {
193:     MatSetValues(Y,1,&i,1,&i,v+i-start,is);
194:   }
195:   VecRestoreArray(D,&v);
196:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
197:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
198:   return(0);
199: }

203: /*@
204:    MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
205:    that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
206:    INSERT_VALUES.

208:    Input Parameters:
209: +  Y - the input matrix
210: .  D - the diagonal matrix, represented as a vector
211: -  i - INSERT_VALUES or ADD_VALUES

213:    Neighbor-wise Collective on Mat and Vec

215:    Level: intermediate

217: .keywords: matrix, add, shift, diagonal

219: .seealso: MatShift()
220: @*/
221: PetscErrorCode  MatDiagonalSet(Mat Y,Vec D,InsertMode is)
222: {

228:   if (Y->ops->diagonalset) {
229:     (*Y->ops->diagonalset)(Y,D,is);
230:   } else {
231:     MatDiagonalSet_Default(Y,D,is);
232:   }
233:   return(0);
234: }

238: /*@
239:    MatAYPX - Computes Y = a*Y + X.

241:    Logically on Mat

243:    Input Parameters:
244: +  a - the PetscScalar multiplier
245: .  Y - the first matrix
246: .  X - the second matrix
247: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN

249:    Level: intermediate

251: .keywords: matrix, add

253: .seealso: MatAXPY()
254:  @*/
255: PetscErrorCode  MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str)
256: {
257:   PetscScalar    one = 1.0;
259:   PetscInt       mX,mY,nX,nY;

265:   MatGetSize(X,&mX,&nX);
266:   MatGetSize(X,&mY,&nY);
267:   if (mX != mY || nX != nY) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrices: %D %D first %D %D second",mX,mY,nX,nY);

269:   MatScale(Y,a);
270:   MatAXPY(Y,one,X,str);
271:   return(0);
272: }

276: /*@
277:     MatComputeExplicitOperator - Computes the explicit matrix

279:     Collective on Mat

281:     Input Parameter:
282: .   inmat - the matrix

284:     Output Parameter:
285: .   mat - the explict preconditioned operator

287:     Notes:
288:     This computation is done by applying the operators to columns of the
289:     identity matrix.

291:     Currently, this routine uses a dense matrix format when 1 processor
292:     is used and a sparse format otherwise.  This routine is costly in general,
293:     and is recommended for use only with relatively small systems.

295:     Level: advanced

297: .keywords: Mat, compute, explicit, operator
298: @*/
299: PetscErrorCode  MatComputeExplicitOperator(Mat inmat,Mat *mat)
300: {
301:   Vec            in,out;
303:   PetscInt       i,m,n,M,N,*rows,start,end;
304:   MPI_Comm       comm;
305:   PetscScalar    *array,zero = 0.0,one = 1.0;
306:   PetscMPIInt    size;


312:   PetscObjectGetComm((PetscObject)inmat,&comm);
313:   MPI_Comm_size(comm,&size);

315:   MatGetLocalSize(inmat,&m,&n);
316:   MatGetSize(inmat,&M,&N);
317:   MatGetVecs(inmat,&in,&out);
318:   VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
319:   VecGetOwnershipRange(out,&start,&end);
320:   PetscMalloc(m*sizeof(PetscInt),&rows);
321:   for (i=0; i<m; i++) rows[i] = start + i;

323:   MatCreate(comm,mat);
324:   MatSetSizes(*mat,m,n,M,N);
325:   if (size == 1) {
326:     MatSetType(*mat,MATSEQDENSE);
327:     MatSeqDenseSetPreallocation(*mat,NULL);
328:   } else {
329:     MatSetType(*mat,MATMPIAIJ);
330:     MatMPIAIJSetPreallocation(*mat,n,NULL,N-n,NULL);
331:   }

333:   for (i=0; i<N; i++) {

335:     VecSet(in,zero);
336:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
337:     VecAssemblyBegin(in);
338:     VecAssemblyEnd(in);

340:     MatMult(inmat,in,out);

342:     VecGetArray(out,&array);
343:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
344:     VecRestoreArray(out,&array);

346:   }
347:   PetscFree(rows);
348:   VecDestroy(&out);
349:   VecDestroy(&in);
350:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
351:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
352:   return(0);
353: }

355: /* Get the map xtoy which is used by MatAXPY() in the case of SUBSET_NONZERO_PATTERN */
358: PetscErrorCode MatAXPYGetxtoy_Private(PetscInt m,PetscInt *xi,PetscInt *xj,PetscInt *xgarray, PetscInt *yi,PetscInt *yj,PetscInt *ygarray, PetscInt **xtoy)
359: {
361:   PetscInt       row,i,nz,xcol,ycol,jx,jy,*x2y;

364:   PetscMalloc(xi[m]*sizeof(PetscInt),&x2y);
365:   i    = 0;
366:   for (row=0; row<m; row++) {
367:     nz = xi[1] - xi[0];
368:     jy = 0;
369:     for (jx=0; jx<nz; jx++,jy++) {
370:       if (xgarray && ygarray) {
371:         xcol = xgarray[xj[*xi + jx]];
372:         ycol = ygarray[yj[*yi + jy]];
373:       } else {
374:         xcol = xj[*xi + jx];
375:         ycol = yj[*yi + jy];  /* col index for y */
376:       }
377:       while (ycol < xcol) {
378:         jy++;
379:         if (ygarray) {
380:           ycol = ygarray[yj[*yi + jy]];
381:         } else {
382:           ycol = yj[*yi + jy];
383:         }
384:       }
385:       if (xcol != ycol) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"X matrix entry (%D,%D) is not in Y matrix",row,ycol);
386:       x2y[i++] = *yi + jy;
387:     }
388:     xi++; yi++;
389:   }
390:   *xtoy = x2y;
391:   return(0);
392: }

396: /*@
397:   MatChop - Set all values in the matrix less than the tolerance to zero

399:   Input Parameters:
400: + A   - The matrix
401: - tol - The zero tolerance

403:   Output Parameters:
404: . A - The chopped matrix

406:   Level: intermediate

408: .seealso: MatCreate(), MatZeroEntries()
409:  @*/
410: PetscErrorCode MatChop(Mat A, PetscReal tol)
411: {
412:   PetscScalar    *newVals;
413:   PetscInt       *newCols;
414:   PetscInt       rStart, rEnd, numRows, maxRows, r, colMax = 0;

418:   MatGetOwnershipRange(A, &rStart, &rEnd);
419:   for (r = rStart; r < rEnd; ++r) {
420:     PetscInt ncols;

422:     MatGetRow(A, r, &ncols, NULL, NULL);
423:     colMax = PetscMax(colMax, ncols);
424:     MatRestoreRow(A, r, &ncols, NULL, NULL);
425:   }
426:   numRows = rEnd - rStart;
427:   MPI_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD);
428:   PetscMalloc2(colMax,PetscInt,&newCols,colMax,PetscScalar,&newVals);
429:   for (r = rStart; r < rStart+maxRows; ++r) {
430:     const PetscScalar *vals;
431:     const PetscInt    *cols;
432:     PetscInt           ncols, newcols, c;

434:     if (r < rEnd) {
435:       MatGetRow(A, r, &ncols, &cols, &vals);
436:       for (c = 0; c < ncols; ++c) {
437:         newCols[c] = cols[c];
438:         newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c];
439:       }
440:       newcols = ncols;
441:       MatRestoreRow(A, r, &ncols, &cols, &vals);
442:       MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
443:     }
444:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
445:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
446:   }
447:   PetscFree2(newCols,newVals);
448:   return(0);
449: }