Actual source code: axpy.c

petsc-main 2021-04-20
Report Typos and Errors

  2: #include <petsc/private/matimpl.h>

  4: static PetscErrorCode MatTransposeAXPY_Private(Mat Y,PetscScalar a,Mat X,MatStructure str,Mat T)
  5: {
  6:   PetscErrorCode ierr,(*f)(Mat,Mat*);
  7:   Mat            A,F;

 10:   PetscObjectQueryFunction((PetscObject)T,"MatTransposeGetMat_C",&f);
 11:   if (f) {
 12:     MatTransposeGetMat(T,&A);
 13:     if (T == X) {
 14:       PetscInfo(NULL,"Explicitly transposing X of type MATTRANSPOSEMAT to perform MatAXPY()\n");
 15:       MatTranspose(A,MAT_INITIAL_MATRIX,&F);
 16:       A = Y;
 17:     } else {
 18:       PetscInfo(NULL,"Transposing X because Y of type MATTRANSPOSEMAT to perform MatAXPY()\n");
 19:       MatTranspose(X,MAT_INITIAL_MATRIX,&F);
 20:     }
 21:   } else {
 22:     MatHermitianTransposeGetMat(T,&A);
 23:     if (T == X) {
 24:       PetscInfo(NULL,"Explicitly Hermitian transposing X of type MATTRANSPOSEMAT to perform MatAXPY()\n");
 25:       MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&F);
 26:       A = Y;
 27:     } else {
 28:       PetscInfo(NULL,"Hermitian transposing X because Y of type MATTRANSPOSEMAT to perform MatAXPY()\n");
 29:       MatHermitianTranspose(X,MAT_INITIAL_MATRIX,&F);
 30:     }
 31:   }
 32:   MatAXPY(A,a,F,str);
 33:   MatDestroy(&F);
 34:   return(0);
 35: }

 37: /*@
 38:    MatAXPY - Computes Y = a*X + Y.

 40:    Logically  Collective on Mat

 42:    Input Parameters:
 43: +  a - the scalar multiplier
 44: .  X - the first matrix
 45: .  Y - the second matrix
 46: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's)

 48:    Notes: No operation is performed when a is zero.

 50:    Level: intermediate

 52: .seealso: MatAYPX()
 53:  @*/
 54: PetscErrorCode MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)
 55: {
 57:   PetscInt       M1,M2,N1,N2;
 58:   PetscInt       m1,m2,n1,n2;
 59:   MatType        t1,t2;
 60:   PetscBool      sametype,transpose;

 67:   MatGetSize(X,&M1,&N1);
 68:   MatGetSize(Y,&M2,&N2);
 69:   MatGetLocalSize(X,&m1,&n1);
 70:   MatGetLocalSize(Y,&m2,&n2);
 71:   if (M1 != M2 || N1 != N2) SETERRQ4(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Non conforming matrix add: global sizes X %D x %D, Y %D x %D",M1,N1,M2,N2);
 72:   if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: local sizes X %D x %D, Y %D x %D",m1,n1,m2,n2);
 73:   if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (Y)");
 74:   if (!X->assembled) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (X)");
 75:   if (a == (PetscScalar)0.0) return(0);
 76:   if (Y == X) {
 77:     MatScale(Y,1.0 + a);
 78:     return(0);
 79:   }
 80:   MatGetType(X,&t1);
 81:   MatGetType(Y,&t2);
 82:   PetscStrcmp(t1,t2,&sametype);
 83:   PetscLogEventBegin(MAT_AXPY,Y,0,0,0);
 84:   if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) {
 85:     (*Y->ops->axpy)(Y,a,X,str);
 86:   } else {
 87:     PetscStrcmp(t1,MATTRANSPOSEMAT,&transpose);
 88:     if (transpose) {
 89:       MatTransposeAXPY_Private(Y,a,X,str,X);
 90:     } else {
 91:       PetscStrcmp(t2,MATTRANSPOSEMAT,&transpose);
 92:       if (transpose) {
 93:         MatTransposeAXPY_Private(Y,a,X,str,Y);
 94:       } else {
 95:         MatAXPY_Basic(Y,a,X,str);
 96:       }
 97:     }
 98:   }
 99:   PetscLogEventEnd(MAT_AXPY,Y,0,0,0);
100:   return(0);
101: }

