Actual source code: dsbasic.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  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: SlepcFinalize()
 34: @*/
 35: PetscErrorCode DSFinalizePackage(void)
 36: {
 37:   PetscFunctionListDestroy(&DSList);
 38:   DSPackageInitialized = PETSC_FALSE;
 39:   DSRegisterAllCalled  = PETSC_FALSE;
 40:   return 0;
 41: }

 43: /*@C
 44:   DSInitializePackage - This function initializes everything in the DS package.
 45:   It is called from PetscDLLibraryRegister() when using dynamic libraries, and
 46:   on the first call to DSCreate() when using static libraries.

 48:   Level: developer

 50: .seealso: SlepcInitialize()
 51: @*/
 52: PetscErrorCode DSInitializePackage(void)
 53: {
 54:   char           logList[256];
 55:   PetscBool      opt,pkg;
 56:   PetscClassId   classids[1];

 58:   if (DSPackageInitialized) return 0;
 59:   DSPackageInitialized = PETSC_TRUE;
 60:   /* Register Classes */
 61:   PetscClassIdRegister("Direct Solver",&DS_CLASSID);
 62:   /* Register Constructors */
 63:   DSRegisterAll();
 64:   /* Register Events */
 65:   PetscLogEventRegister("DSSolve",DS_CLASSID,&DS_Solve);
 66:   PetscLogEventRegister("DSVectors",DS_CLASSID,&DS_Vectors);
 67:   PetscLogEventRegister("DSSynchronize",DS_CLASSID,&DS_Synchronize);
 68:   PetscLogEventRegister("DSOther",DS_CLASSID,&DS_Other);
 69:   /* Process Info */
 70:   classids[0] = DS_CLASSID;
 71:   PetscInfoProcessClass("ds",1,&classids[0]);
 72:   /* Process summary exclusions */
 73:   PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
 74:   if (opt) {
 75:     PetscStrInList("ds",logList,',',&pkg);
 76:     if (pkg) PetscLogEventDeactivateClass(DS_CLASSID);
 77:   }
 78:   /* Register package finalizer */
 79:   PetscRegisterFinalize(DSFinalizePackage);
 80:   return 0;
 81: }

 83: /*@
 84:    DSCreate - Creates a DS context.

 86:    Collective

 88:    Input Parameter:
 89: .  comm - MPI communicator

 91:    Output Parameter:
 92: .  newds - location to put the DS context

 94:    Level: beginner

 96:    Note:
 97:    DS objects are not intended for normal users but only for
 98:    advanced user that for instance implement their own solvers.

100: .seealso: DSDestroy(), DS
101: @*/
102: PetscErrorCode DSCreate(MPI_Comm comm,DS *newds)
103: {
104:   DS             ds;
105:   PetscInt       i;

108:   *newds = NULL;
109:   DSInitializePackage();
110:   SlepcHeaderCreate(ds,DS_CLASSID,"DS","Direct Solver (or Dense System)","DS",comm,DSDestroy,DSView);

112:   ds->state         = DS_STATE_RAW;
113:   ds->method        = 0;
114:   ds->compact       = PETSC_FALSE;
115:   ds->refined       = PETSC_FALSE;
116:   ds->extrarow      = PETSC_FALSE;
117:   ds->ld            = 0;
118:   ds->l             = 0;
119:   ds->n             = 0;
120:   ds->k             = 0;
121:   ds->t             = 0;
122:   ds->bs            = 1;
123:   ds->sc            = NULL;
124:   ds->pmode         = DS_PARALLEL_REDUNDANT;

126:   for (i=0;i<DS_NUM_MAT;i++) ds->omat[i] = NULL;
127:   ds->perm          = NULL;
128:   ds->data          = NULL;
129:   ds->scset         = PETSC_FALSE;
130:   ds->work          = NULL;
131:   ds->rwork         = NULL;
132:   ds->iwork         = NULL;
133:   ds->lwork         = 0;
134:   ds->lrwork        = 0;
135:   ds->liwork        = 0;

137:   *newds = ds;
138:   return 0;
139: }

