Actual source code: dsbasic.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>      /*I "slepcds.h" I*/

 16: PetscFunctionList DSList = 0;
 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_",0};
 23: const char *DSParallelTypes[] = {"REDUNDANT","SYNCHRONIZED","DSParallelType","DS_PARALLEL_",0};
 24: const char *DSMatName[DS_NUM_MAT] = {"A","B","C","T","D","Q","Z","X","Y","U","VT","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: {

 40:   PetscFunctionListDestroy(&DSList);
 41:   DSPackageInitialized = PETSC_FALSE;
 42:   DSRegisterAllCalled  = PETSC_FALSE;
 43:   return(0);
 44: }

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

 51:   Level: developer

 53: .seealso: SlepcInitialize()
 54: @*/
 55: PetscErrorCode DSInitializePackage()
 56: {
 57:   char           logList[256];
 58:   PetscBool      opt,pkg;

 62:   if (DSPackageInitialized) return(0);
 63:   DSPackageInitialized = PETSC_TRUE;
 64:   /* Register Classes */
 65:   PetscClassIdRegister("Direct Solver",&DS_CLASSID);
 66:   /* Register Constructors */
 67:   DSRegisterAll();
 68:   /* Register Events */
 69:   PetscLogEventRegister("DSSolve",DS_CLASSID,&DS_Solve);
 70:   PetscLogEventRegister("DSVectors",DS_CLASSID,&DS_Vectors);
 71:   PetscLogEventRegister("DSSynchronize",DS_CLASSID,&DS_Synchronize);
 72:   PetscLogEventRegister("DSOther",DS_CLASSID,&DS_Other);
 73:   /* Process info exclusions */
 74:   PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
 75:   if (opt) {
 76:     PetscStrInList("ds",logList,',',&pkg);
 77:     if (pkg) { PetscInfoDeactivateClass(DS_CLASSID); }
 78:   }
 79:   /* Process summary exclusions */
 80:   PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
 81:   if (opt) {
 82:     PetscStrInList("ds",logList,',',&pkg);
 83:     if (pkg) { PetscLogEventDeactivateClass(DS_CLASSID); }
 84:   }
 85:   /* Register package finalizer */
 86:   PetscRegisterFinalize(DSFinalizePackage);
 87:   return(0);
 88: }

 90: /*@
 91:    DSCreate - Creates a DS context.

 93:    Collective on MPI_Comm

 95:    Input Parameter:
 96: .  comm - MPI communicator

 98:    Output Parameter:
 99: .  newds - location to put the DS context

101:    Level: beginner

103:    Note:
104:    DS objects are not intended for normal users but only for
105:    advanced user that for instance implement their own solvers.

107: .seealso: DSDestroy(), DS
108: @*/
109: PetscErrorCode DSCreate(MPI_Comm comm,DS *newds)
110: {
111:   DS             ds;
112:   PetscInt       i;

117:   *newds = 0;
118:   DSInitializePackage();
119:   SlepcHeaderCreate(ds,DS_CLASSID,"DS","Direct Solver (or Dense System)","DS",comm,DSDestroy,DSView);

121:   ds->state         = DS_STATE_RAW;
122:   ds->method        = 0;
123:   ds->compact       = PETSC_FALSE;
124:   ds->refined       = PETSC_FALSE;
125:   ds->extrarow      = PETSC_FALSE;
126:   ds->ld            = 0;
127:   ds->l             = 0;
128:   ds->n             = 0;
129:   ds->m             = 0;
130:   ds->k             = 0;
131:   ds->t             = 0;
132:   ds->bs            = 1;
133:   ds->sc            = NULL;
134:   ds->pmode         = DS_PARALLEL_REDUNDANT;

136:   for (i=0;i<DS_NUM_MAT;i++) {
137:     ds->mat[i]      = NULL;
138:     ds->rmat[i]     = NULL;
139:     ds->omat[i]     = NULL;
140:   }
141:   ds->perm          = NULL;
142:   ds->data          = NULL;
143:   ds->scset         = PETSC_FALSE;
144:   ds->work          = NULL;
145:   ds->rwork         = NULL;
146:   ds->iwork         = NULL;
147:   ds->lwork         = 0;
148:   ds->lrwork        = 0;
149:   ds->liwork        = 0;

151:   *newds = ds;
152:   return(0);
153: }

