Actual source code: fdmatrix.c
1: #define PETSCMAT_DLL
3: /*
4: This is where the abstract matrix operations are defined that are
5: used for finite difference computations of Jacobians using coloring.
6: */
8: #include include/private/matimpl.h
10: /* Logging support */
11: PetscCookie MAT_FDCOLORING_COOKIE = 0;
15: PetscErrorCode MatFDColoringSetF(MatFDColoring fd,Vec F)
16: {
18: fd->F = F;
19: return(0);
20: }
24: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
25: {
26: MatFDColoring fd = (MatFDColoring)Aa;
28: PetscInt i,j;
29: PetscReal x,y;
33: /* loop over colors */
34: for (i=0; i<fd->ncolors; i++) {
35: for (j=0; j<fd->nrows[i]; j++) {
36: y = fd->M - fd->rows[i][j] - fd->rstart;
37: x = fd->columnsforrow[i][j];
38: PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
39: }
40: }
41: return(0);
42: }
46: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
47: {
49: PetscTruth isnull;
50: PetscDraw draw;
51: PetscReal xr,yr,xl,yl,h,w;
54: PetscViewerDrawGetDraw(viewer,0,&draw);
55: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
57: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
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: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
63: PetscObjectCompose((PetscObject)fd,"Zoomviewer",PETSC_NULL);
64: return(0);
65: }
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: PetscTruth isdraw,iascii;
103: PetscViewerFormat format;
107: if (!viewer) {
108: PetscViewerASCIIGetStdout(c->comm,&viewer);
109: }
113: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
114: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
115: if (isdraw) {
116: MatFDColoringView_Draw(c,viewer);
117: } else if (iascii) {
118: PetscViewerASCIIPrintf(viewer,"MatFDColoring Object:\n");
119: PetscViewerASCIIPrintf(viewer," Error tolerance=%G\n",c->error_rel);
120: PetscViewerASCIIPrintf(viewer," Umin=%G\n",c->umin);
121: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
123: PetscViewerGetFormat(viewer,&format);
124: if (format != PETSC_VIEWER_ASCII_INFO) {
125: for (i=0; i<c->ncolors; i++) {
126: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
127: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
128: for (j=0; j<c->ncolumns[i]; j++) {
129: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
130: }
131: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
132: for (j=0; j<c->nrows[i]; j++) {
133: PetscViewerASCIIPrintf(viewer," %D %D \n",c->rows[i][j],c->columnsforrow[i][j]);
134: }
135: }
136: }
137: PetscViewerFlush(viewer);
138: } else {
139: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for MatFDColoring",((PetscObject)viewer)->type_name);
140: }
141: return(0);
142: }
146: /*@
147: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
148: a Jacobian matrix using finite differences.
150: Collective on MatFDColoring
152: The Jacobian is estimated with the differencing approximation
153: .vb
154: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
155: h = error_rel*u[i] if abs(u[i]) > umin
156: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
157: dx_{i} = (0, ... 1, .... 0)
158: .ve
160: Input Parameters:
161: + coloring - the coloring context
162: . error_rel - relative error
163: - umin - minimum allowable u-value magnitude
165: Level: advanced
167: .keywords: Mat, finite differences, coloring, set, parameters
169: .seealso: MatFDColoringCreate()
170: @*/
171: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
172: {
176: if (error != PETSC_DEFAULT) matfd->error_rel = error;
177: if (umin != PETSC_DEFAULT) matfd->umin = umin;
178: return(0);
179: }
183: /*@
184: MatFDColoringSetFrequency - Sets the frequency for computing new Jacobian
185: matrices.
187: Collective on MatFDColoring
189: Input Parameters:
190: + coloring - the coloring context
191: - freq - frequency (default is 1)
193: Options Database Keys:
194: . -mat_fd_coloring_freq <freq> - Sets coloring frequency
196: Level: advanced
198: Notes:
199: Using a modified Newton strategy, where the Jacobian remains fixed for several
200: iterations, can be cost effective in terms of overall nonlinear solution
201: efficiency. This parameter indicates that a new Jacobian will be computed every
202: <freq> nonlinear iterations.
204: .keywords: Mat, finite differences, coloring, set, frequency
206: .seealso: MatFDColoringCreate(), MatFDColoringGetFrequency(), MatFDColoringSetRecompute()
207: @*/
208: PetscErrorCode MatFDColoringSetFrequency(MatFDColoring matfd,PetscInt freq)
209: {
213: matfd->freq = freq;
214: return(0);
215: }
219: /*@
220: MatFDColoringGetFrequency - Gets the frequency for computing new Jacobian
221: matrices.
223: Not Collective
225: Input Parameters:
226: . coloring - the coloring context
228: Output Parameters:
229: . freq - frequency (default is 1)
231: Options Database Keys:
232: . -mat_fd_coloring_freq <freq> - Sets coloring frequency
234: Level: advanced
236: Notes:
237: Using a modified Newton strategy, where the Jacobian remains fixed for several
238: iterations, can be cost effective in terms of overall nonlinear solution
239: efficiency. This parameter indicates that a new Jacobian will be computed every
240: <freq> nonlinear iterations.
242: .keywords: Mat, finite differences, coloring, get, frequency
244: .seealso: MatFDColoringSetFrequency()
245: @*/
246: PetscErrorCode MatFDColoringGetFrequency(MatFDColoring matfd,PetscInt *freq)
247: {
250: *freq = matfd->freq;
251: return(0);
252: }
256: /*@C
257: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
259: Collective on MatFDColoring
261: Input Parameters:
262: . coloring - the coloring context
264: Output Parameters:
265: + f - the function
266: - fctx - the optional user-defined function context
268: Level: intermediate
270: .keywords: Mat, Jacobian, finite differences, set, function
271: @*/
272: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
273: {
276: if (f) *f = matfd->f;
277: if (fctx) *fctx = matfd->fctx;
278: return(0);
279: }
283: /*@C
284: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
286: Collective on MatFDColoring
288: Input Parameters:
289: + coloring - the coloring context
290: . f - the function
291: - fctx - the optional user-defined function context
293: Level: intermediate
295: Notes:
296: In Fortran you must call MatFDColoringSetFunctionSNES() for a coloring object to
297: be used with the SNES solvers and MatFDColoringSetFunctionTS() if it is to be used
298: with the TS solvers.
300: .keywords: Mat, Jacobian, finite differences, set, function
301: @*/
302: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
303: {
306: matfd->f = f;
307: matfd->fctx = fctx;
308: return(0);
309: }
313: /*@
314: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
315: the options database.
317: Collective on MatFDColoring
319: The Jacobian, F'(u), is estimated with the differencing approximation
320: .vb
321: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
322: h = error_rel*u[i] if abs(u[i]) > umin
323: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
324: dx_{i} = (0, ... 1, .... 0)
325: .ve
327: Input Parameter:
328: . coloring - the coloring context
330: Options Database Keys:
331: + -mat_fd_coloring_err <err> - Sets <err> (square root
332: of relative error in the function)
333: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
334: . -mat_fd_coloring_freq <freq> - Sets frequency of computing a new Jacobian
335: . -mat_fd_type - "wp" or "ds" (see MATSNESMF_WP or MATSNESMF_DS)
336: . -mat_fd_coloring_view - Activates basic viewing
337: . -mat_fd_coloring_view_info - Activates viewing info
338: - -mat_fd_coloring_view_draw - Activates drawing
340: Level: intermediate
342: .keywords: Mat, finite differences, parameters
344: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
346: @*/
347: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
348: {
350: PetscTruth flg;
351: char value[3];
356: PetscOptionsBegin(matfd->comm,matfd->prefix,"Jacobian computation via finite differences option","MatFD");
357: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
358: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
359: PetscOptionsInt("-mat_fd_coloring_freq","How often Jacobian is recomputed","MatFDColoringSetFrequency",matfd->freq,&matfd->freq,0);
360: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,2,&flg);
361: if (flg) {
362: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
363: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
364: else SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
365: }
366: /* not used here; but so they are presented in the GUI */
367: PetscOptionsName("-mat_fd_coloring_view","Print entire datastructure used for Jacobian","None",0);
368: PetscOptionsName("-mat_fd_coloring_view_info","Print number of colors etc for Jacobian","None",0);
369: PetscOptionsName("-mat_fd_coloring_view_draw","Plot nonzero structure ofJacobian","None",0);
370: PetscOptionsEnd();
371: return(0);
372: }
376: PetscErrorCode MatFDColoringView_Private(MatFDColoring fd)
377: {
379: PetscTruth flg;
380: PetscViewer viewer;
383: PetscViewerASCIIGetStdout(fd->comm,&viewer);
384: PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view",&flg);
385: if (flg) {
386: MatFDColoringView(fd,viewer);
387: }
388: PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_info",&flg);
389: if (flg) {
390: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
391: MatFDColoringView(fd,viewer);
392: PetscViewerPopFormat(viewer);
393: }
394: PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_draw",&flg);
395: if (flg) {
396: MatFDColoringView(fd,PETSC_VIEWER_DRAW_(fd->comm));
397: PetscViewerFlush(PETSC_VIEWER_DRAW_(fd->comm));
398: }
399: return(0);
400: }
404: /*@
405: MatFDColoringCreate - Creates a matrix coloring context for finite difference
406: computation of Jacobians.
408: Collective on Mat
410: Input Parameters:
411: + mat - the matrix containing the nonzero structure of the Jacobian
412: - iscoloring - the coloring of the matrix
414: Output Parameter:
415: . color - the new coloring context
416:
417: Level: intermediate
419: .seealso: MatFDColoringDestroy(),SNESDefaultComputeJacobianColor(), ISColoringCreate(),
420: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
421: MatFDColoringSetFrequency(), MatFDColoringSetRecompute(), MatFDColoringView(),
422: MatFDColoringSetParameters()
423: @*/
424: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
425: {
426: MatFDColoring c;
427: MPI_Comm comm;
429: PetscInt M,N;
430: PetscMPIInt size;
434: MatGetSize(mat,&M,&N);
435: if (M != N) SETERRQ(PETSC_ERR_SUP,"Only for square matrices");
437: PetscObjectGetComm((PetscObject)mat,&comm);
438: PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_COOKIE,0,"MatFDColoring",comm,MatFDColoringDestroy,MatFDColoringView);
440: MPI_Comm_size(comm,&size);
441: c->ctype = iscoloring->ctype;
443: if (mat->ops->fdcoloringcreate) {
444: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
445: } else {
446: SETERRQ(PETSC_ERR_SUP,"Code not yet written for this matrix type");
447: }
449: MatGetVecs(mat,PETSC_NULL,&c->w1);
450: PetscLogObjectParent(c,c->w1);
451: VecDuplicate(c->w1,&c->w2);
452: PetscLogObjectParent(c,c->w2);
454: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
455: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
456: c->freq = 1;
457: c->usersetsrecompute = PETSC_FALSE;
458: c->recompute = PETSC_FALSE;
459: c->currentcolor = -1;
460: c->htype = "wp";
462: *color = c;
464: return(0);
465: }
469: /*@
470: MatFDColoringDestroy - Destroys a matrix coloring context that was created
471: via MatFDColoringCreate().
473: Collective on MatFDColoring
475: Input Parameter:
476: . c - coloring context
478: Level: intermediate
480: .seealso: MatFDColoringCreate()
481: @*/
482: PetscErrorCode MatFDColoringDestroy(MatFDColoring c)
483: {
485: PetscInt i;
488: if (--c->refct > 0) return(0);
490: for (i=0; i<c->ncolors; i++) {
491: PetscFree(c->columns[i]);
492: PetscFree(c->rows[i]);
493: PetscFree(c->columnsforrow[i]);
494: if (c->vscaleforrow) {PetscFree(c->vscaleforrow[i]);}
495: }
496: PetscFree(c->ncolumns);
497: PetscFree(c->columns);
498: PetscFree(c->nrows);
499: PetscFree(c->rows);
500: PetscFree(c->columnsforrow);
501: PetscFree(c->vscaleforrow);
502: if (c->vscale) {VecDestroy(c->vscale);}
503: if (c->w1) {
504: VecDestroy(c->w1);
505: VecDestroy(c->w2);
506: }
507: if (c->w3){
508: VecDestroy(c->w3);
509: }
510: PetscHeaderDestroy(c);
511: return(0);
512: }
516: /*@C
517: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
518: that are currently being perturbed.
520: Not Collective
522: Input Parameters:
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: intermediate
531: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
533: .keywords: coloring, Jacobian, finite differences
534: @*/
535: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
536: {
538: if (coloring->currentcolor >= 0) {
539: *n = coloring->ncolumns[coloring->currentcolor];
540: *cols = coloring->columns[coloring->currentcolor];
541: } else {
542: *n = 0;
543: }
544: return(0);
545: }
547: #include "petscda.h" /*I "petscda.h" I*/
551: /*@
552: MatFDColoringApply - Given a matrix for which a MatFDColoring context
553: has been created, computes the Jacobian for a function via finite differences.
555: Collective on MatFDColoring
557: Input Parameters:
558: + mat - location to store Jacobian
559: . coloring - coloring context created with MatFDColoringCreate()
560: . x1 - location at which Jacobian is to be computed
561: - sctx - optional context required by function (actually a SNES context)
563: Options Database Keys:
564: + -mat_fd_coloring_freq <freq> - Sets coloring frequency
565: . -mat_fd_type - "wp" or "ds" (see MATSNESMF_WP or MATSNESMF_DS)
566: . -mat_fd_coloring_view - Activates basic viewing or coloring
567: . -mat_fd_coloring_view_draw - Activates drawing of coloring
568: - -mat_fd_coloring_view_info - Activates viewing of coloring info
570: Level: intermediate
572: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView()
574: .keywords: coloring, Jacobian, finite differences
575: @*/
576: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
577: {
578: PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f;
580: PetscInt k,start,end,l,row,col,srow,**vscaleforrow,m1,m2;
581: PetscScalar dx,*y,*xx,*w3_array;
582: PetscScalar *vscale_array;
583: PetscReal epsilon = coloring->error_rel,umin = coloring->umin,unorm;
584: Vec w1=coloring->w1,w2=coloring->w2,w3;
585: void *fctx = coloring->fctx;
586: PetscTruth flg;
587: PetscInt ctype=coloring->ctype,N,col_start=0,col_end=0;
588: Vec x1_tmp;
594: if (!f) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
596: if (coloring->usersetsrecompute) {
597: if (!coloring->recompute) {
598: *flag = SAME_PRECONDITIONER;
599: PetscInfo(J,"Skipping Jacobian, since user called MatFDColorSetRecompute()\n");
600: return(0);
601: } else {
602: coloring->recompute = PETSC_FALSE;
603: }
604: }
607: MatSetUnfactored(J);
608: PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);
609: if (flg) {
610: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
611: } else {
612: PetscTruth assembled;
613: MatAssembled(J,&assembled);
614: if (assembled) {
615: MatZeroEntries(J);
616: }
617: }
619: x1_tmp = x1;
620: if (!coloring->vscale){
621: VecDuplicate(x1_tmp,&coloring->vscale);
622: }
623:
624: /*
625: This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets
626: coloring->F for the coarser grids from the finest
627: */
628: if (coloring->F) {
629: VecGetLocalSize(coloring->F,&m1);
630: VecGetLocalSize(w1,&m2);
631: if (m1 != m2) {
632: coloring->F = 0;
633: }
634: }
636: if (coloring->htype[0] == 'w') { /* tacky test; need to make systematic if we add other approaches to computing h*/
637: VecNorm(x1_tmp,NORM_2,&unorm);
638: }
639: VecGetOwnershipRange(w1,&start,&end); /* OwnershipRange is used by ghosted x! */
640:
641: /* Set w1 = F(x1) */
642: if (coloring->F) {
643: w1 = coloring->F; /* use already computed value of function */
644: coloring->F = 0;
645: } else {
647: (*f)(sctx,x1_tmp,w1,fctx);
649: }
650:
651: if (!coloring->w3) {
652: VecDuplicate(x1_tmp,&coloring->w3);
653: PetscLogObjectParent(coloring,coloring->w3);
654: }
655: w3 = coloring->w3;
657: /* Compute all the local scale factors, including ghost points */
658: VecGetLocalSize(x1_tmp,&N);
659: VecGetArray(x1_tmp,&xx);
660: VecGetArray(coloring->vscale,&vscale_array);
661: if (ctype == IS_COLORING_GHOSTED){
662: col_start = 0; col_end = N;
663: } else if (ctype == IS_COLORING_GLOBAL){
664: xx = xx - start;
665: vscale_array = vscale_array - start;
666: col_start = start; col_end = N + start;
667: }
668: for (col=col_start; col<col_end; col++){
669: /* Loop over each local column, vscale[col] = 1./(epsilon*dx[col]) */
670: if (coloring->htype[0] == 'w') {
671: dx = 1.0 + unorm;
672: } else {
673: dx = xx[col];
674: }
675: if (dx == 0.0) dx = 1.0;
676: #if !defined(PETSC_USE_COMPLEX)
677: if (dx < umin && dx >= 0.0) dx = umin;
678: else if (dx < 0.0 && dx > -umin) dx = -umin;
679: #else
680: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
681: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
682: #endif
683: dx *= epsilon;
684: vscale_array[col] = 1.0/dx;
685: }
686: if (ctype == IS_COLORING_GLOBAL) vscale_array = vscale_array + start;
687: VecRestoreArray(coloring->vscale,&vscale_array);
688: if (ctype == IS_COLORING_GLOBAL){
689: VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
690: VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
691: }
692:
693: if (coloring->vscaleforrow) {
694: vscaleforrow = coloring->vscaleforrow;
695: } else {
696: SETERRQ(PETSC_ERR_ARG_NULL,"Null Object: coloring->vscaleforrow");
697: }
699: /*
700: Loop over each color
701: */
702: VecGetArray(coloring->vscale,&vscale_array);
703: for (k=0; k<coloring->ncolors; k++) {
704: coloring->currentcolor = k;
705: VecCopy(x1_tmp,w3);
706: VecGetArray(w3,&w3_array);
707: if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array - start;
708: /*
709: Loop over each column associated with color
710: adding the perturbation to the vector w3.
711: */
712: for (l=0; l<coloring->ncolumns[k]; l++) {
713: col = coloring->columns[k][l]; /* local column of the matrix we are probing for */
714: if (coloring->htype[0] == 'w') {
715: dx = 1.0 + unorm;
716: } else {
717: dx = xx[col];
718: }
719: if (dx == 0.0) dx = 1.0;
720: #if !defined(PETSC_USE_COMPLEX)
721: if (dx < umin && dx >= 0.0) dx = umin;
722: else if (dx < 0.0 && dx > -umin) dx = -umin;
723: #else
724: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
725: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
726: #endif
727: dx *= epsilon;
728: if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter");
729: w3_array[col] += dx;
730: }
731: if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array + start;
732: VecRestoreArray(w3,&w3_array);
734: /*
735: Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
736: w2 = F(x1 + dx) - F(x1)
737: */
739: (*f)(sctx,w3,w2,fctx);
741: VecAXPY(w2,-1.0,w1);
742:
743: /*
744: Loop over rows of vector, putting results into Jacobian matrix
745: */
746: VecGetArray(w2,&y);
747: for (l=0; l<coloring->nrows[k]; l++) {
748: row = coloring->rows[k][l]; /* local row index */
749: col = coloring->columnsforrow[k][l]; /* global column index */
750: y[row] *= vscale_array[vscaleforrow[k][l]];
751: srow = row + start;
752: MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
753: }
754: VecRestoreArray(w2,&y);
755: } /* endof for each color */
756: if (ctype == IS_COLORING_GLOBAL) xx = xx + start;
757: VecRestoreArray(coloring->vscale,&vscale_array);
758: VecRestoreArray(x1_tmp,&xx);
759:
760: coloring->currentcolor = -1;
761: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
762: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
765: PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test",&flg);
766: if (flg) {
767: MatNullSpaceTest(J->nullsp,J);
768: }
769: MatFDColoringView_Private(coloring);
770: return(0);
771: }
775: /*@
776: MatFDColoringApplyTS - Given a matrix for which a MatFDColoring context
777: has been created, computes the Jacobian for a function via finite differences.
779: Collective on Mat, MatFDColoring, and Vec
781: Input Parameters:
782: + mat - location to store Jacobian
783: . coloring - coloring context created with MatFDColoringCreate()
784: . x1 - location at which Jacobian is to be computed
785: - sctx - optional context required by function (actually a SNES context)
787: Options Database Keys:
788: . -mat_fd_coloring_freq <freq> - Sets coloring frequency
790: Level: intermediate
792: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView()
794: .keywords: coloring, Jacobian, finite differences
795: @*/
796: PetscErrorCode MatFDColoringApplyTS(Mat J,MatFDColoring coloring,PetscReal t,Vec x1,MatStructure *flag,void *sctx)
797: {
798: PetscErrorCode (*f)(void*,PetscReal,Vec,Vec,void*)=(PetscErrorCode (*)(void*,PetscReal,Vec,Vec,void *))coloring->f;
800: PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow;
801: PetscScalar dx,*y,*xx,*w3_array;
802: PetscScalar *vscale_array;
803: PetscReal epsilon = coloring->error_rel,umin = coloring->umin;
804: Vec w1=coloring->w1,w2=coloring->w2,w3;
805: void *fctx = coloring->fctx;
806: PetscTruth flg;
814: if (!coloring->w3) {
815: VecDuplicate(x1,&coloring->w3);
816: PetscLogObjectParent(coloring,coloring->w3);
817: }
818: w3 = coloring->w3;
820: MatSetUnfactored(J);
821: PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);
822: if (flg) {
823: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
824: } else {
825: PetscTruth assembled;
826: MatAssembled(J,&assembled);
827: if (assembled) {
828: MatZeroEntries(J);
829: }
830: }
832: VecGetOwnershipRange(x1,&start,&end);
833: VecGetSize(x1,&N);
835: (*f)(sctx,t,x1,w1,fctx);
838: /*
839: Compute all the scale factors and share with other processors
840: */
841: VecGetArray(x1,&xx);xx = xx - start;
842: VecGetArray(coloring->vscale,&vscale_array);vscale_array = vscale_array - start;
843: for (k=0; k<coloring->ncolors; k++) {
844: /*
845: Loop over each column associated with color adding the
846: perturbation to the vector w3.
847: */
848: for (l=0; l<coloring->ncolumns[k]; l++) {
849: col = coloring->columns[k][l]; /* column of the matrix we are probing for */
850: dx = xx[col];
851: if (dx == 0.0) dx = 1.0;
852: #if !defined(PETSC_USE_COMPLEX)
853: if (dx < umin && dx >= 0.0) dx = umin;
854: else if (dx < 0.0 && dx > -umin) dx = -umin;
855: #else
856: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
857: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
858: #endif
859: dx *= epsilon;
860: vscale_array[col] = 1.0/dx;
861: }
862: }
863: vscale_array = vscale_array - start;VecRestoreArray(coloring->vscale,&vscale_array);
864: VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
865: VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
867: if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow;
868: else vscaleforrow = coloring->columnsforrow;
870: VecGetArray(coloring->vscale,&vscale_array);
871: /*
872: Loop over each color
873: */
874: for (k=0; k<coloring->ncolors; k++) {
875: VecCopy(x1,w3);
876: VecGetArray(w3,&w3_array);w3_array = w3_array - start;
877: /*
878: Loop over each column associated with color adding the
879: perturbation to the vector w3.
880: */
881: for (l=0; l<coloring->ncolumns[k]; l++) {
882: col = coloring->columns[k][l]; /* column of the matrix we are probing for */
883: dx = xx[col];
884: if (dx == 0.0) dx = 1.0;
885: #if !defined(PETSC_USE_COMPLEX)
886: if (dx < umin && dx >= 0.0) dx = umin;
887: else if (dx < 0.0 && dx > -umin) dx = -umin;
888: #else
889: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
890: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
891: #endif
892: dx *= epsilon;
893: w3_array[col] += dx;
894: }
895: w3_array = w3_array + start; VecRestoreArray(w3,&w3_array);
897: /*
898: Evaluate function at x1 + dx (here dx is a vector of perturbations)
899: */
901: (*f)(sctx,t,w3,w2,fctx);
903: VecAXPY(w2,-1.0,w1);
905: /*
906: Loop over rows of vector, putting results into Jacobian matrix
907: */
908: VecGetArray(w2,&y);
909: for (l=0; l<coloring->nrows[k]; l++) {
910: row = coloring->rows[k][l];
911: col = coloring->columnsforrow[k][l];
912: y[row] *= vscale_array[vscaleforrow[k][l]];
913: srow = row + start;
914: MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
915: }
916: VecRestoreArray(w2,&y);
917: }
918: VecRestoreArray(coloring->vscale,&vscale_array);
919: xx = xx + start; VecRestoreArray(x1,&xx);
920: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
921: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
923: return(0);
924: }
929: /*@C
930: MatFDColoringSetRecompute - Indicates that the next time a Jacobian preconditioner
931: is needed it sholuld be recomputed. Once this is called and the new Jacobian is computed
932: no additional Jacobian's will be computed (the same one will be used) until this is
933: called again.
935: Collective on MatFDColoring
937: Input Parameters:
938: . c - the coloring context
940: Level: intermediate
942: Notes: The MatFDColoringSetFrequency() is ignored once this is called
944: .seealso: MatFDColoringCreate(), MatFDColoringSetFrequency()
946: .keywords: Mat, finite differences, coloring
947: @*/
948: PetscErrorCode MatFDColoringSetRecompute(MatFDColoring c)
949: {
952: c->usersetsrecompute = PETSC_TRUE;
953: c->recompute = PETSC_TRUE;
954: return(0);
955: }