141: /*@C
142:    DSSetOptionsPrefix - Sets the prefix used for searching for all
143:    DS options in the database.

145:    Logically Collective on ds

147:    Input Parameters:
148: +  ds - the direct solver context
149: -  prefix - the prefix string to prepend to all DS option requests

151:    Notes:
152:    A hyphen (-) must NOT be given at the beginning of the prefix name.
153:    The first character of all runtime options is AUTOMATICALLY the
154:    hyphen.

156:    Level: advanced

158: .seealso: DSAppendOptionsPrefix()
159: @*/
160: PetscErrorCode DSSetOptionsPrefix(DS ds,const char *prefix)
161: {
163:   PetscObjectSetOptionsPrefix((PetscObject)ds,prefix);
164:   return 0;
165: }

167: /*@C
168:    DSAppendOptionsPrefix - Appends to the prefix used for searching for all
169:    DS options in the database.

171:    Logically Collective on ds

173:    Input Parameters:
174: +  ds - the direct solver context
175: -  prefix - the prefix string to prepend to all DS option requests

177:    Notes:
178:    A hyphen (-) must NOT be given at the beginning of the prefix name.
179:    The first character of all runtime options is AUTOMATICALLY the hyphen.

181:    Level: advanced

183: .seealso: DSSetOptionsPrefix()
184: @*/
185: PetscErrorCode DSAppendOptionsPrefix(DS ds,const char *prefix)
186: {
188:   PetscObjectAppendOptionsPrefix((PetscObject)ds,prefix);
189:   return 0;
190: }

192: /*@C
193:    DSGetOptionsPrefix - Gets the prefix used for searching for all
194:    DS options in the database.

196:    Not Collective

198:    Input Parameters:
199: .  ds - the direct solver context

201:    Output Parameters:
202: .  prefix - pointer to the prefix string used is returned

204:    Note:
205:    On the Fortran side, the user should pass in a string 'prefix' of
206:    sufficient length to hold the prefix.

208:    Level: advanced

210: .seealso: DSSetOptionsPrefix(), DSAppendOptionsPrefix()
211: @*/
212: PetscErrorCode DSGetOptionsPrefix(DS ds,const char *prefix[])
213: {
216:   PetscObjectGetOptionsPrefix((PetscObject)ds,prefix);
217:   return 0;
218: }

220: /*@C
221:    DSSetType - Selects the type for the DS object.

223:    Logically Collective on ds

225:    Input Parameters:
226: +  ds   - the direct solver context
227: -  type - a known type

229:    Level: intermediate

231: .seealso: DSGetType()
232: @*/
233: PetscErrorCode DSSetType(DS ds,DSType type)
234: {
235:   PetscErrorCode (*r)(DS);
236:   PetscBool      match;


241:   PetscObjectTypeCompare((PetscObject)ds,type,&match);
242:   if (match) return 0;

244:   PetscFunctionListFind(DSList,type,&r);

247:   PetscMemzero(ds->ops,sizeof(struct _DSOps));

249:   PetscObjectChangeTypeName((PetscObject)ds,type);
250:   (*r)(ds);
251:   return 0;
252: }

254: /*@C
255:    DSGetType - Gets the DS type name (as a string) from the DS context.

257:    Not Collective

259:    Input Parameter:
260: .  ds - the direct solver context

262:    Output Parameter:
263: .  type - name of the direct solver

265:    Level: intermediate

267: .seealso: DSSetType()
268: @*/
269: PetscErrorCode DSGetType(DS ds,DSType *type)
270: {
273:   *type = ((PetscObject)ds)->type_name;
274:   return 0;
275: }

277: /*@
278:    DSDuplicate - Creates a new direct solver object with the same options as
279:    an existing one.

281:    Collective on ds

283:    Input Parameter:
284: .  ds - direct solver context

286:    Output Parameter:
287: .  dsnew - location to put the new DS

289:    Notes:
290:    DSDuplicate() DOES NOT COPY the matrices, and the new DS does not even have
291:    internal arrays allocated. Use DSAllocate() to use the new DS.

293:    The sorting criterion options are not copied, see DSSetSlepcSC().

295:    Level: intermediate

297: .seealso: DSCreate(), DSAllocate(), DSSetSlepcSC()
298: @*/
299: PetscErrorCode DSDuplicate(DS ds,DS *dsnew)
300: {
303:   DSCreate(PetscObjectComm((PetscObject)ds),dsnew);
304:   if (((PetscObject)ds)->type_name) DSSetType(*dsnew,((PetscObject)ds)->type_name);
305:   (*dsnew)->method   = ds->method;
306:   (*dsnew)->compact  = ds->compact;
307:   (*dsnew)->refined  = ds->refined;
308:   (*dsnew)->extrarow = ds->extrarow;
309:   (*dsnew)->bs       = ds->bs;
310:   (*dsnew)->pmode    = ds->pmode;
311:   return 0;
312: }

314: /*@
315:    DSSetMethod - Selects the method to be used to solve the problem.

317:    Logically Collective on ds

319:    Input Parameters:
320: +  ds   - the direct solver context
321: -  meth - an index identifying the method

323:    Options Database Key:
324: .  -ds_method <meth> - Sets the method

326:    Level: intermediate

328: .seealso: DSGetMethod()
329: @*/
330: PetscErrorCode DSSetMethod(DS ds,PetscInt meth)
331: {
336:   ds->method = meth;
337:   return 0;
338: }

340: /*@
341:    DSGetMethod - Gets the method currently used in the DS.

343:    Not Collective

345:    Input Parameter:
346: .  ds - the direct solver context

348:    Output Parameter:
349: .  meth - identifier of the method

351:    Level: intermediate

353: .seealso: DSSetMethod()
354: @*/
355: PetscErrorCode DSGetMethod(DS ds,PetscInt *meth)
356: {
359:   *meth = ds->method;
360:   return 0;
361: }

363: /*@
364:    DSSetParallel - Selects the mode of operation in parallel runs.

366:    Logically Collective on ds

368:    Input Parameters:
369: +  ds    - the direct solver context
370: -  pmode - the parallel mode

372:    Options Database Key:
373: .  -ds_parallel <mode> - Sets the parallel mode, 'redundant', 'synchronized'
374:    or 'distributed'

376:    Notes:
377:    In the 'redundant' parallel mode, all processes will make the computation
378:    redundantly, starting from the same data, and producing the same result.
379:    This result may be slightly different in the different processes if using a
380:    multithreaded BLAS library, which may cause issues in ill-conditioned problems.

382:    In the 'synchronized' parallel mode, only the first MPI process performs the
383:    computation and then the computed quantities are broadcast to the other
384:    processes in the communicator. This communication is not done automatically,
385:    an explicit call to DSSynchronize() is required.

387:    The 'distributed' parallel mode can be used in some DS types only, such
388:    as the contour integral method of DSNEP. In this case, every MPI process
389:    will be in charge of part of the computation.

391:    Level: advanced

393: .seealso: DSSynchronize(), DSGetParallel()
394: @*/
395: PetscErrorCode DSSetParallel(DS ds,DSParallelType pmode)
396: {
399:   ds->pmode = pmode;
400:   return 0;
401: }

403: /*@
404:    DSGetParallel - Gets the mode of operation in parallel runs.

406:    Not Collective

408:    Input Parameter:
409: .  ds - the direct solver context

411:    Output Parameter:
412: .  pmode - the parallel mode

414:    Level: advanced

416: .seealso: DSSetParallel()
417: @*/
418: PetscErrorCode DSGetParallel(DS ds,DSParallelType *pmode)
419: {
422:   *pmode = ds->pmode;
423:   return 0;
424: }

426: /*@
427:    DSSetCompact - Switch to compact storage of matrices.

429:    Logically Collective on ds

431:    Input Parameters:
432: +  ds   - the direct solver context
433: -  comp - a boolean flag

435:    Notes:
436:    Compact storage is used in some DS types such as DSHEP when the matrix
437:    is tridiagonal. This flag can be used to indicate whether the user
438:    provides the matrix entries via the compact form (the tridiagonal DS_MAT_T)
439:    or the non-compact one (DS_MAT_A).

441:    The default is PETSC_FALSE.

443:    Level: advanced

445: .seealso: DSGetCompact()
446: @*/
447: PetscErrorCode DSSetCompact(DS ds,PetscBool comp)
448: {
451:   ds->compact = comp;
452:   return 0;
453: }

455: /*@
456:    DSGetCompact - Gets the compact storage flag.

458:    Not Collective

460:    Input Parameter:
461: .  ds - the direct solver context

463:    Output Parameter:
464: .  comp - the flag

466:    Level: advanced

468: .seealso: DSSetCompact()
469: @*/
470: PetscErrorCode DSGetCompact(DS ds,PetscBool *comp)
471: {
474:   *comp = ds->compact;
475:   return 0;
476: }

478: /*@
479:    DSSetExtraRow - Sets a flag to indicate that the matrix has one extra
480:    row.

482:    Logically Collective on ds

484:    Input Parameters:
485: +  ds  - the direct solver context
486: -  ext - a boolean flag

488:    Notes:
489:    In Krylov methods it is useful that the matrix representing the direct solver
490:    has one extra row, i.e., has dimension (n+1) x n. If this flag is activated, all
491:    transformations applied to the right of the matrix also affect this additional
492:    row. In that case, (n+1) must be less or equal than the leading dimension.

494:    The default is PETSC_FALSE.

496:    Level: advanced

498: .seealso: DSSolve(), DSAllocate(), DSGetExtraRow()
499: @*/
500: PetscErrorCode DSSetExtraRow(DS ds,PetscBool ext)
501: {
505:   ds->extrarow = ext;
506:   return 0;
507: }

509: /*@
510:    DSGetExtraRow - Gets the extra row flag.

512:    Not Collective

514:    Input Parameter:
515: .  ds - the direct solver context

517:    Output Parameter:
518: .  ext - the flag

520:    Level: advanced

522: .seealso: DSSetExtraRow()
523: @*/
524: PetscErrorCode DSGetExtraRow(DS ds,PetscBool *ext)
525: {
528:   *ext = ds->extrarow;
529:   return 0;
530: }

532: /*@
533:    DSSetRefined - Sets a flag to indicate that refined vectors must be
534:    computed.

536:    Logically Collective on ds

538:    Input Parameters:
539: +  ds  - the direct solver context
540: -  ref - a boolean flag

542:    Notes:
543:    Normally the vectors returned in DS_MAT_X are eigenvectors of the
544:    projected matrix. With this flag activated, DSVectors() will return
545:    the right singular vector of the smallest singular value of matrix
546:    \tilde{A}-theta*I, where \tilde{A} is the extended (n+1)xn matrix
547:    and theta is the Ritz value. This is used in the refined Ritz
548:    approximation.

550:    The default is PETSC_FALSE.

552:    Level: advanced

554: .seealso: DSVectors(), DSGetRefined()
555: @*/
556: PetscErrorCode DSSetRefined(DS ds,PetscBool ref)
557: {
560:   ds->refined = ref;
561:   return 0;
562: }