155: /*@C
156:    DSSetOptionsPrefix - Sets the prefix used for searching for all
157:    DS options in the database.

159:    Logically Collective on DS

161:    Input Parameters:
162: +  ds - the direct solver context
163: -  prefix - the prefix string to prepend to all DS option requests

165:    Notes:
166:    A hyphen (-) must NOT be given at the beginning of the prefix name.
167:    The first character of all runtime options is AUTOMATICALLY the
168:    hyphen.

170:    Level: advanced

172: .seealso: DSAppendOptionsPrefix()
173: @*/
174: PetscErrorCode DSSetOptionsPrefix(DS ds,const char *prefix)
175: {

180:   PetscObjectSetOptionsPrefix((PetscObject)ds,prefix);
181:   return(0);
182: }

184: /*@C
185:    DSAppendOptionsPrefix - Appends to the prefix used for searching for all
186:    DS options in the database.

188:    Logically Collective on DS

190:    Input Parameters:
191: +  ds - the direct solver context
192: -  prefix - the prefix string to prepend to all DS option requests

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

198:    Level: advanced

200: .seealso: DSSetOptionsPrefix()
201: @*/
202: PetscErrorCode DSAppendOptionsPrefix(DS ds,const char *prefix)
203: {

208:   PetscObjectAppendOptionsPrefix((PetscObject)ds,prefix);
209:   return(0);
210: }

212: /*@C
213:    DSGetOptionsPrefix - Gets the prefix used for searching for all
214:    DS options in the database.

216:    Not Collective

218:    Input Parameters:
219: .  ds - the direct solver context

221:    Output Parameters:
222: .  prefix - pointer to the prefix string used is returned

224:    Note:
225:    On the Fortran side, the user should pass in a string 'prefix' of
226:    sufficient length to hold the prefix.

228:    Level: advanced

230: .seealso: DSSetOptionsPrefix(), DSAppendOptionsPrefix()
231: @*/
232: PetscErrorCode DSGetOptionsPrefix(DS ds,const char *prefix[])
233: {

239:   PetscObjectGetOptionsPrefix((PetscObject)ds,prefix);
240:   return(0);
241: }

243: /*@C
244:    DSSetType - Selects the type for the DS object.

246:    Logically Collective on DS

248:    Input Parameter:
249: +  ds   - the direct solver context
250: -  type - a known type

252:    Level: intermediate

254: .seealso: DSGetType()
255: @*/
256: PetscErrorCode DSSetType(DS ds,DSType type)
257: {
258:   PetscErrorCode ierr,(*r)(DS);
259:   PetscBool      match;


265:   PetscObjectTypeCompare((PetscObject)ds,type,&match);
266:   if (match) return(0);

268:    PetscFunctionListFind(DSList,type,&r);
269:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DS type %s",type);

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

273:   PetscObjectChangeTypeName((PetscObject)ds,type);
274:   (*r)(ds);
275:   return(0);
276: }

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

281:    Not Collective

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

286:    Output Parameter:
287: .  name - name of the direct solver

289:    Level: intermediate

291: .seealso: DSSetType()
292: @*/
293: PetscErrorCode DSGetType(DS ds,DSType *type)
294: {
298:   *type = ((PetscObject)ds)->type_name;
299:   return(0);
300: }

302: /*@
303:    DSDuplicate - Creates a new direct solver object with the same options as
304:    an existing one.

306:    Collective on DS

308:    Input Parameter:
309: .  ds - direct solver context

311:    Output Parameter:
312: .  dsnew - location to put the new DS

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

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

320:    Level: intermediate

322: .seealso: DSCreate(), DSAllocate(), DSSetSlepcSC()
323: @*/
324: PetscErrorCode DSDuplicate(DS ds,DS *dsnew)
325: {

332:   DSCreate(PetscObjectComm((PetscObject)ds),dsnew);
333:   DSSetType(*dsnew,((PetscObject)ds)->type_name);
334:   (*dsnew)->method   = ds->method;
335:   (*dsnew)->compact  = ds->compact;
336:   (*dsnew)->refined  = ds->refined;
337:   (*dsnew)->extrarow = ds->extrarow;
338:   (*dsnew)->bs       = ds->bs;
339:   (*dsnew)->pmode    = ds->pmode;
340:   return(0);
341: }

