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: Level: intermediate
274: .seealso: [](sec:ds), `DSSetType()`
275: @*/
276: PetscErrorCode DSGetType(DS ds,DSType *type)
277: {
278: PetscFunctionBegin;
280: PetscAssertPointer(type,2);
281: *type = ((PetscObject)ds)->type_name;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@
286: DSDuplicate - Creates a new direct solver object with the same options as
287: an existing one.
289: Collective
291: Input Parameter:
292: . ds - direct solver context
294: Output Parameter:
295: . dsnew - location to put the new `DS`
297: Notes:
298: `DSDuplicate()` DOES NOT COPY the matrices, and the new `DS` does not even have
299: internal arrays allocated. Use `DSAllocate()` to use the new `DS`.
301: The sorting criterion options are not copied, see `DSSetSlepcSC()`.
303: Level: intermediate
305: .seealso: [](sec:ds), `DSCreate()`, `DSAllocate()`, `DSSetSlepcSC()`
306: @*/
307: PetscErrorCode DSDuplicate(DS ds,DS *dsnew)
308: {
309: PetscFunctionBegin;
311: PetscAssertPointer(dsnew,2);
312: PetscCall(DSCreate(PetscObjectComm((PetscObject)ds),dsnew));
313: if (((PetscObject)ds)->type_name) PetscCall(DSSetType(*dsnew,((PetscObject)ds)->type_name));
314: (*dsnew)->method = ds->method;
315: (*dsnew)->compact = ds->compact;
316: (*dsnew)->refined = ds->refined;
317: (*dsnew)->extrarow = ds->extrarow;
318: (*dsnew)->bs = ds->bs;
319: (*dsnew)->pmode = ds->pmode;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@
324: DSSetMethod - Selects the method to be used to solve the problem.
326: Logically Collective
328: Input Parameters:
329: + ds - the direct solver context
330: - meth - an index identifying the method
332: Options Database Key:
333: . -ds_method \<meth\> - sets the method
335: Level: intermediate
337: .seealso: [](sec:ds), `DSGetMethod()`
338: @*/
339: PetscErrorCode DSSetMethod(DS ds,PetscInt meth)
340: {
341: PetscFunctionBegin;
344: PetscCheck(meth>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
345: PetscCheck(meth<=DS_MAX_SOLVE,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
346: ds->method = meth;
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*@
351: DSGetMethod - Gets the method currently used in the `DS`.
353: Not Collective
355: Input Parameter:
356: . ds - the direct solver context
358: Output Parameter:
359: . meth - identifier of the method
361: Level: intermediate
363: .seealso: [](sec:ds), `DSSetMethod()`
364: @*/
365: PetscErrorCode DSGetMethod(DS ds,PetscInt *meth)
366: {
367: PetscFunctionBegin;
369: PetscAssertPointer(meth,2);
370: *meth = ds->method;
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*@
375: DSSetParallel - Selects the mode of operation in parallel runs.
377: Logically Collective
379: Input Parameters:
380: + ds - the direct solver context
381: - pmode - the parallel mode
383: Options Database Key:
384: . -ds_parallel \<pmode\> - sets the parallel mode `redundant`,`synchronized`,`distributed`
386: Notes:
387: In the `redundant` parallel mode, all processes will make the computation
388: redundantly, starting from the same data, and producing the same result.
389: This result may be slightly different in the different processes if using a
390: multithreaded BLAS library, which may cause issues in ill-conditioned problems.
392: In the `synchronized` parallel mode, only the first MPI process performs the
393: computation and then the computed quantities are broadcast to the other
394: processes in the communicator. This communication is not done automatically,
395: an explicit call to `DSSynchronize()` is required.
397: The `distributed` parallel mode can be used in some `DS` types only, such
398: as the contour integral method of `DSNEP`. In this case, every MPI process
399: will be in charge of part of the computation.
401: Level: advanced
403: .seealso: [](sec:ds), `DSSynchronize()`, `DSGetParallel()`
404: @*/
405: PetscErrorCode DSSetParallel(DS ds,DSParallelType pmode)
406: {
407: PetscFunctionBegin;
410: ds->pmode = pmode;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*@
415: DSGetParallel - Gets the mode of operation in parallel runs.
417: Not Collective
419: Input Parameter:
420: . ds - the direct solver context
422: Output Parameter:
423: . pmode - the parallel mode
425: Level: advanced
427: .seealso: [](sec:ds), `DSSetParallel()`
428: @*/
429: PetscErrorCode DSGetParallel(DS ds,DSParallelType *pmode)
430: {
431: PetscFunctionBegin;
433: PetscAssertPointer(pmode,2);
434: *pmode = ds->pmode;
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: DSSetCompact - Switch to compact storage of matrices.
441: Logically Collective
443: Input Parameters:
444: + ds - the direct solver context
445: - comp - a boolean flag
447: Notes:
448: Compact storage is used in some `DS` types such as `DSHEP` when the matrix
449: is tridiagonal. This flag can be used to indicate whether the user
450: provides the matrix entries via the compact form (the tridiagonal `DS_MAT_T`)
451: or the non-compact one (`DS_MAT_A`).
453: The default is `PETSC_FALSE`.
455: Level: advanced
457: .seealso: [](sec:ds), `DSGetCompact()`
458: @*/
459: PetscErrorCode DSSetCompact(DS ds,PetscBool comp)
460: {
461: PetscFunctionBegin;
464: if (ds->compact != comp && ds->ld) PetscTryTypeMethod(ds,setcompact,comp);
465: ds->compact = comp;
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: DSGetCompact - Gets the compact storage flag.
472: Not Collective
474: Input Parameter:
475: . ds - the direct solver context
477: Output Parameter:
478: . comp - the flag
480: Level: advanced
482: .seealso: [](sec:ds), `DSSetCompact()`
483: @*/
484: PetscErrorCode DSGetCompact(DS ds,PetscBool *comp)
485: {
486: PetscFunctionBegin;
488: PetscAssertPointer(comp,2);
489: *comp = ds->compact;
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@
494: DSSetExtraRow - Sets a flag to indicate that the matrix has one extra
495: row.
497: Logically Collective
499: Input Parameters:
500: + ds - the direct solver context
501: - ext - a boolean flag
503: Notes:
504: In Krylov methods it is useful that the matrix representing the direct solver
505: has one extra row, i.e., has dimension $(n+1) \times n$. If this flag is activated, all
506: transformations applied to the right of the matrix also affect this additional
507: row. In that case, $(n+1)$ must be less or equal than the leading dimension.
509: The default is `PETSC_FALSE`.
511: Level: advanced
513: .seealso: [](sec:ds), `DSSolve()`, `DSAllocate()`, `DSGetExtraRow()`
514: @*/
515: PetscErrorCode DSSetExtraRow(DS ds,PetscBool ext)
516: {
517: PetscFunctionBegin;
520: PetscCheck(!ext || ds->n==0 || ds->n!=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot set extra row after setting n=ld");
521: ds->extrarow = ext;
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /*@
526: DSGetExtraRow - Gets the extra row flag.
528: Not Collective
530: Input Parameter:
531: . ds - the direct solver context
533: Output Parameter:
534: . ext - the flag
536: Level: advanced
538: .seealso: [](sec:ds), `DSSetExtraRow()`
539: @*/
540: PetscErrorCode DSGetExtraRow(DS ds,PetscBool *ext)
541: {
542: PetscFunctionBegin;
544: PetscAssertPointer(ext,2);
545: *ext = ds->extrarow;
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: /*@
550: DSSetRefined - Sets a flag to indicate that refined vectors must be
551: computed.
553: Logically Collective
555: Input Parameters:
556: + ds - the direct solver context
557: - ref - a boolean flag
559: Notes:
560: Normally the vectors returned in `DS_MAT_X` are eigenvectors of the
561: projected matrix. With this flag activated, `DSVectors()` will return
562: the right singular vector of the smallest singular value of matrix
563: $\tilde{A}-\theta I$, where $\tilde{A}$ is the extended $(n+1)\times n$
564: matrix and $\theta$ is the Ritz value. This is used in the refined Ritz
565: approximation {cite:p}`Jia97`.
567: The default is `PETSC_FALSE`.
569: Level: advanced
571: .seealso: [](sec:ds), `DSVectors()`, `DSGetRefined()`
572: @*/
573: PetscErrorCode DSSetRefined(DS ds,PetscBool ref)
574: {
575: PetscFunctionBegin;
578: ds->refined = ref;
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*@
583: DSGetRefined - Gets the refined vectors flag.
585: Not Collective
587: Input Parameter:
588: . ds - the direct solver context
590: Output Parameter:
591: . ref - the flag
593: Level: advanced
595: .seealso: [](sec:ds), `DSSetRefined()`
596: @*/
597: PetscErrorCode DSGetRefined(DS ds,PetscBool *ref)
598: {
599: PetscFunctionBegin;
601: PetscAssertPointer(ref,2);
602: *ref = ds->refined;
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /*@
607: DSSetBlockSize - Sets the block size.
609: Logically Collective
611: Input Parameters:
612: + ds - the direct solver context
613: - bs - the block size
615: Options Database Key:
616: . -ds_block_size \<bs\> - sets the block size
618: Level: intermediate
620: .seealso: [](sec:ds), `DSGetBlockSize()`
621: @*/
622: PetscErrorCode DSSetBlockSize(DS ds,PetscInt bs)
623: {
624: PetscFunctionBegin;
627: PetscCheck(bs>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The block size must be at least one");
628: ds->bs = bs;
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /*@
633: DSGetBlockSize - Gets the block size.
635: Not Collective
637: Input Parameter:
638: . ds - the direct solver context
640: Output Parameter:
641: . bs - block size
643: Level: intermediate
645: .seealso: [](sec:ds), `DSSetBlockSize()`
646: @*/
647: PetscErrorCode DSGetBlockSize(DS ds,PetscInt *bs)
648: {
649: PetscFunctionBegin;
651: PetscAssertPointer(bs,2);
652: *bs = ds->bs;
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*@C
657: DSSetSlepcSC - Sets the sorting criterion context.
659: Logically Collective
661: Input Parameters:
662: + ds - the direct solver context
663: - sc - a pointer to the sorting criterion context
665: Note:
666: Not available in Fortran.
668: Level: developer
670: .seealso: [](sec:ds), `DSGetSlepcSC()`, `DSSort()`
671: @*/
672: PetscErrorCode DSSetSlepcSC(DS ds,SlepcSC sc)
673: {
674: PetscFunctionBegin;
676: PetscAssertPointer(sc,2);
677: if (ds->sc && !ds->scset) PetscCall(PetscFree(ds->sc));
678: ds->sc = sc;
679: ds->scset = PETSC_TRUE;
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: /*@C
684: DSGetSlepcSC - Gets the sorting criterion context.
686: Not Collective
688: Input Parameter:
689: . ds - the direct solver context
691: Output Parameter:
692: . sc - a pointer to the sorting criterion context
694: Note:
695: Not available in Fortran.
697: Level: developer
699: .seealso: [](sec:ds), `DSSetSlepcSC()`, `DSSort()`
700: @*/
701: PetscErrorCode DSGetSlepcSC(DS ds,SlepcSC *sc)
702: {
703: PetscFunctionBegin;
705: PetscAssertPointer(sc,2);
706: if (!ds->sc) PetscCall(PetscNew(&ds->sc));
707: *sc = ds->sc;
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: /*@
712: DSSetFromOptions - Sets `DS` options from the options database.
714: Collective
716: Input Parameter:
717: . ds - the direct solver context
719: Notes:
720: To see all options, run your program with the `-help` option.
722: Level: beginner
724: .seealso: [](sec:ds), `DSSetOptionsPrefix()`
725: @*/
726: PetscErrorCode DSSetFromOptions(DS ds)
727: {
728: PetscInt bs,meth;
729: PetscBool flag;
730: DSParallelType pmode;
732: PetscFunctionBegin;
734: PetscCall(DSRegisterAll());
735: /* Set default type (we do not allow changing it with -ds_type) */
736: if (!((PetscObject)ds)->type_name) PetscCall(DSSetType(ds,DSNHEP));
737: PetscObjectOptionsBegin((PetscObject)ds);
739: PetscCall(PetscOptionsInt("-ds_block_size","Block size for the dense system solver","DSSetBlockSize",ds->bs,&bs,&flag));
740: if (flag) PetscCall(DSSetBlockSize(ds,bs));
742: PetscCall(PetscOptionsInt("-ds_method","Method to be used for the dense system","DSSetMethod",ds->method,&meth,&flag));
743: if (flag) PetscCall(DSSetMethod(ds,meth));
745: PetscCall(PetscOptionsEnum("-ds_parallel","Operation mode in parallel runs","DSSetParallel",DSParallelTypes,(PetscEnum)ds->pmode,(PetscEnum*)&pmode,&flag));
746: if (flag) PetscCall(DSSetParallel(ds,pmode));
748: PetscTryTypeMethod(ds,setfromoptions,PetscOptionsObject);
749: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ds,PetscOptionsObject));
750: PetscOptionsEnd();
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: /*@
755: DSView - Prints the `DS` data structure.
757: Collective
759: Input Parameters:
760: + ds - the direct solver context
761: - viewer - optional visualization context
763: Notes:
764: The available visualization contexts include
765: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
766: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
767: first process opens the file; all other processes send their data to the
768: first one to print
770: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
771: to output to a specified file.
773: Use `DSViewFromOptions()` to allow the user to select many different `PetscViewerType`
774: and formats from the options database.
776: Level: beginner
778: .seealso: [](sec:ds), `DSCreate()`, `DSViewFromOptions()`, `DSViewMat()`
779: @*/
780: PetscErrorCode DSView(DS ds,PetscViewer viewer)
781: {
782: PetscBool isascii;
783: PetscViewerFormat format;
784: PetscMPIInt size;
786: PetscFunctionBegin;
788: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer));
790: PetscCheckSameComm(ds,1,viewer,2);
791: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
792: if (isascii) {
793: PetscCall(PetscViewerGetFormat(viewer,&format));
794: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ds,viewer));
795: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
796: if (size>1) PetscCall(PetscViewerASCIIPrintf(viewer," parallel operation mode: %s\n",DSParallelTypes[ds->pmode]));
797: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
798: PetscCall(PetscViewerASCIIPrintf(viewer," current state: %s\n",DSStateTypes[ds->state]));
799: 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: if (ds->state==DS_STATE_TRUNCATED) PetscCall(PetscViewerASCIIPrintf(viewer,", t=%" PetscInt_FMT "\n",ds->t));
801: else PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
802: PetscCall(PetscViewerASCIIPrintf(viewer," flags:%s%s%s\n",ds->compact?" compact":"",ds->extrarow?" extrarow":"",ds->refined?" refined":""));
803: }
804: PetscCall(PetscViewerASCIIPushTab(viewer));
805: PetscTryTypeMethod(ds,view,viewer);
806: PetscCall(PetscViewerASCIIPopTab(viewer));
807: }
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: /*@
812: DSViewFromOptions - View (print) a `DS` object based on values in the options database.
814: Collective
816: Input Parameters:
817: + ds - the direct solver context
818: . obj - optional object that provides the options prefix used to query the options database
819: - name - command line option
821: Level: intermediate
823: .seealso: [](sec:ds), `DSView()`, `DSCreate()`, `PetscObjectViewFromOptions()`
824: @*/
825: PetscErrorCode DSViewFromOptions(DS ds,PetscObject obj,const char name[])
826: {
827: PetscFunctionBegin;
829: PetscCall(PetscObjectViewFromOptions((PetscObject)ds,obj,name));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: /*@
834: DSAllocate - Allocates memory for internal storage or matrices in `DS`.
836: Logically Collective
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)
843: Note:
844: If the leading dimension is different from a previously set value, then
845: all matrices are destroyed with `DSReset()`.
847: Level: intermediate
849: .seealso: [](sec:ds), `DSGetLeadingDimension()`, `DSSetDimensions()`, `DSSetExtraRow()`, `DSReset()`, `DSReallocate()`
850: @*/
851: PetscErrorCode DSAllocate(DS ds,PetscInt ld)
852: {
853: PetscFunctionBegin;
857: PetscCheck(ld>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Leading dimension should be at least one");
858: if (ld!=ds->ld) {
859: PetscCall(PetscInfo(ds,"Allocating memory with leading dimension=%" PetscInt_FMT "\n",ld));
860: PetscCall(DSReset(ds));
861: ds->ld = ld;
862: PetscUseTypeMethod(ds,allocate,ld);
863: }
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: DSReallocate - Reallocates memory for internal storage or matrices in `DS`,
869: keeping the previously set data.
871: Logically Collective
873: Input Parameters:
874: + ds - the direct solver context
875: - ld - new leading dimension
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.
881: This operation is not available in all `DS` types.
883: Level: developer
885: .seealso: [](sec:ds), `DSAllocate()`
886: @*/
887: PetscErrorCode DSReallocate(DS ds,PetscInt ld)
888: {
889: PetscFunctionBegin;
893: PetscCheck(ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"DSAllocate() must be called first");
894: 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: PetscCall(PetscInfo(ds,"Reallocating memory with new leading dimension=%" PetscInt_FMT "\n",ld));
896: PetscUseTypeMethod(ds,reallocate,ld);
897: ds->ld = ld;
898: PetscFunctionReturn(PETSC_SUCCESS);
899: }
901: /*@
902: DSReset - Resets the `DS` context to the initial state.
904: Collective
906: Input Parameter:
907: . ds - the direct solver context
909: Note:
910: All data structures with size depending on the leading dimension
911: of `DSAllocate()` are released.
913: Level: advanced
915: .seealso: [](sec:ds), `DSDestroy()`, `DSAllocate()`
916: @*/
917: PetscErrorCode DSReset(DS ds)
918: {
919: PetscInt i;
921: PetscFunctionBegin;
923: if (!ds) PetscFunctionReturn(PETSC_SUCCESS);
924: ds->state = DS_STATE_RAW;
925: ds->ld = 0;
926: ds->l = 0;
927: ds->n = 0;
928: ds->k = 0;
929: for (i=0;i<DS_NUM_MAT;i++) PetscCall(MatDestroy(&ds->omat[i]));
930: PetscCall(PetscFree(ds->perm));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@
935: DSDestroy - Destroys `DS` context that was created with `DSCreate()`.
937: Collective
939: Input Parameter:
940: . ds - the direct solver context
942: Level: beginner
944: .seealso: [](sec:ds), `DSCreate()`
945: @*/
946: PetscErrorCode DSDestroy(DS *ds)
947: {
948: PetscFunctionBegin;
949: if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
951: if (--((PetscObject)*ds)->refct > 0) { *ds = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
952: PetscCall(DSReset(*ds));
953: PetscTryTypeMethod(*ds,destroy);
954: PetscCall(PetscFree((*ds)->work));
955: PetscCall(PetscFree((*ds)->rwork));
956: PetscCall(PetscFree((*ds)->iwork));
957: if (!(*ds)->scset) PetscCall(PetscFree((*ds)->sc));
958: PetscCall(PetscHeaderDestroy(ds));
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: /*@C
963: DSRegister - Adds a direct solver to the `DS` package.
965: Not Collective
967: Input Parameters:
968: + name - name of a new user-defined `DS`
969: - function - routine to create context
971: Note:
972: `DSRegister()` may be called multiple times to add several user-defined
973: direct solvers.
975: Level: advanced
977: .seealso: [](sec:ds), `DSRegisterAll()`
978: @*/
979: PetscErrorCode DSRegister(const char *name,PetscErrorCode (*function)(DS))
980: {
981: PetscFunctionBegin;
982: PetscCall(DSInitializePackage());
983: PetscCall(PetscFunctionListAdd(&DSList,name,function));
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
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);
999: /*@C
1000: DSRegisterAll - Registers all of the direct solvers in the `DS` package.
1002: Not Collective
1004: Level: advanced
1006: .seealso: [](sec:ds), `DSRegister()`
1007: @*/
1008: PetscErrorCode DSRegisterAll(void)
1009: {
1010: PetscFunctionBegin;
1011: if (DSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1012: DSRegisterAllCalled = PETSC_TRUE;
1013: PetscCall(DSRegister(DSHEP,DSCreate_HEP));
1014: PetscCall(DSRegister(DSNHEP,DSCreate_NHEP));
1015: PetscCall(DSRegister(DSGHEP,DSCreate_GHEP));
1016: PetscCall(DSRegister(DSGHIEP,DSCreate_GHIEP));
1017: PetscCall(DSRegister(DSGNHEP,DSCreate_GNHEP));
1018: PetscCall(DSRegister(DSNHEPTS,DSCreate_NHEPTS));
1019: PetscCall(DSRegister(DSSVD,DSCreate_SVD));
1020: PetscCall(DSRegister(DSHSVD,DSCreate_HSVD));
1021: PetscCall(DSRegister(DSGSVD,DSCreate_GSVD));
1022: PetscCall(DSRegister(DSPEP,DSCreate_PEP));
1023: PetscCall(DSRegister(DSNEP,DSCreate_NEP));
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }