Actual source code: fdmatrix.c
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: {
12: if (F) {
13: VecCopy(F,fd->w1);
14: fd->fset = PETSC_TRUE;
15: } else {
16: fd->fset = PETSC_FALSE;
17: }
18: return 0;
19: }
21: #include <petscdraw.h>
22: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
23: {
24: MatFDColoring fd = (MatFDColoring)Aa;
25: PetscInt i,j,nz,row;
26: PetscReal x,y;
27: MatEntry *Jentry=fd->matentry;
29: /* loop over colors */
30: nz = 0;
31: for (i=0; i<fd->ncolors; i++) {
32: for (j=0; j<fd->nrows[i]; j++) {
33: row = Jentry[nz].row;
34: y = fd->M - row - fd->rstart;
35: x = (PetscReal)Jentry[nz++].col;
36: PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
37: }
38: }
39: return 0;
40: }
42: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
43: {
44: PetscBool isnull;
45: PetscDraw draw;
46: PetscReal xr,yr,xl,yl,h,w;
48: PetscViewerDrawGetDraw(viewer,0,&draw);
49: PetscDrawIsNull(draw,&isnull);
50: if (isnull) return 0;
52: xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0;
53: xr += w; yr += h; xl = -w; yl = -h;
54: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
55: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
56: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
57: PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
58: PetscDrawSave(draw);
59: return 0;
60: }
62: /*@C
63: MatFDColoringView - Views a finite difference coloring context.
65: Collective on MatFDColoring
67: Input Parameters:
68: + c - the coloring context
69: - viewer - visualization context
71: Level: intermediate
73: Notes:
74: The available visualization contexts include
75: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
76: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
77: output where only the first processor opens
78: the file. All other processors send their
79: data to the first processor to print.
80: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
82: Notes:
83: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
84: involves more than 33 then some seemingly identical colors are displayed making it look
85: like an illegal coloring. This is just a graphical artifact.
87: .seealso: MatFDColoringCreate()
89: @*/
90: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
91: {
92: PetscInt i,j;
93: PetscBool isdraw,iascii;
94: PetscViewerFormat format;
97: if (!viewer) {
98: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
99: }
103: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
104: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
105: if (isdraw) {
106: MatFDColoringView_Draw(c,viewer);
107: } else if (iascii) {
108: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
109: PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);
110: PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);
111: PetscViewerASCIIPrintf(viewer," Number of colors=%" PetscInt_FMT "\n",c->ncolors);
113: PetscViewerGetFormat(viewer,&format);
114: if (format != PETSC_VIEWER_ASCII_INFO) {
115: PetscInt row,col,nz;
116: nz = 0;
117: for (i=0; i<c->ncolors; i++) {
118: PetscViewerASCIIPrintf(viewer," Information for color %" PetscInt_FMT "\n",i);
119: PetscViewerASCIIPrintf(viewer," Number of columns %" PetscInt_FMT "\n",c->ncolumns[i]);
120: for (j=0; j<c->ncolumns[i]; j++) {
121: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",c->columns[i][j]);
122: }
123: PetscViewerASCIIPrintf(viewer," Number of rows %" PetscInt_FMT "\n",c->nrows[i]);
124: if (c->matentry) {
125: for (j=0; j<c->nrows[i]; j++) {
126: row = c->matentry[nz].row;
127: col = c->matentry[nz++].col;
128: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " \n",row,col);
129: }
130: }
131: }
132: }
133: PetscViewerFlush(viewer);
134: }
135: return 0;
136: }
138: /*@
139: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
140: a Jacobian matrix using finite differences.
142: Logically Collective on MatFDColoring
144: The Jacobian is estimated with the differencing approximation
145: .vb
146: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
147: htype = 'ds':
148: h = error_rel*u[i] if abs(u[i]) > umin
149: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
150: dx_{i} = (0, ... 1, .... 0)
152: htype = 'wp':
153: h = error_rel * sqrt(1 + ||u||)
154: .ve
156: Input Parameters:
157: + coloring - the coloring context
158: . error_rel - relative error
159: - umin - minimum allowable u-value magnitude
161: Level: advanced
163: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
165: @*/
166: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
167: {
171: if (error != PETSC_DEFAULT) matfd->error_rel = error;
172: if (umin != PETSC_DEFAULT) matfd->umin = umin;
173: return 0;
174: }
176: /*@
177: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
179: Logically Collective on MatFDColoring
181: Input Parameters:
182: + coloring - the coloring context
183: . brows - number of rows in the block
184: - bcols - number of columns in the block
186: Level: intermediate
188: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
190: @*/
191: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
192: {
196: if (brows != PETSC_DEFAULT) matfd->brows = brows;
197: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
198: return 0;
199: }
201: /*@
202: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
204: Collective on Mat
206: Input Parameters:
207: + mat - the matrix containing the nonzero structure of the Jacobian
208: . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
209: - color - the matrix coloring context
211: Level: beginner
213: Notes: When the coloring type is IS_COLORING_LOCAL the coloring is in the local ordering of the unknowns.
215: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
216: @*/
217: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
218: {
219: PetscBool eq;
223: if (color->setupcalled) return 0;
224: PetscObjectCompareId((PetscObject)mat,color->matid,&eq);
227: PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
228: if (mat->ops->fdcoloringsetup) {
229: (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
230: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
232: color->setupcalled = PETSC_TRUE;
233: PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
234: return 0;
235: }
237: /*@C
238: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
240: Not Collective
242: Input Parameter:
243: . coloring - the coloring context
245: Output Parameters:
246: + f - the function
247: - fctx - the optional user-defined function context
249: Level: intermediate
251: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
253: @*/
254: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
255: {
257: if (f) *f = matfd->f;
258: if (fctx) *fctx = matfd->fctx;
259: return 0;
260: }
262: /*@C
263: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
265: Logically Collective on MatFDColoring
267: Input Parameters:
268: + coloring - the coloring context
269: . f - the function
270: - fctx - the optional user-defined function context
272: Calling sequence of (*f) function:
273: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
274: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
276: Level: advanced
278: Notes:
279: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
280: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
281: calling MatFDColoringApply()
283: Fortran Notes:
284: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
285: be used without SNES or within the SNES solvers.
287: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
289: @*/
290: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
291: {
293: matfd->f = f;
294: matfd->fctx = fctx;
295: return 0;
296: }
298: /*@
299: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
300: the options database.
302: Collective on MatFDColoring
304: The Jacobian, F'(u), is estimated with the differencing approximation
305: .vb
306: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
307: h = error_rel*u[i] if abs(u[i]) > umin
308: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
309: dx_{i} = (0, ... 1, .... 0)
310: .ve
312: Input Parameter:
313: . coloring - the coloring context
315: Options Database Keys:
316: + -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
317: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
318: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
319: . -mat_fd_coloring_view - Activates basic viewing
320: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
321: - -mat_fd_coloring_view draw - Activates drawing
323: Level: intermediate
325: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
327: @*/
328: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
329: {
331: PetscBool flg;
332: char value[3];
336: PetscObjectOptionsBegin((PetscObject)matfd);
337: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,NULL);
338: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,NULL);
339: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,sizeof(value),&flg);
340: if (flg) {
341: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
342: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
343: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
344: }
345: PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
346: PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
347: if (flg && matfd->bcols > matfd->ncolors) {
348: /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
349: matfd->bcols = matfd->ncolors;
350: }
352: /* process any options handlers added with PetscObjectAddOptionsHandler() */
353: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)matfd);
354: PetscOptionsEnd();
355: return 0;
356: }
358: /*@C
359: MatFDColoringSetType - Sets the approach for computing the finite difference parameter
361: Collective on MatFDColoring
363: Input Parameters:
364: + coloring - the coloring context
365: - type - either MATMFFD_WP or MATMFFD_DS
367: Options Database Keys:
368: . -mat_fd_type - "wp" or "ds"
370: Note: It is goofy that the argument type is MatMFFDType since the MatFDColoring actually computes the matrix entries
371: but the process of computing the entries is the same as as with the MatMFFD operation so we should reuse the names instead of
372: introducing another one.
374: Level: intermediate
376: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
378: @*/
379: PetscErrorCode MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)
380: {
382: /*
383: It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
384: and this function is being provided as patch for a release so it shouldn't change the implementaton
385: */
386: if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
387: else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
388: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",type);
389: return 0;
390: }
392: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
393: {
394: PetscBool flg;
395: PetscViewer viewer;
396: PetscViewerFormat format;
398: if (prefix) {
399: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,prefix,optionname,&viewer,&format,&flg);
400: } else {
401: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
402: }
403: if (flg) {
404: PetscViewerPushFormat(viewer,format);
405: MatFDColoringView(fd,viewer);
406: PetscViewerPopFormat(viewer);
407: PetscViewerDestroy(&viewer);
408: }
409: return 0;
410: }
412: /*@
413: MatFDColoringCreate - Creates a matrix coloring context for finite difference
414: computation of Jacobians.
416: Collective on Mat
418: Input Parameters:
419: + mat - the matrix containing the nonzero structure of the Jacobian
420: - iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()
422: Output Parameter:
423: . color - the new coloring context
425: Level: intermediate
427: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
428: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
429: MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring(), MatFDColoringSetValues()
430: @*/
431: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
432: {
433: MatFDColoring c;
434: MPI_Comm comm;
435: PetscInt M,N;
439: PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
440: MatGetSize(mat,&M,&N);
442: PetscObjectGetComm((PetscObject)mat,&comm);
443: PetscHeaderCreate(c,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);
445: c->ctype = iscoloring->ctype;
446: PetscObjectGetId((PetscObject)mat,&c->matid);
448: if (mat->ops->fdcoloringcreate) {
449: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
450: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
452: MatCreateVecs(mat,NULL,&c->w1);
453: /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
454: VecBindToCPU(c->w1,PETSC_TRUE);
455: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
456: VecDuplicate(c->w1,&c->w2);
457: /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
458: VecBindToCPU(c->w2,PETSC_TRUE);
459: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);
461: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
462: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
463: c->currentcolor = -1;
464: c->htype = "wp";
465: c->fset = PETSC_FALSE;
466: c->setupcalled = PETSC_FALSE;
468: *color = c;
469: PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
470: PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
471: return 0;
472: }
474: /*@
475: MatFDColoringDestroy - Destroys a matrix coloring context that was created
476: via MatFDColoringCreate().
478: Collective on MatFDColoring
480: Input Parameter:
481: . c - coloring context
483: Level: intermediate
485: .seealso: MatFDColoringCreate()
486: @*/
487: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
488: {
489: PetscInt i;
490: MatFDColoring color = *c;
492: if (!*c) return 0;
493: if (--((PetscObject)color)->refct > 0) {*c = NULL; return 0;}
495: /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
496: for (i=0; i<color->ncolors; i++) {
497: ISDestroy(&color->isa[i]);
498: }
499: PetscFree(color->isa);
500: PetscFree2(color->ncolumns,color->columns);
501: PetscFree(color->nrows);
502: if (color->htype[0] == 'w') {
503: PetscFree(color->matentry2);
504: } else {
505: PetscFree(color->matentry);
506: }
507: PetscFree(color->dy);
508: if (color->vscale) VecDestroy(&color->vscale);
509: VecDestroy(&color->w1);
510: VecDestroy(&color->w2);
511: VecDestroy(&color->w3);
512: PetscHeaderDestroy(c);
513: return 0;
514: }
516: /*@C
517: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
518: that are currently being perturbed.
520: Not Collective
522: Input Parameter:
523: . coloring - coloring context created with MatFDColoringCreate()
525: Output Parameters:
526: + n - the number of local columns being perturbed
527: - cols - the column indices, in global numbering
529: Level: advanced
531: Note: IF the matrix type is BAIJ, then the block column indices are returned
533: Fortran Note:
534: This routine has a different interface for Fortran
535: #include <petsc/finclude/petscmat.h>
536: $ use petscmat
537: $ PetscInt, pointer :: array(:)
538: $ PetscErrorCode ierr
539: $ MatFDColoring i
540: $ call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
541: $ use the entries of array ...
542: $ call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)
544: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
546: @*/
547: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,const PetscInt *cols[])
548: {
549: if (coloring->currentcolor >= 0) {
550: *n = coloring->ncolumns[coloring->currentcolor];
551: *cols = coloring->columns[coloring->currentcolor];
552: } else {
553: *n = 0;
554: }
555: return 0;
556: }
558: /*@
559: MatFDColoringApply - Given a matrix for which a MatFDColoring context
560: has been created, computes the Jacobian for a function via finite differences.
562: Collective on MatFDColoring
564: Input Parameters:
565: + mat - location to store Jacobian
566: . coloring - coloring context created with MatFDColoringCreate()
567: . x1 - location at which Jacobian is to be computed
568: - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
570: Options Database Keys:
571: + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
572: . -mat_fd_coloring_view - Activates basic viewing or coloring
573: . -mat_fd_coloring_view draw - Activates drawing of coloring
574: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
576: Level: intermediate
578: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction(), MatFDColoringSetValues()
580: @*/
581: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
582: {
583: PetscBool eq;
588: PetscObjectCompareId((PetscObject)J,coloring->matid,&eq);
594: MatSetUnfactored(J);
595: PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
596: (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
597: PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
598: if (!coloring->viewed) {
599: MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");
600: coloring->viewed = PETSC_TRUE;
601: }
602: return 0;
603: }