343: /*@
344:    DSSetMethod - Selects the method to be used to solve the problem.

346:    Logically Collective on DS

348:    Input Parameter:
349: +  ds   - the direct solver context
350: -  meth - an index indentifying the method

352:    Options Database Key:
353: .  -ds_method <meth> - Sets the method

355:    Level: intermediate

357: .seealso: DSGetMethod()
358: @*/
359: PetscErrorCode DSSetMethod(DS ds,PetscInt meth)
360: {
364:   if (meth<0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");
365:   if (meth>DS_MAX_SOLVE) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too large value for the method");
366:   ds->method = meth;
367:   return(0);
368: }

370: /*@
371:    DSGetMethod - Gets the method currently used in the DS.

373:    Not Collective

375:    Input Parameter:
376: .  ds - the direct solver context

378:    Output Parameter:
379: .  meth - identifier of the method

381:    Level: intermediate

383: .seealso: DSSetMethod()
384: @*/
385: PetscErrorCode DSGetMethod(DS ds,PetscInt *meth)
386: {
390:   *meth = ds->method;
391:   return(0);
392: }

394: /*@
395:    DSSetParallel - Selects the mode of operation in parallel runs.

397:    Logically Collective on DS

399:    Input Parameter:
400: +  ds    - the direct solver context
401: -  pmode - the parallel mode

403:    Options Database Key:
404: .  -ds_parallel <mode> - Sets the parallel mode, either 'redundant' or 'synchronized'

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

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

417:    Level: advanced

419: .seealso: DSSynchronize(), DSGetParallel()
420: @*/
421: PetscErrorCode DSSetParallel(DS ds,DSParallelType pmode)
422: {
426:   ds->pmode = pmode;
427:   return(0);
428: }

430: /*@
431:    DSGetParallel - Gets the mode of operation in parallel runs.

433:    Not Collective

435:    Input Parameter:
436: .  ds - the direct solver context

438:    Output Parameter:
439: .  pmode - the parallel mode

441:    Level: advanced

443: .seealso: DSSetParallel()
444: @*/
445: PetscErrorCode DSGetParallel(DS ds,DSParallelType *pmode)
446: {
450:   *pmode = ds->pmode;
451:   return(0);
452: }

454: /*@
455:    DSSetCompact - Switch to compact storage of matrices.

457:    Logically Collective on DS

459:    Input Parameter:
460: +  ds   - the direct solver context
461: -  comp - a boolean flag

463:    Notes:
464:    Compact storage is used in some DS types such as DSHEP when the matrix
465:    is tridiagonal. This flag can be used to indicate whether the user
466:    provides the matrix entries via the compact form (the tridiagonal DS_MAT_T)
467:    or the non-compact one (DS_MAT_A).

469:    The default is PETSC_FALSE.

471:    Level: advanced

473: .seealso: DSGetCompact()
474: @*/
475: PetscErrorCode DSSetCompact(DS ds,PetscBool comp)
476: {
480:   ds->compact = comp;
481:   return(0);
482: }

484: /*@
485:    DSGetCompact - Gets the compact storage flag.

487:    Not Collective

489:    Input Parameter:
490: .  ds - the direct solver context

492:    Output Parameter:
493: .  comp - the flag

495:    Level: advanced

497: .seealso: DSSetCompact()
498: @*/
499: PetscErrorCode DSGetCompact(DS ds,PetscBool *comp)
500: {
504:   *comp = ds->compact;
505:   return(0);
506: }

508: /*@
509:    DSSetExtraRow - Sets a flag to indicate that the matrix has one extra
510:    row.

512:    Logically Collective on DS

514:    Input Parameter:
515: +  ds  - the direct solver context
516: -  ext - a boolean flag

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

524:    The default is PETSC_FALSE.

526:    Level: advanced

528: .seealso: DSSolve(), DSAllocate(), DSGetExtraRow()
529: @*/
530: PetscErrorCode DSSetExtraRow(DS ds,PetscBool ext)
531: {
535:   if (ds->n>0 && ds->n==ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot set extra row after setting n=ld");
536:   ds->extrarow = ext;
537:   return(0);
538: }

540: /*@
541:    DSGetExtraRow - Gets the extra row flag.

543:    Not Collective

545:    Input Parameter:
546: .  ds - the direct solver context

548:    Output Parameter:
549: .  ext - the flag

551:    Level: advanced

553: .seealso: DSSetExtraRow()
554: @*/
555: PetscErrorCode DSGetExtraRow(DS ds,PetscBool *ext)
556: {
560:   *ext = ds->extrarow;
561:   return(0);
562: }

564: /*@
565:    DSSetRefined - Sets a flag to indicate that refined vectors must be
566:    computed.

568:    Logically Collective on DS

570:    Input Parameter:
571: +  ds  - the direct solver context
572: -  ref - a boolean flag

574:    Notes:
575:    Normally the vectors returned in DS_MAT_X are eigenvectors of the
576:    projected matrix. With this flag activated, DSVectors() will return
577:    the right singular vector of the smallest singular value of matrix
578:    \tilde{A}-theta*I, where \tilde{A} is the extended (n+1)xn matrix
579:    and theta is the Ritz value. This is used in the refined Ritz
580:    approximation.

582:    The default is PETSC_FALSE.

584:    Level: advanced

586: .seealso: DSVectors(), DSGetRefined()
587: @*/
588: PetscErrorCode DSSetRefined(DS ds,PetscBool ref)
589: {
593:   ds->refined = ref;
594:   return(0);
595: }

597: /*@
598:    DSGetRefined - Gets the refined vectors flag.

600:    Not Collective

602:    Input Parameter:
603: .  ds - the direct solver context

605:    Output Parameter:
606: .  ref - the flag

608:    Level: advanced

610: .seealso: DSSetRefined()
611: @*/
612: PetscErrorCode DSGetRefined(DS ds,PetscBool *ref)
613: {
617:   *ref = ds->refined;
618:   return(0);
619: }

621: /*@
622:    DSSetBlockSize - Sets the block size.

624:    Logically Collective on DS

626:    Input Parameter:
627: +  ds - the direct solver context
628: -  bs - the block size

630:    Options Database Key:
631: .  -ds_block_size <bs> - Sets the block size

633:    Level: intermediate

635: .seealso: DSGetBlockSize()
636: @*/
637: PetscErrorCode DSSetBlockSize(DS ds,PetscInt bs)
638: {
642:   if (bs<1) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The block size must be at least one");
643:   ds->bs = bs;
644:   return(0);
645: }

647: /*@
648:    DSGetBlockSize - Gets the block size.

650:    Not Collective

652:    Input Parameter:
653: .  ds - the direct solver context

655:    Output Parameter:
656: .  bs - block size

658:    Level: intermediate

660: .seealso: DSSetBlockSize()
661: @*/
662: PetscErrorCode DSGetBlockSize(DS ds,PetscInt *bs)
663: {
667:   *bs = ds->bs;
668:   return(0);
669: }

671: /*@C
672:    DSSetSlepcSC - Sets the sorting criterion context.

674:    Not Collective

676:    Input Parameters:
677: +  ds - the direct solver context
678: -  sc - a pointer to the sorting criterion context

680:    Level: developer

682: .seealso: DSGetSlepcSC(), DSSort()
683: @*/
684: PetscErrorCode DSSetSlepcSC(DS ds,SlepcSC sc)
685: {

691:   if (ds->sc && !ds->scset) {
692:     PetscFree(ds->sc);
693:   }
694:   ds->sc    = sc;
695:   ds->scset = PETSC_TRUE;
696:   return(0);
697: }

699: /*@C
700:    DSGetSlepcSC - Gets the sorting criterion context.

702:    Not Collective

704:    Input Parameter:
705: .  ds - the direct solver context

707:    Output Parameters:
708: .  sc - a pointer to the sorting criterion context

710:    Level: developer

712: .seealso: DSSetSlepcSC(), DSSort()
713: @*/
714: PetscErrorCode DSGetSlepcSC(DS ds,SlepcSC *sc)
715: {

721:   if (!ds->sc) {
722:     PetscNewLog(ds,&ds->sc);
723:   }
724:   *sc = ds->sc;
725:   return(0);
726: }

728: /*@
729:    DSSetFromOptions - Sets DS options from the options database.

731:    Collective on DS

733:    Input Parameters:
734: .  ds - the direct solver context

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

739:    Level: beginner
740: @*/
741: PetscErrorCode DSSetFromOptions(DS ds)
742: {
744:   PetscInt       bs,meth;
745:   PetscBool      flag;
746:   DSParallelType pmode;

750:   DSRegisterAll();
751:   /* Set default type (we do not allow changing it with -ds_type) */
752:   if (!((PetscObject)ds)->type_name) {
753:     DSSetType(ds,DSNHEP);
754:   }
755:   PetscObjectOptionsBegin((PetscObject)ds);

757:     PetscOptionsInt("-ds_block_size","Block size for the dense system solver","DSSetBlockSize",ds->bs,&bs,&flag);
758:     if (flag) { DSSetBlockSize(ds,bs); }

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

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

766:     PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)ds);
767:   PetscOptionsEnd();
768:   return(0);
769: }

