1
2 #include <petsc/private/matimpl.h>
3 #include <../src/mat/impls/mffd/mffdimpl.h> /*I "petscmat.h" I*/
4
5 PetscFunctionList MatMFFDList = NULL;
6 PetscBool MatMFFDRegisterAllCalled = PETSC_FALSE;
7
8 PetscClassId MATMFFD_CLASSID;
9 PetscLogEvent MATMFFD_Mult;
10
11 static PetscBool MatMFFDPackageInitialized = PETSC_FALSE;
12 /*@C
13 MatMFFDFinalizePackage - This function destroys everything in the MatMFFD package. It is
14 called from PetscFinalize().
15
16 Level: developer
17
18 .seealso: PetscFinalize(), MatCreateMFFD(), MatCreateSNESMF()
19 @*/
MatMFFDFinalizePackage(void)20 PetscErrorCode MatMFFDFinalizePackage(void)
21 {
22 PetscErrorCode ierr;
23
24 PetscFunctionBegin;
25 ierr = PetscFunctionListDestroy(&MatMFFDList);CHKERRQ(ierr);
26 MatMFFDPackageInitialized = PETSC_FALSE;
27 MatMFFDRegisterAllCalled = PETSC_FALSE;
28 PetscFunctionReturn(0);
29 }
30
31 /*@C
32 MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
33 from MatInitializePackage().
34
35 Level: developer
36
37 .seealso: PetscInitialize()
38 @*/
MatMFFDInitializePackage(void)39 PetscErrorCode MatMFFDInitializePackage(void)
40 {
41 char logList[256];
42 PetscBool opt,pkg;
43 PetscErrorCode ierr;
44
45 PetscFunctionBegin;
46 if (MatMFFDPackageInitialized) PetscFunctionReturn(0);
47 MatMFFDPackageInitialized = PETSC_TRUE;
48 /* Register Classes */
49 ierr = PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);CHKERRQ(ierr);
50 /* Register Constructors */
51 ierr = MatMFFDRegisterAll();CHKERRQ(ierr);
52 /* Register Events */
53 ierr = PetscLogEventRegister("MatMult MF",MATMFFD_CLASSID,&MATMFFD_Mult);CHKERRQ(ierr);
54 /* Process Info */
55 {
56 PetscClassId classids[1];
57
58 classids[0] = MATMFFD_CLASSID;
59 ierr = PetscInfoProcessClass("matmffd", 1, classids);CHKERRQ(ierr);
60 }
61 /* Process summary exclusions */
62 ierr = PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);CHKERRQ(ierr);
63 if (opt) {
64 ierr = PetscStrInList("matmffd",logList,',',&pkg);CHKERRQ(ierr);
65 if (pkg) {ierr = PetscLogEventExcludeClass(MATMFFD_CLASSID);CHKERRQ(ierr);}
66 }
67 /* Register package finalizer */
68 ierr = PetscRegisterFinalize(MatMFFDFinalizePackage);CHKERRQ(ierr);
69 PetscFunctionReturn(0);
70 }
71
MatMFFDSetType_MFFD(Mat mat,MatMFFDType ftype)72 static PetscErrorCode MatMFFDSetType_MFFD(Mat mat,MatMFFDType ftype)
73 {
74 PetscErrorCode ierr,(*r)(MatMFFD);
75 MatMFFD ctx;
76 PetscBool match;
77
78 PetscFunctionBegin;
79 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
80 PetscValidCharPointer(ftype,2);
81 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
82
83 /* already set, so just return */
84 ierr = PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);CHKERRQ(ierr);
85 if (match) PetscFunctionReturn(0);
86
87 /* destroy the old one if it exists */
88 if (ctx->ops->destroy) {
89 ierr = (*ctx->ops->destroy)(ctx);CHKERRQ(ierr);
90 }
91
92 ierr = PetscFunctionListFind(MatMFFDList,ftype,&r);CHKERRQ(ierr);
93 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
94 ierr = (*r)(ctx);CHKERRQ(ierr);
95 ierr = PetscObjectChangeTypeName((PetscObject)ctx,ftype);CHKERRQ(ierr);
96 PetscFunctionReturn(0);
97 }
98
99 /*@C
100 MatMFFDSetType - Sets the method that is used to compute the
101 differencing parameter for finite differene matrix-free formulations.
102
103 Input Parameters:
104 + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
105 or MatSetType(mat,MATMFFD);
106 - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
107
108 Level: advanced
109
110 Notes:
111 For example, such routines can compute h for use in
112 Jacobian-vector products of the form
113
114 F(x+ha) - F(x)
115 F'(u)a ~= ----------------
116 h
117
118 .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
119 @*/
MatMFFDSetType(Mat mat,MatMFFDType ftype)120 PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
121 {
122 PetscErrorCode ierr;
123
124 PetscFunctionBegin;
125 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
126 PetscValidCharPointer(ftype,2);
127 ierr = PetscTryMethod(mat,"MatMFFDSetType_C",(Mat,MatMFFDType),(mat,ftype));CHKERRQ(ierr);
128 PetscFunctionReturn(0);
129 }
130
131 static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);
132
133 typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)134 static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
135 {
136 MatMFFD ctx;
137 PetscErrorCode ierr;
138
139 PetscFunctionBegin;
140 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
141 ctx->funcisetbase = func;
142 PetscFunctionReturn(0);
143 }
144
145 typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)146 static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
147 {
148 MatMFFD ctx;
149 PetscErrorCode ierr;
150
151 PetscFunctionBegin;
152 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
153 ctx->funci = funci;
154 ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_MFFD);CHKERRQ(ierr);
155 PetscFunctionReturn(0);
156 }
157
MatMFFDGetH_MFFD(Mat mat,PetscScalar * h)158 static PetscErrorCode MatMFFDGetH_MFFD(Mat mat,PetscScalar *h)
159 {
160 MatMFFD ctx;
161 PetscErrorCode ierr;
162
163 PetscFunctionBegin;
164 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
165 *h = ctx->currenth;
166 PetscFunctionReturn(0);
167 }
168
MatMFFDResetHHistory_MFFD(Mat J)169 static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
170 {
171 MatMFFD ctx;
172 PetscErrorCode ierr;
173
174 PetscFunctionBegin;
175 ierr = MatShellGetContext(J,&ctx);CHKERRQ(ierr);
176 ctx->ncurrenth = 0;
177 PetscFunctionReturn(0);
178 }
179
180 /*@C
181 MatMFFDRegister - Adds a method to the MatMFFD registry.
182
183 Not Collective
184
185 Input Parameters:
186 + name_solver - name of a new user-defined compute-h module
187 - routine_create - routine to create method context
188
189 Level: developer
190
191 Notes:
192 MatMFFDRegister() may be called multiple times to add several user-defined solvers.
193
194 Sample usage:
195 .vb
196 MatMFFDRegister("my_h",MyHCreate);
197 .ve
198
199 Then, your solver can be chosen with the procedural interface via
200 $ MatMFFDSetType(mfctx,"my_h")
201 or at runtime via the option
202 $ -mat_mffd_type my_h
203
204 .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
205 @*/
MatMFFDRegister(const char sname[],PetscErrorCode (* function)(MatMFFD))206 PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
207 {
208 PetscErrorCode ierr;
209
210 PetscFunctionBegin;
211 ierr = MatInitializePackage();CHKERRQ(ierr);
212 ierr = PetscFunctionListAdd(&MatMFFDList,sname,function);CHKERRQ(ierr);
213 PetscFunctionReturn(0);
214 }
215
216 /* ----------------------------------------------------------------------------------------*/
MatDestroy_MFFD(Mat mat)217 static PetscErrorCode MatDestroy_MFFD(Mat mat)
218 {
219 PetscErrorCode ierr;
220 MatMFFD ctx;
221
222 PetscFunctionBegin;
223 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
224 ierr = VecDestroy(&ctx->w);CHKERRQ(ierr);
225 ierr = VecDestroy(&ctx->current_u);CHKERRQ(ierr);
226 if (ctx->current_f_allocated) {
227 ierr = VecDestroy(&ctx->current_f);CHKERRQ(ierr);
228 }
229 if (ctx->ops->destroy) {ierr = (*ctx->ops->destroy)(ctx);CHKERRQ(ierr);}
230 ierr = PetscHeaderDestroy(&ctx);CHKERRQ(ierr);
231
232 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);CHKERRQ(ierr);
233 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);CHKERRQ(ierr);
234 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);CHKERRQ(ierr);
235 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);CHKERRQ(ierr);
236 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);CHKERRQ(ierr);
237 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);CHKERRQ(ierr);
238 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);CHKERRQ(ierr);
239 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);CHKERRQ(ierr);
240 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetHHistory_C",NULL);CHKERRQ(ierr);
241 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetType_C",NULL);CHKERRQ(ierr);
242 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMFFDGetH_C",NULL);CHKERRQ(ierr);
243 PetscFunctionReturn(0);
244 }
245
246 /*
247 MatMFFDView_MFFD - Views matrix-free parameters.
248
249 */
MatView_MFFD(Mat J,PetscViewer viewer)250 static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
251 {
252 PetscErrorCode ierr;
253 MatMFFD ctx;
254 PetscBool iascii, viewbase, viewfunction;
255 const char *prefix;
256
257 PetscFunctionBegin;
258 ierr = MatShellGetContext(J,&ctx);CHKERRQ(ierr);
259 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
260 if (iascii) {
261 ierr = PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");CHKERRQ(ierr);
262 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
263 ierr = PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);CHKERRQ(ierr);
264 if (!((PetscObject)ctx)->type_name) {
265 ierr = PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");CHKERRQ(ierr);
266 } else {
267 ierr = PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);CHKERRQ(ierr);
268 }
269 #if defined(PETSC_USE_COMPLEX)
270 if (ctx->usecomplex) {
271 ierr = PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");CHKERRQ(ierr);
272 }
273 #endif
274 if (ctx->ops->view) {
275 ierr = (*ctx->ops->view)(ctx,viewer);CHKERRQ(ierr);
276 }
277 ierr = PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);CHKERRQ(ierr);
278
279 ierr = PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);CHKERRQ(ierr);
280 if (viewbase) {
281 ierr = PetscViewerASCIIPrintf(viewer, "Base:\n");CHKERRQ(ierr);
282 ierr = VecView(ctx->current_u, viewer);CHKERRQ(ierr);
283 }
284 ierr = PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);CHKERRQ(ierr);
285 if (viewfunction) {
286 ierr = PetscViewerASCIIPrintf(viewer, "Function:\n");CHKERRQ(ierr);
287 ierr = VecView(ctx->current_f, viewer);CHKERRQ(ierr);
288 }
289 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
290 }
291 PetscFunctionReturn(0);
292 }
293
294 /*
295 MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
296 allows the user to indicate the beginning of a new linear solve by calling
297 MatAssemblyXXX() on the matrix free matrix. This then allows the
298 MatCreateMFFD_WP() to properly compute ||U|| only the first time
299 in the linear solver rather than every time.
300
301 This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
302 it must be labeled as PETSC_EXTERN
303 */
MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)304 PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
305 {
306 PetscErrorCode ierr;
307 MatMFFD j;
308
309 PetscFunctionBegin;
310 ierr = MatShellGetContext(J,&j);CHKERRQ(ierr);
311 ierr = MatMFFDResetHHistory(J);CHKERRQ(ierr);
312 PetscFunctionReturn(0);
313 }
314
315 /*
316 MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
317
318 y ~= (F(u + ha) - F(u))/h,
319 where F = nonlinear function, as set by SNESSetFunction()
320 u = current iterate
321 h = difference interval
322 */
MatMult_MFFD(Mat mat,Vec a,Vec y)323 static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
324 {
325 MatMFFD ctx;
326 PetscScalar h;
327 Vec w,U,F;
328 PetscErrorCode ierr;
329 PetscBool zeroa;
330
331 PetscFunctionBegin;
332 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
333 if (!ctx->current_u) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MatMFFDSetBase() has not been called, this is often caused by forgetting to call \n\t\tMatAssemblyBegin/End on the first Mat in the SNES compute function");
334 /* We log matrix-free matrix-vector products separately, so that we can
335 separate the performance monitoring from the cases that use conventional
336 storage. We may eventually modify event logging to associate events
337 with particular objects, hence alleviating the more general problem. */
338 ierr = PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr);
339
340 w = ctx->w;
341 U = ctx->current_u;
342 F = ctx->current_f;
343 /*
344 Compute differencing parameter
345 */
346 if (!((PetscObject)ctx)->type_name) {
347 ierr = MatMFFDSetType(mat,MATMFFD_WP);CHKERRQ(ierr);
348 ierr = MatSetFromOptions(mat);CHKERRQ(ierr);
349 }
350 ierr = (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);CHKERRQ(ierr);
351 if (zeroa) {
352 ierr = VecSet(y,0.0);CHKERRQ(ierr);
353 PetscFunctionReturn(0);
354 }
355
356 if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
357 if (ctx->checkh) {
358 ierr = (*ctx->checkh)(ctx->checkhctx,U,a,&h);CHKERRQ(ierr);
359 }
360
361 /* keep a record of the current differencing parameter h */
362 ctx->currenth = h;
363 #if defined(PETSC_USE_COMPLEX)
364 ierr = PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));CHKERRQ(ierr);
365 #else
366 ierr = PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);CHKERRQ(ierr);
367 #endif
368 if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
369 ctx->historyh[ctx->ncurrenth] = h;
370 }
371 ctx->ncurrenth++;
372
373 #if defined(PETSC_USE_COMPLEX)
374 if (ctx->usecomplex) h = PETSC_i*h;
375 #endif
376
377 /* w = u + ha */
378 ierr = VecWAXPY(w,h,a,U);CHKERRQ(ierr);
379
380 /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
381 if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
382 ierr = (*ctx->func)(ctx->funcctx,U,F);CHKERRQ(ierr);
383 }
384 ierr = (*ctx->func)(ctx->funcctx,w,y);CHKERRQ(ierr);
385
386 #if defined(PETSC_USE_COMPLEX)
387 if (ctx->usecomplex) {
388 ierr = VecImaginaryPart(y);CHKERRQ(ierr);
389 h = PetscImaginaryPart(h);
390 } else {
391 ierr = VecAXPY(y,-1.0,F);CHKERRQ(ierr);
392 }
393 #else
394 ierr = VecAXPY(y,-1.0,F);CHKERRQ(ierr);
395 #endif
396 ierr = VecScale(y,1.0/h);CHKERRQ(ierr);
397 if (mat->nullsp) {ierr = MatNullSpaceRemove(mat->nullsp,y);CHKERRQ(ierr);}
398
399 ierr = PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr);
400 PetscFunctionReturn(0);
401 }
402
403 /*
404 MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
405
406 y ~= (F(u + ha) - F(u))/h,
407 where F = nonlinear function, as set by SNESSetFunction()
408 u = current iterate
409 h = difference interval
410 */
MatGetDiagonal_MFFD(Mat mat,Vec a)411 PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
412 {
413 MatMFFD ctx;
414 PetscScalar h,*aa,*ww,v;
415 PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
416 Vec w,U;
417 PetscErrorCode ierr;
418 PetscInt i,rstart,rend;
419
420 PetscFunctionBegin;
421 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
422 if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
423 if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
424 w = ctx->w;
425 U = ctx->current_u;
426 ierr = (*ctx->func)(ctx->funcctx,U,a);CHKERRQ(ierr);
427 if (ctx->funcisetbase) {
428 ierr = (*ctx->funcisetbase)(ctx->funcctx,U);CHKERRQ(ierr);
429 }
430 ierr = VecCopy(U,w);CHKERRQ(ierr);
431
432 ierr = VecGetOwnershipRange(a,&rstart,&rend);CHKERRQ(ierr);
433 ierr = VecGetArray(a,&aa);CHKERRQ(ierr);
434 for (i=rstart; i<rend; i++) {
435 ierr = VecGetArray(w,&ww);CHKERRQ(ierr);
436 h = ww[i-rstart];
437 if (h == 0.0) h = 1.0;
438 if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
439 else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
440 h *= epsilon;
441
442 ww[i-rstart] += h;
443 ierr = VecRestoreArray(w,&ww);CHKERRQ(ierr);
444 ierr = (*ctx->funci)(ctx->funcctx,i,w,&v);CHKERRQ(ierr);
445 aa[i-rstart] = (v - aa[i-rstart])/h;
446
447 ierr = VecGetArray(w,&ww);CHKERRQ(ierr);
448 ww[i-rstart] -= h;
449 ierr = VecRestoreArray(w,&ww);CHKERRQ(ierr);
450 }
451 ierr = VecRestoreArray(a,&aa);CHKERRQ(ierr);
452 PetscFunctionReturn(0);
453 }
454
MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)455 PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
456 {
457 PetscErrorCode ierr;
458 MatMFFD ctx;
459
460 PetscFunctionBegin;
461 ierr = MatShellGetContext(J,&ctx);CHKERRQ(ierr);
462 ierr = MatMFFDResetHHistory(J);CHKERRQ(ierr);
463 if (!ctx->current_u) {
464 ierr = VecDuplicate(U,&ctx->current_u);CHKERRQ(ierr);
465 ierr = VecLockReadPush(ctx->current_u);CHKERRQ(ierr);
466 }
467 ierr = VecLockReadPop(ctx->current_u);CHKERRQ(ierr);
468 ierr = VecCopy(U,ctx->current_u);CHKERRQ(ierr);
469 ierr = VecLockReadPush(ctx->current_u);CHKERRQ(ierr);
470 if (F) {
471 if (ctx->current_f_allocated) {ierr = VecDestroy(&ctx->current_f);CHKERRQ(ierr);}
472 ctx->current_f = F;
473 ctx->current_f_allocated = PETSC_FALSE;
474 } else if (!ctx->current_f_allocated) {
475 ierr = MatCreateVecs(J,NULL,&ctx->current_f);CHKERRQ(ierr);
476
477 ctx->current_f_allocated = PETSC_TRUE;
478 }
479 if (!ctx->w) {
480 ierr = VecDuplicate(ctx->current_u,&ctx->w);CHKERRQ(ierr);
481 }
482 J->assembled = PETSC_TRUE;
483 PetscFunctionReturn(0);
484 }
485
486 typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
487
MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void * ectx)488 static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
489 {
490 MatMFFD ctx;
491 PetscErrorCode ierr;
492
493 PetscFunctionBegin;
494 ierr = MatShellGetContext(J,&ctx);CHKERRQ(ierr);
495 ctx->checkh = fun;
496 ctx->checkhctx = ectx;
497 PetscFunctionReturn(0);
498 }
499
500 /*@C
501 MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
502 MatMFFD options in the database.
503
504 Collective on Mat
505
506 Input Parameter:
507 + A - the Mat context
508 - prefix - the prefix to prepend to all option names
509
510 Notes:
511 A hyphen (-) must NOT be given at the beginning of the prefix name.
512 The first character of all runtime options is AUTOMATICALLY the hyphen.
513
514 Level: advanced
515
516 .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
517 @*/
MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])518 PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
519
520 {
521 MatMFFD mfctx;
522 PetscErrorCode ierr;
523
524 PetscFunctionBegin;
525 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
526 ierr = MatShellGetContext(mat,&mfctx);CHKERRQ(ierr);
527 PetscValidHeaderSpecific(mfctx,MATMFFD_CLASSID,1);
528 ierr = PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);CHKERRQ(ierr);
529 PetscFunctionReturn(0);
530 }
531
MatSetFromOptions_MFFD(PetscOptionItems * PetscOptionsObject,Mat mat)532 static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
533 {
534 MatMFFD mfctx;
535 PetscErrorCode ierr;
536 PetscBool flg;
537 char ftype[256];
538
539 PetscFunctionBegin;
540 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
541 ierr = MatShellGetContext(mat,&mfctx);CHKERRQ(ierr);
542 PetscValidHeaderSpecific(mfctx,MATMFFD_CLASSID,1);
543 ierr = PetscObjectOptionsBegin((PetscObject)mfctx);CHKERRQ(ierr);
544 ierr = PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);CHKERRQ(ierr);
545 if (flg) {
546 ierr = MatMFFDSetType(mat,ftype);CHKERRQ(ierr);
547 }
548
549 ierr = PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,NULL);CHKERRQ(ierr);
550 ierr = PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,NULL);CHKERRQ(ierr);
551
552 flg = PETSC_FALSE;
553 ierr = PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);CHKERRQ(ierr);
554 if (flg) {
555 ierr = MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,NULL);CHKERRQ(ierr);
556 }
557 #if defined(PETSC_USE_COMPLEX)
558 ierr = PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);CHKERRQ(ierr);
559 #endif
560 if (mfctx->ops->setfromoptions) {
561 ierr = (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);CHKERRQ(ierr);
562 }
563 ierr = PetscOptionsEnd();CHKERRQ(ierr);
564 PetscFunctionReturn(0);
565 }
566
MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)567 static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
568 {
569 MatMFFD ctx;
570 PetscErrorCode ierr;
571
572 PetscFunctionBegin;
573 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
574 ctx->recomputeperiod = period;
575 PetscFunctionReturn(0);
576 }
577
MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (* func)(void *,Vec,Vec),void * funcctx)578 static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
579 {
580 MatMFFD ctx;
581 PetscErrorCode ierr;
582
583 PetscFunctionBegin;
584 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
585 ctx->func = func;
586 ctx->funcctx = funcctx;
587 PetscFunctionReturn(0);
588 }
589
MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)590 static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
591 {
592 MatMFFD ctx;
593 PetscErrorCode ierr;
594
595 PetscFunctionBegin;
596 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
597 if (error != PETSC_DEFAULT) ctx->error_rel = error;
598 PetscFunctionReturn(0);
599 }
600
MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)601 PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)
602 {
603 MatMFFD ctx;
604 PetscErrorCode ierr;
605
606 PetscFunctionBegin;
607 ierr = MatShellGetContext(J,&ctx);CHKERRQ(ierr);
608 ctx->historyh = history;
609 ctx->maxcurrenth = nhistory;
610 ctx->currenth = 0.;
611 PetscFunctionReturn(0);
612 }
613
614 /*MC
615 MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
616
617 Level: advanced
618
619 Developers Note: This is implemented on top of MATSHELL to get support for scaling and shifting without requiring duplicate code
620
621 .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
622 MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
623 MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
624 MatMFFDGetH(),
625 M*/
MatCreate_MFFD(Mat A)626 PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
627 {
628 MatMFFD mfctx;
629 PetscErrorCode ierr;
630
631 PetscFunctionBegin;
632 ierr = MatMFFDInitializePackage();CHKERRQ(ierr);
633
634 ierr = PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),NULL,NULL);CHKERRQ(ierr);
635
636 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
637 mfctx->recomputeperiod = 1;
638 mfctx->count = 0;
639 mfctx->currenth = 0.0;
640 mfctx->historyh = NULL;
641 mfctx->ncurrenth = 0;
642 mfctx->maxcurrenth = 0;
643 ((PetscObject)mfctx)->type_name = NULL;
644
645 /*
646 Create the empty data structure to contain compute-h routines.
647 These will be filled in below from the command line options or
648 a later call with MatMFFDSetType() or if that is not called
649 then it will default in the first use of MatMult_MFFD()
650 */
651 mfctx->ops->compute = NULL;
652 mfctx->ops->destroy = NULL;
653 mfctx->ops->view = NULL;
654 mfctx->ops->setfromoptions = NULL;
655 mfctx->hctx = NULL;
656
657 mfctx->func = NULL;
658 mfctx->funcctx = NULL;
659 mfctx->w = NULL;
660 mfctx->mat = A;
661
662 ierr = MatSetType(A,MATSHELL);CHKERRQ(ierr);
663 ierr = MatShellSetContext(A,mfctx);CHKERRQ(ierr);
664 ierr = MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_MFFD);CHKERRQ(ierr);
665 ierr = MatShellSetOperation(A,MATOP_DESTROY,(void (*)(void))MatDestroy_MFFD);CHKERRQ(ierr);
666 ierr = MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MFFD);CHKERRQ(ierr);
667 ierr = MatShellSetOperation(A,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MFFD);CHKERRQ(ierr);
668 ierr = MatShellSetOperation(A,MATOP_SET_FROM_OPTIONS,(void (*)(void))MatSetFromOptions_MFFD);CHKERRQ(ierr);
669
670 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);CHKERRQ(ierr);
671 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);CHKERRQ(ierr);
672 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);CHKERRQ(ierr);
673 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);CHKERRQ(ierr);
674 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);CHKERRQ(ierr);
675 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);CHKERRQ(ierr);
676 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);CHKERRQ(ierr);
677 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);CHKERRQ(ierr);
678 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetHHistory_C",MatMFFDSetHHistory_MFFD);CHKERRQ(ierr);
679 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetType_C",MatMFFDSetType_MFFD);CHKERRQ(ierr);
680 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDGetH_C",MatMFFDGetH_MFFD);CHKERRQ(ierr);
681 ierr = PetscObjectChangeTypeName((PetscObject)A,MATMFFD);CHKERRQ(ierr);
682 PetscFunctionReturn(0);
683 }
684
685 /*@
686 MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
687
688 Collective on Vec
689
690 Input Parameters:
691 + comm - MPI communicator
692 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
693 This value should be the same as the local size used in creating the
694 y vector for the matrix-vector product y = Ax.
695 . n - This value should be the same as the local size used in creating the
696 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
697 calculated if N is given) For square matrices n is almost always m.
698 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
699 - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
700
701
702 Output Parameter:
703 . J - the matrix-free matrix
704
705 Options Database Keys: call MatSetFromOptions() to trigger these
706 + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
707 . -mat_mffd_err - square root of estimated relative error in function evaluation
708 . -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
709 . -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
710 - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
711 (requires real valued functions but that PETSc be configured for complex numbers)
712
713
714 Level: advanced
715
716 Notes:
717 The matrix-free matrix context merely contains the function pointers
718 and work space for performing finite difference approximations of
719 Jacobian-vector products, F'(u)*a,
720
721 The default code uses the following approach to compute h
722
723 .vb
724 F'(u)*a = [F(u+h*a) - F(u)]/h where
725 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
726 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
727 where
728 error_rel = square root of relative error in function evaluation
729 umin = minimum iterate parameter
730 .ve
731
732 You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
733 preconditioner matrix
734
735 The user can set the error_rel via MatMFFDSetFunctionError() and
736 umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
737
738 The user should call MatDestroy() when finished with the matrix-free
739 matrix context.
740
741 Options Database Keys:
742 + -mat_mffd_err <error_rel> - Sets error_rel
743 . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
744 - -mat_mffd_check_positivity
745
746 .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
747 MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
748 MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
749
750 @*/
MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat * J)751 PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
752 {
753 PetscErrorCode ierr;
754
755 PetscFunctionBegin;
756 ierr = MatCreate(comm,J);CHKERRQ(ierr);
757 ierr = MatSetSizes(*J,m,n,M,N);CHKERRQ(ierr);
758 ierr = MatSetType(*J,MATMFFD);CHKERRQ(ierr);
759 ierr = MatSetUp(*J);CHKERRQ(ierr);
760 PetscFunctionReturn(0);
761 }
762
763 /*@
764 MatMFFDGetH - Gets the last value that was used as the differencing
765 parameter.
766
767 Not Collective
768
769 Input Parameters:
770 . mat - the matrix obtained with MatCreateSNESMF()
771
772 Output Parameter:
773 . h - the differencing step size
774
775 Level: advanced
776
777 .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
778 @*/
MatMFFDGetH(Mat mat,PetscScalar * h)779 PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
780 {
781 PetscErrorCode ierr;
782
783 PetscFunctionBegin;
784 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
785 PetscValidPointer(h,2);
786 ierr = PetscUseMethod(mat,"MatMFFDGetH_C",(Mat,PetscScalar*),(mat,h));CHKERRQ(ierr);
787 PetscFunctionReturn(0);
788 }
789
790 /*@C
791 MatMFFDSetFunction - Sets the function used in applying the matrix free.
792
793 Logically Collective on Mat
794
795 Input Parameters:
796 + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
797 . func - the function to use
798 - funcctx - optional function context passed to function
799
800 Calling Sequence of func:
801 $ func (void *funcctx, Vec x, Vec f)
802
803 + funcctx - user provided context
804 . x - input vector
805 - f - computed output function
806
807 Level: advanced
808
809 Notes:
810 If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
811 matrix inside your compute Jacobian routine
812
813 If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
814
815 .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
816 MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
817 @*/
MatMFFDSetFunction(Mat mat,PetscErrorCode (* func)(void *,Vec,Vec),void * funcctx)818 PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
819 {
820 PetscErrorCode ierr;
821
822 PetscFunctionBegin;
823 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
824 ierr = PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));CHKERRQ(ierr);
825 PetscFunctionReturn(0);
826 }
827
828 /*@C
829 MatMFFDSetFunctioni - Sets the function for a single component
830
831 Logically Collective on Mat
832
833 Input Parameters:
834 + mat - the matrix free matrix created via MatCreateSNESMF()
835 - funci - the function to use
836
837 Level: advanced
838
839 Notes:
840 If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
841 matrix inside your compute Jacobian routine.
842 This function is necessary to compute the diagonal of the matrix.
843 funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
844
845 .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
846
847 @*/
MatMFFDSetFunctioni(Mat mat,PetscErrorCode (* funci)(void *,PetscInt,Vec,PetscScalar *))848 PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
849 {
850 PetscErrorCode ierr;
851
852 PetscFunctionBegin;
853 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
854 ierr = PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));CHKERRQ(ierr);
855 PetscFunctionReturn(0);
856 }
857
858 /*@C
859 MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
860
861 Logically Collective on Mat
862
863 Input Parameters:
864 + mat - the matrix free matrix created via MatCreateSNESMF()
865 - func - the function to use
866
867 Level: advanced
868
869 Notes:
870 If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
871 matrix inside your compute Jacobian routine.
872 This function is necessary to compute the diagonal of the matrix.
873
874
875 .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
876 MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
877 @*/
MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (* func)(void *,Vec))878 PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
879 {
880 PetscErrorCode ierr;
881
882 PetscFunctionBegin;
883 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
884 ierr = PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));CHKERRQ(ierr);
885 PetscFunctionReturn(0);
886 }
887
888 /*@
889 MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
890
891 Logically Collective on Mat
892
893 Input Parameters:
894 + mat - the matrix free matrix created via MatCreateSNESMF()
895 - period - 1 for everytime, 2 for every second etc
896
897 Options Database Keys:
898 . -mat_mffd_period <period>
899
900 Level: advanced
901
902
903 .seealso: MatCreateSNESMF(),MatMFFDGetH(),
904 MatMFFDSetHHistory(), MatMFFDResetHHistory()
905 @*/
MatMFFDSetPeriod(Mat mat,PetscInt period)906 PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
907 {
908 PetscErrorCode ierr;
909
910 PetscFunctionBegin;
911 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
912 PetscValidLogicalCollectiveInt(mat,period,2);
913 ierr = PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));CHKERRQ(ierr);
914 PetscFunctionReturn(0);
915 }
916
917 /*@
918 MatMFFDSetFunctionError - Sets the error_rel for the approximation of
919 matrix-vector products using finite differences.
920
921 Logically Collective on Mat
922
923 Input Parameters:
924 + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
925 - error_rel - relative error (should be set to the square root of
926 the relative error in the function evaluations)
927
928 Options Database Keys:
929 . -mat_mffd_err <error_rel> - Sets error_rel
930
931 Level: advanced
932
933 Notes:
934 The default matrix-free matrix-vector product routine computes
935 .vb
936 F'(u)*a = [F(u+h*a) - F(u)]/h where
937 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
938 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
939 .ve
940
941 .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
942 MatMFFDSetHHistory(), MatMFFDResetHHistory()
943 @*/
MatMFFDSetFunctionError(Mat mat,PetscReal error)944 PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
945 {
946 PetscErrorCode ierr;
947
948 PetscFunctionBegin;
949 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
950 PetscValidLogicalCollectiveReal(mat,error,2);
951 ierr = PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));CHKERRQ(ierr);
952 PetscFunctionReturn(0);
953 }
954
955 /*@
956 MatMFFDSetHHistory - Sets an array to collect a history of the
957 differencing values (h) computed for the matrix-free product.
958
959 Logically Collective on Mat
960
961 Input Parameters:
962 + J - the matrix-free matrix context
963 . histroy - space to hold the history
964 - nhistory - number of entries in history, if more entries are generated than
965 nhistory, then the later ones are discarded
966
967 Level: advanced
968
969 Notes:
970 Use MatMFFDResetHHistory() to reset the history counter and collect
971 a new batch of differencing parameters, h.
972
973 .seealso: MatMFFDGetH(), MatCreateSNESMF(),
974 MatMFFDResetHHistory(), MatMFFDSetFunctionError()
975
976 @*/
MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)977 PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
978 {
979 PetscErrorCode ierr;
980
981 PetscFunctionBegin;
982 PetscValidHeaderSpecific(J,MAT_CLASSID,1);
983 if (history) PetscValidPointer(history,2);
984 PetscValidLogicalCollectiveInt(J,nhistory,3);
985 ierr = PetscUseMethod(J,"MatMFFDSetHHistory_C",(Mat,PetscScalar[],PetscInt),(J,history,nhistory));CHKERRQ(ierr);
986 PetscFunctionReturn(0);
987 }
988
989 /*@
990 MatMFFDResetHHistory - Resets the counter to zero to begin
991 collecting a new set of differencing histories.
992
993 Logically Collective on Mat
994
995 Input Parameters:
996 . J - the matrix-free matrix context
997
998 Level: advanced
999
1000 Notes:
1001 Use MatMFFDSetHHistory() to create the original history counter.
1002
1003 .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1004 MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1005
1006 @*/
MatMFFDResetHHistory(Mat J)1007 PetscErrorCode MatMFFDResetHHistory(Mat J)
1008 {
1009 PetscErrorCode ierr;
1010
1011 PetscFunctionBegin;
1012 PetscValidHeaderSpecific(J,MAT_CLASSID,1);
1013 ierr = PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));CHKERRQ(ierr);
1014 PetscFunctionReturn(0);
1015 }
1016
1017 /*@
1018 MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1019 Jacobian are computed
1020
1021 Logically Collective on Mat
1022
1023 Input Parameters:
1024 + J - the MatMFFD matrix
1025 . U - the vector
1026 - F - (optional) vector that contains F(u) if it has been already computed
1027
1028 Notes:
1029 This is rarely used directly
1030
1031 If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1032 point during the first MatMult() after each call to MatMFFDSetBase().
1033
1034 Level: advanced
1035
1036 @*/
MatMFFDSetBase(Mat J,Vec U,Vec F)1037 PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1038 {
1039 PetscErrorCode ierr;
1040
1041 PetscFunctionBegin;
1042 PetscValidHeaderSpecific(J,MAT_CLASSID,1);
1043 PetscValidHeaderSpecific(U,VEC_CLASSID,2);
1044 if (F) PetscValidHeaderSpecific(F,VEC_CLASSID,3);
1045 ierr = PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));CHKERRQ(ierr);
1046 PetscFunctionReturn(0);
1047 }
1048
1049 /*@C
1050 MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1051 it to satisfy some criteria
1052
1053 Logically Collective on Mat
1054
1055 Input Parameters:
1056 + J - the MatMFFD matrix
1057 . fun - the function that checks h
1058 - ctx - any context needed by the function
1059
1060 Options Database Keys:
1061 . -mat_mffd_check_positivity
1062
1063 Level: advanced
1064
1065 Notes:
1066 For example, MatMFFDCheckPositivity() insures that all entries
1067 of U + h*a are non-negative
1068
1069 The function you provide is called after the default h has been computed and allows you to
1070 modify it.
1071
1072 .seealso: MatMFFDCheckPositivity()
1073 @*/
MatMFFDSetCheckh(Mat J,PetscErrorCode (* fun)(void *,Vec,Vec,PetscScalar *),void * ctx)1074 PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1075 {
1076 PetscErrorCode ierr;
1077
1078 PetscFunctionBegin;
1079 PetscValidHeaderSpecific(J,MAT_CLASSID,1);
1080 ierr = PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));CHKERRQ(ierr);
1081 PetscFunctionReturn(0);
1082 }
1083
1084 /*@
1085 MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1086 zero, decreases h until this is satisfied.
1087
1088 Logically Collective on Vec
1089
1090 Input Parameters:
1091 + U - base vector that is added to
1092 . a - vector that is added
1093 . h - scaling factor on a
1094 - dummy - context variable (unused)
1095
1096 Options Database Keys:
1097 . -mat_mffd_check_positivity
1098
1099 Level: advanced
1100
1101 Notes:
1102 This is rarely used directly, rather it is passed as an argument to
1103 MatMFFDSetCheckh()
1104
1105 .seealso: MatMFFDSetCheckh()
1106 @*/
MatMFFDCheckPositivity(void * dummy,Vec U,Vec a,PetscScalar * h)1107 PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1108 {
1109 PetscReal val, minval;
1110 PetscScalar *u_vec, *a_vec;
1111 PetscErrorCode ierr;
1112 PetscInt i,n;
1113 MPI_Comm comm;
1114
1115 PetscFunctionBegin;
1116 PetscValidHeaderSpecific(U,VEC_CLASSID,2);
1117 PetscValidHeaderSpecific(a,VEC_CLASSID,3);
1118 PetscValidPointer(h,4);
1119 ierr = PetscObjectGetComm((PetscObject)U,&comm);CHKERRQ(ierr);
1120 ierr = VecGetArray(U,&u_vec);CHKERRQ(ierr);
1121 ierr = VecGetArray(a,&a_vec);CHKERRQ(ierr);
1122 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr);
1123 minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1124 for (i=0; i<n; i++) {
1125 if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1126 val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1127 if (val < minval) minval = val;
1128 }
1129 }
1130 ierr = VecRestoreArray(U,&u_vec);CHKERRQ(ierr);
1131 ierr = VecRestoreArray(a,&a_vec);CHKERRQ(ierr);
1132 ierr = MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
1133 if (val <= PetscAbsScalar(*h)) {
1134 ierr = PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));CHKERRQ(ierr);
1135 if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1136 else *h = -0.99*val;
1137 }
1138 PetscFunctionReturn(0);
1139 }
1140