103: PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B)
104: {
106:   PetscErrorCode (*preall)(Mat,Mat,Mat*) = NULL;

109:   /* look for any available faster alternative to the general preallocator */
110:   PetscObjectQueryFunction((PetscObject)Y,"MatAXPYGetPreallocation_C",&preall);
111:   if (preall) {
112:     (*preall)(Y,X,B);
113:   } else { /* Use MatPrellocator, assumes same row-col distribution */
114:     Mat      preallocator;
115:     PetscInt r,rstart,rend;
116:     PetscInt m,n,M,N;

118:     MatGetRowUpperTriangular(Y);
119:     MatGetRowUpperTriangular(X);
120:     MatGetSize(Y,&M,&N);
121:     MatGetLocalSize(Y,&m,&n);
122:     MatCreate(PetscObjectComm((PetscObject)Y),&preallocator);
123:     MatSetType(preallocator,MATPREALLOCATOR);
124:     MatSetLayouts(preallocator,Y->rmap,Y->cmap);
125:     MatSetUp(preallocator);
126:     MatGetOwnershipRange(preallocator,&rstart,&rend);
127:     for (r = rstart; r < rend; ++r) {
128:       PetscInt          ncols;
129:       const PetscInt    *row;
130:       const PetscScalar *vals;

132:       MatGetRow(Y,r,&ncols,&row,&vals);
133:       MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES);
134:       MatRestoreRow(Y,r,&ncols,&row,&vals);
135:       MatGetRow(X,r,&ncols,&row,&vals);
136:       MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES);
137:       MatRestoreRow(X,r,&ncols,&row,&vals);
138:     }
139:     MatSetOption(preallocator,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
140:     MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
141:     MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);
142:     MatRestoreRowUpperTriangular(Y);
143:     MatRestoreRowUpperTriangular(X);

145:     MatCreate(PetscObjectComm((PetscObject)Y),B);
146:     MatSetType(*B,((PetscObject)Y)->type_name);
147:     MatSetLayouts(*B,Y->rmap,Y->cmap);
148:     MatPreallocatorPreallocate(preallocator,PETSC_FALSE,*B);
149:     MatDestroy(&preallocator);
150:   }
151:   return(0);
152: }

154: PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
155: {
157:   PetscBool      isshell,isdense,isnest;

160:   MatIsShell(Y,&isshell);
161:   if (isshell) { /* MatShell has special support for AXPY */
162:     PetscErrorCode (*f)(Mat,PetscScalar,Mat,MatStructure);

164:     MatGetOperation(Y,MATOP_AXPY,(void (**)(void))&f);
165:     if (f) {
166:       (*f)(Y,a,X,str);
167:       return(0);
168:     }
169:   }
170:   /* no need to preallocate if Y is dense */
171:   PetscObjectBaseTypeCompareAny((PetscObject)Y,&isdense,MATSEQDENSE,MATMPIDENSE,"");
172:   if (isdense) {
173:     PetscObjectTypeCompare((PetscObject)X,MATNEST,&isnest);
174:     if (isnest) {
175:       MatAXPY_Dense_Nest(Y,a,X);
176:       return(0);
177:     }
178:     if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) str = SUBSET_NONZERO_PATTERN;
179:   }
180:   if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) {
181:     PetscInt          i,start,end,j,ncols,m,n;
182:     const PetscInt    *row;
183:     PetscScalar       *val;
184:     const PetscScalar *vals;

186:     MatGetSize(X,&m,&n);
187:     MatGetOwnershipRange(X,&start,&end);
188:     MatGetRowUpperTriangular(X);
189:     if (a == 1.0) {
190:       for (i = start; i < end; i++) {
191:         MatGetRow(X,i,&ncols,&row,&vals);
192:         MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);
193:         MatRestoreRow(X,i,&ncols,&row,&vals);
194:       }
195:     } else {
196:       PetscInt vs = 100;
197:       /* realloc if needed, as this function may be used in parallel */
198:       PetscMalloc1(vs,&val);
199:       for (i=start; i<end; i++) {
200:         MatGetRow(X,i,&ncols,&row,&vals);
201:         if (vs < ncols) {
202:           vs   = PetscMin(2*ncols,n);
203:           PetscRealloc(vs*sizeof(*val),&val);
204:         }
205:         for (j=0; j<ncols; j++) val[j] = a*vals[j];
206:         MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);
207:         MatRestoreRow(X,i,&ncols,&row,&vals);
208:       }
209:       PetscFree(val);
210:     }
211:     MatRestoreRowUpperTriangular(X);
212:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
213:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
214:   } else {
215:     Mat B;

217:     MatAXPY_Basic_Preallocate(Y,X,&B);
218:     MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
219:     /* TODO mat_tests-ex37_nsize-1_mat_type-baij_mat_block_size-2 fails with MatHeaderMerge */
220:     MatHeaderReplace(Y,&B);
221:   }
222:   return(0);
223: }