771: /*@C
772:    DSView - Prints the DS data structure.

774:    Collective on DS

776:    Input Parameters:
777: +  ds - the direct solver context
778: -  viewer - optional visualization context

780:    Note:
781:    The available visualization contexts include
782: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
783: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
784:          output where only the first processor opens
785:          the file.  All other processors send their
786:          data to the first processor to print.

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

791:    Level: beginner

793: .seealso: DSViewMat()
794: @*/
795: PetscErrorCode DSView(DS ds,PetscViewer viewer)
796: {
797:   PetscBool         isascii,issvd;
798:   PetscViewerFormat format;
799:   PetscErrorCode    ierr;
800:   PetscMPIInt       size;

804:   if (!viewer) {
805:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer);
806:   }
809:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
810:   if (isascii) {
811:     PetscViewerGetFormat(viewer,&format);
812:     PetscObjectPrintClassNamePrefixType((PetscObject)ds,viewer);
813:     MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
814:     if (size>1) {
815:       PetscViewerASCIIPrintf(viewer,"  parallel operation mode: %s\n",DSParallelTypes[ds->pmode]);
816:     }
817:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
818:       PetscViewerASCIIPrintf(viewer,"  current state: %s\n",DSStateTypes[ds->state]);
819:       PetscObjectTypeCompare((PetscObject)ds,DSSVD,&issvd);
820:       if (issvd) {
821:         PetscViewerASCIIPrintf(viewer,"  dimensions: ld=%D, n=%D, m=%D, l=%D, k=%D",ds->ld,ds->n,ds->m,ds->l,ds->k);
822:       } else {
823:         PetscViewerASCIIPrintf(viewer,"  dimensions: ld=%D, n=%D, l=%D, k=%D",ds->ld,ds->n,ds->l,ds->k);
824:       }
825:       if (ds->state==DS_STATE_TRUNCATED) {
826:         PetscViewerASCIIPrintf(viewer,", t=%D\n",ds->t);
827:       } else {
828:         PetscViewerASCIIPrintf(viewer,"\n");
829:       }
830:       PetscViewerASCIIPrintf(viewer,"  flags:%s%s%s\n",ds->compact?" compact":"",ds->extrarow?" extrarow":"",ds->refined?" refined":"");
831:     }
832:     if (ds->ops->view) {
833:       PetscViewerASCIIPushTab(viewer);
834:       (*ds->ops->view)(ds,viewer);
835:       PetscViewerASCIIPopTab(viewer);
836:     }
837:   }
838:   return(0);
839: }

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

