Actual source code: fdmatrix.c
petsc-3.5.4 2015-05-23
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*/
11: PetscErrorCode MatFDColoringSetF(MatFDColoring fd,Vec F)
12: {
16: if (F) {
17: VecCopy(F,fd->w1);
18: fd->fset = PETSC_TRUE;
19: } else {
20: fd->fset = PETSC_FALSE;
21: }
22: return(0);
23: }
25: #include <petscdraw.h>
28: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
29: {
30: MatFDColoring fd = (MatFDColoring)Aa;
32: PetscInt i,j,nz,row;
33: PetscReal x,y;
34: MatEntry *Jentry=fd->matentry;
37: /* loop over colors */
38: nz = 0;
39: for (i=0; i<fd->ncolors; i++) {
40: for (j=0; j<fd->nrows[i]; j++) {
41: row = Jentry[nz].row;
42: y = fd->M - row - fd->rstart;
43: x = (PetscReal)Jentry[nz++].col;
44: PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
45: }
46: }
47: return(0);
48: }
52: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
53: {
55: PetscBool isnull;
56: PetscDraw draw;
57: PetscReal xr,yr,xl,yl,h,w;
60: PetscViewerDrawGetDraw(viewer,0,&draw);
61: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
63: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
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: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
69: PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
70: return(0);
71: }
75: /*@C
76: MatFDColoringView - Views a finite difference coloring context.
78: Collective on MatFDColoring
80: Input Parameters:
81: + c - the coloring context
82: - viewer - visualization context
84: Level: intermediate
86: Notes:
87: The available visualization contexts include
88: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
89: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
90: output where only the first processor opens
91: the file. All other processors send their
92: data to the first processor to print.
93: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
95: Notes:
96: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
97: involves more than 33 then some seemingly identical colors are displayed making it look
98: like an illegal coloring. This is just a graphical artifact.
100: .seealso: MatFDColoringCreate()
102: .keywords: Mat, finite differences, coloring, view
103: @*/
104: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
105: {
106: PetscErrorCode ierr;
107: PetscInt i,j;
108: PetscBool isdraw,iascii;
109: PetscViewerFormat format;
113: if (!viewer) {
114: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
115: }
119: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
120: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
121: if (isdraw) {
122: MatFDColoringView_Draw(c,viewer);
123: } else if (iascii) {
124: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
125: PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);
126: PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);
127: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
129: PetscViewerGetFormat(viewer,&format);
130: if (format != PETSC_VIEWER_ASCII_INFO) {
131: PetscInt row,col,nz;
132: nz = 0;
133: for (i=0; i<c->ncolors; i++) {
134: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
135: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
136: for (j=0; j<c->ncolumns[i]; j++) {
137: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
138: }
139: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
140: for (j=0; j<c->nrows[i]; j++) {
141: row = c->matentry[nz].row;
142: col = c->matentry[nz++].col;
143: PetscViewerASCIIPrintf(viewer," %D %D \n",row,col);
144: }
145: }
146: }
147: PetscViewerFlush(viewer);
148: }
149: return(0);
150: }
154: /*@
155: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
156: a Jacobian matrix using finite differences.
158: Logically Collective on MatFDColoring
160: The Jacobian is estimated with the differencing approximation
161: .vb
162: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
163: h = error_rel*u[i] if abs(u[i]) > umin
164: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
165: dx_{i} = (0, ... 1, .... 0)
166: .ve
168: Input Parameters:
169: + coloring - the coloring context
170: . error_rel - relative error
171: - umin - minimum allowable u-value magnitude
173: Level: advanced
175: .keywords: Mat, finite differences, coloring, set, parameters
177: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
179: @*/
180: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
181: {
186: if (error != PETSC_DEFAULT) matfd->error_rel = error;
187: if (umin != PETSC_DEFAULT) matfd->umin = umin;
188: return(0);
189: }
193: /*@
194: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
196: Logically Collective on MatFDColoring
198: Input Parameters:
199: + coloring - the coloring context
200: . brows - number of rows in the block
201: - bcols - number of columns in the block
203: Level: intermediate
205: .keywords: Mat, coloring
207: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
209: @*/
210: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
211: {
216: if (brows != PETSC_DEFAULT) matfd->brows = brows;
217: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
218: return(0);
219: }
223: /*@
224: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
226: Collective on Mat
228: Input Parameters:
229: + mat - the matrix containing the nonzero structure of the Jacobian
230: . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
231: - color - the matrix coloring context
233: Level: beginner
235: .keywords: MatFDColoring, setup
237: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
238: @*/
239: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
240: {
246: if (color->setupcalled) return(0);
248: PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
249: if (mat->ops->fdcoloringsetup) {
250: (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
251: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
253: color->setupcalled = PETSC_TRUE;
254: PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
255: return(0);
256: }
260: /*@C
261: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
263: Not Collective
265: Input Parameters:
266: . coloring - the coloring context
268: Output Parameters:
269: + f - the function
270: - fctx - the optional user-defined function context
272: Level: intermediate
274: .keywords: Mat, Jacobian, finite differences, set, function
276: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
278: @*/
279: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
280: {
283: if (f) *f = matfd->f;
284: if (fctx) *fctx = matfd->fctx;
285: return(0);
286: }
290: /*@C
291: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
293: Logically Collective on MatFDColoring
295: Input Parameters:
296: + coloring - the coloring context
297: . f - the function
298: - fctx - the optional user-defined function context
300: Calling sequence of (*f) function:
301: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
302: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
304: Level: advanced
306: Notes: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
307: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
308: calling MatFDColoringApply()
310: Fortran Notes:
311: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
312: be used without SNES or within the SNES solvers.
314: .keywords: Mat, Jacobian, finite differences, set, function
316: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
318: @*/
319: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
320: {
323: matfd->f = f;
324: matfd->fctx = fctx;
325: return(0);
326: }
330: /*@
331: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
332: the options database.
334: Collective on MatFDColoring
336: The Jacobian, F'(u), is estimated with the differencing approximation
337: .vb
338: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
339: h = error_rel*u[i] if abs(u[i]) > umin
340: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
341: dx_{i} = (0, ... 1, .... 0)
342: .ve
344: Input Parameter:
345: . coloring - the coloring context
347: Options Database Keys:
348: + -mat_fd_coloring_err <err> - Sets <err> (square root
349: of relative error in the function)
350: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
351: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
352: . -mat_fd_coloring_view - Activates basic viewing
353: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
354: - -mat_fd_coloring_view draw - Activates drawing
356: Level: intermediate
358: .keywords: Mat, finite differences, parameters
360: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
362: @*/
363: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
364: {
366: PetscBool flg;
367: char value[3];
372: PetscObjectOptionsBegin((PetscObject)matfd);
373: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
374: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
375: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
376: if (flg) {
377: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
378: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
379: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
380: }
381: PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
382: PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
383: if (flg && matfd->bcols > matfd->ncolors) {
384: /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
385: matfd->bcols = matfd->ncolors;
386: }
388: /* process any options handlers added with PetscObjectAddOptionsHandler() */
389: PetscObjectProcessOptionsHandlers((PetscObject)matfd);
390: PetscOptionsEnd();
391: return(0);
392: }
396: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
397: {
398: PetscErrorCode ierr;
399: PetscBool flg;
400: PetscViewer viewer;
401: PetscViewerFormat format;
404: if (prefix) {
405: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),prefix,optionname,&viewer,&format,&flg);
406: } else {
407: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
408: }
409: if (flg) {
410: PetscViewerPushFormat(viewer,format);
411: MatFDColoringView(fd,viewer);
412: PetscViewerPopFormat(viewer);
413: PetscViewerDestroy(&viewer);
414: }
415: return(0);
416: }
420: /*@
421: MatFDColoringCreate - Creates a matrix coloring context for finite difference
422: computation of Jacobians.
424: Collective on Mat
426: Input Parameters:
427: + mat - the matrix containing the nonzero structure of the Jacobian
428: - iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()
430: Output Parameter:
431: . color - the new coloring context
433: Level: intermediate
435: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
436: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
437: MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring()
438: @*/
439: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
440: {
441: MatFDColoring c;
442: MPI_Comm comm;
444: PetscInt M,N;
448: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
449: PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
450: MatGetSize(mat,&M,&N);
451: if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
452: PetscObjectGetComm((PetscObject)mat,&comm);
453: PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);
455: c->ctype = iscoloring->ctype;
457: if (mat->ops->fdcoloringcreate) {
458: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
459: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
461: MatGetVecs(mat,NULL,&c->w1);
462: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
463: VecDuplicate(c->w1,&c->w2);
464: PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);
466: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
467: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
468: c->currentcolor = -1;
469: c->htype = "wp";
470: c->fset = PETSC_FALSE;
471: c->setupcalled = PETSC_FALSE;
473: *color = c;
474: PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
475: PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
476: return(0);
477: }
481: /*@
482: MatFDColoringDestroy - Destroys a matrix coloring context that was created
483: via MatFDColoringCreate().
485: Collective on MatFDColoring
487: Input Parameter:
488: . c - coloring context
490: Level: intermediate
492: .seealso: MatFDColoringCreate()
493: @*/
494: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
495: {
497: PetscInt i;
498: MatFDColoring color = *c;
501: if (!*c) return(0);
502: if (--((PetscObject)color)->refct > 0) {*c = 0; return(0);}
504: for (i=0; i<color->ncolors; i++) {
505: PetscFree(color->columns[i]);
506: }
507: PetscFree(color->ncolumns);
508: PetscFree(color->columns);
509: PetscFree(color->nrows);
510: if (color->htype[0] == 'w') {
511: PetscFree(color->matentry2);
512: } else {
513: PetscFree(color->matentry);
514: }
515: PetscFree(color->dy);
516: if (color->vscale) {VecDestroy(&color->vscale);}
517: VecDestroy(&color->w1);
518: VecDestroy(&color->w2);
519: VecDestroy(&color->w3);
520: PetscHeaderDestroy(c);
521: return(0);
522: }
526: /*@C
527: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
528: that are currently being perturbed.
530: Not Collective
532: Input Parameters:
533: . coloring - coloring context created with MatFDColoringCreate()
535: Output Parameters:
536: + n - the number of local columns being perturbed
537: - cols - the column indices, in global numbering
539: Level: intermediate
541: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
543: .keywords: coloring, Jacobian, finite differences
544: @*/
545: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
546: {
548: if (coloring->currentcolor >= 0) {
549: *n = coloring->ncolumns[coloring->currentcolor];
550: *cols = coloring->columns[coloring->currentcolor];
551: } else {
552: *n = 0;
553: }
554: return(0);
555: }
559: /*@
560: MatFDColoringApply - Given a matrix for which a MatFDColoring context
561: has been created, computes the Jacobian for a function via finite differences.
563: Collective on MatFDColoring
565: Input Parameters:
566: + mat - location to store Jacobian
567: . coloring - coloring context created with MatFDColoringCreate()
568: . x1 - location at which Jacobian is to be computed
569: - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
571: Options Database Keys:
572: + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
573: . -mat_fd_coloring_view - Activates basic viewing or coloring
574: . -mat_fd_coloring_view draw - Activates drawing of coloring
575: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
577: Level: intermediate
579: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()
581: .keywords: coloring, Jacobian, finite differences
582: @*/
583: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
584: {
586: PetscBool flg = PETSC_FALSE;
592: if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
593: if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
594: if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");
596: MatSetUnfactored(J);
597: PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);
598: if (flg) {
599: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
600: } else {
601: PetscBool assembled;
602: MatAssembled(J,&assembled);
603: if (assembled) {
604: MatZeroEntries(J);
605: }
606: }
608: PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
609: (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
610: PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
611: return(0);
612: }