564: /*@
565:    DSGetRefined - Gets the refined vectors flag.

567:    Not Collective

569:    Input Parameter:
570: .  ds - the direct solver context

572:    Output Parameter:
573: .  ref - the flag

575:    Level: advanced

577: .seealso: DSSetRefined()
578: @*/
579: PetscErrorCode DSGetRefined(DS ds,PetscBool *ref)
580: {
583:   *ref = ds->refined;
584:   return 0;
585: }

587: /*@
588:    DSSetBlockSize - Sets the block size.

590:    Logically Collective on ds

592:    Input Parameters:
593: +  ds - the direct solver context
594: -  bs - the block size

596:    Options Database Key:
597: .  -ds_block_size <bs> - Sets the block size

599:    Level: intermediate

601: .seealso: DSGetBlockSize()
602: @*/
603: PetscErrorCode DSSetBlockSize(DS ds,PetscInt bs)
604: {
608:   ds->bs = bs;
609:   return 0;
610: }

612: /*@
613:    DSGetBlockSize - Gets the block size.

615:    Not Collective

617:    Input Parameter:
618: .  ds - the direct solver context

620:    Output Parameter:
621: .  bs - block size

623:    Level: intermediate

625: .seealso: DSSetBlockSize()
626: @*/
627: PetscErrorCode DSGetBlockSize(DS ds,PetscInt *bs)
628: {
631:   *bs = ds->bs;
632:   return 0;
633: }

635: /*@C
636:    DSSetSlepcSC - Sets the sorting criterion context.

638:    Not Collective

640:    Input Parameters:
641: +  ds - the direct solver context
642: -  sc - a pointer to the sorting criterion context

644:    Level: developer

646: .seealso: DSGetSlepcSC(), DSSort()
647: @*/
648: PetscErrorCode DSSetSlepcSC(DS ds,SlepcSC sc)
649: {
652:   if (ds->sc && !ds->scset) PetscFree(ds->sc);
653:   ds->sc    = sc;
654:   ds->scset = PETSC_TRUE;
655:   return 0;
656: }

658: /*@C
659:    DSGetSlepcSC - Gets the sorting criterion context.

661:    Not Collective

663:    Input Parameter:
664: .  ds - the direct solver context

666:    Output Parameters:
667: .  sc - a pointer to the sorting criterion context

669:    Level: developer

671: .seealso: DSSetSlepcSC(), DSSort()
672: @*/
673: PetscErrorCode DSGetSlepcSC(DS ds,SlepcSC *sc)
674: {
677:   if (!ds->sc) PetscNew(&ds->sc);
678:   *sc = ds->sc;
679:   return 0;
680: }

682: /*@
683:    DSSetFromOptions - Sets DS options from the options database.

685:    Collective on ds

687:    Input Parameters:
688: .  ds - the direct solver context

690:    Notes:
691:    To see all options, run your program with the -help option.

693:    Level: beginner

695: .seealso: DSSetOptionsPrefix()
696: @*/
697: PetscErrorCode DSSetFromOptions(DS ds)
698: {
699:   PetscInt       bs,meth;
700:   PetscBool      flag;
701:   DSParallelType pmode;

704:   DSRegisterAll();
705:   /* Set default type (we do not allow changing it with -ds_type) */
706:   if (!((PetscObject)ds)->type_name) DSSetType(ds,DSNHEP);
707:   PetscObjectOptionsBegin((PetscObject)ds);

709:     PetscOptionsInt("-ds_block_size","Block size for the dense system solver","DSSetBlockSize",ds->bs,&bs,&flag);
710:     if (flag) DSSetBlockSize(ds,bs);

712:     PetscOptionsInt("-ds_method","Method to be used for the dense system","DSSetMethod",ds->method,&meth,&flag);
713:     if (flag) DSSetMethod(ds,meth);

715:     PetscOptionsEnum("-ds_parallel","Operation mode in parallel runs","DSSetParallel",DSParallelTypes,(PetscEnum)ds->pmode,(PetscEnum*)&pmode,&flag);
716:     if (flag) DSSetParallel(ds,pmode);

718:     PetscTryTypeMethod(ds,setfromoptions,PetscOptionsObject);
719:     PetscObjectProcessOptionsHandlers((PetscObject)ds,PetscOptionsObject);
720:   PetscOptionsEnd();
721:   return 0;
722: }

