Actual source code: fdmatrix.c
petsc-3.13.6 2020-09-29
2: /*
3: This is where the abstract matrix operations are defined that are
4: used for finite difference computations of Jacobians using coloring.
5: */
7: #include <petsc/private/matimpl.h>
8: #include <petsc/private/isimpl.h>
10: PetscErrorCode MatFDColoringSetF(MatFDColoring fd,Vec F)
11: {
15: if (F) {
16: VecCopy(F,fd->w1);
17: fd->fset = PETSC_TRUE;
18: } else {
19: fd->fset = PETSC_FALSE;
20: }
21: return(0);
22: }
24: #include <petscdraw.h>
25: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
26: {
27: MatFDColoring fd = (MatFDColoring)Aa;
29: PetscInt i,j,nz,row;
30: PetscReal x,y;
31: MatEntry *Jentry=fd->matentry;
34: /* loop over colors */
35: nz = 0;
36: for (i=0; i<fd->ncolors; i++) {
37: for (j=0; j<fd->nrows[i]; j++) {
38: row = Jentry[nz].row;
39: y = fd->M - row - fd->rstart;
40: x = (PetscReal)Jentry[nz++].col;
41: PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
42: }
43: }
44: return(0);
45: }
47: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
48: {
50: PetscBool isnull;
51: PetscDraw draw;
52: PetscReal xr,yr,xl,yl,h,w;
55: PetscViewerDrawGetDraw(viewer,0,&draw);
56: PetscDrawIsNull(draw,&isnull);
57: if (isnull) return(0);
59: xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0;
60: xr += w; yr += h; xl = -w; yl = -h;
61: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
62: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
63: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
64: PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
65: PetscDrawSave(draw);
66: return(0);
67: }
69: /*@C
70: MatFDColoringView - Views a finite difference coloring context.
72: Collective on MatFDColoring
74: Input Parameters:
75: + c - the coloring context
76: - viewer - visualization context
78: Level: intermediate
80: Notes:
81: The available visualization contexts include
82: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
83: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
84: output where only the first processor opens
85: the file. All other processors send their
86: data to the first processor to print.
87: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
89: Notes:
90: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
91: involves more than 33 then some seemingly identical colors are displayed making it look
92: like an illegal coloring. This is just a graphical artifact.
94: .seealso: MatFDColoringCreate()
96: @*/
97: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
98: {
99: PetscErrorCode ierr;
100: PetscInt i,j;
101: PetscBool isdraw,iascii;
102: PetscViewerFormat format;
106: if (!viewer) {
107: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
108: }
112: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
113: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
114: if (isdraw) {
115: MatFDColoringView_Draw(c,viewer);
116: } else if (iascii) {
117: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
118: PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);
119: PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);
120: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
122: PetscViewerGetFormat(viewer,&format);
123: if (format != PETSC_VIEWER_ASCII_INFO) {
124: PetscInt row,col,nz;
125: nz = 0;
126: for (i=0; i<c->ncolors; i++) {
127: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
128: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
129: for (j=0; j<c->ncolumns[i]; j++) {
130: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
131: }
132: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
133: if (c->matentry) {
134: for (j=0; j<c->nrows[i]; j++) {
135: row = c->matentry[nz].row;
136: col = c->matentry[nz++].col;
137: PetscViewerASCIIPrintf(viewer," %D %D \n",row,col);
138: }
139: }
140: }
141: }
142: PetscViewerFlush(viewer);
143: }
144: return(0);
145: }
147: /*@
148: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
149: a Jacobian matrix using finite differences.
151: Logically Collective on MatFDColoring
153: The Jacobian is estimated with the differencing approximation
154: .vb
155: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
156: htype = 'ds':
157: h = error_rel*u[i] if abs(u[i]) > umin
158: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
159: dx_{i} = (0, ... 1, .... 0)
161: htype = 'wp':
162: h = error_rel * sqrt(1 + ||u||)
163: .ve
165: Input Parameters:
166: + coloring - the coloring context
167: . error_rel - relative error
168: - umin - minimum allowable u-value magnitude
170: Level: advanced
172: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
174: @*/
175: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
176: {
181: if (error != PETSC_DEFAULT) matfd->error_rel = error;
182: if (umin != PETSC_DEFAULT) matfd->umin = umin;
183: return(0);
184: }
186: /*@
187: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
189: Logically Collective on MatFDColoring
191: Input Parameters:
192: + coloring - the coloring context
193: . brows - number of rows in the block
194: - bcols - number of columns in the block
196: Level: intermediate
198: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
200: @*/
201: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
202: {
207: if (brows != PETSC_DEFAULT) matfd->brows = brows;
208: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
209: return(0);
210: }
212: /*@
213: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
215: Collective on Mat
217: Input Parameters:
218: + mat - the matrix containing the nonzero structure of the Jacobian
219: . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
220: - color - the matrix coloring context
222: Level: beginner
224: Notes: When the coloring type is IS_COLORING_LOCAL the coloring is in the local ordering of the unknowns.
226: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
227: @*/
228: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
229: {
231: PetscBool eq;
236: if (color->setupcalled) return(0);
237: PetscObjectCompareId((PetscObject)mat,color->matid,&eq);
238: if (!eq) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Matrix used with MatFDColoringSetUp() must be that used with MatFDColoringCreate()");
240: PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
241: if (mat->ops->fdcoloringsetup) {
242: (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
243: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
245: color->setupcalled = PETSC_TRUE;
246: PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
247: return(0);
248: }
250: /*@C
251: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
253: Not Collective
255: Input Parameters:
256: . coloring - the coloring context
258: Output Parameters:
259: + f - the function
260: - fctx - the optional user-defined function context
262: Level: intermediate
264: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
266: @*/
267: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
268: {
271: if (f) *f = matfd->f;
272: if (fctx) *fctx = matfd->fctx;
273: return(0);
274: }
276: /*@C
277: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
279: Logically Collective on MatFDColoring
281: Input Parameters:
282: + coloring - the coloring context
283: . f - the function
284: - fctx - the optional user-defined function context
286: Calling sequence of (*f) function:
287: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
288: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
290: Level: advanced
292: Notes:
293: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
294: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
295: calling MatFDColoringApply()
297: Fortran Notes:
298: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
299: be used without SNES or within the SNES solvers.
301: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
303: @*/
304: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
305: {
308: matfd->f = f;
309: matfd->fctx = fctx;
310: return(0);
311: }
313: /*@
314: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
315: the options database.
317: Collective on MatFDColoring
319: The Jacobian, F'(u), is estimated with the differencing approximation
320: .vb
321: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
322: h = error_rel*u[i] if abs(u[i]) > umin
323: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
324: dx_{i} = (0, ... 1, .... 0)
325: .ve
327: Input Parameter:
328: . coloring - the coloring context
330: Options Database Keys:
331: + -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
332: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
333: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
334: . -mat_fd_coloring_view - Activates basic viewing
335: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
336: - -mat_fd_coloring_view draw - Activates drawing
338: Level: intermediate
340: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
342: @*/
343: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
344: {
346: PetscBool flg;
347: char value[3];
352: PetscObjectOptionsBegin((PetscObject)matfd);
353: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
354: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
355: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
356: if (flg) {
357: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
358: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
359: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
360: }
361: PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
362: PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
363: if (flg && matfd->bcols > matfd->ncolors) {
364: /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
365: matfd->bcols = matfd->ncolors;
366: }
368: /* process any options handlers added with PetscObjectAddOptionsHandler() */
369: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)matfd);
370: PetscOptionsEnd();
371: return(0);
372: }
374: /*@C
375: MatFDColoringSetType - Sets the approach for computing the finite difference parameter
377: Collective on MatFDColoring
379: Input Parameters:
380: + coloring - the coloring context
381: - type - either MATMFFD_WP or MATMFFD_DS
383: Options Database Keys:
384: . -mat_fd_type - "wp" or "ds"
386: Note: It is goofy that the argument type is MatMFFDType since the MatFDColoring actually computes the matrix entries
387: but the process of computing the entries is the same as as with the MatMFFD operation so we should reuse the names instead of
388: introducing another one.
390: Level: intermediate
392: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
394: @*/
395: PetscErrorCode MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)
396: {
399: /*
400: It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
401: and this function is being provided as patch for a release so it shouldn't change the implementaton
402: */
403: if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
404: else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
405: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",type);
406: return(0);
407: }
409: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
410: {
411: PetscErrorCode ierr;
412: PetscBool flg;
413: PetscViewer viewer;
414: PetscViewerFormat format;
417: if (prefix) {
418: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,prefix,optionname,&viewer,&format,&flg);
419: } else {
420: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
421: }
422: if (flg) {
423: PetscViewerPushFormat(viewer,format);
424: MatFDColoringView(fd,viewer);
425: PetscViewerPopFormat(viewer);
426: PetscViewerDestroy(&viewer);
427: }
428: return(0);
429: }
431: /*@
432: MatFDColoringCreate - Creates a matrix coloring context for finite difference
433: computation of Jacobians.
435: Collective on Mat
437: Input Parameters:
438: + mat - the matrix containing the nonzero structure of the Jacobian
439: - iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()
441: Output Parameter:
442: . color - the new coloring context
444: Level: intermediate
446: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
447: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
448: MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring()
449: @*/
450: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
451: {
452: MatFDColoring c;
453: MPI_Comm comm;
455: PetscInt M,N;
459: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
460: PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
461: MatGetSize(mat,&M,&N);
462: if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
463: PetscObjectGetComm((PetscObject)mat,&comm);
464: PetscHeaderCreate(c,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);
466: c->ctype = iscoloring->ctype;
467: PetscObjectGetId((PetscObject)mat,&c->matid);
469: if (mat->ops->fdcoloringcreate) {
470: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
471: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
473: MatCreateVecs(mat,NULL,&c->w1);
474: /* Vec is used instensively in particular piece of scalar CPU code; won't benifit from bouncing back and forth to the GPU */
475: VecBindToCPU(c->w1,PETSC_TRUE);
476: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
477: VecDuplicate(c->w1,&c->w2);
478: /* Vec is used instensively in particular piece of scalar CPU code; won't benifit from bouncing back and forth to the GPU */
479: VecBindToCPU(c->w2,PETSC_TRUE);
480: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);
482: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
483: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
484: c->currentcolor = -1;
485: c->htype = "wp";
486: c->fset = PETSC_FALSE;
487: c->setupcalled = PETSC_FALSE;
489: *color = c;
490: PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
491: PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
492: return(0);
493: }
495: /*@
496: MatFDColoringDestroy - Destroys a matrix coloring context that was created
497: via MatFDColoringCreate().
499: Collective on MatFDColoring
501: Input Parameter:
502: . c - coloring context
504: Level: intermediate
506: .seealso: MatFDColoringCreate()
507: @*/
508: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
509: {
511: PetscInt i;
512: MatFDColoring color = *c;
515: if (!*c) return(0);
516: if (--((PetscObject)color)->refct > 0) {*c = 0; return(0);}
518: /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
519: for (i=0; i<color->ncolors; i++) {
520: ISDestroy(&color->isa[i]);
521: }
522: PetscFree(color->isa);
523: PetscFree2(color->ncolumns,color->columns);
524: PetscFree(color->nrows);
525: if (color->htype[0] == 'w') {
526: PetscFree(color->matentry2);
527: } else {
528: PetscFree(color->matentry);
529: }
530: PetscFree(color->dy);
531: if (color->vscale) {VecDestroy(&color->vscale);}
532: VecDestroy(&color->w1);
533: VecDestroy(&color->w2);
534: VecDestroy(&color->w3);
535: PetscHeaderDestroy(c);
536: return(0);
537: }
539: /*@C
540: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
541: that are currently being perturbed.
543: Not Collective
545: Input Parameters:
546: . coloring - coloring context created with MatFDColoringCreate()
548: Output Parameters:
549: + n - the number of local columns being perturbed
550: - cols - the column indices, in global numbering
552: Level: advanced
554: Note: IF the matrix type is BAIJ, then the block column indices are returned
556: Fortran Note:
557: This routine has a different interface for Fortran
558: #include <petsc/finclude/petscmat.h>
559: $ use petscmat
560: $ PetscInt, pointer :: array(:)
561: $ PetscErrorCode ierr
562: $ MatFDColoring i
563: $ call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
564: $ use the entries of array ...
565: $ call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)
567: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
569: @*/
570: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,const PetscInt *cols[])
571: {
573: if (coloring->currentcolor >= 0) {
574: *n = coloring->ncolumns[coloring->currentcolor];
575: *cols = coloring->columns[coloring->currentcolor];
576: } else {
577: *n = 0;
578: }
579: return(0);
580: }
582: /*@
583: MatFDColoringApply - Given a matrix for which a MatFDColoring context
584: has been created, computes the Jacobian for a function via finite differences.
586: Collective on MatFDColoring
588: Input Parameters:
589: + mat - location to store Jacobian
590: . coloring - coloring context created with MatFDColoringCreate()
591: . x1 - location at which Jacobian is to be computed
592: - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
594: Options Database Keys:
595: + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
596: . -mat_fd_coloring_view - Activates basic viewing or coloring
597: . -mat_fd_coloring_view draw - Activates drawing of coloring
598: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
600: Level: intermediate
602: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()
604: @*/
605: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
606: {
608: PetscBool eq;
614: PetscObjectCompareId((PetscObject)J,coloring->matid,&eq);
615: if (!eq) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Matrix used with MatFDColoringApply() must be that used with MatFDColoringCreate()");
616: if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
617: if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
618: if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");
620: MatSetUnfactored(J);
621: PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
622: (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
623: PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
624: if (!coloring->viewed) {
625: MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");
626: coloring->viewed = PETSC_TRUE;
627: }
628: return(0);
629: }