225: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str)
226: {
227:   PetscInt          i,start,end,j,ncols,m,n;
228:   PetscErrorCode    ierr;
229:   const PetscInt    *row;
230:   PetscScalar       *val;
231:   const PetscScalar *vals;

234:   MatGetSize(X,&m,&n);
235:   MatGetOwnershipRange(X,&start,&end);
236:   MatGetRowUpperTriangular(Y);
237:   MatGetRowUpperTriangular(X);
238:   if (a == 1.0) {
239:     for (i = start; i < end; i++) {
240:       MatGetRow(Y,i,&ncols,&row,&vals);
241:       MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
242:       MatRestoreRow(Y,i,&ncols,&row,&vals);

244:       MatGetRow(X,i,&ncols,&row,&vals);
245:       MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
246:       MatRestoreRow(X,i,&ncols,&row,&vals);
247:     }
248:   } else {
249:     PetscInt vs = 100;
250:     /* realloc if needed, as this function may be used in parallel */
251:     PetscMalloc1(vs,&val);
252:     for (i=start; i<end; i++) {
253:       MatGetRow(Y,i,&ncols,&row,&vals);
254:       MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
255:       MatRestoreRow(Y,i,&ncols,&row,&vals);

257:       MatGetRow(X,i,&ncols,&row,&vals);
258:       if (vs < ncols) {
259:         vs   = PetscMin(2*ncols,n);
260:         PetscRealloc(vs*sizeof(*val),&val);
261:       }
262:       for (j=0; j<ncols; j++) val[j] = a*vals[j];
263:       MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);
264:       MatRestoreRow(X,i,&ncols,&row,&vals);
265:     }
266:     PetscFree(val);
267:   }
268:   MatRestoreRowUpperTriangular(Y);
269:   MatRestoreRowUpperTriangular(X);
270:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
271:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
272:   return(0);
273: }

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

278:    Neighbor-wise Collective on Mat

280:    Input Parameters:
281: +  Y - the matrices
282: -  a - the PetscScalar

284:    Level: intermediate

286:    Notes:
287:     If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially
288:    fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
289:    entry). No operation is performed when a is zero.

291:    To form Y = Y + diag(V) use MatDiagonalSet()

293: .seealso: MatDiagonalSet(), MatScale(), MatDiagonalScale()
294:  @*/
295: PetscErrorCode  MatShift(Mat Y,PetscScalar a)
296: {

301:   if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
302:   if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
303:   MatCheckPreallocated(Y,1);
304:   if (a == 0.0) return(0);

306:   if (Y->ops->shift) {
307:     (*Y->ops->shift)(Y,a);
308:   } else {
309:     MatShift_Basic(Y,a);
310:   }

312:   PetscObjectStateIncrease((PetscObject)Y);
313:   return(0);
314: }

316: PetscErrorCode  MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is)
317: {
318:   PetscErrorCode    ierr;
319:   PetscInt          i,start,end;
320:   const PetscScalar *v;

323:   MatGetOwnershipRange(Y,&start,&end);
324:   VecGetArrayRead(D,&v);
325:   for (i=start; i<end; i++) {
326:     MatSetValues(Y,1,&i,1,&i,v+i-start,is);
327:   }
328:   VecRestoreArrayRead(D,&v);
329:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
330:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
331:   return(0);
332: }

334: /*@
335:    MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
336:    that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
337:    INSERT_VALUES.

339:    Neighbor-wise Collective on Mat

341:    Input Parameters:
342: +  Y - the input matrix
343: .  D - the diagonal matrix, represented as a vector
344: -  i - INSERT_VALUES or ADD_VALUES

346:    Notes:
347:     If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially
348:    fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
349:    entry).

351:    Level: intermediate

353: .seealso: MatShift(), MatScale(), MatDiagonalScale()
354: @*/
355: PetscErrorCode  MatDiagonalSet(Mat Y,Vec D,InsertMode is)
356: {
358:   PetscInt       matlocal,veclocal;

363:   MatGetLocalSize(Y,&matlocal,NULL);
364:   VecGetLocalSize(D,&veclocal);
365:   if (matlocal != veclocal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number local rows of matrix %D does not match that of vector for diagonal %D",matlocal,veclocal);
366:   if (Y->ops->diagonalset) {
367:     (*Y->ops->diagonalset)(Y,D,is);
368:   } else {
369:     MatDiagonalSet_Default(Y,D,is);
370:   }
371:   PetscObjectStateIncrease((PetscObject)Y);
372:   return(0);
373: }

375: /*@
376:    MatAYPX - Computes Y = a*Y + X.

378:    Logically on Mat

380:    Input Parameters:
381: +  a - the PetscScalar multiplier
382: .  Y - the first matrix
383: .  X - the second matrix
384: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's)

386:    Level: intermediate

388: .seealso: MatAXPY()
389:  @*/
390: PetscErrorCode  MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str)
391: {

395:   MatScale(Y,a);
396:   MatAXPY(Y,1.0,X,str);
397:   return(0);
398: }