724: /*@C
725:    DSView - Prints the DS data structure.

727:    Collective on ds

729:    Input Parameters:
730: +  ds - the direct solver context
731: -  viewer - optional visualization context

733:    Note:
734:    The available visualization contexts include
735: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
736: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
737:          output where only the first processor opens
738:          the file.  All other processors send their
739:          data to the first processor to print.

741:    The user can open an alternative visualization context with
742:    PetscViewerASCIIOpen() - output to a specified file.

744:    Level: beginner

746: .seealso: DSViewMat()
747: @*/
748: PetscErrorCode DSView(DS ds,PetscViewer viewer)
749: {
750:   PetscBool         isascii;
751:   PetscViewerFormat format;
752:   PetscMPIInt       size;

755:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer);
758:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
759:   if (isascii) {
760:     PetscViewerGetFormat(viewer,&format);
761:     PetscObjectPrintClassNamePrefixType((PetscObject)ds,viewer);
762:     MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
763:     if (size>1) PetscViewerASCIIPrintf(viewer,"  parallel operation mode: %s\n",DSParallelTypes[ds->pmode]);
764:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
765:       PetscViewerASCIIPrintf(viewer,"  current state: %s\n",DSStateTypes[ds->state]);
766:       PetscViewerASCIIPrintf(viewer,"  dimensions: ld=%" PetscInt_FMT ", n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT,ds->ld,ds->n,ds->l,ds->k);
767:       if (ds->state==DS_STATE_TRUNCATED) PetscViewerASCIIPrintf(viewer,", t=%" PetscInt_FMT "\n",ds->t);
768:       else PetscViewerASCIIPrintf(viewer,"\n");
769:       PetscViewerASCIIPrintf(viewer,"  flags:%s%s%s\n",ds->compact?" compact":"",ds->extrarow?" extrarow":"",ds->refined?" refined":"");
770:     }
771:     PetscViewerASCIIPushTab(viewer);
772:     PetscTryTypeMethod(ds,view,viewer);
773:     PetscViewerASCIIPopTab(viewer);
774:   }
775:   return 0;
776: }

778: /*@C
779:    DSViewFromOptions - View from options

781:    Collective on DS

783:    Input Parameters:
784: +  ds   - the direct solver context
785: .  obj  - optional object
786: -  name - command line option

788:    Level: intermediate

790: .seealso: DSView(), DSCreate()
791: @*/
792: PetscErrorCode DSViewFromOptions(DS ds,PetscObject obj,const char name[])
793: {
795:   PetscObjectViewFromOptions((PetscObject)ds,obj,name);
796:   return 0;
797: }

799: /*@
800:    DSAllocate - Allocates memory for internal storage or matrices in DS.

802:    Logically Collective on ds

804:    Input Parameters:
805: +  ds - the direct solver context
806: -  ld - leading dimension (maximum allowed dimension for the matrices, including
807:         the extra row if present)

809:    Note:
810:    If the leading dimension is different from a previously set value, then
811:    all matrices are destroyed with DSReset().

813:    Level: intermediate

815: .seealso: DSGetLeadingDimension(), DSSetDimensions(), DSSetExtraRow(), DSReset()
816: @*/
817: PetscErrorCode DSAllocate(DS ds,PetscInt ld)
818: {
823:   if (ld!=ds->ld) {
824:     PetscInfo(ds,"Allocating memory with leading dimension=%" PetscInt_FMT "\n",ld);
825:     DSReset(ds);
826:     ds->ld = ld;
827:     PetscUseTypeMethod(ds,allocate,ld);
828:   }
829:   return 0;
830: }

