Actual source code: fdmatrix.c
petsc-3.6.1 2015-08-06
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); if (isnull) return(0);
64: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
66: xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0;
67: xr += w; yr += h; xl = -w; yl = -h;
68: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
69: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
70: PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
71: return(0);
72: }
76: /*@C
77: MatFDColoringView - Views a finite difference coloring context.
79: Collective on MatFDColoring
81: Input Parameters:
82: + c - the coloring context
83: - viewer - visualization context
85: Level: intermediate
87: Notes:
88: The available visualization contexts include
89: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
90: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
91: output where only the first processor opens
92: the file. All other processors send their
93: data to the first processor to print.
94: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
96: Notes:
97: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
98: involves more than 33 then some seemingly identical colors are displayed making it look
99: like an illegal coloring. This is just a graphical artifact.
101: .seealso: MatFDColoringCreate()
103: .keywords: Mat, finite differences, coloring, view
104: @*/
105: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
106: {
107: PetscErrorCode ierr;
108: PetscInt i,j;
109: PetscBool isdraw,iascii;
110: PetscViewerFormat format;
114: if (!viewer) {
115: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
116: }
120: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
121: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
122: if (isdraw) {
123: MatFDColoringView_Draw(c,viewer);
124: } else if (iascii) {
125: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
126: PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);
127: PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);
128: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
130: PetscViewerGetFormat(viewer,&format);
131: if (format != PETSC_VIEWER_ASCII_INFO) {
132: PetscInt row,col,nz;
133: nz = 0;
134: for (i=0; i<c->ncolors; i++) {
135: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
136: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
137: for (j=0; j<c->ncolumns[i]; j++) {
138: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
139: }
140: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
141: for (j=0; j<c->nrows[i]; j++) {
142: row = c->matentry[nz].row;
143: col = c->matentry[nz++].col;
144: PetscViewerASCIIPrintf(viewer," %D %D \n",row,col);
145: }
146: }
147: }
148: PetscViewerFlush(viewer);
149: }
150: return(0);
151: }
155: /*@
156: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
157: a Jacobian matrix using finite differences.
159: Logically Collective on MatFDColoring
161: The Jacobian is estimated with the differencing approximation
162: .vb
163: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
164: h = error_rel*u[i] if abs(u[i]) > umin
165: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
166: dx_{i} = (0, ... 1, .... 0)
167: .ve
169: Input Parameters:
170: + coloring - the coloring context
171: . error_rel - relative error
172: - umin - minimum allowable u-value magnitude
174: Level: advanced
176: .keywords: Mat, finite differences, coloring, set, parameters
178: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
180: @*/
181: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
182: {
187: if (error != PETSC_DEFAULT) matfd->error_rel = error;
188: if (umin != PETSC_DEFAULT) matfd->umin = umin;
189: return(0);
190: }
194: /*@
195: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
197: Logically Collective on MatFDColoring
199: Input Parameters:
200: + coloring - the coloring context
201: . brows - number of rows in the block
202: - bcols - number of columns in the block
204: Level: intermediate
206: .keywords: Mat, coloring
208: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
210: @*/
211: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
212: {
217: if (brows != PETSC_DEFAULT) matfd->brows = brows;
218: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
219: return(0);
220: }
224: /*@
225: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
227: Collective on Mat
229: Input Parameters:
230: + mat - the matrix containing the nonzero structure of the Jacobian
231: . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
232: - color - the matrix coloring context
234: Level: beginner
236: .keywords: MatFDColoring, setup
238: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
239: @*/
240: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
241: {
247: if (color->setupcalled) return(0);
249: PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
250: if (mat->ops->fdcoloringsetup) {
251: (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
252: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
254: color->setupcalled = PETSC_TRUE;
255: PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
256: return(0);
257: }
261: /*@C
262: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
264: Not Collective
266: Input Parameters:
267: . coloring - the coloring context
269: Output Parameters:
270: + f - the function
271: - fctx - the optional user-defined function context
273: Level: intermediate
275: .keywords: Mat, Jacobian, finite differences, set, function
277: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
279: @*/
280: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
281: {
284: if (f) *f = matfd->f;
285: if (fctx) *fctx = matfd->fctx;
286: return(0);
287: }
291: /*@C
292: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
294: Logically Collective on MatFDColoring
296: Input Parameters:
297: + coloring - the coloring context
298: . f - the function
299: - fctx - the optional user-defined function context
301: Calling sequence of (*f) function:
302: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
303: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
305: Level: advanced
307: Notes: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
308: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
309: calling MatFDColoringApply()
311: Fortran Notes:
312: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
313: be used without SNES or within the SNES solvers.
315: .keywords: Mat, Jacobian, finite differences, set, function
317: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
319: @*/
320: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
321: {
324: matfd->f = f;
325: matfd->fctx = fctx;
326: return(0);
327: }
331: /*@
332: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
333: the options database.
335: Collective on MatFDColoring
337: The Jacobian, F'(u), is estimated with the differencing approximation
338: .vb
339: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
340: h = error_rel*u[i] if abs(u[i]) > umin
341: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
342: dx_{i} = (0, ... 1, .... 0)
343: .ve
345: Input Parameter:
346: . coloring - the coloring context
348: Options Database Keys:
349: + -mat_fd_coloring_err <err> - Sets <err> (square root 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,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: MatCreateVecs(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: }