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 999 : PetscErrorCode DSFinalizePackage(void)
36 : {
37 999 : PetscFunctionBegin;
38 999 : PetscCall(PetscFunctionListDestroy(&DSList));
39 999 : DSPackageInitialized = PETSC_FALSE;
40 999 : DSRegisterAllCalled = PETSC_FALSE;
41 999 : 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 12181 : PetscErrorCode DSInitializePackage(void)
54 : {
55 12181 : char logList[256];
56 12181 : PetscBool opt,pkg;
57 12181 : PetscClassId classids[1];
58 :
59 12181 : PetscFunctionBegin;
60 12181 : if (DSPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
61 999 : DSPackageInitialized = PETSC_TRUE;
62 : /* Register Classes */
63 999 : PetscCall(PetscClassIdRegister("Direct Solver",&DS_CLASSID));
64 : /* Register Constructors */
65 999 : PetscCall(DSRegisterAll());
66 : /* Register Events */
67 999 : PetscCall(PetscLogEventRegister("DSSolve",DS_CLASSID,&DS_Solve));
68 999 : PetscCall(PetscLogEventRegister("DSVectors",DS_CLASSID,&DS_Vectors));
69 999 : PetscCall(PetscLogEventRegister("DSSynchronize",DS_CLASSID,&DS_Synchronize));
70 999 : PetscCall(PetscLogEventRegister("DSOther",DS_CLASSID,&DS_Other));
71 : /* Process Info */
72 999 : classids[0] = DS_CLASSID;
73 999 : PetscCall(PetscInfoProcessClass("ds",1,&classids[0]));
74 : /* Process summary exclusions */
75 999 : PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
76 999 : if (opt) {
77 8 : PetscCall(PetscStrInList("ds",logList,',',&pkg));
78 8 : if (pkg) PetscCall(PetscLogEventDeactivateClass(DS_CLASSID));
79 : }
80 : /* Register package finalizer */
81 999 : PetscCall(PetscRegisterFinalize(DSFinalizePackage));
82 999 : 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 1191 : PetscErrorCode DSCreate(MPI_Comm comm,DS *newds)
105 : {
106 1191 : DS ds;
107 1191 : PetscInt i;
108 :
109 1191 : PetscFunctionBegin;
110 1191 : PetscAssertPointer(newds,2);
111 1191 : PetscCall(DSInitializePackage());
112 1191 : PetscCall(SlepcHeaderCreate(ds,DS_CLASSID,"DS","Direct Solver (or Dense System)","DS",comm,DSDestroy,DSView));
113 :
114 1191 : ds->state = DS_STATE_RAW;
115 1191 : ds->method = 0;
116 1191 : ds->compact = PETSC_FALSE;
117 1191 : ds->refined = PETSC_FALSE;
118 1191 : ds->extrarow = PETSC_FALSE;
119 1191 : ds->ld = 0;
120 1191 : ds->l = 0;
121 1191 : ds->n = 0;
122 1191 : ds->k = 0;
123 1191 : ds->t = 0;
124 1191 : ds->bs = 1;
125 1191 : ds->sc = NULL;
126 1191 : ds->pmode = DS_PARALLEL_REDUNDANT;
127 :
128 27393 : for (i=0;i<DS_NUM_MAT;i++) ds->omat[i] = NULL;
129 1191 : ds->perm = NULL;
130 1191 : ds->data = NULL;
131 1191 : ds->scset = PETSC_FALSE;
132 1191 : ds->work = NULL;
133 1191 : ds->rwork = NULL;
134 1191 : ds->iwork = NULL;
135 1191 : ds->lwork = 0;
136 1191 : ds->lrwork = 0;
137 1191 : ds->liwork = 0;
138 :
139 1191 : *newds = ds;
140 1191 : 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 214 : PetscErrorCode DSSetOptionsPrefix(DS ds,const char *prefix)
163 : {
164 214 : PetscFunctionBegin;
165 214 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
166 214 : PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ds,prefix));
167 214 : 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 178 : PetscErrorCode DSAppendOptionsPrefix(DS ds,const char *prefix)
189 : {
190 178 : PetscFunctionBegin;
191 178 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
192 178 : PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ds,prefix));
193 178 : 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 2382 : PetscErrorCode DSSetType(DS ds,DSType type)
239 : {
240 2382 : PetscErrorCode (*r)(DS);
241 2382 : PetscBool match;
242 :
243 2382 : PetscFunctionBegin;
244 2382 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
245 2382 : PetscAssertPointer(type,2);
246 :
247 2382 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,type,&match));
248 2382 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
249 :
250 1567 : PetscCall(PetscFunctionListFind(DSList,type,&r));
251 1567 : PetscCheck(r,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DS type %s",type);
252 :
253 1567 : PetscTryTypeMethod(ds,destroy);
254 1567 : PetscCall(DSReset(ds));
255 1567 : PetscCall(PetscMemzero(ds->ops,sizeof(struct _DSOps)));
256 :
257 1567 : PetscCall(PetscObjectChangeTypeName((PetscObject)ds,type));
258 1567 : PetscCall((*r)(ds));
259 1567 : 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 510 : PetscErrorCode DSSetCompact(DS ds,PetscBool comp)
462 : {
463 510 : PetscFunctionBegin;
464 510 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
465 1530 : PetscValidLogicalCollectiveBool(ds,comp,2);
466 510 : if (ds->compact != comp && ds->ld) PetscTryTypeMethod(ds,setcompact,comp);
467 510 : ds->compact = comp;
468 510 : 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 807 : PetscErrorCode DSSetExtraRow(DS ds,PetscBool ext)
518 : {
519 807 : PetscFunctionBegin;
520 807 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
521 2421 : PetscValidLogicalCollectiveBool(ds,ext,2);
522 807 : PetscCheck(!ext || ds->n==0 || ds->n!=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot set extra row after setting n=ld");
523 807 : ds->extrarow = ext;
524 807 : 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 2117 : PetscErrorCode DSGetExtraRow(DS ds,PetscBool *ext)
543 : {
544 2117 : PetscFunctionBegin;
545 2117 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
546 2117 : PetscAssertPointer(ext,2);
547 2117 : *ext = ds->extrarow;
548 2117 : 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 3921 : PetscErrorCode DSGetRefined(DS ds,PetscBool *ref)
600 : {
601 3921 : PetscFunctionBegin;
602 3921 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
603 3921 : PetscAssertPointer(ref,2);
604 3921 : *ref = ds->refined;
605 3921 : 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 1578 : PetscErrorCode DSGetSlepcSC(DS ds,SlepcSC *sc)
704 : {
705 1578 : PetscFunctionBegin;
706 1578 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
707 1578 : PetscAssertPointer(sc,2);
708 1578 : if (!ds->sc) PetscCall(PetscNew(&ds->sc));
709 1578 : *sc = ds->sc;
710 1578 : 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 1116 : PetscErrorCode DSSetFromOptions(DS ds)
729 : {
730 1116 : PetscInt bs,meth;
731 1116 : PetscBool flag;
732 1116 : DSParallelType pmode;
733 :
734 1116 : PetscFunctionBegin;
735 1116 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
736 1116 : PetscCall(DSRegisterAll());
737 : /* Set default type (we do not allow changing it with -ds_type) */
738 1116 : if (!((PetscObject)ds)->type_name) PetscCall(DSSetType(ds,DSNHEP));
739 3348 : PetscObjectOptionsBegin((PetscObject)ds);
740 :
741 1116 : PetscCall(PetscOptionsInt("-ds_block_size","Block size for the dense system solver","DSSetBlockSize",ds->bs,&bs,&flag));
742 1116 : if (flag) PetscCall(DSSetBlockSize(ds,bs));
743 :
744 1116 : PetscCall(PetscOptionsInt("-ds_method","Method to be used for the dense system","DSSetMethod",ds->method,&meth,&flag));
745 1116 : if (flag) PetscCall(DSSetMethod(ds,meth));
746 :
747 1116 : PetscCall(PetscOptionsEnum("-ds_parallel","Operation mode in parallel runs","DSSetParallel",DSParallelTypes,(PetscEnum)ds->pmode,(PetscEnum*)&pmode,&flag));
748 1116 : if (flag) PetscCall(DSSetParallel(ds,pmode));
749 :
750 1116 : PetscTryTypeMethod(ds,setfromoptions,PetscOptionsObject);
751 1116 : PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ds,PetscOptionsObject));
752 1116 : PetscOptionsEnd();
753 1116 : 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(), DSReallocate()
850 : @*/
851 1394 : PetscErrorCode DSAllocate(DS ds,PetscInt ld)
852 : {
853 1394 : PetscFunctionBegin;
854 1394 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
855 4182 : PetscValidLogicalCollectiveInt(ds,ld,2);
856 1394 : PetscValidType(ds,1);
857 1394 : PetscCheck(ld>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Leading dimension should be at least one");
858 1394 : if (ld!=ds->ld) {
859 1098 : PetscCall(PetscInfo(ds,"Allocating memory with leading dimension=%" PetscInt_FMT "\n",ld));
860 1098 : PetscCall(DSReset(ds));
861 1098 : ds->ld = ld;
862 1098 : PetscUseTypeMethod(ds,allocate,ld);
863 : }
864 1394 : PetscFunctionReturn(PETSC_SUCCESS);
865 : }
866 :
867 : /*@
868 : DSReallocate - Reallocates memory for internal storage or matrices in DS,
869 : keeping the previously set data.
870 :
871 : Logically Collective
872 :
873 : Input Parameters:
874 : + ds - the direct solver context
875 : - ld - new leading dimension
876 :
877 : Notes:
878 : The new leading dimension must be larger than the previous one. The relevant
879 : data previously set is copied over to the new data structures.
880 :
881 : This operation is not available in all DS types.
882 :
883 : Level: developer
884 :
885 : .seealso: DSAllocate()
886 : @*/
887 12 : PetscErrorCode DSReallocate(DS ds,PetscInt ld)
888 : {
889 12 : PetscFunctionBegin;
890 12 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
891 36 : PetscValidLogicalCollectiveInt(ds,ld,2);
892 12 : PetscValidType(ds,1);
893 12 : PetscCheck(ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"DSAllocate() must be called first");
894 12 : PetscCheck(ld>ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"New leading dimension %" PetscInt_FMT " must be larger than the previous one %" PetscInt_FMT,ld,ds->ld);
895 12 : PetscCall(PetscInfo(ds,"Reallocating memory with new leading dimension=%" PetscInt_FMT "\n",ld));
896 12 : PetscUseTypeMethod(ds,reallocate,ld);
897 12 : ds->ld = ld;
898 12 : PetscFunctionReturn(PETSC_SUCCESS);
899 : }
900 :
901 : /*@
902 : DSReset - Resets the DS context to the initial state.
903 :
904 : Collective
905 :
906 : Input Parameter:
907 : . ds - the direct solver context
908 :
909 : Note:
910 : All data structures with size depending on the leading dimension
911 : of DSAllocate() are released.
912 :
913 : Level: advanced
914 :
915 : .seealso: DSDestroy(), DSAllocate()
916 : @*/
917 3876 : PetscErrorCode DSReset(DS ds)
918 : {
919 3876 : PetscInt i;
920 :
921 3876 : PetscFunctionBegin;
922 3876 : if (ds) PetscValidHeaderSpecific(ds,DS_CLASSID,1);
923 3876 : if (!ds) PetscFunctionReturn(PETSC_SUCCESS);
924 3876 : ds->state = DS_STATE_RAW;
925 3876 : ds->ld = 0;
926 3876 : ds->l = 0;
927 3876 : ds->n = 0;
928 3876 : ds->k = 0;
929 89148 : for (i=0;i<DS_NUM_MAT;i++) PetscCall(MatDestroy(&ds->omat[i]));
930 3876 : PetscCall(PetscFree(ds->perm));
931 3876 : PetscFunctionReturn(PETSC_SUCCESS);
932 : }
933 :
934 : /*@
935 : DSDestroy - Destroys DS context that was created with DSCreate().
936 :
937 : Collective
938 :
939 : Input Parameter:
940 : . ds - the direct solver context
941 :
942 : Level: beginner
943 :
944 : .seealso: DSCreate()
945 : @*/
946 1252 : PetscErrorCode DSDestroy(DS *ds)
947 : {
948 1252 : PetscFunctionBegin;
949 1252 : if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
950 1191 : PetscValidHeaderSpecific(*ds,DS_CLASSID,1);
951 1191 : if (--((PetscObject)*ds)->refct > 0) { *ds = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
952 1191 : PetscCall(DSReset(*ds));
953 1191 : PetscTryTypeMethod(*ds,destroy);
954 1191 : PetscCall(PetscFree((*ds)->work));
955 1191 : PetscCall(PetscFree((*ds)->rwork));
956 1191 : PetscCall(PetscFree((*ds)->iwork));
957 1191 : if (!(*ds)->scset) PetscCall(PetscFree((*ds)->sc));
958 1191 : PetscCall(PetscHeaderDestroy(ds));
959 1191 : PetscFunctionReturn(PETSC_SUCCESS);
960 : }
961 :
962 : /*@C
963 : DSRegister - Adds a direct solver to the DS package.
964 :
965 : Not Collective
966 :
967 : Input Parameters:
968 : + name - name of a new user-defined DS
969 : - function - routine to create context
970 :
971 : Note:
972 : DSRegister() may be called multiple times to add several user-defined
973 : direct solvers.
974 :
975 : Level: advanced
976 :
977 : .seealso: DSRegisterAll()
978 : @*/
979 10989 : PetscErrorCode DSRegister(const char *name,PetscErrorCode (*function)(DS))
980 : {
981 10989 : PetscFunctionBegin;
982 10989 : PetscCall(DSInitializePackage());
983 10989 : PetscCall(PetscFunctionListAdd(&DSList,name,function));
984 10989 : PetscFunctionReturn(PETSC_SUCCESS);
985 : }
986 :
987 : SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS);
988 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS);
989 : SLEPC_EXTERN PetscErrorCode DSCreate_GHEP(DS);
990 : SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS);
991 : SLEPC_EXTERN PetscErrorCode DSCreate_GNHEP(DS);
992 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEPTS(DS);
993 : SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS);
994 : SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS);
995 : SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS);
996 : SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS);
997 : SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS);
998 :
999 : /*@C
1000 : DSRegisterAll - Registers all of the direct solvers in the DS package.
1001 :
1002 : Not Collective
1003 :
1004 : Level: advanced
1005 :
1006 : .seealso: DSRegister()
1007 : @*/
1008 2115 : PetscErrorCode DSRegisterAll(void)
1009 : {
1010 2115 : PetscFunctionBegin;
1011 2115 : if (DSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1012 999 : DSRegisterAllCalled = PETSC_TRUE;
1013 999 : PetscCall(DSRegister(DSHEP,DSCreate_HEP));
1014 999 : PetscCall(DSRegister(DSNHEP,DSCreate_NHEP));
1015 999 : PetscCall(DSRegister(DSGHEP,DSCreate_GHEP));
1016 999 : PetscCall(DSRegister(DSGHIEP,DSCreate_GHIEP));
1017 999 : PetscCall(DSRegister(DSGNHEP,DSCreate_GNHEP));
1018 999 : PetscCall(DSRegister(DSNHEPTS,DSCreate_NHEPTS));
1019 999 : PetscCall(DSRegister(DSSVD,DSCreate_SVD));
1020 999 : PetscCall(DSRegister(DSHSVD,DSCreate_HSVD));
1021 999 : PetscCall(DSRegister(DSGSVD,DSCreate_GSVD));
1022 999 : PetscCall(DSRegister(DSPEP,DSCreate_PEP));
1023 999 : PetscCall(DSRegister(DSNEP,DSCreate_NEP));
1024 999 : PetscFunctionReturn(PETSC_SUCCESS);
1025 : }
|
