1
2 /*
3 This is where the abstract matrix operations are defined that are
4 used for finite difference computations of Jacobians using coloring.
5 */
6
7 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
8 #include <petsc/private/isimpl.h>
9
MatFDColoringSetF(MatFDColoring fd,Vec F)10 PetscErrorCode MatFDColoringSetF(MatFDColoring fd,Vec F)
11 {
12 PetscErrorCode ierr;
13
14 PetscFunctionBegin;
15 if (F) {
16 ierr = VecCopy(F,fd->w1);CHKERRQ(ierr);
17 fd->fset = PETSC_TRUE;
18 } else {
19 fd->fset = PETSC_FALSE;
20 }
21 PetscFunctionReturn(0);
22 }
23
24 #include <petscdraw.h>
MatFDColoringView_Draw_Zoom(PetscDraw draw,void * Aa)25 static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
26 {
27 MatFDColoring fd = (MatFDColoring)Aa;
28 PetscErrorCode ierr;
29 PetscInt i,j,nz,row;
30 PetscReal x,y;
31 MatEntry *Jentry=fd->matentry;
32
33 PetscFunctionBegin;
34 /* loop over colors */
35 nz = 0;
36 for (i=0; i<fd->ncolors; i++) {
37 for (j=0; j<fd->nrows[i]; j++) {
38 row = Jentry[nz].row;
39 y = fd->M - row - fd->rstart;
40 x = (PetscReal)Jentry[nz++].col;
41 ierr = PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);CHKERRQ(ierr);
42 }
43 }
44 PetscFunctionReturn(0);
45 }
46
MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)47 static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
48 {
49 PetscErrorCode ierr;
50 PetscBool isnull;
51 PetscDraw draw;
52 PetscReal xr,yr,xl,yl,h,w;
53
54 PetscFunctionBegin;
55 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
56 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr);
57 if (isnull) PetscFunctionReturn(0);
58
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 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr);
62 ierr = PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr);
63 ierr = PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);CHKERRQ(ierr);
64 ierr = PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);CHKERRQ(ierr);
65 ierr = PetscDrawSave(draw);CHKERRQ(ierr);
66 PetscFunctionReturn(0);
67 }
68
69 /*@C
70 MatFDColoringView - Views a finite difference coloring context.
71
72 Collective on MatFDColoring
73
74 Input Parameters:
75 + c - the coloring context
76 - viewer - visualization context
77
78 Level: intermediate
79
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
88
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.
93
94 .seealso: MatFDColoringCreate()
95
96 @*/
MatFDColoringView(MatFDColoring c,PetscViewer viewer)97 PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
98 {
99 PetscErrorCode ierr;
100 PetscInt i,j;
101 PetscBool isdraw,iascii;
102 PetscViewerFormat format;
103
104 PetscFunctionBegin;
105 PetscValidHeaderSpecific(c,MAT_FDCOLORING_CLASSID,1);
106 if (!viewer) {
107 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);CHKERRQ(ierr);
108 }
109 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
110 PetscCheckSameComm(c,1,viewer,2);
111
112 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
113 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
114 if (isdraw) {
115 ierr = MatFDColoringView_Draw(c,viewer);CHKERRQ(ierr);
116 } else if (iascii) {
117 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);CHKERRQ(ierr);
118 ierr = PetscViewerASCIIPrintf(viewer," Error tolerance=%g\n",(double)c->error_rel);CHKERRQ(ierr);
119 ierr = PetscViewerASCIIPrintf(viewer," Umin=%g\n",(double)c->umin);CHKERRQ(ierr);
120 ierr = PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);CHKERRQ(ierr);
121
122 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
123 if (format != PETSC_VIEWER_ASCII_INFO) {
124 PetscInt row,col,nz;
125 nz = 0;
126 for (i=0; i<c->ncolors; i++) {
127 ierr = PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);CHKERRQ(ierr);
128 ierr = PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);CHKERRQ(ierr);
129 for (j=0; j<c->ncolumns[i]; j++) {
130 ierr = PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);CHKERRQ(ierr);
131 }
132 ierr = PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);CHKERRQ(ierr);
133 if (c->matentry) {
134 for (j=0; j<c->nrows[i]; j++) {
135 row = c->matentry[nz].row;
136 col = c->matentry[nz++].col;
137 ierr = PetscViewerASCIIPrintf(viewer," %D %D \n",row,col);CHKERRQ(ierr);
138 }
139 }
140 }
141 }
142 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
143 }
144 PetscFunctionReturn(0);
145 }
146
147 /*@
148 MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
149 a Jacobian matrix using finite differences.
150
151 Logically Collective on MatFDColoring
152
153 The Jacobian is estimated with the differencing approximation
154 .vb
155 F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
156 htype = 'ds':
157 h = error_rel*u[i] if abs(u[i]) > umin
158 = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
159 dx_{i} = (0, ... 1, .... 0)
160
161 htype = 'wp':
162 h = error_rel * sqrt(1 + ||u||)
163 .ve
164
165 Input Parameters:
166 + coloring - the coloring context
167 . error_rel - relative error
168 - umin - minimum allowable u-value magnitude
169
170 Level: advanced
171
172 .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
173
174 @*/
MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)175 PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
176 {
177 PetscFunctionBegin;
178 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);
179 PetscValidLogicalCollectiveReal(matfd,error,2);
180 PetscValidLogicalCollectiveReal(matfd,umin,3);
181 if (error != PETSC_DEFAULT) matfd->error_rel = error;
182 if (umin != PETSC_DEFAULT) matfd->umin = umin;
183 PetscFunctionReturn(0);
184 }
185
186 /*@
187 MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
188
189 Logically Collective on MatFDColoring
190
191 Input Parameters:
192 + coloring - the coloring context
193 . brows - number of rows in the block
194 - bcols - number of columns in the block
195
196 Level: intermediate
197
198 .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
199
200 @*/
MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)201 PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
202 {
203 PetscFunctionBegin;
204 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);
205 PetscValidLogicalCollectiveInt(matfd,brows,2);
206 PetscValidLogicalCollectiveInt(matfd,bcols,3);
207 if (brows != PETSC_DEFAULT) matfd->brows = brows;
208 if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
209 PetscFunctionReturn(0);
210 }
211
212 /*@
213 MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
214
215 Collective on Mat
216
217 Input Parameters:
218 + mat - the matrix containing the nonzero structure of the Jacobian
219 . iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
220 - color - the matrix coloring context
221
222 Level: beginner
223
224 Notes: When the coloring type is IS_COLORING_LOCAL the coloring is in the local ordering of the unknowns.
225
226 .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
227 @*/
MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)228 PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
229 {
230 PetscErrorCode ierr;
231 PetscBool eq;
232
233 PetscFunctionBegin;
234 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
235 PetscValidHeaderSpecific(color,MAT_FDCOLORING_CLASSID,3);
236 if (color->setupcalled) PetscFunctionReturn(0);
237 ierr = PetscObjectCompareId((PetscObject)mat,color->matid,&eq);CHKERRQ(ierr);
238 if (!eq) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Matrix used with MatFDColoringSetUp() must be that used with MatFDColoringCreate()");
239
240 ierr = PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);CHKERRQ(ierr);
241 if (mat->ops->fdcoloringsetup) {
242 ierr = (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);CHKERRQ(ierr);
243 } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
244
245 color->setupcalled = PETSC_TRUE;
246 ierr = PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);CHKERRQ(ierr);
247 PetscFunctionReturn(0);
248 }
249
250 /*@C
251 MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
252
253 Not Collective
254
255 Input Parameters:
256 . coloring - the coloring context
257
258 Output Parameters:
259 + f - the function
260 - fctx - the optional user-defined function context
261
262 Level: intermediate
263
264 .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
265
266 @*/
MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (** f)(void),void ** fctx)267 PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
268 {
269 PetscFunctionBegin;
270 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);
271 if (f) *f = matfd->f;
272 if (fctx) *fctx = matfd->fctx;
273 PetscFunctionReturn(0);
274 }
275
276 /*@C
277 MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
278
279 Logically Collective on MatFDColoring
280
281 Input Parameters:
282 + coloring - the coloring context
283 . f - the function
284 - fctx - the optional user-defined function context
285
286 Calling sequence of (*f) function:
287 For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
288 If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
289
290 Level: advanced
291
292 Notes:
293 This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
294 SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
295 calling MatFDColoringApply()
296
297 Fortran Notes:
298 In Fortran you must call MatFDColoringSetFunction() for a coloring object to
299 be used without SNES or within the SNES solvers.
300
301 .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
302
303 @*/
MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (* f)(void),void * fctx)304 PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
305 {
306 PetscFunctionBegin;
307 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);
308 matfd->f = f;
309 matfd->fctx = fctx;
310 PetscFunctionReturn(0);
311 }
312
313 /*@
314 MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
315 the options database.
316
317 Collective on MatFDColoring
318
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
326
327 Input Parameter:
328 . coloring - the coloring context
329
330 Options Database Keys:
331 + -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
332 . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
333 . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
334 . -mat_fd_coloring_view - Activates basic viewing
335 . -mat_fd_coloring_view ::ascii_info - Activates viewing info
336 - -mat_fd_coloring_view draw - Activates drawing
337
338 Level: intermediate
339
340 .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
341
342 @*/
MatFDColoringSetFromOptions(MatFDColoring matfd)343 PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
344 {
345 PetscErrorCode ierr;
346 PetscBool flg;
347 char value[3];
348
349 PetscFunctionBegin;
350 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);
351
352 ierr = PetscObjectOptionsBegin((PetscObject)matfd);CHKERRQ(ierr);
353 ierr = PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,NULL);CHKERRQ(ierr);
354 ierr = PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,NULL);CHKERRQ(ierr);
355 ierr = PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,sizeof(value),&flg);CHKERRQ(ierr);
356 if (flg) {
357 if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
358 else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
359 else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
360 }
361 ierr = PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);CHKERRQ(ierr);
362 ierr = PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);CHKERRQ(ierr);
363 if (flg && matfd->bcols > matfd->ncolors) {
364 /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
365 matfd->bcols = matfd->ncolors;
366 }
367
368 /* process any options handlers added with PetscObjectAddOptionsHandler() */
369 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)matfd);CHKERRQ(ierr);
370 PetscOptionsEnd();CHKERRQ(ierr);
371 PetscFunctionReturn(0);
372 }
373
374 /*@C
375 MatFDColoringSetType - Sets the approach for computing the finite difference parameter
376
377 Collective on MatFDColoring
378
379 Input Parameters:
380 + coloring - the coloring context
381 - type - either MATMFFD_WP or MATMFFD_DS
382
383 Options Database Keys:
384 . -mat_fd_type - "wp" or "ds"
385
386 Note: It is goofy that the argument type is MatMFFDType since the MatFDColoring actually computes the matrix entries
387 but the process of computing the entries is the same as as with the MatMFFD operation so we should reuse the names instead of
388 introducing another one.
389
390 Level: intermediate
391
392 .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
393
394 @*/
MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)395 PetscErrorCode MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)
396 {
397 PetscFunctionBegin;
398 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);
399 /*
400 It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
401 and this function is being provided as patch for a release so it shouldn't change the implementaton
402 */
403 if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
404 else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
405 else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",type);
406 PetscFunctionReturn(0);
407 }
408
MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])409 PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
410 {
411 PetscErrorCode ierr;
412 PetscBool flg;
413 PetscViewer viewer;
414 PetscViewerFormat format;
415
416 PetscFunctionBegin;
417 if (prefix) {
418 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,prefix,optionname,&viewer,&format,&flg);CHKERRQ(ierr);
419 } else {
420 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);CHKERRQ(ierr);
421 }
422 if (flg) {
423 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr);
424 ierr = MatFDColoringView(fd,viewer);CHKERRQ(ierr);
425 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
426 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
427 }
428 PetscFunctionReturn(0);
429 }
430
431 /*@
432 MatFDColoringCreate - Creates a matrix coloring context for finite difference
433 computation of Jacobians.
434
435 Collective on Mat
436
437 Input Parameters:
438 + mat - the matrix containing the nonzero structure of the Jacobian
439 - iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()
440
441 Output Parameter:
442 . color - the new coloring context
443
444 Level: intermediate
445
446 .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
447 MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
448 MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring(), MatFDColoringSetValues()
449 @*/
MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring * color)450 PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
451 {
452 MatFDColoring c;
453 MPI_Comm comm;
454 PetscErrorCode ierr;
455 PetscInt M,N;
456
457 PetscFunctionBegin;
458 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
459 if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
460 ierr = PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);CHKERRQ(ierr);
461 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
462 if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
463 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr);
464 ierr = PetscHeaderCreate(c,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);CHKERRQ(ierr);
465
466 c->ctype = iscoloring->ctype;
467 ierr = PetscObjectGetId((PetscObject)mat,&c->matid);CHKERRQ(ierr);
468
469 if (mat->ops->fdcoloringcreate) {
470 ierr = (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);CHKERRQ(ierr);
471 } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
472
473 ierr = MatCreateVecs(mat,NULL,&c->w1);CHKERRQ(ierr);
474 /* Vec is used instensively in particular piece of scalar CPU code; won't benifit from bouncing back and forth to the GPU */
475 ierr = VecBindToCPU(c->w1,PETSC_TRUE);CHKERRQ(ierr);
476 ierr = PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);CHKERRQ(ierr);
477 ierr = VecDuplicate(c->w1,&c->w2);CHKERRQ(ierr);
478 /* Vec is used instensively in particular piece of scalar CPU code; won't benifit from bouncing back and forth to the GPU */
479 ierr = VecBindToCPU(c->w2,PETSC_TRUE);CHKERRQ(ierr);
480 ierr = PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);CHKERRQ(ierr);
481
482 c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
483 c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
484 c->currentcolor = -1;
485 c->htype = "wp";
486 c->fset = PETSC_FALSE;
487 c->setupcalled = PETSC_FALSE;
488
489 *color = c;
490 ierr = PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);CHKERRQ(ierr);
491 ierr = PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);CHKERRQ(ierr);
492 PetscFunctionReturn(0);
493 }
494
495 /*@
496 MatFDColoringDestroy - Destroys a matrix coloring context that was created
497 via MatFDColoringCreate().
498
499 Collective on MatFDColoring
500
501 Input Parameter:
502 . c - coloring context
503
504 Level: intermediate
505
506 .seealso: MatFDColoringCreate()
507 @*/
MatFDColoringDestroy(MatFDColoring * c)508 PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
509 {
510 PetscErrorCode ierr;
511 PetscInt i;
512 MatFDColoring color = *c;
513
514 PetscFunctionBegin;
515 if (!*c) PetscFunctionReturn(0);
516 if (--((PetscObject)color)->refct > 0) {*c = NULL; PetscFunctionReturn(0);}
517
518 /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
519 for (i=0; i<color->ncolors; i++) {
520 ierr = ISDestroy(&color->isa[i]);CHKERRQ(ierr);
521 }
522 ierr = PetscFree(color->isa);CHKERRQ(ierr);
523 ierr = PetscFree2(color->ncolumns,color->columns);CHKERRQ(ierr);
524 ierr = PetscFree(color->nrows);CHKERRQ(ierr);
525 if (color->htype[0] == 'w') {
526 ierr = PetscFree(color->matentry2);CHKERRQ(ierr);
527 } else {
528 ierr = PetscFree(color->matentry);CHKERRQ(ierr);
529 }
530 ierr = PetscFree(color->dy);CHKERRQ(ierr);
531 if (color->vscale) {ierr = VecDestroy(&color->vscale);CHKERRQ(ierr);}
532 ierr = VecDestroy(&color->w1);CHKERRQ(ierr);
533 ierr = VecDestroy(&color->w2);CHKERRQ(ierr);
534 ierr = VecDestroy(&color->w3);CHKERRQ(ierr);
535 ierr = PetscHeaderDestroy(c);CHKERRQ(ierr);
536 PetscFunctionReturn(0);
537 }
538
539 /*@C
540 MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
541 that are currently being perturbed.
542
543 Not Collective
544
545 Input Parameters:
546 . coloring - coloring context created with MatFDColoringCreate()
547
548 Output Parameters:
549 + n - the number of local columns being perturbed
550 - cols - the column indices, in global numbering
551
552 Level: advanced
553
554 Note: IF the matrix type is BAIJ, then the block column indices are returned
555
556 Fortran Note:
557 This routine has a different interface for Fortran
558 $ #include <petsc/finclude/petscmat.h>
559 $ use petscmat
560 $ PetscInt, pointer :: array(:)
561 $ PetscErrorCode ierr
562 $ MatFDColoring i
563 $ call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
564 $ use the entries of array ...
565 $ call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)
566
567 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
568
569 @*/
MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt * n,const PetscInt * cols[])570 PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,const PetscInt *cols[])
571 {
572 PetscFunctionBegin;
573 if (coloring->currentcolor >= 0) {
574 *n = coloring->ncolumns[coloring->currentcolor];
575 *cols = coloring->columns[coloring->currentcolor];
576 } else {
577 *n = 0;
578 }
579 PetscFunctionReturn(0);
580 }
581
582 /*@
583 MatFDColoringApply - Given a matrix for which a MatFDColoring context
584 has been created, computes the Jacobian for a function via finite differences.
585
586 Collective on MatFDColoring
587
588 Input Parameters:
589 + mat - location to store Jacobian
590 . coloring - coloring context created with MatFDColoringCreate()
591 . x1 - location at which Jacobian is to be computed
592 - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
593
594 Options Database Keys:
595 + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
596 . -mat_fd_coloring_view - Activates basic viewing or coloring
597 . -mat_fd_coloring_view draw - Activates drawing of coloring
598 - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
599
600 Level: intermediate
601
602 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction(), MatFDColoringSetValues()
603
604 @*/
MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void * sctx)605 PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
606 {
607 PetscErrorCode ierr;
608 PetscBool eq;
609
610 PetscFunctionBegin;
611 PetscValidHeaderSpecific(J,MAT_CLASSID,1);
612 PetscValidHeaderSpecific(coloring,MAT_FDCOLORING_CLASSID,2);
613 PetscValidHeaderSpecific(x1,VEC_CLASSID,3);
614 ierr = PetscObjectCompareId((PetscObject)J,coloring->matid,&eq);CHKERRQ(ierr);
615 if (!eq) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Matrix used with MatFDColoringApply() must be that used with MatFDColoringCreate()");
616 if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
617 if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
618 if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");
619
620 ierr = MatSetUnfactored(J);CHKERRQ(ierr);
621 ierr = PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);CHKERRQ(ierr);
622 ierr = (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);CHKERRQ(ierr);
623 ierr = PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);CHKERRQ(ierr);
624 if (!coloring->viewed) {
625 ierr = MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");CHKERRQ(ierr);
626 coloring->viewed = PETSC_TRUE;
627 }
628 PetscFunctionReturn(0);
629 }
630