Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : Basic RG routines
12 : */
13 :
14 : #include <slepc/private/rgimpl.h> /*I "slepcrg.h" I*/
15 :
16 : PetscFunctionList RGList = NULL;
17 : PetscBool RGRegisterAllCalled = PETSC_FALSE;
18 : PetscClassId RG_CLASSID = 0;
19 : static PetscBool RGPackageInitialized = PETSC_FALSE;
20 :
21 : /*@C
22 : RGFinalizePackage - This function destroys everything in the Slepc interface
23 : to the RG package. It is called from SlepcFinalize().
24 :
25 : Level: developer
26 :
27 : .seealso: SlepcFinalize()
28 : @*/
29 772 : PetscErrorCode RGFinalizePackage(void)
30 : {
31 772 : PetscFunctionBegin;
32 772 : PetscCall(PetscFunctionListDestroy(&RGList));
33 772 : RGPackageInitialized = PETSC_FALSE;
34 772 : RGRegisterAllCalled = PETSC_FALSE;
35 772 : PetscFunctionReturn(PETSC_SUCCESS);
36 : }
37 :
38 : /*@C
39 : RGInitializePackage - This function initializes everything in the RG package.
40 : It is called from PetscDLLibraryRegister() when using dynamic libraries, and
41 : on the first call to RGCreate() when using static libraries.
42 :
43 : Level: developer
44 :
45 : .seealso: SlepcInitialize()
46 : @*/
47 3995 : PetscErrorCode RGInitializePackage(void)
48 : {
49 3995 : char logList[256];
50 3995 : PetscBool opt,pkg;
51 3995 : PetscClassId classids[1];
52 :
53 3995 : PetscFunctionBegin;
54 3995 : if (RGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
55 772 : RGPackageInitialized = PETSC_TRUE;
56 : /* Register Classes */
57 772 : PetscCall(PetscClassIdRegister("Region",&RG_CLASSID));
58 : /* Register Constructors */
59 772 : PetscCall(RGRegisterAll());
60 : /* Process Info */
61 772 : classids[0] = RG_CLASSID;
62 772 : PetscCall(PetscInfoProcessClass("rg",1,&classids[0]));
63 : /* Process summary exclusions */
64 772 : PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
65 772 : if (opt) {
66 8 : PetscCall(PetscStrInList("rg",logList,',',&pkg));
67 8 : if (pkg) PetscCall(PetscLogEventDeactivateClass(RG_CLASSID));
68 : }
69 : /* Register package finalizer */
70 772 : PetscCall(PetscRegisterFinalize(RGFinalizePackage));
71 772 : PetscFunctionReturn(PETSC_SUCCESS);
72 : }
73 :
74 : /*@
75 : RGCreate - Creates an RG context.
76 :
77 : Collective
78 :
79 : Input Parameter:
80 : . comm - MPI communicator
81 :
82 : Output Parameter:
83 : . newrg - location to put the RG context
84 :
85 : Level: beginner
86 :
87 : .seealso: RGDestroy(), RG
88 : @*/
89 906 : PetscErrorCode RGCreate(MPI_Comm comm,RG *newrg)
90 : {
91 906 : RG rg;
92 :
93 906 : PetscFunctionBegin;
94 906 : PetscAssertPointer(newrg,2);
95 906 : *newrg = NULL;
96 906 : PetscCall(RGInitializePackage());
97 906 : PetscCall(SlepcHeaderCreate(rg,RG_CLASSID,"RG","Region","RG",comm,RGDestroy,RGView));
98 906 : rg->complement = PETSC_FALSE;
99 906 : rg->sfactor = 1.0;
100 906 : rg->osfactor = 0.0;
101 906 : rg->data = NULL;
102 :
103 906 : *newrg = rg;
104 906 : PetscFunctionReturn(PETSC_SUCCESS);
105 : }
106 :
107 : /*@C
108 : RGSetOptionsPrefix - Sets the prefix used for searching for all
109 : RG options in the database.
110 :
111 : Logically Collective
112 :
113 : Input Parameters:
114 : + rg - the region context
115 : - prefix - the prefix string to prepend to all RG option requests
116 :
117 : Notes:
118 : A hyphen (-) must NOT be given at the beginning of the prefix name.
119 : The first character of all runtime options is AUTOMATICALLY the
120 : hyphen.
121 :
122 : Level: advanced
123 :
124 : .seealso: RGAppendOptionsPrefix()
125 : @*/
126 227 : PetscErrorCode RGSetOptionsPrefix(RG rg,const char *prefix)
127 : {
128 227 : PetscFunctionBegin;
129 227 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
130 227 : PetscCall(PetscObjectSetOptionsPrefix((PetscObject)rg,prefix));
131 227 : PetscFunctionReturn(PETSC_SUCCESS);
132 : }
133 :
134 : /*@C
135 : RGAppendOptionsPrefix - Appends to the prefix used for searching for all
136 : RG options in the database.
137 :
138 : Logically Collective
139 :
140 : Input Parameters:
141 : + rg - the region context
142 : - prefix - the prefix string to prepend to all RG option requests
143 :
144 : Notes:
145 : A hyphen (-) must NOT be given at the beginning of the prefix name.
146 : The first character of all runtime options is AUTOMATICALLY the hyphen.
147 :
148 : Level: advanced
149 :
150 : .seealso: RGSetOptionsPrefix()
151 : @*/
152 187 : PetscErrorCode RGAppendOptionsPrefix(RG rg,const char *prefix)
153 : {
154 187 : PetscFunctionBegin;
155 187 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
156 187 : PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)rg,prefix));
157 187 : PetscFunctionReturn(PETSC_SUCCESS);
158 : }
159 :
160 : /*@C
161 : RGGetOptionsPrefix - Gets the prefix used for searching for all
162 : RG options in the database.
163 :
164 : Not Collective
165 :
166 : Input Parameters:
167 : . rg - the region context
168 :
169 : Output Parameters:
170 : . prefix - pointer to the prefix string used is returned
171 :
172 : Note:
173 : On the Fortran side, the user should pass in a string 'prefix' of
174 : sufficient length to hold the prefix.
175 :
176 : Level: advanced
177 :
178 : .seealso: RGSetOptionsPrefix(), RGAppendOptionsPrefix()
179 : @*/
180 0 : PetscErrorCode RGGetOptionsPrefix(RG rg,const char *prefix[])
181 : {
182 0 : PetscFunctionBegin;
183 0 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
184 0 : PetscAssertPointer(prefix,2);
185 0 : PetscCall(PetscObjectGetOptionsPrefix((PetscObject)rg,prefix));
186 0 : PetscFunctionReturn(PETSC_SUCCESS);
187 : }
188 :
189 : /*@C
190 : RGSetType - Selects the type for the RG object.
191 :
192 : Logically Collective
193 :
194 : Input Parameters:
195 : + rg - the region context
196 : - type - a known type
197 :
198 : Level: intermediate
199 :
200 : .seealso: RGGetType()
201 : @*/
202 952 : PetscErrorCode RGSetType(RG rg,RGType type)
203 : {
204 952 : PetscErrorCode (*r)(RG);
205 952 : PetscBool match;
206 :
207 952 : PetscFunctionBegin;
208 952 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
209 952 : PetscAssertPointer(type,2);
210 :
211 952 : PetscCall(PetscObjectTypeCompare((PetscObject)rg,type,&match));
212 952 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
213 :
214 910 : PetscCall(PetscFunctionListFind(RGList,type,&r));
215 910 : PetscCheck(r,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested RG type %s",type);
216 :
217 910 : PetscTryTypeMethod(rg,destroy);
218 910 : PetscCall(PetscMemzero(rg->ops,sizeof(struct _RGOps)));
219 :
220 910 : PetscCall(PetscObjectChangeTypeName((PetscObject)rg,type));
221 910 : PetscCall((*r)(rg));
222 910 : PetscFunctionReturn(PETSC_SUCCESS);
223 : }
224 :
225 : /*@C
226 : RGGetType - Gets the RG type name (as a string) from the RG context.
227 :
228 : Not Collective
229 :
230 : Input Parameter:
231 : . rg - the region context
232 :
233 : Output Parameter:
234 : . type - name of the region
235 :
236 : Level: intermediate
237 :
238 : .seealso: RGSetType()
239 : @*/
240 5 : PetscErrorCode RGGetType(RG rg,RGType *type)
241 : {
242 5 : PetscFunctionBegin;
243 5 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
244 5 : PetscAssertPointer(type,2);
245 5 : *type = ((PetscObject)rg)->type_name;
246 5 : PetscFunctionReturn(PETSC_SUCCESS);
247 : }
248 :
249 : /*@
250 : RGSetFromOptions - Sets RG options from the options database.
251 :
252 : Collective
253 :
254 : Input Parameters:
255 : . rg - the region context
256 :
257 : Notes:
258 : To see all options, run your program with the -help option.
259 :
260 : Level: beginner
261 :
262 : .seealso: RGSetOptionsPrefix()
263 : @*/
264 884 : PetscErrorCode RGSetFromOptions(RG rg)
265 : {
266 884 : char type[256];
267 884 : PetscBool flg;
268 884 : PetscReal sfactor;
269 :
270 884 : PetscFunctionBegin;
271 884 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
272 884 : PetscCall(RGRegisterAll());
273 2652 : PetscObjectOptionsBegin((PetscObject)rg);
274 953 : PetscCall(PetscOptionsFList("-rg_type","Region type","RGSetType",RGList,(char*)(((PetscObject)rg)->type_name?((PetscObject)rg)->type_name:RGINTERVAL),type,sizeof(type),&flg));
275 884 : if (flg) PetscCall(RGSetType(rg,type));
276 812 : else if (!((PetscObject)rg)->type_name) PetscCall(RGSetType(rg,RGINTERVAL));
277 :
278 884 : PetscCall(PetscOptionsBool("-rg_complement","Whether region is complemented or not","RGSetComplement",rg->complement,&rg->complement,NULL));
279 :
280 884 : PetscCall(PetscOptionsReal("-rg_scale","Scaling factor","RGSetScale",1.0,&sfactor,&flg));
281 884 : if (flg) PetscCall(RGSetScale(rg,sfactor));
282 :
283 884 : PetscTryTypeMethod(rg,setfromoptions,PetscOptionsObject);
284 884 : PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)rg,PetscOptionsObject));
285 884 : PetscOptionsEnd();
286 884 : PetscFunctionReturn(PETSC_SUCCESS);
287 : }
288 :
289 : /*@C
290 : RGView - Prints the RG data structure.
291 :
292 : Collective
293 :
294 : Input Parameters:
295 : + rg - the region context
296 : - viewer - optional visualization context
297 :
298 : Note:
299 : The available visualization contexts include
300 : + PETSC_VIEWER_STDOUT_SELF - standard output (default)
301 : - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
302 : output where only the first processor opens
303 : the file. All other processors send their
304 : data to the first processor to print.
305 :
306 : The user can open an alternative visualization context with
307 : PetscViewerASCIIOpen() - output to a specified file.
308 :
309 : Level: beginner
310 :
311 : .seealso: RGCreate()
312 : @*/
313 17 : PetscErrorCode RGView(RG rg,PetscViewer viewer)
314 : {
315 17 : PetscBool isdraw,isascii;
316 :
317 17 : PetscFunctionBegin;
318 17 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
319 17 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rg),&viewer));
320 17 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
321 17 : PetscCheckSameComm(rg,1,viewer,2);
322 17 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
323 17 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
324 17 : if (isascii) {
325 13 : PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rg,viewer));
326 13 : PetscCall(PetscViewerASCIIPushTab(viewer));
327 13 : PetscTryTypeMethod(rg,view,viewer);
328 13 : PetscCall(PetscViewerASCIIPopTab(viewer));
329 13 : if (rg->complement) PetscCall(PetscViewerASCIIPrintf(viewer," selected region is the complement of the specified one\n"));
330 13 : if (rg->sfactor!=1.0) PetscCall(PetscViewerASCIIPrintf(viewer," scaling factor = %g\n",(double)rg->sfactor));
331 4 : } else if (isdraw) PetscTryTypeMethod(rg,view,viewer);
332 17 : PetscFunctionReturn(PETSC_SUCCESS);
333 : }
334 :
335 : /*@C
336 : RGViewFromOptions - View from options
337 :
338 : Collective
339 :
340 : Input Parameters:
341 : + rg - the region context
342 : . obj - optional object
343 : - name - command line option
344 :
345 : Level: intermediate
346 :
347 : .seealso: RGView(), RGCreate()
348 : @*/
349 1168 : PetscErrorCode RGViewFromOptions(RG rg,PetscObject obj,const char name[])
350 : {
351 1168 : PetscFunctionBegin;
352 1168 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
353 1168 : PetscCall(PetscObjectViewFromOptions((PetscObject)rg,obj,name));
354 1168 : PetscFunctionReturn(PETSC_SUCCESS);
355 : }
356 :
357 : /*@
358 : RGIsTrivial - Whether it is the trivial region (whole complex plane).
359 :
360 : Not Collective
361 :
362 : Input Parameter:
363 : . rg - the region context
364 :
365 : Output Parameter:
366 : . trivial - true if the region is equal to the whole complex plane, e.g.,
367 : an interval region with all four endpoints unbounded or an
368 : ellipse with infinite radius.
369 :
370 : Level: beginner
371 :
372 : .seealso: RGCheckInside()
373 : @*/
374 9757 : PetscErrorCode RGIsTrivial(RG rg,PetscBool *trivial)
375 : {
376 9757 : PetscFunctionBegin;
377 9757 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
378 9757 : PetscValidType(rg,1);
379 9757 : PetscAssertPointer(trivial,2);
380 9757 : *trivial = PETSC_FALSE;
381 9757 : PetscTryTypeMethod(rg,istrivial,trivial);
382 9757 : PetscFunctionReturn(PETSC_SUCCESS);
383 : }
384 :
385 : /*@
386 : RGCheckInside - Determines if a set of given points are inside the region or not.
387 :
388 : Not Collective
389 :
390 : Input Parameters:
391 : + rg - the region context
392 : . n - number of points to check
393 : . ar - array of real parts
394 : - ai - array of imaginary parts
395 :
396 : Output Parameter:
397 : . inside - array of results (1=inside, 0=on the contour, -1=outside)
398 :
399 : Note:
400 : The point a is expressed as a couple of PetscScalar variables ar,ai.
401 : If built with complex scalars, the point is supposed to be stored in ar,
402 : otherwise ar,ai contain the real and imaginary parts, respectively.
403 :
404 : If a scaling factor was set, the points are scaled before checking.
405 :
406 : Level: intermediate
407 :
408 : .seealso: RGSetScale(), RGSetComplement()
409 : @*/
410 66194 : PetscErrorCode RGCheckInside(RG rg,PetscInt n,PetscScalar *ar,PetscScalar *ai,PetscInt *inside)
411 : {
412 66194 : PetscReal px,py;
413 66194 : PetscInt i;
414 :
415 66194 : PetscFunctionBegin;
416 66194 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
417 66194 : PetscValidType(rg,1);
418 66194 : PetscAssertPointer(ar,3);
419 : #if !defined(PETSC_USE_COMPLEX)
420 66194 : PetscAssertPointer(ai,4);
421 : #endif
422 66194 : PetscAssertPointer(inside,5);
423 :
424 198300 : for (i=0;i<n;i++) {
425 : #if defined(PETSC_USE_COMPLEX)
426 : px = PetscRealPart(ar[i]);
427 : py = PetscImaginaryPart(ar[i]);
428 : #else
429 132106 : px = ar[i];
430 132106 : py = ai[i];
431 : #endif
432 132106 : if (PetscUnlikely(rg->sfactor != 1.0)) {
433 2 : px /= rg->sfactor;
434 2 : py /= rg->sfactor;
435 : }
436 132106 : PetscUseTypeMethod(rg,checkinside,px,py,inside+i);
437 132106 : if (PetscUnlikely(rg->complement)) inside[i] = -inside[i];
438 : }
439 66194 : PetscFunctionReturn(PETSC_SUCCESS);
440 : }
441 :
442 : /*@
443 : RGIsAxisymmetric - Determines if the region is symmetric with respect
444 : to the real or imaginary axis.
445 :
446 : Not Collective
447 :
448 : Input Parameters:
449 : + rg - the region context
450 : - vertical - true if symmetry must be checked against the vertical axis
451 :
452 : Output Parameter:
453 : . symm - true if the region is axisymmetric
454 :
455 : Note:
456 : If the vertical argument is true, symmetry is checked with respect to
457 : the vertical axis, otherwise with respect to the horizontal axis.
458 :
459 : Level: intermediate
460 :
461 : .seealso: RGCanUseConjugates()
462 : @*/
463 9 : PetscErrorCode RGIsAxisymmetric(RG rg,PetscBool vertical,PetscBool *symm)
464 : {
465 9 : PetscFunctionBegin;
466 9 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
467 9 : PetscValidType(rg,1);
468 9 : PetscAssertPointer(symm,3);
469 9 : *symm = PETSC_FALSE;
470 9 : PetscTryTypeMethod(rg,isaxisymmetric,vertical,symm);
471 9 : PetscFunctionReturn(PETSC_SUCCESS);
472 : }
473 :
474 : /*@
475 : RGCanUseConjugates - Used in contour integral methods to determine whether
476 : half of integration points can be avoided (use their conjugates).
477 :
478 : Not Collective
479 :
480 : Input Parameters:
481 : + rg - the region context
482 : - realmats - true if the problem matrices are real
483 :
484 : Output Parameter:
485 : . useconj - whether it is possible to use conjugates
486 :
487 : Notes:
488 : If some integration points are the conjugates of other points, then the
489 : associated computational cost can be saved. This depends on the problem
490 : matrices being real and also the region being symmetric with respect to
491 : the horizontal axis. The result is false if using real arithmetic or
492 : in the case of a flat region (height equal to zero).
493 :
494 : Level: developer
495 :
496 : .seealso: RGIsAxisymmetric()
497 : @*/
498 26 : PetscErrorCode RGCanUseConjugates(RG rg,PetscBool realmats,PetscBool *useconj)
499 : {
500 : #if defined(PETSC_USE_COMPLEX)
501 : PetscReal c,d;
502 : PetscBool isaxisymm;
503 : #endif
504 :
505 26 : PetscFunctionBegin;
506 26 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
507 26 : PetscValidType(rg,1);
508 26 : PetscAssertPointer(useconj,3);
509 26 : *useconj = PETSC_FALSE;
510 : #if defined(PETSC_USE_COMPLEX)
511 : if (realmats) {
512 : PetscCall(RGIsAxisymmetric(rg,PETSC_FALSE,&isaxisymm));
513 : if (isaxisymm) {
514 : PetscCall(RGComputeBoundingBox(rg,NULL,NULL,&c,&d));
515 : if (c!=d) *useconj = PETSC_TRUE;
516 : }
517 : }
518 : #endif
519 26 : PetscFunctionReturn(PETSC_SUCCESS);
520 : }
521 :
522 : /*@
523 : RGComputeContour - Computes the coordinates of several points lying on the
524 : contour of the region.
525 :
526 : Not Collective
527 :
528 : Input Parameters:
529 : + rg - the region context
530 : - n - number of points to compute
531 :
532 : Output Parameters:
533 : + cr - location to store real parts
534 : - ci - location to store imaginary parts
535 :
536 : Notes:
537 : In real scalars, either cr or ci can be NULL (but not both). In complex
538 : scalars, the coordinates are stored in cr, which cannot be NULL (ci is
539 : not referenced).
540 :
541 : Level: intermediate
542 :
543 : .seealso: RGComputeBoundingBox()
544 : @*/
545 64 : PetscErrorCode RGComputeContour(RG rg,PetscInt n,PetscScalar cr[],PetscScalar ci[])
546 : {
547 64 : PetscInt i;
548 :
549 64 : PetscFunctionBegin;
550 64 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
551 64 : PetscValidType(rg,1);
552 : #if defined(PETSC_USE_COMPLEX)
553 : PetscAssertPointer(cr,3);
554 : #else
555 64 : PetscCheck(cr || ci,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"cr and ci cannot be NULL simultaneously");
556 : #endif
557 64 : PetscCheck(!rg->complement,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Cannot compute contour of region with complement flag set");
558 64 : PetscUseTypeMethod(rg,computecontour,n,cr,ci);
559 152159 : for (i=0;i<n;i++) {
560 152095 : if (cr) cr[i] *= rg->sfactor;
561 152095 : if (ci) ci[i] *= rg->sfactor;
562 : }
563 64 : PetscFunctionReturn(PETSC_SUCCESS);
564 : }
565 :
566 : /*@
567 : RGComputeBoundingBox - Determines the endpoints of a rectangle in the complex plane that
568 : contains the region.
569 :
570 : Not Collective
571 :
572 : Input Parameter:
573 : . rg - the region context
574 :
575 : Output Parameters:
576 : + a - left endpoint of the bounding box in the real axis
577 : . b - right endpoint of the bounding box in the real axis
578 : . c - bottom endpoint of the bounding box in the imaginary axis
579 : - d - top endpoint of the bounding box in the imaginary axis
580 :
581 : Notes:
582 : The bounding box is defined as [a,b]x[c,d]. In regions that are not bounded (e.g. an
583 : open interval) or with the complement flag set, it makes no sense to compute a bounding
584 : box, so the return values are infinite.
585 :
586 : Level: intermediate
587 :
588 : .seealso: RGComputeContour()
589 : @*/
590 13 : PetscErrorCode RGComputeBoundingBox(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
591 : {
592 13 : PetscFunctionBegin;
593 13 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
594 13 : PetscValidType(rg,1);
595 :
596 13 : if (rg->complement) { /* cannot compute bounding box */
597 0 : if (a) *a = -PETSC_MAX_REAL;
598 0 : if (b) *b = PETSC_MAX_REAL;
599 0 : if (c) *c = -PETSC_MAX_REAL;
600 0 : if (d) *d = PETSC_MAX_REAL;
601 : } else {
602 13 : PetscUseTypeMethod(rg,computebbox,a,b,c,d);
603 13 : if (a && *a!=-PETSC_MAX_REAL) *a *= rg->sfactor;
604 13 : if (b && *b!= PETSC_MAX_REAL) *b *= rg->sfactor;
605 13 : if (c && *c!=-PETSC_MAX_REAL) *c *= rg->sfactor;
606 13 : if (d && *d!= PETSC_MAX_REAL) *d *= rg->sfactor;
607 : }
608 13 : PetscFunctionReturn(PETSC_SUCCESS);
609 : }
610 :
611 : /*@
612 : RGComputeQuadrature - Computes the values of the parameters used in a
613 : quadrature rule for a contour integral around the boundary of the region.
614 :
615 : Not Collective
616 :
617 : Input Parameters:
618 : + rg - the region context
619 : . quad - the type of quadrature
620 : - n - number of quadrature points to compute
621 :
622 : Output Parameters:
623 : + z - quadrature points
624 : . zn - normalized quadrature points
625 : - w - quadrature weights
626 :
627 : Notes:
628 : In complex scalars, the values returned in z are often the same as those
629 : computed by RGComputeContour(), but this is not the case in real scalars
630 : where all output arguments are real.
631 :
632 : The computed values change for different quadrature rules (trapezoidal
633 : or Chebyshev).
634 :
635 : Level: intermediate
636 :
637 : .seealso: RGComputeContour()
638 : @*/
639 26 : PetscErrorCode RGComputeQuadrature(RG rg,RGQuadRule quad,PetscInt n,PetscScalar z[],PetscScalar zn[],PetscScalar w[])
640 : {
641 26 : PetscFunctionBegin;
642 26 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
643 26 : PetscValidType(rg,1);
644 26 : PetscAssertPointer(z,4);
645 26 : PetscAssertPointer(zn,5);
646 26 : PetscAssertPointer(w,6);
647 :
648 26 : PetscCall(RGComputeContour(rg,n,z,NULL));
649 26 : PetscUseTypeMethod(rg,computequadrature,quad,n,z,zn,w);
650 26 : PetscFunctionReturn(PETSC_SUCCESS);
651 : }
652 :
653 : /*@
654 : RGSetComplement - Sets a flag to indicate that the region is the complement
655 : of the specified one.
656 :
657 : Logically Collective
658 :
659 : Input Parameters:
660 : + rg - the region context
661 : - flg - the boolean flag
662 :
663 : Options Database Key:
664 : . -rg_complement <bool> - Activate/deactivate the complementation of the region
665 :
666 : Level: intermediate
667 :
668 : .seealso: RGGetComplement()
669 : @*/
670 1 : PetscErrorCode RGSetComplement(RG rg,PetscBool flg)
671 : {
672 1 : PetscFunctionBegin;
673 1 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
674 4 : PetscValidLogicalCollectiveBool(rg,flg,2);
675 1 : rg->complement = flg;
676 1 : PetscFunctionReturn(PETSC_SUCCESS);
677 : }
678 :
679 : /*@
680 : RGGetComplement - Gets a flag that indicates whether the region
681 : is complemented or not.
682 :
683 : Not Collective
684 :
685 : Input Parameter:
686 : . rg - the region context
687 :
688 : Output Parameter:
689 : . flg - the flag
690 :
691 : Level: intermediate
692 :
693 : .seealso: RGSetComplement()
694 : @*/
695 26 : PetscErrorCode RGGetComplement(RG rg,PetscBool *flg)
696 : {
697 26 : PetscFunctionBegin;
698 26 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
699 26 : PetscAssertPointer(flg,2);
700 26 : *flg = rg->complement;
701 26 : PetscFunctionReturn(PETSC_SUCCESS);
702 : }
703 :
704 : /*@
705 : RGSetScale - Sets the scaling factor to be used when checking that a
706 : point is inside the region and when computing the contour.
707 :
708 : Logically Collective
709 :
710 : Input Parameters:
711 : + rg - the region context
712 : - sfactor - the scaling factor
713 :
714 : Options Database Key:
715 : . -rg_scale <real> - Sets the scaling factor
716 :
717 : Level: advanced
718 :
719 : .seealso: RGGetScale(), RGCheckInside()
720 : @*/
721 1 : PetscErrorCode RGSetScale(RG rg,PetscReal sfactor)
722 : {
723 1 : PetscFunctionBegin;
724 1 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
725 4 : PetscValidLogicalCollectiveReal(rg,sfactor,2);
726 1 : if (sfactor == (PetscReal)PETSC_DEFAULT || sfactor == (PetscReal)PETSC_DECIDE) sfactor = 1.0;
727 1 : PetscCheck(sfactor>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of scaling factor. Must be > 0");
728 1 : rg->sfactor = sfactor;
729 1 : PetscFunctionReturn(PETSC_SUCCESS);
730 : }
731 :
732 : /*@
733 : RGGetScale - Gets the scaling factor.
734 :
735 : Not Collective
736 :
737 : Input Parameter:
738 : . rg - the region context
739 :
740 : Output Parameter:
741 : . sfactor - the scaling factor
742 :
743 : Level: advanced
744 :
745 : .seealso: RGSetScale()
746 : @*/
747 52 : PetscErrorCode RGGetScale(RG rg,PetscReal *sfactor)
748 : {
749 52 : PetscFunctionBegin;
750 52 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
751 52 : PetscAssertPointer(sfactor,2);
752 52 : *sfactor = rg->sfactor;
753 52 : PetscFunctionReturn(PETSC_SUCCESS);
754 : }
755 :
756 : /*@
757 : RGPushScale - Sets an additional scaling factor, that will multiply the
758 : user-defined scaling factor.
759 :
760 : Logically Collective
761 :
762 : Input Parameters:
763 : + rg - the region context
764 : - sfactor - the scaling factor
765 :
766 : Notes:
767 : The current implementation does not allow pushing several scaling factors.
768 :
769 : This is intended for internal use, for instance in polynomial eigensolvers
770 : that use parameter scaling.
771 :
772 : Level: developer
773 :
774 : .seealso: RGPopScale(), RGSetScale()
775 : @*/
776 79 : PetscErrorCode RGPushScale(RG rg,PetscReal sfactor)
777 : {
778 79 : PetscFunctionBegin;
779 79 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
780 316 : PetscValidLogicalCollectiveReal(rg,sfactor,2);
781 79 : PetscCheck(sfactor>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of scaling factor. Must be > 0");
782 79 : PetscCheck(!rg->osfactor,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Current implementation does not allow pushing several scaling factors");
783 79 : rg->osfactor = rg->sfactor;
784 79 : rg->sfactor *= sfactor;
785 79 : PetscFunctionReturn(PETSC_SUCCESS);
786 : }
787 :
788 : /*@
789 : RGPopScale - Pops the scaling factor set with RGPushScale().
790 :
791 : Logically Collective
792 :
793 : Input Parameter:
794 : . rg - the region context
795 :
796 : Level: developer
797 :
798 : .seealso: RGPushScale()
799 : @*/
800 79 : PetscErrorCode RGPopScale(RG rg)
801 : {
802 79 : PetscFunctionBegin;
803 79 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
804 79 : PetscCheck(rg->osfactor,PetscObjectComm((PetscObject)rg),PETSC_ERR_ORDER,"Must call RGPushScale first");
805 79 : rg->sfactor = rg->osfactor;
806 79 : rg->osfactor = 0.0;
807 79 : PetscFunctionReturn(PETSC_SUCCESS);
808 : }
809 :
810 : /*@C
811 : RGDestroy - Destroys RG context that was created with RGCreate().
812 :
813 : Collective
814 :
815 : Input Parameter:
816 : . rg - the region context
817 :
818 : Level: beginner
819 :
820 : .seealso: RGCreate()
821 : @*/
822 940 : PetscErrorCode RGDestroy(RG *rg)
823 : {
824 940 : PetscFunctionBegin;
825 940 : if (!*rg) PetscFunctionReturn(PETSC_SUCCESS);
826 916 : PetscValidHeaderSpecific(*rg,RG_CLASSID,1);
827 916 : if (--((PetscObject)*rg)->refct > 0) { *rg = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
828 906 : PetscTryTypeMethod(*rg,destroy);
829 906 : PetscCall(PetscHeaderDestroy(rg));
830 906 : PetscFunctionReturn(PETSC_SUCCESS);
831 : }
832 :
833 : /*@C
834 : RGRegister - Adds a region to the RG package.
835 :
836 : Not Collective
837 :
838 : Input Parameters:
839 : + name - name of a new user-defined RG
840 : - function - routine to create context
841 :
842 : Notes:
843 : RGRegister() may be called multiple times to add several user-defined regions.
844 :
845 : Level: advanced
846 :
847 : .seealso: RGRegisterAll()
848 : @*/
849 3088 : PetscErrorCode RGRegister(const char *name,PetscErrorCode (*function)(RG))
850 : {
851 3088 : PetscFunctionBegin;
852 3088 : PetscCall(RGInitializePackage());
853 3088 : PetscCall(PetscFunctionListAdd(&RGList,name,function));
854 3088 : PetscFunctionReturn(PETSC_SUCCESS);
855 : }
|