832: /*@
833:    DSReset - Resets the DS context to the initial state.

835:    Collective on ds

837:    Input Parameter:
838: .  ds - the direct solver context

840:    Note:
841:    All data structures with size depending on the leading dimension
842:    of DSAllocate() are released.

844:    Level: advanced

846: .seealso: DSDestroy(), DSAllocate()
847: @*/
848: PetscErrorCode DSReset(DS ds)
849: {
850:   PetscInt       i;

853:   if (!ds) return 0;
854:   ds->state    = DS_STATE_RAW;
855:   ds->ld       = 0;
856:   ds->l        = 0;
857:   ds->n        = 0;
858:   ds->k        = 0;
859:   for (i=0;i<DS_NUM_MAT;i++) MatDestroy(&ds->omat[i]);
860:   PetscFree(ds->perm);
861:   return 0;
862: }

864: /*@C
865:    DSDestroy - Destroys DS context that was created with DSCreate().

867:    Collective on ds

869:    Input Parameter:
870: .  ds - the direct solver context

872:    Level: beginner

874: .seealso: DSCreate()
875: @*/
876: PetscErrorCode DSDestroy(DS *ds)
877: {
878:   if (!*ds) return 0;
880:   if (--((PetscObject)(*ds))->refct > 0) { *ds = NULL; return 0; }
881:   DSReset(*ds);
882:   PetscTryTypeMethod(*ds,destroy);
883:   PetscFree((*ds)->work);
884:   PetscFree((*ds)->rwork);
885:   PetscFree((*ds)->iwork);
886:   if (!(*ds)->scset) PetscFree((*ds)->sc);
887:   PetscHeaderDestroy(ds);
888:   return 0;
889: }

891: /*@C
892:    DSRegister - Adds a direct solver to the DS package.

894:    Not collective

896:    Input Parameters:
897: +  name - name of a new user-defined DS
898: -  function - routine to create context

900:    Notes:
901:    DSRegister() may be called multiple times to add several user-defined
902:    direct solvers.

904:    Level: advanced

906: .seealso: DSRegisterAll()
907: @*/
908: PetscErrorCode DSRegister(const char *name,PetscErrorCode (*function)(DS))
909: {
910:   DSInitializePackage();
911:   PetscFunctionListAdd(&DSList,name,function);
912:   return 0;
913: }

915: SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS);
916: SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS);
917: SLEPC_EXTERN PetscErrorCode DSCreate_GHEP(DS);
918: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS);
919: SLEPC_EXTERN PetscErrorCode DSCreate_GNHEP(DS);
920: SLEPC_EXTERN PetscErrorCode DSCreate_NHEPTS(DS);
921: SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS);
922: SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS);
923: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS);
924: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS);

926: /*@C
927:    DSRegisterAll - Registers all of the direct solvers in the DS package.

929:    Not Collective

931:    Level: advanced

933: .seealso: DSRegister()
934: @*/
935: PetscErrorCode DSRegisterAll(void)
936: {
937:   if (DSRegisterAllCalled) return 0;
938:   DSRegisterAllCalled = PETSC_TRUE;
939:   DSRegister(DSHEP,DSCreate_HEP);
940:   DSRegister(DSNHEP,DSCreate_NHEP);
941:   DSRegister(DSGHEP,DSCreate_GHEP);
942:   DSRegister(DSGHIEP,DSCreate_GHIEP);
943:   DSRegister(DSGNHEP,DSCreate_GNHEP);
944:   DSRegister(DSNHEPTS,DSCreate_NHEPTS);
945:   DSRegister(DSSVD,DSCreate_SVD);
946:   DSRegister(DSGSVD,DSCreate_GSVD);
947:   DSRegister(DSPEP,DSCreate_PEP);
948:   DSRegister(DSNEP,DSCreate_NEP);
949:   return 0;
950: }