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 DS routines
12 : */
13 :
14 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
15 :
16 : PetscFunctionList DSList = NULL;
17 : PetscBool DSRegisterAllCalled = PETSC_FALSE;
18 : PetscClassId DS_CLASSID = 0;
19 : PetscLogEvent DS_Solve = 0,DS_Vectors = 0,DS_Synchronize = 0,DS_Other = 0;
20 : static PetscBool DSPackageInitialized = PETSC_FALSE;
21 :
22 : const char *DSStateTypes[] = {"RAW","INTERMEDIATE","CONDENSED","TRUNCATED","DSStateType","DS_STATE_",NULL};
23 : const char *DSParallelTypes[] = {"REDUNDANT","SYNCHRONIZED","DISTRIBUTED","DSParallelType","DS_PARALLEL_",NULL};
24 : const char *DSMatName[DS_NUM_MAT] = {"A","B","C","T","D","Q","Z","X","Y","U","V","W","E0","E1","E2","E3","E4","E5","E6","E7","E8","E9"};
25 : DSMatType DSMatExtra[DS_NUM_EXTRA] = {DS_MAT_E0,DS_MAT_E1,DS_MAT_E2,DS_MAT_E3,DS_MAT_E4,DS_MAT_E5,DS_MAT_E6,DS_MAT_E7,DS_MAT_E8,DS_MAT_E9};
26 :
27 : /*@C
28 : DSFinalizePackage - This function destroys everything in the SLEPc interface
29 : to the DS package. It is called from SlepcFinalize().
30 :
31 : Level: developer
32 :
33 : .seealso: SlepcFinalize()
34 : @*/
35 976 : PetscErrorCode DSFinalizePackage(void)
36 : {
37 976 : PetscFunctionBegin;
38 976 : PetscCall(PetscFunctionListDestroy(&DSList));
39 976 : DSPackageInitialized = PETSC_FALSE;
40 976 : DSRegisterAllCalled = PETSC_FALSE;
41 976 : PetscFunctionReturn(PETSC_SUCCESS);
42 : }
43 :
44 : /*@C
45 : DSInitializePackage - This function initializes everything in the DS package.
46 : It is called from PetscDLLibraryRegister() when using dynamic libraries, and
47 : on the first call to DSCreate() when using static libraries.
48 :
49 : Level: developer
50 :
51 : .seealso: SlepcInitialize()
52 : @*/
53 11901 : PetscErrorCode DSInitializePackage(void)
54 : {
55 11901 : char logList[256];
56 11901 : PetscBool opt,pkg;
57 11901 : PetscClassId classids[1];
58 :
59 11901 : PetscFunctionBegin;
60 11901 : if (DSPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
61 976 : DSPackageInitialized = PETSC_TRUE;
62 : /* Register Classes */
63 976 : PetscCall(PetscClassIdRegister("Direct Solver",&DS_CLASSID));
64 : /* Register Constructors */
65 976 : PetscCall(DSRegisterAll());
66 : /* Register Events */
67 976 : PetscCall(PetscLogEventRegister("DSSolve",DS_CLASSID,&DS_Solve));
68 976 : PetscCall(PetscLogEventRegister("DSVectors",DS_CLASSID,&DS_Vectors));
69 976 : PetscCall(PetscLogEventRegister("DSSynchronize",DS_CLASSID,&DS_Synchronize));
70 976 : PetscCall(PetscLogEventRegister("DSOther",DS_CLASSID,&DS_Other));
71 : /* Process Info */
72 976 : classids[0] = DS_CLASSID;
73 976 : PetscCall(PetscInfoProcessClass("ds",1,&classids[0]));
74 : /* Process summary exclusions */
75 976 : PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
76 976 : if (opt) {
77 8 : PetscCall(PetscStrInList("ds",logList,',',&pkg));
78 8 : if (pkg) PetscCall(PetscLogEventDeactivateClass(DS_CLASSID));
79 : }
80 : /* Register package finalizer */
81 976 : PetscCall(PetscRegisterFinalize(DSFinalizePackage));
82 976 : PetscFunctionReturn(PETSC_SUCCESS);
83 : }
84 :
85 : /*@
86 : DSCreate - Creates a DS context.
87 :
88 : Collective
89 :
90 : Input Parameter:
91 : . comm - MPI communicator
92 :
93 : Output Parameter:
94 : . newds - location to put the DS context
95 :
96 : Level: beginner
97 :
98 : Note:
99 : DS objects are not intended for normal users but only for
100 : advanced user that for instance implement their own solvers.
101 :
102 : .seealso: DSDestroy(), DS
103 : @*/
104 1164 : PetscErrorCode DSCreate(MPI_Comm comm,DS *newds)
105 : {
106 1164 : DS ds;
107 1164 : PetscInt i;
108 :
109 1164 : PetscFunctionBegin;
110 1164 : PetscAssertPointer(newds,2);
111 1164 : PetscCall(DSInitializePackage());
112 1164 : PetscCall(SlepcHeaderCreate(ds,DS_CLASSID,"DS","Direct Solver (or Dense System)","DS",comm,DSDestroy,DSView));
113 :
114 1164 : ds->state = DS_STATE_RAW;
115 1164 : ds->method = 0;
116 1164 : ds->compact = PETSC_FALSE;
117 1164 : ds->refined = PETSC_FALSE;
118 1164 : ds->extrarow = PETSC_FALSE;
119 1164 : ds->ld = 0;
120 1164 : ds->l = 0;
121 1164 : ds->n = 0;
122 1164 : ds->k = 0;
123 1164 : ds->t = 0;
124 1164 : ds->bs = 1;
125 1164 : ds->sc = NULL;
126 1164 : ds->pmode = DS_PARALLEL_REDUNDANT;
127 :
128 26772 : for (i=0;i<DS_NUM_MAT;i++) ds->omat[i] = NULL;
129 1164 : ds->perm = NULL;
130 1164 : ds->data = NULL;
131 1164 : ds->scset = PETSC_FALSE;
132 1164 : ds->work = NULL;
133 1164 : ds->rwork = NULL;
134 1164 : ds->iwork = NULL;
135 1164 : ds->lwork = 0;
136 1164 : ds->lrwork = 0;
137 1164 : ds->liwork = 0;
138 :
139 1164 : *newds = ds;
140 1164 : PetscFunctionReturn(PETSC_SUCCESS);
141 : }
142 :
143 : /*@
144 : DSSetOptionsPrefix - Sets the prefix used for searching for all
145 : DS options in the database.
146 :
147 : Logically Collective
148 :
149 : Input Parameters:
150 : + ds - the direct solver context
151 : - prefix - the prefix string to prepend to all DS option requests
152 :
153 : Notes:
154 : A hyphen (-) must NOT be given at the beginning of the prefix name.
155 : The first character of all runtime options is AUTOMATICALLY the
156 : hyphen.
157 :
158 : Level: advanced
159 :
160 : .seealso: DSAppendOptionsPrefix()
161 : @*/
162 210 : PetscErrorCode DSSetOptionsPrefix(DS ds,const char *prefix)
163 : {
164 210 : PetscFunctionBegin;
165 210 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
166 210 : PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ds,prefix));
167 210 : PetscFunctionReturn(PETSC_SUCCESS);
168 : }
169 :
170 : /*@
171 : DSAppendOptionsPrefix - Appends to the prefix used for searching for all
172 : DS options in the database.
173 :
174 : Logically Collective
175 :
176 : Input Parameters:
177 : + ds - the direct solver context
178 : - prefix - the prefix string to prepend to all DS option requests
179 :
180 : Notes:
181 : A hyphen (-) must NOT be given at the beginning of the prefix name.
182 : The first character of all runtime options is AUTOMATICALLY the hyphen.
183 :
184 : Level: advanced
185 :
186 : .seealso: DSSetOptionsPrefix()
187 : @*/
188 174 : PetscErrorCode DSAppendOptionsPrefix(DS ds,const char *prefix)
189 : {
190 174 : PetscFunctionBegin;
191 174 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
192 174 : PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ds,prefix));
193 174 : PetscFunctionReturn(PETSC_SUCCESS);
194 : }
195 :
196 : /*@
197 : DSGetOptionsPrefix - Gets the prefix used for searching for all
198 : DS options in the database.
199 :
200 : Not Collective
201 :
202 : Input Parameters:
203 : . ds - the direct solver context
204 :
205 : Output Parameters:
206 : . prefix - pointer to the prefix string used is returned
207 :
208 : Note:
209 : On the Fortran side, the user should pass in a string 'prefix' of
210 : sufficient length to hold the prefix.
211 :
212 : Level: advanced
213 :
214 : .seealso: DSSetOptionsPrefix(), DSAppendOptionsPrefix()
215 : @*/
216 0 : PetscErrorCode DSGetOptionsPrefix(DS ds,const char *prefix[])
217 : {
218 0 : PetscFunctionBegin;
219 0 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
220 0 : PetscAssertPointer(prefix,2);
221 0 : PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ds,prefix));
222 0 : PetscFunctionReturn(PETSC_SUCCESS);
223 : }
224 :
225 : /*@
226 : DSSetType - Selects the type for the DS object.
227 :
228 : Logically Collective
229 :
230 : Input Parameters:
231 : + ds - the direct solver context
232 : - type - a known type
233 :
234 : Level: intermediate
235 :
236 : .seealso: DSGetType()
237 : @*/
238 2332 : PetscErrorCode DSSetType(DS ds,DSType type)
239 : {
240 2332 : PetscErrorCode (*r)(DS);
241 2332 : PetscBool match;
242 :
243 2332 : PetscFunctionBegin;
244 2332 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
245 2332 : PetscAssertPointer(type,2);
246 :
247 2332 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,type,&match));
248 2332 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
249 :
250 1524 : PetscCall(PetscFunctionListFind(DSList,type,&r));
251 1524 : PetscCheck(r,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DS type %s",type);
252 :
253 1524 : PetscTryTypeMethod(ds,destroy);
254 1524 : PetscCall(DSReset(ds));
255 1524 : PetscCall(PetscMemzero(ds->ops,sizeof(struct _DSOps)));
256 :
257 1524 : PetscCall(PetscObjectChangeTypeName((PetscObject)ds,type));
258 1524 : PetscCall((*r)(ds));
259 1524 : PetscFunctionReturn(PETSC_SUCCESS);
260 : }
261 :
262 : /*@
263 : DSGetType - Gets the DS type name (as a string) from the DS context.
264 :
265 : Not Collective
266 :
267 : Input Parameter:
268 : . ds - the direct solver context
269 :
270 : Output Parameter:
271 : . type - name of the direct solver
272 :
273 : Level: intermediate
274 :
275 : .seealso: DSSetType()
276 : @*/
277 2 : PetscErrorCode DSGetType(DS ds,DSType *type)
278 : {
279 2 : PetscFunctionBegin;
280 2 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
281 2 : PetscAssertPointer(type,2);
282 2 : *type = ((PetscObject)ds)->type_name;
283 2 : PetscFunctionReturn(PETSC_SUCCESS);
284 : }
285 :
286 : /*@
287 : DSDuplicate - Creates a new direct solver object with the same options as
288 : an existing one.
289 :
290 : Collective
291 :
292 : Input Parameter:
293 : . ds - direct solver context
294 :
295 : Output Parameter:
296 : . dsnew - location to put the new DS
297 :
298 : Notes:
299 : DSDuplicate() DOES NOT COPY the matrices, and the new DS does not even have
300 : internal arrays allocated. Use DSAllocate() to use the new DS.
301 :
302 : The sorting criterion options are not copied, see DSSetSlepcSC().
303 :
304 : Level: intermediate
305 :
306 : .seealso: DSCreate(), DSAllocate(), DSSetSlepcSC()
307 : @*/
308 0 : PetscErrorCode DSDuplicate(DS ds,DS *dsnew)
309 : {
310 0 : PetscFunctionBegin;
311 0 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
312 0 : PetscAssertPointer(dsnew,2);
313 0 : PetscCall(DSCreate(PetscObjectComm((PetscObject)ds),dsnew));
314 0 : if (((PetscObject)ds)->type_name) PetscCall(DSSetType(*dsnew,((PetscObject)ds)->type_name));
315 0 : (*dsnew)->method = ds->method;
316 0 : (*dsnew)->compact = ds->compact;
317 0 : (*dsnew)->refined = ds->refined;
318 0 : (*dsnew)->extrarow = ds->extrarow;
319 0 : (*dsnew)->bs = ds->bs;
320 0 : (*dsnew)->pmode = ds->pmode;
321 0 : PetscFunctionReturn(PETSC_SUCCESS);
322 : }
323 :
324 : /*@
325 : DSSetMethod - Selects the method to be used to solve the problem.
326 :
327 : Logically Collective
328 :
329 : Input Parameters:
330 : + ds - the direct solver context
331 : - meth - an index identifying the method
332 :
333 : Options Database Key:
334 : . -ds_method <meth> - Sets the method
335 :
336 : Level: intermediate
337 :
338 : .seealso: DSGetMethod()
339 : @*/
340 79 : PetscErrorCode DSSetMethod(DS ds,PetscInt meth)
341 : {
342 79 : PetscFunctionBegin;
343 79 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
344 237 : PetscValidLogicalCollectiveInt(ds,meth,2);
345 79 : PetscCheck(meth>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
346 79 : PetscCheck(meth<=DS_MAX_SOLVE,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
347 79 : ds->method = meth;
348 79 : PetscFunctionReturn(PETSC_SUCCESS);
349 : }
350 :
351 : /*@
352 : DSGetMethod - Gets the method currently used in the DS.
353 :
354 : Not Collective
355 :
356 : Input Parameter:
357 : . ds - the direct solver context
358 :
359 : Output Parameter:
360 : . meth - identifier of the method
361 :
362 : Level: intermediate
363 :
364 : .seealso: DSSetMethod()
365 : @*/
366 21 : PetscErrorCode DSGetMethod(DS ds,PetscInt *meth)
367 : {
368 21 : PetscFunctionBegin;
369 21 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
370 21 : PetscAssertPointer(meth,2);
371 21 : *meth = ds->method;
372 21 : PetscFunctionReturn(PETSC_SUCCESS);
373 : }
374 :
375 : /*@
376 : DSSetParallel - Selects the mode of operation in parallel runs.
377 :
378 : Logically Collective
379 :
380 : Input Parameters:
381 : + ds - the direct solver context
382 : - pmode - the parallel mode
383 :
384 : Options Database Key:
385 : . -ds_parallel <mode> - Sets the parallel mode, 'redundant', 'synchronized'
386 : or 'distributed'
387 :
388 : Notes:
389 : In the 'redundant' parallel mode, all processes will make the computation
390 : redundantly, starting from the same data, and producing the same result.
391 : This result may be slightly different in the different processes if using a
392 : multithreaded BLAS library, which may cause issues in ill-conditioned problems.
393 :
394 : In the 'synchronized' parallel mode, only the first MPI process performs the
395 : computation and then the computed quantities are broadcast to the other
396 : processes in the communicator. This communication is not done automatically,
397 : an explicit call to DSSynchronize() is required.
398 :
399 : The 'distributed' parallel mode can be used in some DS types only, such
400 : as the contour integral method of DSNEP. In this case, every MPI process
401 : will be in charge of part of the computation.
402 :
403 : Level: advanced
404 :
405 : .seealso: DSSynchronize(), DSGetParallel()
406 : @*/
407 84 : PetscErrorCode DSSetParallel(DS ds,DSParallelType pmode)
408 : {
409 84 : PetscFunctionBegin;
410 84 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
411 252 : PetscValidLogicalCollectiveEnum(ds,pmode,2);
412 84 : ds->pmode = pmode;
413 84 : PetscFunctionReturn(PETSC_SUCCESS);
414 : }
415 :
416 : /*@
417 : DSGetParallel - Gets the mode of operation in parallel runs.
418 :
419 : Not Collective
420 :
421 : Input Parameter:
422 : . ds - the direct solver context
423 :
424 : Output Parameter:
425 : . pmode - the parallel mode
426 :
427 : Level: advanced
428 :
429 : .seealso: DSSetParallel()
430 : @*/
431 18 : PetscErrorCode DSGetParallel(DS ds,DSParallelType *pmode)
432 : {
433 18 : PetscFunctionBegin;
434 18 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
435 18 : PetscAssertPointer(pmode,2);
436 18 : *pmode = ds->pmode;
437 18 : PetscFunctionReturn(PETSC_SUCCESS);
438 : }
439 :
440 : /*@
441 : DSSetCompact - Switch to compact storage of matrices.
442 :
443 : Logically Collective
444 :
445 : Input Parameters:
446 : + ds - the direct solver context
447 : - comp - a boolean flag
448 :
449 : Notes:
450 : Compact storage is used in some DS types such as DSHEP when the matrix
451 : is tridiagonal. This flag can be used to indicate whether the user
452 : provides the matrix entries via the compact form (the tridiagonal DS_MAT_T)
453 : or the non-compact one (DS_MAT_A).
454 :
455 : The default is PETSC_FALSE.
456 :
457 : Level: advanced
458 :
459 : .seealso: DSGetCompact()
460 : @*/
461 488 : PetscErrorCode DSSetCompact(DS ds,PetscBool comp)
462 : {
463 488 : PetscFunctionBegin;
464 488 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
465 1464 : PetscValidLogicalCollectiveBool(ds,comp,2);
466 488 : if (ds->compact != comp && ds->ld) PetscTryTypeMethod(ds,setcompact,comp);
467 488 : ds->compact = comp;
468 488 : PetscFunctionReturn(PETSC_SUCCESS);
469 : }
470 :
471 : /*@
472 : DSGetCompact - Gets the compact storage flag.
473 :
474 : Not Collective
475 :
476 : Input Parameter:
477 : . ds - the direct solver context
478 :
479 : Output Parameter:
480 : . comp - the flag
481 :
482 : Level: advanced
483 :
484 : .seealso: DSSetCompact()
485 : @*/
486 10 : PetscErrorCode DSGetCompact(DS ds,PetscBool *comp)
487 : {
488 10 : PetscFunctionBegin;
489 10 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
490 10 : PetscAssertPointer(comp,2);
491 10 : *comp = ds->compact;
492 10 : PetscFunctionReturn(PETSC_SUCCESS);
493 : }
494 :
495 : /*@
496 : DSSetExtraRow - Sets a flag to indicate that the matrix has one extra
497 : row.
498 :
499 : Logically Collective
500 :
501 : Input Parameters:
502 : + ds - the direct solver context
503 : - ext - a boolean flag
504 :
505 : Notes:
506 : In Krylov methods it is useful that the matrix representing the direct solver
507 : has one extra row, i.e., has dimension (n+1) x n. If this flag is activated, all
508 : transformations applied to the right of the matrix also affect this additional
509 : row. In that case, (n+1) must be less or equal than the leading dimension.
510 :
511 : The default is PETSC_FALSE.
512 :
513 : Level: advanced
514 :
515 : .seealso: DSSolve(), DSAllocate(), DSGetExtraRow()
516 : @*/
517 784 : PetscErrorCode DSSetExtraRow(DS ds,PetscBool ext)
518 : {
519 784 : PetscFunctionBegin;
520 784 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
521 2352 : PetscValidLogicalCollectiveBool(ds,ext,2);
522 784 : PetscCheck(!ext || ds->n==0 || ds->n!=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot set extra row after setting n=ld");
523 784 : ds->extrarow = ext;
524 784 : PetscFunctionReturn(PETSC_SUCCESS);
525 : }
526 :
527 : /*@
528 : DSGetExtraRow - Gets the extra row flag.
529 :
530 : Not Collective
531 :
532 : Input Parameter:
533 : . ds - the direct solver context
534 :
535 : Output Parameter:
536 : . ext - the flag
537 :
538 : Level: advanced
539 :
540 : .seealso: DSSetExtraRow()
541 : @*/
542 1833 : PetscErrorCode DSGetExtraRow(DS ds,PetscBool *ext)
543 : {
544 1833 : PetscFunctionBegin;
545 1833 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
546 1833 : PetscAssertPointer(ext,2);
547 1833 : *ext = ds->extrarow;
548 1833 : PetscFunctionReturn(PETSC_SUCCESS);
549 : }
550 :
551 : /*@
552 : DSSetRefined - Sets a flag to indicate that refined vectors must be
553 : computed.
554 :
555 : Logically Collective
556 :
557 : Input Parameters:
558 : + ds - the direct solver context
559 : - ref - a boolean flag
560 :
561 : Notes:
562 : Normally the vectors returned in DS_MAT_X are eigenvectors of the
563 : projected matrix. With this flag activated, DSVectors() will return
564 : the right singular vector of the smallest singular value of matrix
565 : \tilde{A}-theta*I, where \tilde{A} is the extended (n+1)xn matrix
566 : and theta is the Ritz value. This is used in the refined Ritz
567 : approximation.
568 :
569 : The default is PETSC_FALSE.
570 :
571 : Level: advanced
572 :
573 : .seealso: DSVectors(), DSGetRefined()
574 : @*/
575 1 : PetscErrorCode DSSetRefined(DS ds,PetscBool ref)
576 : {
577 1 : PetscFunctionBegin;
578 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
579 3 : PetscValidLogicalCollectiveBool(ds,ref,2);
580 1 : ds->refined = ref;
581 1 : PetscFunctionReturn(PETSC_SUCCESS);
582 : }
583 :
584 : /*@
585 : DSGetRefined - Gets the refined vectors flag.
586 :
587 : Not Collective
588 :
589 : Input Parameter:
590 : . ds - the direct solver context
591 :
592 : Output Parameter:
593 : . ref - the flag
594 :
595 : Level: advanced
596 :
597 : .seealso: DSSetRefined()
598 : @*/
599 3192 : PetscErrorCode DSGetRefined(DS ds,PetscBool *ref)
600 : {
601 3192 : PetscFunctionBegin;
602 3192 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
603 3192 : PetscAssertPointer(ref,2);
604 3192 : *ref = ds->refined;
605 3192 : PetscFunctionReturn(PETSC_SUCCESS);
606 : }
607 :
608 : /*@
609 : DSSetBlockSize - Sets the block size.
610 :
611 : Logically Collective
612 :
613 : Input Parameters:
614 : + ds - the direct solver context
615 : - bs - the block size
616 :
617 : Options Database Key:
618 : . -ds_block_size <bs> - Sets the block size
619 :
620 : Level: intermediate
621 :
622 : .seealso: DSGetBlockSize()
623 : @*/
624 0 : PetscErrorCode DSSetBlockSize(DS ds,PetscInt bs)
625 : {
626 0 : PetscFunctionBegin;
627 0 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
628 0 : PetscValidLogicalCollectiveInt(ds,bs,2);
629 0 : PetscCheck(bs>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The block size must be at least one");
630 0 : ds->bs = bs;
631 0 : PetscFunctionReturn(PETSC_SUCCESS);
632 : }
633 :
634 : /*@
635 : DSGetBlockSize - Gets the block size.
636 :
637 : Not Collective
638 :
639 : Input Parameter:
640 : . ds - the direct solver context
641 :
642 : Output Parameter:
643 : . bs - block size
644 :
645 : Level: intermediate
646 :
647 : .seealso: DSSetBlockSize()
648 : @*/
649 0 : PetscErrorCode DSGetBlockSize(DS ds,PetscInt *bs)
650 : {
651 0 : PetscFunctionBegin;
652 0 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
653 0 : PetscAssertPointer(bs,2);
654 0 : *bs = ds->bs;
655 0 : PetscFunctionReturn(PETSC_SUCCESS);
656 : }
657 :
658 : /*@C
659 : DSSetSlepcSC - Sets the sorting criterion context.
660 :
661 : Logically Collective
662 :
663 : Input Parameters:
664 : + ds - the direct solver context
665 : - sc - a pointer to the sorting criterion context
666 :
667 : Note:
668 : Not available in Fortran.
669 :
670 : Level: developer
671 :
672 : .seealso: DSGetSlepcSC(), DSSort()
673 : @*/
674 1 : PetscErrorCode DSSetSlepcSC(DS ds,SlepcSC sc)
675 : {
676 1 : PetscFunctionBegin;
677 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
678 1 : PetscAssertPointer(sc,2);
679 1 : if (ds->sc && !ds->scset) PetscCall(PetscFree(ds->sc));
680 1 : ds->sc = sc;
681 1 : ds->scset = PETSC_TRUE;
682 1 : PetscFunctionReturn(PETSC_SUCCESS);
683 : }
684 :
685 : /*@C
686 : DSGetSlepcSC - Gets the sorting criterion context.
687 :
688 : Not Collective
689 :
690 : Input Parameter:
691 : . ds - the direct solver context
692 :
693 : Output Parameter:
694 : . sc - a pointer to the sorting criterion context
695 :
696 : Note:
697 : Not available in Fortran.
698 :
699 : Level: developer
700 :
701 : .seealso: DSSetSlepcSC(), DSSort()
702 : @*/
703 1551 : PetscErrorCode DSGetSlepcSC(DS ds,SlepcSC *sc)
704 : {
705 1551 : PetscFunctionBegin;
706 1551 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
707 1551 : PetscAssertPointer(sc,2);
708 1551 : if (!ds->sc) PetscCall(PetscNew(&ds->sc));
709 1551 : *sc = ds->sc;
710 1551 : PetscFunctionReturn(PETSC_SUCCESS);
711 : }
712 :
713 : /*@
714 : DSSetFromOptions - Sets DS options from the options database.
715 :
716 : Collective
717 :
718 : Input Parameters:
719 : . ds - the direct solver context
720 :
721 : Notes:
722 : To see all options, run your program with the -help option.
723 :
724 : Level: beginner
725 :
726 : .seealso: DSSetOptionsPrefix()
727 : @*/
728 1089 : PetscErrorCode DSSetFromOptions(DS ds)
729 : {
730 1089 : PetscInt bs,meth;
731 1089 : PetscBool flag;
732 1089 : DSParallelType pmode;
733 :
734 1089 : PetscFunctionBegin;
735 1089 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
736 1089 : PetscCall(DSRegisterAll());
737 : /* Set default type (we do not allow changing it with -ds_type) */
738 1089 : if (!((PetscObject)ds)->type_name) PetscCall(DSSetType(ds,DSNHEP));
739 3267 : PetscObjectOptionsBegin((PetscObject)ds);
740 :
741 1089 : PetscCall(PetscOptionsInt("-ds_block_size","Block size for the dense system solver","DSSetBlockSize",ds->bs,&bs,&flag));
742 1089 : if (flag) PetscCall(DSSetBlockSize(ds,bs));
743 :
744 1089 : PetscCall(PetscOptionsInt("-ds_method","Method to be used for the dense system","DSSetMethod",ds->method,&meth,&flag));
745 1089 : if (flag) PetscCall(DSSetMethod(ds,meth));
746 :
747 1089 : PetscCall(PetscOptionsEnum("-ds_parallel","Operation mode in parallel runs","DSSetParallel",DSParallelTypes,(PetscEnum)ds->pmode,(PetscEnum*)&pmode,&flag));
748 1089 : if (flag) PetscCall(DSSetParallel(ds,pmode));
749 :
750 1089 : PetscTryTypeMethod(ds,setfromoptions,PetscOptionsObject);
751 1089 : PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ds,PetscOptionsObject));
752 1089 : PetscOptionsEnd();
753 1089 : PetscFunctionReturn(PETSC_SUCCESS);
754 : }
755 :
756 : /*@
757 : DSView - Prints the DS data structure.
758 :
759 : Collective
760 :
761 : Input Parameters:
762 : + ds - the direct solver context
763 : - viewer - optional visualization context
764 :
765 : Note:
766 : The available visualization contexts include
767 : + PETSC_VIEWER_STDOUT_SELF - standard output (default)
768 : - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
769 : output where only the first processor opens
770 : the file. All other processors send their
771 : data to the first processor to print.
772 :
773 : The user can open an alternative visualization context with
774 : PetscViewerASCIIOpen() - output to a specified file.
775 :
776 : Level: beginner
777 :
778 : .seealso: DSViewMat()
779 : @*/
780 122 : PetscErrorCode DSView(DS ds,PetscViewer viewer)
781 : {
782 122 : PetscBool isascii;
783 122 : PetscViewerFormat format;
784 122 : PetscMPIInt size;
785 :
786 122 : PetscFunctionBegin;
787 122 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
788 122 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer));
789 122 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
790 122 : PetscCheckSameComm(ds,1,viewer,2);
791 122 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
792 122 : if (isascii) {
793 122 : PetscCall(PetscViewerGetFormat(viewer,&format));
794 122 : PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ds,viewer));
795 122 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
796 122 : if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",DSParallelTypes[ds->pmode]));
797 122 : if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
798 65 : PetscCall(PetscViewerASCIIPrintf(viewer," current state: %s\n",DSStateTypes[ds->state]));
799 65 : PetscCall(PetscViewerASCIIPrintf(viewer," dimensions: ld=%" PetscInt_FMT ", n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT,ds->ld,ds->n,ds->l,ds->k));
800 65 : if (ds->state==DS_STATE_TRUNCATED) PetscCall(PetscViewerASCIIPrintf(viewer,", t=%" PetscInt_FMT "\n",ds->t));
801 65 : else PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
802 215 : PetscCall(PetscViewerASCIIPrintf(viewer," flags:%s%s%s\n",ds->compact?" compact":"",ds->extrarow?" extrarow":"",ds->refined?" refined":""));
803 : }
804 122 : PetscCall(PetscViewerASCIIPushTab(viewer));
805 122 : PetscTryTypeMethod(ds,view,viewer);
806 122 : PetscCall(PetscViewerASCIIPopTab(viewer));
807 : }
808 122 : PetscFunctionReturn(PETSC_SUCCESS);
809 : }
810 :
811 : /*@
812 : DSViewFromOptions - View from options
813 :
814 : Collective
815 :
816 : Input Parameters:
817 : + ds - the direct solver context
818 : . obj - optional object
819 : - name - command line option
820 :
821 : Level: intermediate
822 :
823 : .seealso: DSView(), DSCreate()
824 : @*/
825 2 : PetscErrorCode DSViewFromOptions(DS ds,PetscObject obj,const char name[])
826 : {
827 2 : PetscFunctionBegin;
828 2 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
829 2 : PetscCall(PetscObjectViewFromOptions((PetscObject)ds,obj,name));
830 2 : PetscFunctionReturn(PETSC_SUCCESS);
831 : }
832 :
833 : /*@
834 : DSAllocate - Allocates memory for internal storage or matrices in DS.
835 :
836 : Logically Collective
837 :
838 : Input Parameters:
839 : + ds - the direct solver context
840 : - ld - leading dimension (maximum allowed dimension for the matrices, including
841 : the extra row if present)
842 :
843 : Note:
844 : If the leading dimension is different from a previously set value, then
845 : all matrices are destroyed with DSReset().
846 :
847 : Level: intermediate
848 :
849 : .seealso: DSGetLeadingDimension(), DSSetDimensions(), DSSetExtraRow(), DSReset()
850 : @*/
851 1371 : PetscErrorCode DSAllocate(DS ds,PetscInt ld)
852 : {
853 1371 : PetscFunctionBegin;
854 1371 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
855 4113 : PetscValidLogicalCollectiveInt(ds,ld,2);
856 1371 : PetscValidType(ds,1);
857 1371 : PetscCheck(ld>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Leading dimension should be at least one");
858 1371 : if (ld!=ds->ld) {
859 1075 : PetscCall(PetscInfo(ds,"Allocating memory with leading dimension=%" PetscInt_FMT "\n",ld));
860 1075 : PetscCall(DSReset(ds));
861 1075 : ds->ld = ld;
862 1075 : PetscUseTypeMethod(ds,allocate,ld);
863 : }
864 1371 : PetscFunctionReturn(PETSC_SUCCESS);
865 : }
866 :
867 : /*@
868 : DSReset - Resets the DS context to the initial state.
869 :
870 : Collective
871 :
872 : Input Parameter:
873 : . ds - the direct solver context
874 :
875 : Note:
876 : All data structures with size depending on the leading dimension
877 : of DSAllocate() are released.
878 :
879 : Level: advanced
880 :
881 : .seealso: DSDestroy(), DSAllocate()
882 : @*/
883 3783 : PetscErrorCode DSReset(DS ds)
884 : {
885 3783 : PetscInt i;
886 :
887 3783 : PetscFunctionBegin;
888 3783 : if (ds) PetscValidHeaderSpecific(ds,DS_CLASSID,1);
889 3783 : if (!ds) PetscFunctionReturn(PETSC_SUCCESS);
890 3783 : ds->state = DS_STATE_RAW;
891 3783 : ds->ld = 0;
892 3783 : ds->l = 0;
893 3783 : ds->n = 0;
894 3783 : ds->k = 0;
895 87009 : for (i=0;i<DS_NUM_MAT;i++) PetscCall(MatDestroy(&ds->omat[i]));
896 3783 : PetscCall(PetscFree(ds->perm));
897 3783 : PetscFunctionReturn(PETSC_SUCCESS);
898 : }
899 :
900 : /*@
901 : DSDestroy - Destroys DS context that was created with DSCreate().
902 :
903 : Collective
904 :
905 : Input Parameter:
906 : . ds - the direct solver context
907 :
908 : Level: beginner
909 :
910 : .seealso: DSCreate()
911 : @*/
912 1225 : PetscErrorCode DSDestroy(DS *ds)
913 : {
914 1225 : PetscFunctionBegin;
915 1225 : if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
916 1164 : PetscValidHeaderSpecific(*ds,DS_CLASSID,1);
917 1164 : if (--((PetscObject)*ds)->refct > 0) { *ds = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
918 1164 : PetscCall(DSReset(*ds));
919 1164 : PetscTryTypeMethod(*ds,destroy);
920 1164 : PetscCall(PetscFree((*ds)->work));
921 1164 : PetscCall(PetscFree((*ds)->rwork));
922 1164 : PetscCall(PetscFree((*ds)->iwork));
923 1164 : if (!(*ds)->scset) PetscCall(PetscFree((*ds)->sc));
924 1164 : PetscCall(PetscHeaderDestroy(ds));
925 1164 : PetscFunctionReturn(PETSC_SUCCESS);
926 : }
927 :
928 : /*@C
929 : DSRegister - Adds a direct solver to the DS package.
930 :
931 : Not Collective
932 :
933 : Input Parameters:
934 : + name - name of a new user-defined DS
935 : - function - routine to create context
936 :
937 : Note:
938 : DSRegister() may be called multiple times to add several user-defined
939 : direct solvers.
940 :
941 : Level: advanced
942 :
943 : .seealso: DSRegisterAll()
944 : @*/
945 10736 : PetscErrorCode DSRegister(const char *name,PetscErrorCode (*function)(DS))
946 : {
947 10736 : PetscFunctionBegin;
948 10736 : PetscCall(DSInitializePackage());
949 10736 : PetscCall(PetscFunctionListAdd(&DSList,name,function));
950 10736 : PetscFunctionReturn(PETSC_SUCCESS);
951 : }
952 :
953 : SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS);
954 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS);
955 : SLEPC_EXTERN PetscErrorCode DSCreate_GHEP(DS);
956 : SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS);
957 : SLEPC_EXTERN PetscErrorCode DSCreate_GNHEP(DS);
958 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEPTS(DS);
959 : SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS);
960 : SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS);
961 : SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS);
962 : SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS);
963 : SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS);
964 :
965 : /*@C
966 : DSRegisterAll - Registers all of the direct solvers in the DS package.
967 :
968 : Not Collective
969 :
970 : Level: advanced
971 :
972 : .seealso: DSRegister()
973 : @*/
974 2065 : PetscErrorCode DSRegisterAll(void)
975 : {
976 2065 : PetscFunctionBegin;
977 2065 : if (DSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
978 976 : DSRegisterAllCalled = PETSC_TRUE;
979 976 : PetscCall(DSRegister(DSHEP,DSCreate_HEP));
980 976 : PetscCall(DSRegister(DSNHEP,DSCreate_NHEP));
981 976 : PetscCall(DSRegister(DSGHEP,DSCreate_GHEP));
982 976 : PetscCall(DSRegister(DSGHIEP,DSCreate_GHIEP));
983 976 : PetscCall(DSRegister(DSGNHEP,DSCreate_GNHEP));
984 976 : PetscCall(DSRegister(DSNHEPTS,DSCreate_NHEPTS));
985 976 : PetscCall(DSRegister(DSSVD,DSCreate_SVD));
986 976 : PetscCall(DSRegister(DSHSVD,DSCreate_HSVD));
987 976 : PetscCall(DSRegister(DSGSVD,DSCreate_GSVD));
988 976 : PetscCall(DSRegister(DSPEP,DSCreate_PEP));
989 976 : PetscCall(DSRegister(DSNEP,DSCreate_NEP));
990 976 : PetscFunctionReturn(PETSC_SUCCESS);
991 : }
|