844:    Logically Collective on DS

846:    Input Parameters:
847: +  ds - the direct solver context
848: -  ld - leading dimension (maximum allowed dimension for the matrices, including
849:         the extra row if present)

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

855:    Level: intermediate

857: .seealso: DSGetLeadingDimension(), DSSetDimensions(), DSSetExtraRow(), DSReset()
858: @*/
859: PetscErrorCode DSAllocate(DS ds,PetscInt ld)
860: {

867:   if (ld<1) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Leading dimension should be at least one");
868:   if (ld!=ds->ld) {
869:     DSReset(ds);
870:     ds->ld = ld;
871:     (*ds->ops->allocate)(ds,ld);
872:   }
873:   return(0);
874: }

876: /*@
877:    DSReset - Resets the DS context to the initial state.

879:    Collective on DS

881:    Input Parameter:
882: .  ds - the direct solver context

884:    Note:
885:    All data structures with size depending on the leading dimension
886:    of DSAllocate() are released.

888:    Level: advanced

890: .seealso: DSDestroy(), DSAllocate()
891: @*/
892: PetscErrorCode DSReset(DS ds)
893: {
894:   PetscInt       i;

899:   if (!ds) return(0);
900:   ds->state    = DS_STATE_RAW;
901:   ds->ld       = 0;
902:   ds->l        = 0;
903:   ds->n        = 0;
904:   ds->m        = 0;
905:   ds->k        = 0;
906:   for (i=0;i<DS_NUM_MAT;i++) {
907:     PetscFree(ds->mat[i]);
908:     PetscFree(ds->rmat[i]);
909:     MatDestroy(&ds->omat[i]);
910:   }
911:   PetscFree(ds->perm);
912:   return(0);
913: }

