Actual source code: fdmatrix.c
petsc-3.7.5 2017-01-01
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> /*I "petscmat.h" I*/
8: #include <petsc/private/isimpl.h>
12: PetscErrorCode MatFDColoringSetF(MatFDColoring fd,Vec F)
13: {
17: if (F) {
18: VecCopy(F,fd->w1);
19: fd->fset = PETSC_TRUE;
20: } else {
21: fd->fset = PETSC_FALSE;
22: }
23: return(0);
24: }
26: #include <petscdraw.h>
29: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
30: {
31: MatFDColoring fd = (MatFDColoring)Aa;
33: PetscInt i,j,nz,row;
34: PetscReal x,y;
35: MatEntry *Jentry=fd->matentry;
38: /* loop over colors */
39: nz = 0;
40: for (i=0; i<fd->ncolors; i++) {
41: for (j=0; j<fd->nrows[i]; j++) {
42: row = Jentry[nz].row;
43: y = fd->M - row - fd->rstart;
44: x = (PetscReal)Jentry[nz++].col;
45: PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
46: }
47: }
48: return(0);
49: }
53: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
54: {
56: PetscBool isnull;
57: PetscDraw draw;
58: PetscReal xr,yr,xl,yl,h,w;
61: PetscViewerDrawGetDraw(viewer,0,&draw);
62: PetscDrawIsNull(draw,&isnull);
63: if (isnull) return(0);
65: xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0;
66: xr += w; yr += h; xl = -w; yl = -h;
67: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
68: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
69: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
70: PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
71: PetscDrawSave(draw);
72: return(0);
73: }
77: /*@C
78: MatFDColoringView - Views a finite difference coloring context.
80: Collective on MatFDColoring
82: Input Parameters:
83: + c - the coloring context
84: - viewer - visualization context
86: Level: intermediate
88: Notes:
89: The available visualization contexts include
90: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
91: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
92: output where only the first processor opens
93: the file. All other processors send their
94: data to the first processor to print.
95: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
97: Notes:
98: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
99: involves more than 33 then some seemingly identical colors are displayed making it look
100: like an illegal coloring. This is just a graphical artifact.
102: .seealso: MatFDColoringCreate()
104: .keywords: Mat, finite differences, coloring, view
105: @*/
106: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
107: {
108: PetscErrorCode ierr;
109: PetscInt i,j;
110: PetscBool isdraw,iascii;
111: PetscViewerFormat format;
115: if (!viewer) {
116: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
117: }
121: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
122: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
123: if (isdraw) {
124: MatFDColoringView_Draw(c,viewer);
125: } else if (iascii) {
126: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
127: PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);
128: PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);
129: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
131: PetscViewerGetFormat(viewer,&format);
132: if (format != PETSC_VIEWER_ASCII_INFO) {
133: PetscInt row,col,nz;
134: nz = 0;
135: for (i=0; i<c->ncolors; i++) {
136: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
137: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
138: for (j=0; j<c->ncolumns[i]; j++) {
139: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
140: }
141: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
142: for (j=0; j<c->nrows[i]; j++) {
143: row = c->matentry[nz].row;
144: col = c->matentry[nz++].col;
145: PetscViewerASCIIPrintf(viewer," %D %D \n",row,col);
146: }
147: }
148: }
149: PetscViewerFlush(viewer);
150: }
151: return(0);
152: }
156: /*@
157: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
158: a Jacobian matrix using finite differences.
160: Logically Collective on MatFDColoring
162: The Jacobian is estimated with the differencing approximation
163: .vb
164: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
165: h = error_rel*u[i] if abs(u[i]) > umin
166: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
167: dx_{i} = (0, ... 1, .... 0)
168: .ve
170: Input Parameters:
171: + coloring - the coloring context
172: . error_rel - relative error
173: - umin - minimum allowable u-value magnitude
175: Level: advanced
177: .keywords: Mat, finite differences, coloring, set, parameters
179: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
181: @*/
182: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
183: {
188: if (error != PETSC_DEFAULT) matfd->error_rel = error;
189: if (umin != PETSC_DEFAULT) matfd->umin = umin;
190: return(0);
191: }
195: /*@
196: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
198: Logically Collective on MatFDColoring
200: Input Parameters:
201: + coloring - the coloring context
202: . brows - number of rows in the block
203: - bcols - number of columns in the block
205: Level: intermediate
207: .keywords: Mat, coloring
209: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
211: @*/
212: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
213: {
218: if (brows != PETSC_DEFAULT) matfd->brows = brows;
219: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
220: return(0);
221: }
225: /*@
226: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
228: Collective on Mat
230: Input Parameters:
231: + mat - the matrix containing the nonzero structure of the Jacobian
232: . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
233: - color - the matrix coloring context
235: Level: beginner
237: .keywords: MatFDColoring, setup
239: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
240: @*/
241: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
242: {
248: if (color->setupcalled) return(0);
250: PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
251: if (mat->ops->fdcoloringsetup) {
252: (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
253: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
255: color->setupcalled = PETSC_TRUE;
256: PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
257: return(0);
258: }
262: /*@C
263: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
265: Not Collective
267: Input Parameters:
268: . coloring - the coloring context
270: Output Parameters:
271: + f - the function
272: - fctx - the optional user-defined function context
274: Level: intermediate
276: .keywords: Mat, Jacobian, finite differences, set, function
278: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
280: @*/
281: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
282: {
285: if (f) *f = matfd->f;
286: if (fctx) *fctx = matfd->fctx;
287: return(0);
288: }
292: /*@C
293: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
295: Logically Collective on MatFDColoring
297: Input Parameters:
298: + coloring - the coloring context
299: . f - the function
300: - fctx - the optional user-defined function context
302: Calling sequence of (*f) function:
303: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
304: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
306: Level: advanced
308: Notes: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
309: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
310: calling MatFDColoringApply()
312: Fortran Notes:
313: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
314: be used without SNES or within the SNES solvers.
316: .keywords: Mat, Jacobian, finite differences, set, function
318: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
320: @*/
321: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
322: {
325: matfd->f = f;
326: matfd->fctx = fctx;
327: return(0);
328: }
332: /*@
333: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
334: the options database.
336: Collective on MatFDColoring
338: The Jacobian, F'(u), is estimated with the differencing approximation
339: .vb
340: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
341: h = error_rel*u[i] if abs(u[i]) > umin
342: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
343: dx_{i} = (0, ... 1, .... 0)
344: .ve
346: Input Parameter:
347: . coloring - the coloring context
349: Options Database Keys:
350: + -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
351: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
352: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
353: . -mat_fd_coloring_view - Activates basic viewing
354: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
355: - -mat_fd_coloring_view draw - Activates drawing
357: Level: intermediate
359: .keywords: Mat, finite differences, parameters
361: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
363: @*/
364: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
365: {
367: PetscBool flg;
368: char value[3];
373: PetscObjectOptionsBegin((PetscObject)matfd);
374: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
375: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
376: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
377: if (flg) {
378: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
379: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
380: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
381: }
382: PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
383: PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
384: if (flg && matfd->bcols > matfd->ncolors) {
385: /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
386: matfd->bcols = matfd->ncolors;
387: }
389: /* process any options handlers added with PetscObjectAddOptionsHandler() */
390: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)matfd);
391: PetscOptionsEnd();
392: return(0);
393: }
397: /*@C
398: MatFDColoringSetType - Sets the approach for computing the finite difference parameter
400: Collective on MatFDColoring
402: Input Parameters:
403: + coloring - the coloring context
404: - type - either MATMFFD_WP or MATMFFD_DS
406: Options Database Keys:
407: . -mat_fd_type - "wp" or "ds"
409: Note: It is goofy that the argument type is MatMFFDType since the MatFDColoring actually computes the matrix entries
410: but the process of computing the entries is the same as as with the MatMFFD operation so we should reuse the names instead of
411: introducing another one.
413: Level: intermediate
415: .keywords: Mat, finite differences, parameters
417: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
419: @*/
420: PetscErrorCode MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)
421: {
424: /*
425: It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
426: and this function is being provided as patch for a release so it shouldn't change the implementaton
427: */
428: if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
429: else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
430: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",type);
431: return(0);
432: }
436: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
437: {
438: PetscErrorCode ierr;
439: PetscBool flg;
440: PetscViewer viewer;
441: PetscViewerFormat format;
444: if (prefix) {
445: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),prefix,optionname,&viewer,&format,&flg);
446: } else {
447: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
448: }
449: if (flg) {
450: PetscViewerPushFormat(viewer,format);
451: MatFDColoringView(fd,viewer);
452: PetscViewerPopFormat(viewer);
453: PetscViewerDestroy(&viewer);
454: }
455: return(0);
456: }
460: /*@
461: MatFDColoringCreate - Creates a matrix coloring context for finite difference
462: computation of Jacobians.
464: Collective on Mat
466: Input Parameters:
467: + mat - the matrix containing the nonzero structure of the Jacobian
468: - iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()
470: Output Parameter:
471: . color - the new coloring context
473: Level: intermediate
475: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
476: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
477: MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring()
478: @*/
479: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
480: {
481: MatFDColoring c;
482: MPI_Comm comm;
484: PetscInt M,N;
488: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
489: PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
490: MatGetSize(mat,&M,&N);
491: if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
492: PetscObjectGetComm((PetscObject)mat,&comm);
493: PetscHeaderCreate(c,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);
495: c->ctype = iscoloring->ctype;
497: if (mat->ops->fdcoloringcreate) {
498: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
499: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
501: MatCreateVecs(mat,NULL,&c->w1);
502: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
503: VecDuplicate(c->w1,&c->w2);
504: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);
506: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
507: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
508: c->currentcolor = -1;
509: c->htype = "wp";
510: c->fset = PETSC_FALSE;
511: c->setupcalled = PETSC_FALSE;
513: *color = c;
514: PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
515: PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
516: return(0);
517: }
521: /*@
522: MatFDColoringDestroy - Destroys a matrix coloring context that was created
523: via MatFDColoringCreate().
525: Collective on MatFDColoring
527: Input Parameter:
528: . c - coloring context
530: Level: intermediate
532: .seealso: MatFDColoringCreate()
533: @*/
534: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
535: {
537: PetscInt i;
538: MatFDColoring color = *c;
541: if (!*c) return(0);
542: if (--((PetscObject)color)->refct > 0) {*c = 0; return(0);}
544: for (i=0; i<color->ncolors; i++) {
545: PetscFree(color->columns[i]);
546: }
547: PetscFree(color->ncolumns);
548: PetscFree(color->columns);
549: PetscFree(color->nrows);
550: if (color->htype[0] == 'w') {
551: PetscFree(color->matentry2);
552: } else {
553: PetscFree(color->matentry);
554: }
555: PetscFree(color->dy);
556: if (color->vscale) {VecDestroy(&color->vscale);}
557: VecDestroy(&color->w1);
558: VecDestroy(&color->w2);
559: VecDestroy(&color->w3);
560: PetscHeaderDestroy(c);
561: return(0);
562: }
566: /*@C
567: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
568: that are currently being perturbed.
570: Not Collective
572: Input Parameters:
573: . coloring - coloring context created with MatFDColoringCreate()
575: Output Parameters:
576: + n - the number of local columns being perturbed
577: - cols - the column indices, in global numbering
579: Level: intermediate
581: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
583: .keywords: coloring, Jacobian, finite differences
584: @*/
585: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
586: {
588: if (coloring->currentcolor >= 0) {
589: *n = coloring->ncolumns[coloring->currentcolor];
590: *cols = coloring->columns[coloring->currentcolor];
591: } else {
592: *n = 0;
593: }
594: return(0);
595: }
599: /*@
600: MatFDColoringApply - Given a matrix for which a MatFDColoring context
601: has been created, computes the Jacobian for a function via finite differences.
603: Collective on MatFDColoring
605: Input Parameters:
606: + mat - location to store Jacobian
607: . coloring - coloring context created with MatFDColoringCreate()
608: . x1 - location at which Jacobian is to be computed
609: - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
611: Options Database Keys:
612: + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
613: . -mat_fd_coloring_view - Activates basic viewing or coloring
614: . -mat_fd_coloring_view draw - Activates drawing of coloring
615: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
617: Level: intermediate
619: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()
621: .keywords: coloring, Jacobian, finite differences
622: @*/
623: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
624: {
626: PetscBool flg = PETSC_FALSE;
632: if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
633: if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
634: if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");
636: MatSetUnfactored(J);
637: PetscOptionsGetBool(((PetscObject)coloring)->options,NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);
638: if (flg) {
639: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
640: } else {
641: PetscBool assembled;
642: MatAssembled(J,&assembled);
643: if (assembled) {
644: MatZeroEntries(J);
645: }
646: }
648: PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
649: (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
650: PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
651: return(0);
652: }