400: /*@
401:     MatComputeOperator - Computes the explicit matrix

403:     Collective on Mat

405:     Input Parameter:
406: +   inmat - the matrix
407: -   mattype - the matrix type for the explicit operator

409:     Output Parameter:
410: .   mat - the explict  operator

412:     Notes:
413:     This computation is done by applying the operators to columns of the identity matrix.
414:     This routine is costly in general, and is recommended for use only with relatively small systems.
415:     Currently, this routine uses a dense matrix format if mattype == NULL.

417:     Level: advanced

419: @*/
420: PetscErrorCode  MatComputeOperator(Mat inmat,MatType mattype,Mat *mat)
421: {

427:   MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);
428:   return(0);
429: }

431: /*@
432:     MatComputeOperatorTranspose - Computes the explicit matrix representation of
433:         a give matrix that can apply MatMultTranspose()

435:     Collective on Mat

437:     Input Parameter:
438: .   inmat - the matrix

440:     Output Parameter:
441: .   mat - the explict  operator transposed

443:     Notes:
444:     This computation is done by applying the transpose of the operator to columns of the identity matrix.
445:     This routine is costly in general, and is recommended for use only with relatively small systems.
446:     Currently, this routine uses a dense matrix format if mattype == NULL.

448:     Level: advanced

450: @*/
451: PetscErrorCode  MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat *mat)
452: {
453:   Mat            A;

459:   MatCreateTranspose(inmat,&A);
460:   MatConvert_Shell(A,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);
461:   MatDestroy(&A);
462:   return(0);
463: }

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

468:   Input Parameters:
469: + A   - The matrix
470: - tol - The zero tolerance

472:   Output Parameters:
473: . A - The chopped matrix

475:   Level: intermediate

477: .seealso: MatCreate(), MatZeroEntries()
478:  @*/
479: PetscErrorCode MatChop(Mat A, PetscReal tol)
480: {
481:   Mat            a;
482:   PetscScalar    *newVals;
483:   PetscInt       *newCols;
484:   PetscInt       rStart, rEnd, numRows, maxRows, r, colMax = 0;
485:   PetscBool      flg;

489:   PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, "");
490:   if (flg) {
491:     MatDenseGetLocalMatrix(A, &a);
492:     MatDenseGetLDA(a, &r);
493:     MatGetSize(a, &rStart, &rEnd);
494:     MatDenseGetArray(a, &newVals);
495:     for (; colMax < rEnd; ++colMax) {
496:       for (maxRows = 0; maxRows < rStart; ++maxRows) {
497:         newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) < tol ? 0.0 : newVals[maxRows + colMax * r];
498:       }
499:     }
500:     MatDenseRestoreArray(a, &newVals);
501:   } else {
502:     MatGetOwnershipRange(A, &rStart, &rEnd);
503:     MatGetRowUpperTriangular(A);
504:     for (r = rStart; r < rEnd; ++r) {
505:       PetscInt ncols;

507:       MatGetRow(A, r, &ncols, NULL, NULL);
508:       colMax = PetscMax(colMax, ncols);
509:       MatRestoreRow(A, r, &ncols, NULL, NULL);
510:     }
511:     numRows = rEnd - rStart;
512:     MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));
513:     PetscMalloc2(colMax,&newCols,colMax,&newVals);
514:     MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
515:     for (r = rStart; r < rStart+maxRows; ++r) {
516:       const PetscScalar *vals;
517:       const PetscInt    *cols;
518:       PetscInt           ncols, newcols, c;

520:       if (r < rEnd) {
521:         MatGetRow(A, r, &ncols, &cols, &vals);
522:         for (c = 0; c < ncols; ++c) {
523:           newCols[c] = cols[c];
524:           newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c];
525:         }
526:         newcols = ncols;
527:         MatRestoreRow(A, r, &ncols, &cols, &vals);
528:         MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
529:       }
530:       MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
531:       MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
532:     }
533:     MatRestoreRowUpperTriangular(A);
534:     PetscFree2(newCols,newVals);
535:   }
536:   return(0);
537: }