915: /*@
916:    DSDestroy - Destroys DS context that was created with DSCreate().

918:    Collective on DS

920:    Input Parameter:
921: .  ds - the direct solver context

923:    Level: beginner

925: .seealso: DSCreate()
926: @*/
927: PetscErrorCode DSDestroy(DS *ds)
928: {

932:   if (!*ds) return(0);
934:   if (--((PetscObject)(*ds))->refct > 0) { *ds = 0; return(0); }
935:   DSReset(*ds);
936:   if ((*ds)->ops->destroy) { (*(*ds)->ops->destroy)(*ds); }
937:   PetscFree((*ds)->work);
938:   PetscFree((*ds)->rwork);
939:   PetscFree((*ds)->iwork);
940:   if (!(*ds)->scset) { PetscFree((*ds)->sc); }
941:   PetscHeaderDestroy(ds);
942:   return(0);
943: }

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

948:    Not collective

950:    Input Parameters:
951: +  name - name of a new user-defined DS
952: -  routine_create - routine to create context

954:    Notes:
955:    DSRegister() may be called multiple times to add several user-defined
956:    direct solvers.

958:    Level: advanced

960: .seealso: DSRegisterAll()
961: @*/
962: PetscErrorCode DSRegister(const char *name,PetscErrorCode (*function)(DS))
963: {

967:   DSInitializePackage();
968:   PetscFunctionListAdd(&DSList,name,function);
969:   return(0);
970: }

972: SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS);
973: SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS);
974: SLEPC_EXTERN PetscErrorCode DSCreate_GHEP(DS);
975: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS);
976: SLEPC_EXTERN PetscErrorCode DSCreate_GNHEP(DS);
977: SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS);
978: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS);
979: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS);

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

984:    Not Collective

986:    Level: advanced
987: @*/
988: PetscErrorCode DSRegisterAll(void)
989: {

993:   if (DSRegisterAllCalled) return(0);
994:   DSRegisterAllCalled = PETSC_TRUE;
995:   DSRegister(DSHEP,DSCreate_HEP);
996:   DSRegister(DSNHEP,DSCreate_NHEP);
997:   DSRegister(DSGHEP,DSCreate_GHEP);
998:   DSRegister(DSGHIEP,DSCreate_GHIEP);
999:   DSRegister(DSGNHEP,DSCreate_GNHEP);
1000:   DSRegister(DSSVD,DSCreate_SVD);
1001:   DSRegister(DSPEP,DSCreate_PEP);
1002:   DSRegister(DSNEP,DSCreate_NEP);
1003:   return(0);
1004: }