Actual source code: fdmatrix.c
petsc-3.11.4 2019-09-28
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: .keywords: Mat, finite differences, coloring, view
97: @*/
98: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
99: {
100: PetscErrorCode ierr;
101: PetscInt i,j;
102: PetscBool isdraw,iascii;
103: PetscViewerFormat format;
107: if (!viewer) {
108: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
109: }
113: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
114: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
115: if (isdraw) {
116: MatFDColoringView_Draw(c,viewer);
117: } else if (iascii) {
118: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
119: PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);
120: PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);
121: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
123: PetscViewerGetFormat(viewer,&format);
124: if (format != PETSC_VIEWER_ASCII_INFO) {
125: PetscInt row,col,nz;
126: nz = 0;
127: for (i=0; i<c->ncolors; i++) {
128: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
129: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
130: for (j=0; j<c->ncolumns[i]; j++) {
131: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
132: }
133: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
134: if (c->matentry) {
135: for (j=0; j<c->nrows[i]; j++) {
136: row = c->matentry[nz].row;
137: col = c->matentry[nz++].col;
138: PetscViewerASCIIPrintf(viewer," %D %D \n",row,col);
139: }
140: }
141: }
142: }
143: PetscViewerFlush(viewer);
144: }
145: return(0);
146: }
148: /*@
149: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
150: a Jacobian matrix using finite differences.
152: Logically Collective on MatFDColoring
154: The Jacobian is estimated with the differencing approximation
155: .vb
156: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
157: htype = 'ds':
158: h = error_rel*u[i] if abs(u[i]) > umin
159: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
160: dx_{i} = (0, ... 1, .... 0)
162: htype = 'wp':
163: h = error_rel * sqrt(1 + ||u||)
164: .ve
166: Input Parameters:
167: + coloring - the coloring context
168: . error_rel - relative error
169: - umin - minimum allowable u-value magnitude
171: Level: advanced
173: .keywords: Mat, finite differences, coloring, set, parameters
175: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
177: @*/
178: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
179: {
184: if (error != PETSC_DEFAULT) matfd->error_rel = error;
185: if (umin != PETSC_DEFAULT) matfd->umin = umin;
186: return(0);
187: }
189: /*@
190: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
192: Logically Collective on MatFDColoring
194: Input Parameters:
195: + coloring - the coloring context
196: . brows - number of rows in the block
197: - bcols - number of columns in the block
199: Level: intermediate
201: .keywords: Mat, coloring
203: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
205: @*/
206: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
207: {
212: if (brows != PETSC_DEFAULT) matfd->brows = brows;
213: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
214: return(0);
215: }
217: /*@
218: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
220: Collective on Mat
222: Input Parameters:
223: + mat - the matrix containing the nonzero structure of the Jacobian
224: . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
225: - color - the matrix coloring context
227: Level: beginner
229: .keywords: MatFDColoring, setup
231: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
232: @*/
233: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
234: {
240: if (color->setupcalled) return(0);
242: PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
243: if (mat->ops->fdcoloringsetup) {
244: (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
245: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
247: color->setupcalled = PETSC_TRUE;
248: PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
249: return(0);
250: }
252: /*@C
253: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
255: Not Collective
257: Input Parameters:
258: . coloring - the coloring context
260: Output Parameters:
261: + f - the function
262: - fctx - the optional user-defined function context
264: Level: intermediate
266: .keywords: Mat, Jacobian, finite differences, set, function
268: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
270: @*/
271: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
272: {
275: if (f) *f = matfd->f;
276: if (fctx) *fctx = matfd->fctx;
277: return(0);
278: }
280: /*@C
281: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
283: Logically Collective on MatFDColoring
285: Input Parameters:
286: + coloring - the coloring context
287: . f - the function
288: - fctx - the optional user-defined function context
290: Calling sequence of (*f) function:
291: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
292: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
294: Level: advanced
296: Notes:
297: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
298: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
299: calling MatFDColoringApply()
301: Fortran Notes:
302: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
303: be used without SNES or within the SNES solvers.
305: .keywords: Mat, Jacobian, finite differences, set, function
307: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
309: @*/
310: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
311: {
314: matfd->f = f;
315: matfd->fctx = fctx;
316: return(0);
317: }
319: /*@
320: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
321: the options database.
323: Collective on MatFDColoring
325: The Jacobian, F'(u), is estimated with the differencing approximation
326: .vb
327: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
328: h = error_rel*u[i] if abs(u[i]) > umin
329: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
330: dx_{i} = (0, ... 1, .... 0)
331: .ve
333: Input Parameter:
334: . coloring - the coloring context
336: Options Database Keys:
337: + -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
338: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
339: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
340: . -mat_fd_coloring_view - Activates basic viewing
341: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
342: - -mat_fd_coloring_view draw - Activates drawing
344: Level: intermediate
346: .keywords: Mat, finite differences, parameters
348: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
350: @*/
351: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
352: {
354: PetscBool flg;
355: char value[3];
360: PetscObjectOptionsBegin((PetscObject)matfd);
361: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
362: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
363: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
364: if (flg) {
365: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
366: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
367: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
368: }
369: PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
370: PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
371: if (flg && matfd->bcols > matfd->ncolors) {
372: /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
373: matfd->bcols = matfd->ncolors;
374: }
376: /* process any options handlers added with PetscObjectAddOptionsHandler() */
377: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)matfd);
378: PetscOptionsEnd();
379: return(0);
380: }
382: /*@C
383: MatFDColoringSetType - Sets the approach for computing the finite difference parameter
385: Collective on MatFDColoring
387: Input Parameters:
388: + coloring - the coloring context
389: - type - either MATMFFD_WP or MATMFFD_DS
391: Options Database Keys:
392: . -mat_fd_type - "wp" or "ds"
394: Note: It is goofy that the argument type is MatMFFDType since the MatFDColoring actually computes the matrix entries
395: but the process of computing the entries is the same as as with the MatMFFD operation so we should reuse the names instead of
396: introducing another one.
398: Level: intermediate
400: .keywords: Mat, finite differences, parameters
402: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
404: @*/
405: PetscErrorCode MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)
406: {
409: /*
410: It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
411: and this function is being provided as patch for a release so it shouldn't change the implementaton
412: */
413: if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
414: else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
415: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",type);
416: return(0);
417: }
419: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
420: {
421: PetscErrorCode ierr;
422: PetscBool flg;
423: PetscViewer viewer;
424: PetscViewerFormat format;
427: if (prefix) {
428: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,prefix,optionname,&viewer,&format,&flg);
429: } else {
430: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
431: }
432: if (flg) {
433: PetscViewerPushFormat(viewer,format);
434: MatFDColoringView(fd,viewer);
435: PetscViewerPopFormat(viewer);
436: PetscViewerDestroy(&viewer);
437: }
438: return(0);
439: }
441: /*@
442: MatFDColoringCreate - Creates a matrix coloring context for finite difference
443: computation of Jacobians.
445: Collective on Mat
447: Input Parameters:
448: + mat - the matrix containing the nonzero structure of the Jacobian
449: - iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()
451: Output Parameter:
452: . color - the new coloring context
454: Level: intermediate
456: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
457: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
458: MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring()
459: @*/
460: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
461: {
462: MatFDColoring c;
463: MPI_Comm comm;
465: PetscInt M,N;
469: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
470: PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
471: MatGetSize(mat,&M,&N);
472: if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
473: PetscObjectGetComm((PetscObject)mat,&comm);
474: PetscHeaderCreate(c,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);
476: c->ctype = iscoloring->ctype;
478: if (mat->ops->fdcoloringcreate) {
479: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
480: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
482: MatCreateVecs(mat,NULL,&c->w1);
483: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
484: VecDuplicate(c->w1,&c->w2);
485: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);
487: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
488: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
489: c->currentcolor = -1;
490: c->htype = "wp";
491: c->fset = PETSC_FALSE;
492: c->setupcalled = PETSC_FALSE;
494: *color = c;
495: PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
496: PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
497: return(0);
498: }
500: /*@
501: MatFDColoringDestroy - Destroys a matrix coloring context that was created
502: via MatFDColoringCreate().
504: Collective on MatFDColoring
506: Input Parameter:
507: . c - coloring context
509: Level: intermediate
511: .seealso: MatFDColoringCreate()
512: @*/
513: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
514: {
516: PetscInt i;
517: MatFDColoring color = *c;
520: if (!*c) return(0);
521: if (--((PetscObject)color)->refct > 0) {*c = 0; return(0);}
523: for (i=0; i<color->ncolors; i++) {
524: PetscFree(color->columns[i]);
525: }
526: PetscFree(color->ncolumns);
527: PetscFree(color->columns);
528: PetscFree(color->nrows);
529: if (color->htype[0] == 'w') {
530: PetscFree(color->matentry2);
531: } else {
532: PetscFree(color->matentry);
533: }
534: PetscFree(color->dy);
535: if (color->vscale) {VecDestroy(&color->vscale);}
536: VecDestroy(&color->w1);
537: VecDestroy(&color->w2);
538: VecDestroy(&color->w3);
539: PetscHeaderDestroy(c);
540: return(0);
541: }
543: /*@C
544: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
545: that are currently being perturbed.
547: Not Collective
549: Input Parameters:
550: . coloring - coloring context created with MatFDColoringCreate()
552: Output Parameters:
553: + n - the number of local columns being perturbed
554: - cols - the column indices, in global numbering
556: Level: advanced
558: Fortran Note:
559: This routine has a different interface for Fortran
560: $ #include <petsc/finclude/petscmat.h>
561: $ use petscmat
562: $ PetscInt, pointer :: array(:)
563: $ PetscErrorCode ierr
564: $ MatFDColoring i
565: $ call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
566: $ use the entries of array ...
567: $ call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)
569: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
571: .keywords: coloring, Jacobian, finite differences
572: @*/
573: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,const PetscInt *cols[])
574: {
576: if (coloring->currentcolor >= 0) {
577: *n = coloring->ncolumns[coloring->currentcolor];
578: *cols = coloring->columns[coloring->currentcolor];
579: } else {
580: *n = 0;
581: }
582: return(0);
583: }
585: /*@
586: MatFDColoringApply - Given a matrix for which a MatFDColoring context
587: has been created, computes the Jacobian for a function via finite differences.
589: Collective on MatFDColoring
591: Input Parameters:
592: + mat - location to store Jacobian
593: . coloring - coloring context created with MatFDColoringCreate()
594: . x1 - location at which Jacobian is to be computed
595: - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
597: Options Database Keys:
598: + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
599: . -mat_fd_coloring_view - Activates basic viewing or coloring
600: . -mat_fd_coloring_view draw - Activates drawing of coloring
601: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
603: Level: intermediate
605: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()
607: .keywords: coloring, Jacobian, finite differences
608: @*/
609: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
610: {
612: PetscBool flg = PETSC_FALSE;
618: if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
619: if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
620: if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");
622: MatSetUnfactored(J);
623: PetscOptionsGetBool(((PetscObject)coloring)->options,NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);
624: if (flg) {
625: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
626: } else {
627: PetscBool assembled;
628: MatAssembled(J,&assembled);
629: if (assembled) {
630: MatZeroEntries(J);
631: }
632: }
634: PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
635: (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
636: PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
637: if (!coloring->viewed) {
638: MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");
639: coloring->viewed = PETSC_TRUE;
640: }
641: return(0);
642: }