Actual source code: dsbasic.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
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: */
14: #include <slepc/private/dsimpl.h>
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;
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};
27: /*@C
28: DSFinalizePackage - This function destroys everything in the SLEPc interface
29: to the `DS` package. It is called from `SlepcFinalize()`.
31: Level: developer
33: .seealso: [](sec:ds), `SlepcFinalize()`, `DSInitializePackage()`
34: @*/
35: PetscErrorCode DSFinalizePackage(void)
36: {
37: PetscFunctionBegin;
38: PetscCall(PetscFunctionListDestroy(&DSList));
39: DSPackageInitialized = PETSC_FALSE;
40: DSRegisterAllCalled = PETSC_FALSE;
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*@C
45: DSInitializePackage - This function initializes everything in the `DS` package.
46: It is called from `PetscDLLibraryRegister_slepc()` when using dynamic libraries, and
47: on the first call to `DSCreate()` when using shared or static libraries.
49: Note:
50: This function never needs to be called by SLEPc users.
52: Level: developer
54: .seealso: [](sec:ds), `DS`, `SlepcInitialize()`, `DSFinalizePackage()`
55: @*/
56: PetscErrorCode DSInitializePackage(void)
57: {
58: char logList[256];
59: PetscBool opt,pkg;
60: PetscClassId classids[1];
62: PetscFunctionBegin;
63: if (DSPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
64: DSPackageInitialized = PETSC_TRUE;
65: /* Register Classes */
66: PetscCall(PetscClassIdRegister("Direct Solver",&DS_CLASSID));
67: /* Register Constructors */
68: PetscCall(DSRegisterAll());
69: /* Register Events */
70: PetscCall(PetscLogEventRegister("DSSolve",DS_CLASSID,&DS_Solve));
71: PetscCall(PetscLogEventRegister("DSVectors",DS_CLASSID,&DS_Vectors));
72: PetscCall(PetscLogEventRegister("DSSynchronize",DS_CLASSID,&DS_Synchronize));
73: PetscCall(PetscLogEventRegister("DSOther",DS_CLASSID,&DS_Other));
74: /* Process Info */
75: classids[0] = DS_CLASSID;
76: PetscCall(PetscInfoProcessClass("ds",1,&classids[0]));
77: /* Process summary exclusions */
78: PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
79: if (opt) {
80: PetscCall(PetscStrInList("ds",logList,',',&pkg));
81: if (pkg) PetscCall(PetscLogEventDeactivateClass(DS_CLASSID));
82: }
83: /* Register package finalizer */
84: PetscCall(PetscRegisterFinalize(DSFinalizePackage));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: /*@
89: DSCreate - Creates a `DS` context.
91: Collective
93: Input Parameter:
94: . comm - MPI communicator
96: Output Parameter:
97: . newds - location to put the `DS` context
99: Level: beginner
101: Note:
102: `DS` objects are not intended for normal users but only for
103: advanced user that for instance implement their own solvers.
105: .seealso: [](sec:ds), `DSDestroy()`, `DS`
106: @*/
107: PetscErrorCode DSCreate(MPI_Comm comm,DS *newds)
108: {
109: DS ds;
110: PetscInt i;
112: PetscFunctionBegin;
113: PetscAssertPointer(newds,2);
114: PetscCall(DSInitializePackage());
115: PetscCall(SlepcHeaderCreate(ds,DS_CLASSID,"DS","Direct Solver (or Dense System)","DS",comm,DSDestroy,DSView));
117: ds->state = DS_STATE_RAW;
118: ds->method = 0;
119: ds->compact = PETSC_FALSE;
120: ds->refined = PETSC_FALSE;
121: ds->extrarow = PETSC_FALSE;
122: ds->ld = 0;
123: ds->l = 0;
124: ds->n = 0;
125: ds->k = 0;
126: ds->t = 0;
127: ds->bs = 1;
128: ds->sc = NULL;
129: ds->pmode = DS_PARALLEL_REDUNDANT;
131: for (i=0;i<DS_NUM_MAT;i++) ds->omat[i] = NULL;
132: ds->perm = NULL;
133: ds->data = NULL;
134: ds->scset = PETSC_FALSE;
135: ds->work = NULL;
136: ds->rwork = NULL;
137: ds->iwork = NULL;
138: ds->lwork = 0;
139: ds->lrwork = 0;
140: ds->liwork = 0;
142: *newds = ds;
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*@
147: DSSetOptionsPrefix - Sets the prefix used for searching for all
148: `DS` options in the database.
150: Logically Collective
152: Input Parameters:
153: + ds - the direct solver context
154: - prefix - the prefix string to prepend to all `DS` option requests
156: Notes:
157: A hyphen (-) must NOT be given at the beginning of the prefix name.
158: The first character of all runtime options is AUTOMATICALLY the
159: hyphen.
161: Level: advanced
163: .seealso: [](sec:ds), `DSAppendOptionsPrefix()`
164: @*/
165: PetscErrorCode DSSetOptionsPrefix(DS ds,const char prefix[])
166: {
167: PetscFunctionBegin;
169: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ds,prefix));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@
174: DSAppendOptionsPrefix - Appends to the prefix used for searching for all
175: `DS` options in the database.
177: Logically Collective
179: Input Parameters:
180: + ds - the direct solver context
181: - prefix - the prefix string to prepend to all DS option requests
183: Notes:
184: A hyphen (-) must NOT be given at the beginning of the prefix name.
185: The first character of all runtime options is AUTOMATICALLY the hyphen.
187: Level: advanced
189: .seealso: [](sec:ds), `DSSetOptionsPrefix()`
190: @*/
191: PetscErrorCode DSAppendOptionsPrefix(DS ds,const char prefix[])
192: {
193: PetscFunctionBegin;
195: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ds,prefix));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: DSGetOptionsPrefix - Gets the prefix used for searching for all
201: `DS` options in the database.
203: Not Collective
205: Input Parameter:
206: . ds - the direct solver context
208: Output Parameter:
209: . prefix - pointer to the prefix string used is returned
211: Level: advanced
213: .seealso: [](sec:ds), `DSSetOptionsPrefix()`, `DSAppendOptionsPrefix()`
214: @*/
215: PetscErrorCode DSGetOptionsPrefix(DS ds,const char *prefix[])
216: {
217: PetscFunctionBegin;
219: PetscAssertPointer(prefix,2);
220: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ds,prefix));
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*@
225: DSSetType - Selects the type for the DS object.
227: Logically Collective
229: Input Parameters:
230: + ds - the direct solver context
231: - type - a known type
233: Level: intermediate
235: .seealso: [](sec:ds), `DSGetType()`
236: @*/
237: PetscErrorCode DSSetType(DS ds,DSType type)
238: {
239: PetscErrorCode (*r)(DS);
240: PetscBool match;
242: PetscFunctionBegin;
244: PetscAssertPointer(type,2);
246: PetscCall(PetscObjectTypeCompare((PetscObject)ds,type,&match));
247: if (match) PetscFunctionReturn(PETSC_SUCCESS);
249: PetscCall(PetscFunctionListFind(DSList,type,&r));
250: PetscCheck(r,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DS type %s",type);
252: PetscTryTypeMethod(ds,destroy);
253: PetscCall(DSReset(ds));
254: PetscCall(PetscMemzero(ds->ops,sizeof(struct _DSOps)));
256: PetscCall(PetscObjectChangeTypeName((PetscObject)ds,type));
257: PetscCall((*r)(ds));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*@
262: DSGetType - Gets the `DS` type name (as a string) from the `DS` context.
264: Not Collective
266: Input Parameter:
267: . ds - the direct solver context
269: Output Parameter:
270: . type - name of the direct solver
272: Note:
273: `type` should not be retained for later use as it will be an invalid pointer
274: if the `DSType` of `ds` is changed.
276: Level: intermediate
278: .seealso: [](sec:ds), `DSSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
279: @*/
280: PetscErrorCode DSGetType(DS ds,DSType *type)
281: {
282: PetscFunctionBegin;
284: PetscAssertPointer(type,2);
285: *type = ((PetscObject)ds)->type_name;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: DSDuplicate - Creates a new direct solver object with the same options as
291: an existing one.
293: Collective
295: Input Parameter:
296: . ds - direct solver context
298: Output Parameter:
299: . dsnew - location to put the new `DS`
301: Notes:
302: `DSDuplicate()` DOES NOT COPY the matrices, and the new `DS` does not even have
303: internal arrays allocated. Use `DSAllocate()` to use the new `DS`.
305: The sorting criterion options are not copied, see `DSSetSlepcSC()`.
307: Level: intermediate
309: .seealso: [](sec:ds), `DSCreate()`, `DSAllocate()`, `DSSetSlepcSC()`
310: @*/
311: PetscErrorCode DSDuplicate(DS ds,DS *dsnew)
312: {
313: PetscFunctionBegin;
315: PetscAssertPointer(dsnew,2);
316: PetscCall(DSCreate(PetscObjectComm((PetscObject)ds),dsnew));
317: if (((PetscObject)ds)->type_name) PetscCall(DSSetType(*dsnew,((PetscObject)ds)->type_name));
318: (*dsnew)->method = ds->method;
319: (*dsnew)->compact = ds->compact;
320: (*dsnew)->refined = ds->refined;
321: (*dsnew)->extrarow = ds->extrarow;
322: (*dsnew)->bs = ds->bs;
323: (*dsnew)->pmode = ds->pmode;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: DSSetMethod - Selects the method to be used to solve the problem.
330: Logically Collective
332: Input Parameters:
333: + ds - the direct solver context
334: - meth - an index identifying the method
336: Options Database Key:
337: . -ds_method meth - sets the method
339: Level: intermediate
341: .seealso: [](sec:ds), `DSGetMethod()`
342: @*/
343: PetscErrorCode DSSetMethod(DS ds,PetscInt meth)
344: {
345: PetscFunctionBegin;
348: PetscCheck(meth>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
349: PetscCheck(meth<=DS_MAX_SOLVE,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
350: ds->method = meth;
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /*@
355: DSGetMethod - Gets the method currently used in the `DS`.
357: Not Collective
359: Input Parameter:
360: . ds - the direct solver context
362: Output Parameter:
363: . meth - identifier of the method
365: Level: intermediate
367: .seealso: [](sec:ds), `DSSetMethod()`
368: @*/
369: PetscErrorCode DSGetMethod(DS ds,PetscInt *meth)
370: {
371: PetscFunctionBegin;
373: PetscAssertPointer(meth,2);
374: *meth = ds->method;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: DSSetParallel - Selects the mode of operation in parallel runs.
381: Logically Collective
383: Input Parameters:
384: + ds - the direct solver context
385: - pmode - the parallel mode
387: Options Database Key:
388: . -ds_parallel (redundant|synchronized|distributed) - sets the parallel mode
390: Notes:
391: In the `redundant` parallel mode, all processes will make the computation
392: redundantly, starting from the same data, and producing the same result.
393: This result may be slightly different in the different processes if using a
394: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
396: In the `synchronized` parallel mode, only the first MPI process performs the
397: computation and then the computed quantities are broadcast to the other
398: processes in the communicator. This communication is not done automatically,
399: an explicit call to `DSSynchronize()` is required.
401: The `distributed` parallel mode can be used in some `DS` types only, such
402: as the contour integral method of `DSNEP`. In this case, every MPI process
403: will be in charge of part of the computation.
405: Level: advanced
407: .seealso: [](sec:ds), `DSSynchronize()`, `DSGetParallel()`
408: @*/
409: PetscErrorCode DSSetParallel(DS ds,DSParallelType pmode)
410: {
411: PetscFunctionBegin;
414: ds->pmode = pmode;
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /*@
419: DSGetParallel - Gets the mode of operation in parallel runs.
421: Not Collective
423: Input Parameter:
424: . ds - the direct solver context
426: Output Parameter:
427: . pmode - the parallel mode
429: Level: advanced
431: .seealso: [](sec:ds), `DSSetParallel()`
432: @*/
433: PetscErrorCode DSGetParallel(DS ds,DSParallelType *pmode)
434: {
435: PetscFunctionBegin;
437: PetscAssertPointer(pmode,2);
438: *pmode = ds->pmode;
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: DSSetCompact - Switch to compact storage of matrices.
445: Logically Collective
447: Input Parameters:
448: + ds - the direct solver context
449: - comp - a boolean flag
451: Notes:
452: Compact storage is used in some `DS` types such as `DSHEP` when the matrix
453: is tridiagonal. This flag can be used to indicate whether the user
454: provides the matrix entries via the compact form (the tridiagonal `DS_MAT_T`)
455: or the non-compact one (`DS_MAT_A`).
457: The default is `PETSC_FALSE`.
459: Level: advanced
461: .seealso: [](sec:ds), `DSGetCompact()`
462: @*/
463: PetscErrorCode DSSetCompact(DS ds,PetscBool comp)
464: {
465: PetscFunctionBegin;
468: if (ds->compact != comp && ds->ld) PetscTryTypeMethod(ds,setcompact,comp);
469: ds->compact = comp;
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: /*@
474: DSGetCompact - Gets the compact storage flag.
476: Not Collective
478: Input Parameter:
479: . ds - the direct solver context
481: Output Parameter:
482: . comp - the flag
484: Level: advanced
486: .seealso: [](sec:ds), `DSSetCompact()`
487: @*/
488: PetscErrorCode DSGetCompact(DS ds,PetscBool *comp)
489: {
490: PetscFunctionBegin;
492: PetscAssertPointer(comp,2);
493: *comp = ds->compact;
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*@
498: DSSetExtraRow - Sets a flag to indicate that the matrix has one extra
499: row.
501: Logically Collective
503: Input Parameters:
504: + ds - the direct solver context
505: - ext - a boolean flag
507: Notes:
508: In Krylov methods it is useful that the matrix representing the direct solver
509: has one extra row, i.e., has dimension $(n+1) \times n$. If this flag is activated, all
510: transformations applied to the right of the matrix also affect this additional
511: row. In that case, $(n+1)$ must be less or equal than the leading dimension.
513: The default is `PETSC_FALSE`.
515: Level: advanced
517: .seealso: [](sec:ds), `DSSolve()`, `DSAllocate()`, `DSGetExtraRow()`
518: @*/
519: PetscErrorCode DSSetExtraRow(DS ds,PetscBool ext)
520: {
521: PetscFunctionBegin;
524: PetscCheck(!ext || ds->n==0 || ds->n!=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot set extra row after setting n=ld");
525: ds->extrarow = ext;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: DSGetExtraRow - Gets the extra row flag.
532: Not Collective
534: Input Parameter:
535: . ds - the direct solver context
537: Output Parameter:
538: . ext - the flag
540: Level: advanced
542: .seealso: [](sec:ds), `DSSetExtraRow()`
543: @*/
544: PetscErrorCode DSGetExtraRow(DS ds,PetscBool *ext)
545: {
546: PetscFunctionBegin;
548: PetscAssertPointer(ext,2);
549: *ext = ds->extrarow;
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: /*@
554: DSSetRefined - Sets a flag to indicate that refined vectors must be
555: computed.
557: Logically Collective
559: Input Parameters:
560: + ds - the direct solver context
561: - ref - a boolean flag
563: Notes:
564: Normally the vectors returned in `DS_MAT_X` are eigenvectors of the
565: projected matrix. With this flag activated, `DSVectors()` will return
566: the right singular vector of the smallest singular value of matrix
567: $\tilde{A}-\theta I$, where $\tilde{A}$ is the extended $(n+1)\times n$
568: matrix and $\theta$ is the Ritz value. This is used in the refined Ritz
569: approximation {cite:p}`Jia97`.
571: The default is `PETSC_FALSE`.
573: Level: advanced
575: .seealso: [](sec:ds), `DSVectors()`, `DSGetRefined()`
576: @*/
577: PetscErrorCode DSSetRefined(DS ds,PetscBool ref)
578: {
579: PetscFunctionBegin;
582: ds->refined = ref;
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: DSGetRefined - Gets the refined vectors flag.
589: Not Collective
591: Input Parameter:
592: . ds - the direct solver context
594: Output Parameter:
595: . ref - the flag
597: Level: advanced
599: .seealso: [](sec:ds), `DSSetRefined()`
600: @*/
601: PetscErrorCode DSGetRefined(DS ds,PetscBool *ref)
602: {
603: PetscFunctionBegin;
605: PetscAssertPointer(ref,2);
606: *ref = ds->refined;
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /*@
611: DSSetBlockSize - Sets the block size.
613: Logically Collective
615: Input Parameters:
616: + ds - the direct solver context
617: - bs - the block size
619: Options Database Key:
620: . -ds_block_size bs - sets the block size
622: Level: intermediate
624: .seealso: [](sec:ds), `DSGetBlockSize()`
625: @*/
626: PetscErrorCode DSSetBlockSize(DS ds,PetscInt bs)
627: {
628: PetscFunctionBegin;
631: PetscCheck(bs>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The block size must be at least one");
632: ds->bs = bs;
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@
637: DSGetBlockSize - Gets the block size.
639: Not Collective
641: Input Parameter:
642: . ds - the direct solver context
644: Output Parameter:
645: . bs - block size
647: Level: intermediate
649: .seealso: [](sec:ds), `DSSetBlockSize()`
650: @*/
651: PetscErrorCode DSGetBlockSize(DS ds,PetscInt *bs)
652: {
653: PetscFunctionBegin;
655: PetscAssertPointer(bs,2);
656: *bs = ds->bs;
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /*@C
661: DSSetSlepcSC - Sets the sorting criterion context.
663: Logically Collective
665: Input Parameters:
666: + ds - the direct solver context
667: - sc - a pointer to the sorting criterion context
669: Note:
670: Not available in Fortran.
672: Level: developer
674: .seealso: [](sec:ds), `DSGetSlepcSC()`, `DSSort()`
675: @*/
676: PetscErrorCode DSSetSlepcSC(DS ds,SlepcSC sc)
677: {
678: PetscFunctionBegin;
680: PetscAssertPointer(sc,2);
681: if (ds->sc && !ds->scset) PetscCall(PetscFree(ds->sc));
682: ds->sc = sc;
683: ds->scset = PETSC_TRUE;
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@C
688: DSGetSlepcSC - Gets the sorting criterion context.
690: Not Collective
692: Input Parameter:
693: . ds - the direct solver context
695: Output Parameter:
696: . sc - a pointer to the sorting criterion context
698: Note:
699: Not available in Fortran.
701: Level: developer
703: .seealso: [](sec:ds), `DSSetSlepcSC()`, `DSSort()`
704: @*/
705: PetscErrorCode DSGetSlepcSC(DS ds,SlepcSC *sc)
706: {
707: PetscFunctionBegin;
709: PetscAssertPointer(sc,2);
710: if (!ds->sc) PetscCall(PetscNew(&ds->sc));
711: *sc = ds->sc;
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: /*@
716: DSSetFromOptions - Sets `DS` options from the options database.
718: Collective
720: Input Parameter:
721: . ds - the direct solver context
723: Notes:
724: To see all options, run your program with the `-help` option.
726: Level: beginner
728: .seealso: [](sec:ds), `DSSetOptionsPrefix()`
729: @*/
730: PetscErrorCode DSSetFromOptions(DS ds)
731: {
732: PetscInt bs,meth;
733: PetscBool flag;
734: DSParallelType pmode;
736: PetscFunctionBegin;
738: PetscCall(DSRegisterAll());
739: /* Set default type (we do not allow changing it with -ds_type) */
740: if (!((PetscObject)ds)->type_name) PetscCall(DSSetType(ds,DSNHEP));
741: PetscObjectOptionsBegin((PetscObject)ds);
743: PetscCall(PetscOptionsInt("-ds_block_size","Block size for the dense system solver","DSSetBlockSize",ds->bs,&bs,&flag));
744: if (flag) PetscCall(DSSetBlockSize(ds,bs));
746: PetscCall(PetscOptionsInt("-ds_method","Method to be used for the dense system","DSSetMethod",ds->method,&meth,&flag));
747: if (flag) PetscCall(DSSetMethod(ds,meth));
749: PetscCall(PetscOptionsEnum("-ds_parallel","Operation mode in parallel runs","DSSetParallel",DSParallelTypes,(PetscEnum)ds->pmode,(PetscEnum*)&pmode,&flag));
750: if (flag) PetscCall(DSSetParallel(ds,pmode));
752: PetscTryTypeMethod(ds,setfromoptions,PetscOptionsObject);
753: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ds,PetscOptionsObject));
754: PetscOptionsEnd();
755: PetscFunctionReturn(PETSC_SUCCESS);
756: }
758: /*@
759: DSView - Prints the `DS` data structure.
761: Collective
763: Input Parameters:
764: + ds - the direct solver context
765: - viewer - optional visualization context
767: Notes:
768: The available visualization contexts include
769: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
770: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
771: first process opens the file; all other processes send their data to the
772: first one to print
774: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
775: to output to a specified file.
777: Use `DSViewFromOptions()` to allow the user to select many different `PetscViewerType`
778: and formats from the options database.
780: Level: beginner
782: .seealso: [](sec:ds), `DSCreate()`, `DSViewFromOptions()`, `DSViewMat()`
783: @*/
784: PetscErrorCode DSView(DS ds,PetscViewer viewer)
785: {
786: PetscBool isascii;
787: PetscViewerFormat format;
788: PetscMPIInt size;
790: PetscFunctionBegin;
792: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer));
794: PetscCheckSameComm(ds,1,viewer,2);
795: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
796: if (isascii) {
797: PetscCall(PetscViewerGetFormat(viewer,&format));
798: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ds,viewer));
799: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
800: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",DSParallelTypes[ds->pmode]));
801: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
802: PetscCall(PetscViewerASCIIPrintf(viewer," current state: %s\n",DSStateTypes[ds->state]));
803: PetscCall(PetscViewerASCIIPrintf(viewer," dimensions: ld=%" PetscInt_FMT ", n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT,ds->ld,ds->n,ds->l,ds->k));
804: if (ds->state==DS_STATE_TRUNCATED) PetscCall(PetscViewerASCIIPrintf(viewer,", t=%" PetscInt_FMT "\n",ds->t));
805: else PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
806: PetscCall(PetscViewerASCIIPrintf(viewer," flags:%s%s%s\n",ds->compact?" compact":"",ds->extrarow?" extrarow":"",ds->refined?" refined":""));
807: }
808: PetscCall(PetscViewerASCIIPushTab(viewer));
809: PetscTryTypeMethod(ds,view,viewer);
810: PetscCall(PetscViewerASCIIPopTab(viewer));
811: }
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: /*@
816: DSViewFromOptions - View (print) a `DS` object based on values in the options database.
818: Collective
820: Input Parameters:
821: + ds - the direct solver context
822: . obj - optional object that provides the options prefix used to query the options database
823: - name - command line option
825: Level: intermediate
827: .seealso: [](sec:ds), `DSView()`, `DSCreate()`, `PetscObjectViewFromOptions()`
828: @*/
829: PetscErrorCode DSViewFromOptions(DS ds,PetscObject obj,const char name[])
830: {
831: PetscFunctionBegin;
833: PetscCall(PetscObjectViewFromOptions((PetscObject)ds,obj,name));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: /*@
838: DSAllocate - Allocates memory for internal storage or matrices in `DS`.
840: Logically Collective
842: Input Parameters:
843: + ds - the direct solver context
844: - ld - leading dimension (maximum allowed dimension for the matrices, including
845: the extra row if present)
847: Note:
848: If the leading dimension is different from a previously set value, then
849: all matrices are destroyed with `DSReset()`.
851: Level: intermediate
853: .seealso: [](sec:ds), `DSGetLeadingDimension()`, `DSSetDimensions()`, `DSSetExtraRow()`, `DSReset()`, `DSReallocate()`
854: @*/
855: PetscErrorCode DSAllocate(DS ds,PetscInt ld)
856: {
857: PetscFunctionBegin;
861: PetscCheck(ld>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Leading dimension should be at least one");
862: if (ld!=ds->ld) {
863: PetscCall(PetscInfo(ds,"Allocating memory with leading dimension=%" PetscInt_FMT "\n",ld));
864: PetscCall(DSReset(ds));
865: ds->ld = ld;
866: PetscUseTypeMethod(ds,allocate,ld);
867: }
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: /*@
872: DSReallocate - Reallocates memory for internal storage or matrices in `DS`,
873: keeping the previously set data.
875: Logically Collective
877: Input Parameters:
878: + ds - the direct solver context
879: - ld - new leading dimension
881: Notes:
882: The new leading dimension must be larger than the previous one. The relevant
883: data previously set is copied over to the new data structures.
885: This operation is not available in all `DS` types.
887: Level: developer
889: .seealso: [](sec:ds), `DSAllocate()`
890: @*/
891: PetscErrorCode DSReallocate(DS ds,PetscInt ld)
892: {
893: PetscInt i;
895: PetscFunctionBegin;
899: PetscCheck(ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"DSAllocate() must be called first");
900: 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);
901: PetscCall(PetscInfo(ds,"Reallocating memory with new leading dimension=%" PetscInt_FMT "\n",ld));
902: PetscUseTypeMethod(ds,reallocate,ld);
903: for (i=ds->ld;i<ld;i++) ds->perm[i] = i;
904: ds->ld = ld;
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: /*@
909: DSReset - Resets the `DS` context to the initial state.
911: Collective
913: Input Parameter:
914: . ds - the direct solver context
916: Note:
917: All data structures with size depending on the leading dimension
918: of `DSAllocate()` are released.
920: Level: advanced
922: .seealso: [](sec:ds), `DSDestroy()`, `DSAllocate()`
923: @*/
924: PetscErrorCode DSReset(DS ds)
925: {
926: PetscInt i;
928: PetscFunctionBegin;
930: if (!ds) PetscFunctionReturn(PETSC_SUCCESS);
931: ds->state = DS_STATE_RAW;
932: ds->ld = 0;
933: ds->l = 0;
934: ds->n = 0;
935: ds->k = 0;
936: for (i=0;i<DS_NUM_MAT;i++) PetscCall(MatDestroy(&ds->omat[i]));
937: PetscCall(PetscFree(ds->perm));
938: PetscFunctionReturn(PETSC_SUCCESS);
939: }
941: /*@
942: DSDestroy - Destroys `DS` context that was created with `DSCreate()`.
944: Collective
946: Input Parameter:
947: . ds - the direct solver context
949: Level: beginner
951: .seealso: [](sec:ds), `DSCreate()`
952: @*/
953: PetscErrorCode DSDestroy(DS *ds)
954: {
955: PetscFunctionBegin;
956: if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
958: if (--((PetscObject)*ds)->refct > 0) { *ds = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
959: PetscCall(DSReset(*ds));
960: PetscTryTypeMethod(*ds,destroy);
961: PetscCall(PetscFree((*ds)->work));
962: PetscCall(PetscFree((*ds)->rwork));
963: PetscCall(PetscFree((*ds)->iwork));
964: if (!(*ds)->scset) PetscCall(PetscFree((*ds)->sc));
965: PetscCall(PetscHeaderDestroy(ds));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }
969: /*@C
970: DSRegister - Adds a direct solver to the `DS` package.
972: Not Collective
974: Input Parameters:
975: + name - name of a new user-defined `DS`
976: - function - routine to create context
978: Note:
979: `DSRegister()` may be called multiple times to add several user-defined
980: direct solvers.
982: Level: advanced
984: .seealso: [](sec:ds), `DSRegisterAll()`
985: @*/
986: PetscErrorCode DSRegister(const char *name,PetscErrorCode (*function)(DS))
987: {
988: PetscFunctionBegin;
989: PetscCall(DSInitializePackage());
990: PetscCall(PetscFunctionListAdd(&DSList,name,function));
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS);
995: SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS);
996: SLEPC_EXTERN PetscErrorCode DSCreate_GHEP(DS);
997: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS);
998: SLEPC_EXTERN PetscErrorCode DSCreate_GNHEP(DS);
999: SLEPC_EXTERN PetscErrorCode DSCreate_NHEPTS(DS);
1000: SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS);
1001: SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS);
1002: SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS);
1003: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS);
1004: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS);
1006: /*@C
1007: DSRegisterAll - Registers all of the direct solvers in the `DS` package.
1009: Not Collective
1011: Level: advanced
1013: .seealso: [](sec:ds), `DSRegister()`
1014: @*/
1015: PetscErrorCode DSRegisterAll(void)
1016: {
1017: PetscFunctionBegin;
1018: if (DSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1019: DSRegisterAllCalled = PETSC_TRUE;
1020: PetscCall(DSRegister(DSHEP,DSCreate_HEP));
1021: PetscCall(DSRegister(DSNHEP,DSCreate_NHEP));
1022: PetscCall(DSRegister(DSGHEP,DSCreate_GHEP));
1023: PetscCall(DSRegister(DSGHIEP,DSCreate_GHIEP));
1024: PetscCall(DSRegister(DSGNHEP,DSCreate_GNHEP));
1025: PetscCall(DSRegister(DSNHEPTS,DSCreate_NHEPTS));
1026: PetscCall(DSRegister(DSSVD,DSCreate_SVD));
1027: PetscCall(DSRegister(DSHSVD,DSCreate_HSVD));
1028: PetscCall(DSRegister(DSGSVD,DSCreate_GSVD));
1029: PetscCall(DSRegister(DSPEP,DSCreate_PEP));
1030: PetscCall(DSRegister(DSNEP,DSCreate_NEP));
1031: PetscFunctionReturn(PETSC_SUCCESS);
1032: }