Actual source code: dsops.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:    DS operations: DSSolve(), DSVectors(), etc
 12: */

 14: #include <slepc/private/dsimpl.h>

 16: /*@
 17:    DSGetLeadingDimension - Returns the leading dimension of the allocated
 18:    matrices.

 20:    Not Collective

 22:    Input Parameter:
 23: .  ds - the direct solver context

 25:    Output Parameter:
 26: .  ld - leading dimension (maximum allowed dimension for the matrices)

 28:    Level: advanced

 30: .seealso: DSAllocate(), DSSetDimensions()
 31: @*/
 32: PetscErrorCode DSGetLeadingDimension(DS ds,PetscInt *ld)
 33: {
 36:   *ld = ds->ld;
 37:   return 0;
 38: }

 40: /*@
 41:    DSSetState - Change the state of the DS object.

 43:    Logically Collective on ds

 45:    Input Parameters:
 46: +  ds    - the direct solver context
 47: -  state - the new state

 49:    Notes:
 50:    The state indicates that the dense system is in an initial state (raw),
 51:    in an intermediate state (such as tridiagonal, Hessenberg or
 52:    Hessenberg-triangular), in a condensed state (such as diagonal, Schur or
 53:    generalized Schur), or in a truncated state.

 55:    This function is normally used to return to the raw state when the
 56:    condensed structure is destroyed.

 58:    Level: advanced

 60: .seealso: DSGetState()
 61: @*/
 62: PetscErrorCode DSSetState(DS ds,DSStateType state)
 63: {
 66:   switch (state) {
 67:     case DS_STATE_RAW:
 68:     case DS_STATE_INTERMEDIATE:
 69:     case DS_STATE_CONDENSED:
 70:     case DS_STATE_TRUNCATED:
 71:       if (ds->state!=state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]);
 72:       ds->state = state;
 73:       break;
 74:     default:
 75:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
 76:   }
 77:   return 0;
 78: }

 80: /*@
 81:    DSGetState - Returns the current state.

 83:    Not Collective

 85:    Input Parameter:
 86: .  ds - the direct solver context

 88:    Output Parameter:
 89: .  state - current state

 91:    Level: advanced

 93: .seealso: DSSetState()
 94: @*/
 95: PetscErrorCode DSGetState(DS ds,DSStateType *state)
 96: {
 99:   *state = ds->state;
100:   return 0;
101: }

103: /*@
104:    DSSetDimensions - Resize the matrices in the DS object.

106:    Logically Collective on ds

108:    Input Parameters:
109: +  ds - the direct solver context
110: .  n  - the new size
111: .  l  - number of locked (inactive) leading columns
112: -  k  - intermediate dimension (e.g., position of arrow)

114:    Notes:
115:    The internal arrays are not reallocated.

117:    Some DS types have additional dimensions, e.g. the number of columns
118:    in DSSVD. For these, you should call a specific interface function.

120:    Level: intermediate

122: .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
123: @*/
124: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
125: {
126:   PetscInt       on,ol,ok;
127:   PetscBool      issvd;

130:   DSCheckAlloc(ds,1);
134:   on = ds->n; ol = ds->l; ok = ds->k;
135:   if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
136:     ds->n = ds->extrarow? ds->ld-1: ds->ld;
137:   } else {
139:     PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,"");  /* SVD and GSVD have extra column instead of extra row */
141:     ds->n = n;
142:   }
143:   ds->t = ds->n;   /* truncated length equal to the new dimension */
144:   if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
145:     ds->l = 0;
146:   } else {
148:     ds->l = l;
149:   }
150:   if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
151:     ds->k = ds->n/2;
152:   } else {
154:     ds->k = k;
155:   }
156:   if (on!=ds->n || ol!=ds->l || ok!=ds->k) PetscInfo(ds,"New dimensions are: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
157:   return 0;
158: }

160: /*@
161:    DSGetDimensions - Returns the current dimensions.

163:    Not Collective

165:    Input Parameter:
166: .  ds - the direct solver context

168:    Output Parameters:
169: +  n  - the current size
170: .  l  - number of locked (inactive) leading columns
171: .  k  - intermediate dimension (e.g., position of arrow)
172: -  t  - truncated length

174:    Note:
175:    The t parameter makes sense only if DSTruncate() has been called.
176:    Otherwise its value equals n.

178:    Some DS types have additional dimensions, e.g. the number of columns
179:    in DSSVD. For these, you should call a specific interface function.

181:    Level: intermediate

183: .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
184: @*/
185: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
186: {
188:   DSCheckAlloc(ds,1);
189:   if (n) *n = ds->n;
190:   if (l) *l = ds->l;
191:   if (k) *k = ds->k;
192:   if (t) *t = ds->t;
193:   return 0;
194: }

196: /*@
197:    DSTruncate - Truncates the system represented in the DS object.

199:    Logically Collective on ds

201:    Input Parameters:
202: +  ds   - the direct solver context
203: .  n    - the new size
204: -  trim - a flag to indicate if the factorization must be trimmed

206:    Note:
207:    The new size is set to n. Note that in some cases the new size could
208:    be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
209:    Schur form). In cases where the extra row is meaningful, the first
210:    n elements are kept as the extra row for the new system.

212:    If the flag trim is turned on, it resets the locked and intermediate
213:    dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
214:    It also cleans the extra row if being used.

216:    The typical usage of trim=true is to truncate the Schur decomposition
217:    at the end of a Krylov iteration. In this case, the state must be
218:    changed to RAW so that DSVectors() computes eigenvectors from scratch.

220:    Level: advanced

222: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
223: @*/
224: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
225: {
226:   DSStateType    old;

230:   DSCheckAlloc(ds,1);
234:   PetscLogEventBegin(DS_Other,ds,0,0,0);
235:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
236:   PetscUseTypeMethod(ds,truncate,n,trim);
237:   PetscFPTrapPop();
238:   PetscLogEventEnd(DS_Other,ds,0,0,0);
239:   PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n);
240:   old = ds->state;
241:   ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
242:   if (old!=ds->state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]);
243:   PetscObjectStateIncrease((PetscObject)ds);
244:   return 0;
245: }

247: /*@
248:    DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.

250:    Not Collective

252:    Input Parameters:
253: +  ds - the direct solver context
254: -  t  - the requested matrix

256:    Output Parameters:
257: +  n  - the number of rows
258: -  m  - the number of columns

260:    Note:
261:    This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().

263:    Level: developer

265: .seealso: DSSetDimensions(), DSGetMat()
266: @*/
267: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
268: {
269:   PetscInt       rows,cols;

273:   DSCheckValidMat(ds,t,2);
274:   if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
275:   else {
276:     if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
277:     else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
278:     if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
279:     else if (t==DS_MAT_D) cols = 1;
280:     else cols = ds->n;
281:   }
282:   if (m) *m = rows;
283:   if (n) *n = cols;
284:   return 0;
285: }

287: /*@
288:    DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.

290:    Not Collective

292:    Input Parameters:
293: +  ds - the direct solver context
294: -  t  - the requested matrix

296:    Output Parameter:
297: .  flg - the Hermitian flag

299:    Note:
300:    Does not check the matrix values directly. The flag is set according to the
301:    problem structure. For instance, in DSHEP matrix A is Hermitian.

303:    Level: developer

305: .seealso: DSGetMat()
306: @*/
307: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
308: {
311:   DSCheckValidMat(ds,t,2);
313:   *flg = PETSC_FALSE;
314:   PetscTryTypeMethod(ds,hermitian,t,flg);
315:   return 0;
316: }

318: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
319: {
320: #if !defined(PETSC_USE_COMPLEX)
321:   PetscScalar val;
322: #endif

324: #if !defined(PETSC_USE_COMPLEX)
325:   MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val);
326:   if (val != 0.0) {
327:     if (l+(*k)<n-1) (*k)++;
328:     else (*k)--;
329:   }
330: #endif
331:   return 0;
332: }

334: /*@
335:    DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
336:    to avoid breaking a 2x2 block.

338:    Not Collective

340:    Input Parameters:
341: +  ds - the direct solver context
342: .  l  - the size of the locked part (set to 0 to use ds->l)
343: .  n  - the total matrix size (set to 0 to use ds->n)
344: -  k  - the wanted truncation size

346:    Output Parameter:
347: .  k  - the possibly modified value of the truncation size

349:    Notes:
350:    This should be called before DSTruncate() to make sure that the truncation
351:    does not break a 2x2 block corresponding to a complex conjugate eigenvalue.

353:    The total size is n (either user-provided or ds->n if 0 is passed). The
354:    size where the truncation is intended is equal to l+k (where l can be
355:    equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
356:    at the l+k limit, the value of k is increased (or decreased) by 1.

358:    Level: advanced

360: .seealso: DSTruncate(), DSSetDimensions()
361: @*/
362: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
363: {
365:   DSCheckAlloc(ds,1);
369:   PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
370:   return 0;
371: }

373: /*@
374:    DSGetMat - Returns a sequential dense Mat object containing the requested
375:    matrix.

377:    Not Collective

379:    Input Parameters:
380: +  ds - the direct solver context
381: -  m  - the requested matrix

383:    Output Parameter:
384: .  A  - Mat object

386:    Notes:
387:    The returned Mat has sizes equal to the current DS dimensions (nxm),
388:    and contains the values that would be obtained with DSGetArray()
389:    (not DSGetArrayReal()). If the DS was truncated, then the number of rows
390:    is equal to the dimension prior to truncation, see DSTruncate().
391:    The communicator is always PETSC_COMM_SELF.

393:    It is implemented with MatDenseGetSubMatrix(), and when no longer needed
394:    the user must call DSRestoreMat() which will invoke MatDenseRestoreSubMatrix().

396:    For matrices DS_MAT_T and DS_MAT_D, this function will return a Mat object
397:    that cannot be used directly for computations, since it uses compact storage
398:    (three and one diagonals for T and D, respectively). In complex scalars, the
399:    internal array stores real values, so it is sufficient with 2 columns for T.

401:    Level: advanced

403: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
404: @*/
405: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
406: {
407:   PetscInt  rows,cols;
408:   PetscBool flg;

411:   DSCheckAlloc(ds,1);
412:   DSCheckValidMat(ds,m,2);

415:   DSMatGetSize(ds,m,&rows,&cols);
417:   MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A);

419:   /* set Hermitian flag */
420:   DSMatIsHermitian(ds,m,&flg);
421:   MatSetOption(*A,MAT_SYMMETRIC,flg);
422:   MatSetOption(*A,MAT_HERMITIAN,flg);
423:   return 0;
424: }

426: /*@
427:    DSRestoreMat - Restores the matrix after DSGetMat() was called.

429:    Not Collective

431:    Input Parameters:
432: +  ds - the direct solver context
433: .  m  - the requested matrix
434: -  A  - the fetched Mat object

436:    Note:
437:    A call to this function must match a previous call of DSGetMat().

439:    Level: advanced

441: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
442: @*/
443: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
444: {
446:   DSCheckAlloc(ds,1);
447:   DSCheckValidMat(ds,m,2);

450:   MatDenseRestoreSubMatrix(ds->omat[m],A);
451:   PetscObjectStateIncrease((PetscObject)ds);
452:   return 0;
453: }

455: /*@
456:    DSGetMatAndColumn - Returns a sequential dense Mat object containing the requested
457:    matrix and one of its columns as a Vec.

459:    Not Collective

461:    Input Parameters:
462: +  ds  - the direct solver context
463: .  m   - the requested matrix
464: -  col - the requested column

466:    Output Parameters:
467: +  A   - Mat object
468: -  v   - Vec object (the column)

470:    Notes:
471:    This calls DSGetMat() and then it extracts the selected column.
472:    The user must call DSRestoreMatAndColumn() to recover the original state.
473:    For matrices DS_MAT_T and DS_MAT_D, in complex scalars this function implies
474:    copying from real values stored internally to scalar values in the Vec.

476:    Level: advanced

478: .seealso: DSRestoreMatAndColumn(), DSGetMat()
479: @*/
480: PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
481: {
483:   DSCheckAlloc(ds,1);
484:   DSCheckValidMat(ds,m,2);

488:   DSGetMat(ds,m,A);
489:   if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
490:     const PetscScalar *as;
491:     PetscScalar       *vs;
492:     PetscReal         *ar;
493:     PetscInt          i,n,lda;
494:     MatCreateVecs(*A,NULL,v);
495:     VecGetSize(*v,&n);
496:     MatDenseGetLDA(*A,&lda);
497:     MatDenseGetArrayRead(*A,&as);
498:     VecGetArrayWrite(*v,&vs);
499:     ar = (PetscReal*)as;
500:     for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
501:     VecRestoreArrayWrite(*v,&vs);
502:     MatDenseRestoreArrayRead(*A,&as);
503:   } else MatDenseGetColumnVec(*A,col,v);
504:   return 0;
505: }

507: /*@
508:    DSRestoreMatAndColumn - Restores the matrix and vector after DSGetMatAndColumn()
509:    was called.

511:    Not Collective

513:    Input Parameters:
514: +  ds  - the direct solver context
515: .  m   - the requested matrix
516: .  col - the requested column
517: .  A   - the fetched Mat object
518: -  v   - the fetched Vec object

520:    Note:
521:    A call to this function must match a previous call of DSGetMatAndColumn().

523:    Level: advanced

525: .seealso: DSGetMatAndColumn(), DSRestoreMat()
526: @*/
527: PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
528: {
530:   DSCheckAlloc(ds,1);
531:   DSCheckValidMat(ds,m,2);

535:   if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
536:     const PetscScalar *vs;
537:     PetscScalar       *as;
538:     PetscReal         *ar;
539:     PetscInt          i,n,lda;
540:     VecGetSize(*v,&n);
541:     MatDenseGetLDA(*A,&lda);
542:     MatDenseGetArray(*A,&as);
543:     VecGetArrayRead(*v,&vs);
544:     ar = (PetscReal*)as;
545:     for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
546:     VecRestoreArrayRead(*v,&vs);
547:     MatDenseRestoreArray(*A,&as);
548:     VecDestroy(v);
549:   } else MatDenseRestoreColumnVec(*A,col,v);
550:   DSRestoreMat(ds,m,A);
551:   return 0;
552: }

554: /*@C
555:    DSGetArray - Returns a pointer to the internal array of one of the
556:    matrices. You MUST call DSRestoreArray() when you no longer
557:    need to access the array.

559:    Not Collective

561:    Input Parameters:
562: +  ds - the direct solver context
563: -  m  - the requested matrix

565:    Output Parameter:
566: .  a  - pointer to the values

568:    Note:
569:    To get read-only access, use DSGetMat() followed by MatDenseGetArrayRead().

571:    Level: advanced

573: .seealso: DSRestoreArray(), DSGetArrayReal()
574: @*/
575: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
576: {
578:   DSCheckAlloc(ds,1);
579:   DSCheckValidMat(ds,m,2);
581:   MatDenseGetArray(ds->omat[m],a);
582:   return 0;
583: }

585: /*@C
586:    DSRestoreArray - Restores the matrix after DSGetArray() was called.

588:    Not Collective

590:    Input Parameters:
591: +  ds - the direct solver context
592: .  m  - the requested matrix
593: -  a  - pointer to the values

595:    Level: advanced

597: .seealso: DSGetArray(), DSGetArrayReal()
598: @*/
599: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
600: {
602:   DSCheckAlloc(ds,1);
603:   DSCheckValidMat(ds,m,2);
605:   MatDenseRestoreArray(ds->omat[m],a);
606:   PetscObjectStateIncrease((PetscObject)ds);
607:   return 0;
608: }

610: /*@C
611:    DSGetArrayReal - Returns a real pointer to the internal array of T or D.
612:    You MUST call DSRestoreArrayReal() when you no longer need to access the array.

614:    Not Collective

616:    Input Parameters:
617: +  ds - the direct solver context
618: -  m  - the requested matrix

620:    Output Parameter:
621: .  a  - pointer to the values

623:    Note:
624:    This function can be used only for DS_MAT_T and DS_MAT_D. These matrices always
625:    store real values, even in complex scalars, that is why the returned pointer is
626:    PetscReal.

628:    Level: advanced

630: .seealso: DSRestoreArrayReal(), DSGetArray()
631: @*/
632: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
633: {
634: #if defined(PETSC_USE_COMPLEX)
635:   PetscScalar *as;
636: #endif

639:   DSCheckAlloc(ds,1);
640:   DSCheckValidMatReal(ds,m,2);
642: #if defined(PETSC_USE_COMPLEX)
643:   MatDenseGetArray(ds->omat[m],&as);
644:   *a = (PetscReal*)as;
645: #else
646:   MatDenseGetArray(ds->omat[m],a);
647: #endif
648:   return 0;
649: }

651: /*@C
652:    DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.

654:    Not Collective

656:    Input Parameters:
657: +  ds - the direct solver context
658: .  m  - the requested matrix
659: -  a  - pointer to the values

661:    Level: advanced

663: .seealso: DSGetArrayReal(), DSGetArray()
664: @*/
665: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
666: {
667: #if defined(PETSC_USE_COMPLEX)
668:   PetscScalar *as;
669: #endif

672:   DSCheckAlloc(ds,1);
673:   DSCheckValidMatReal(ds,m,2);
675: #if defined(PETSC_USE_COMPLEX)
676:   MatDenseRestoreArray(ds->omat[m],&as);
677:   *a = NULL;
678: #else
679:   MatDenseRestoreArray(ds->omat[m],a);
680: #endif
681:   PetscObjectStateIncrease((PetscObject)ds);
682:   return 0;
683: }

685: /*@
686:    DSSolve - Solves the problem.

688:    Logically Collective on ds

690:    Input Parameters:
691: +  ds   - the direct solver context
692: .  eigr - array to store the computed eigenvalues (real part)
693: -  eigi - array to store the computed eigenvalues (imaginary part)

695:    Note:
696:    This call brings the dense system to condensed form. No ordering
697:    of the eigenvalues is enforced (for this, call DSSort() afterwards).

699:    Level: intermediate

701: .seealso: DSSort(), DSStateType
702: @*/
703: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
704: {
707:   DSCheckAlloc(ds,1);
709:   if (ds->state>=DS_STATE_CONDENSED) return 0;
711:   PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
712:   PetscLogEventBegin(DS_Solve,ds,0,0,0);
713:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
714:   PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
715:   PetscFPTrapPop();
716:   PetscLogEventEnd(DS_Solve,ds,0,0,0);
717:   PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
718:   ds->state = DS_STATE_CONDENSED;
719:   PetscObjectStateIncrease((PetscObject)ds);
720:   return 0;
721: }

723: /*@C
724:    DSSort - Sorts the result of DSSolve() according to a given sorting
725:    criterion.

727:    Logically Collective on ds

729:    Input Parameters:
730: +  ds   - the direct solver context
731: .  rr   - (optional) array containing auxiliary values (real part)
732: -  ri   - (optional) array containing auxiliary values (imaginary part)

734:    Input/Output Parameters:
735: +  eigr - array containing the computed eigenvalues (real part)
736: .  eigi - array containing the computed eigenvalues (imaginary part)
737: -  k    - (optional) number of elements in the leading group

739:    Notes:
740:    This routine sorts the arrays provided in eigr and eigi, and also
741:    sorts the dense system stored inside ds (assumed to be in condensed form).
742:    The sorting criterion is specified with DSSetSlepcSC().

744:    If arrays rr and ri are provided, then a (partial) reordering based on these
745:    values rather than on the eigenvalues is performed. In symmetric problems
746:    a total order is obtained (parameter k is ignored), but otherwise the result
747:    is sorted only partially. In this latter case, it is only guaranteed that
748:    all the first k elements satisfy the comparison with any of the last n-k
749:    elements. The output value of parameter k is the final number of elements in
750:    the first set.

752:    Level: intermediate

754: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
755: @*/
756: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
757: {
758:   PetscInt       i;

762:   DSCheckSolved(ds,1);

769:   for (i=0;i<ds->n;i++) ds->perm[i] = i;   /* initialize to trivial permutation */
770:   PetscLogEventBegin(DS_Other,ds,0,0,0);
771:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
772:   PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
773:   PetscFPTrapPop();
774:   PetscLogEventEnd(DS_Other,ds,0,0,0);
775:   PetscObjectStateIncrease((PetscObject)ds);
776:   PetscInfo(ds,"Finished sorting\n");
777:   return 0;
778: }

780: /*@C
781:    DSSortWithPermutation - Reorders the result of DSSolve() according to a given
782:    permutation.

784:    Logically Collective on ds

786:    Input Parameters:
787: +  ds   - the direct solver context
788: -  perm - permutation that indicates the new ordering

790:    Input/Output Parameters:
791: +  eigr - array with the reordered eigenvalues (real part)
792: -  eigi - array with the reordered eigenvalues (imaginary part)

794:    Notes:
795:    This routine reorders the arrays provided in eigr and eigi, and also the dense
796:    system stored inside ds (assumed to be in condensed form). There is no sorting
797:    criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.

799:    Level: advanced

801: .seealso: DSSolve(), DSSort()
802: @*/
803: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)
804: {
807:   DSCheckSolved(ds,1);

812:   PetscLogEventBegin(DS_Other,ds,0,0,0);
813:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
814:   PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
815:   PetscFPTrapPop();
816:   PetscLogEventEnd(DS_Other,ds,0,0,0);
817:   PetscObjectStateIncrease((PetscObject)ds);
818:   PetscInfo(ds,"Finished sorting\n");
819:   return 0;
820: }

822: /*@
823:    DSSynchronize - Make sure that all processes have the same data, performing
824:    communication if necessary.

826:    Collective on ds

828:    Input Parameter:
829: .  ds   - the direct solver context

831:    Input/Output Parameters:
832: +  eigr - (optional) array with the computed eigenvalues (real part)
833: -  eigi - (optional) array with the computed eigenvalues (imaginary part)

835:    Notes:
836:    When the DS has been created with a communicator with more than one process,
837:    the internal data, especially the computed matrices, may diverge in the
838:    different processes. This happens when using multithreaded BLAS and may
839:    cause numerical issues in some ill-conditioned problems. This function
840:    performs the necessary communication among the processes so that the
841:    internal data is exactly equal in all of them.

843:    Depending on the parallel mode as set with DSSetParallel(), this function
844:    will either do nothing or synchronize the matrices computed by DSSolve()
845:    and DSSort(). The arguments eigr and eigi are typically those used in the
846:    calls to DSSolve() and DSSort().

848:    Level: developer

850: .seealso: DSSetParallel(), DSSolve(), DSSort()
851: @*/
852: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
853: {
854:   PetscMPIInt    size;

858:   DSCheckAlloc(ds,1);
859:   MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
860:   if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
861:     PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
862:     PetscUseTypeMethod(ds,synchronize,eigr,eigi);
863:     PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
864:     PetscObjectStateIncrease((PetscObject)ds);
865:     PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
866:   }
867:   return 0;
868: }

870: /*@C
871:    DSVectors - Compute vectors associated to the dense system such
872:    as eigenvectors.

874:    Logically Collective on ds

876:    Input Parameters:
877: +  ds  - the direct solver context
878: -  mat - the matrix, used to indicate which vectors are required

880:    Input/Output Parameter:
881: .  j   - (optional) index of vector to be computed

883:    Output Parameter:
884: .  rnorm - (optional) computed residual norm

886:    Notes:
887:    Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
888:    compute right or left eigenvectors, or left or right singular vectors,
889:    respectively.

891:    If NULL is passed in argument j then all vectors are computed,
892:    otherwise j indicates which vector must be computed. In real non-symmetric
893:    problems, on exit the index j will be incremented when a complex conjugate
894:    pair is found.

896:    This function can be invoked after the dense problem has been solved,
897:    to get the residual norm estimate of the associated Ritz pair. In that
898:    case, the relevant information is returned in rnorm.

900:    For computing eigenvectors, LAPACK's _trevc is used so the matrix must
901:    be in (quasi-)triangular form, or call DSSolve() first.

903:    Level: intermediate

905: .seealso: DSSolve()
906: @*/
907: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
908: {
911:   DSCheckAlloc(ds,1);
915:   if (!ds->omat[mat]) DSAllocateMat_Private(ds,mat);
916:   if (!j) PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]);
917:   PetscLogEventBegin(DS_Vectors,ds,0,0,0);
918:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
919:   PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
920:   PetscFPTrapPop();
921:   PetscLogEventEnd(DS_Vectors,ds,0,0,0);
922:   PetscObjectStateIncrease((PetscObject)ds);
923:   return 0;
924: }

926: /*@
927:    DSUpdateExtraRow - Performs all necessary operations so that the extra
928:    row gets up-to-date after a call to DSSolve().

930:    Not Collective

932:    Input Parameters:
933: .  ds - the direct solver context

935:    Level: advanced

937: .seealso: DSSolve(), DSSetExtraRow()
938: @*/
939: PetscErrorCode DSUpdateExtraRow(DS ds)
940: {
943:   DSCheckAlloc(ds,1);
945:   PetscInfo(ds,"Updating extra row\n");
946:   PetscLogEventBegin(DS_Other,ds,0,0,0);
947:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
948:   PetscUseTypeMethod(ds,update);
949:   PetscFPTrapPop();
950:   PetscLogEventEnd(DS_Other,ds,0,0,0);
951:   return 0;
952: }

954: /*@
955:    DSCond - Compute the condition number.

957:    Not Collective

959:    Input Parameters:
960: +  ds - the direct solver context
961: -  cond - the computed condition number

963:    Notes:
964:    In standard eigenvalue problems, returns the inf-norm condition number of the first
965:    matrix, computed as cond(A) = norm(A)*norm(inv(A)).

967:    In GSVD problems, returns the maximum of cond(A) and cond(B), where cond(.) is
968:    computed as the ratio of the largest and smallest singular values.

970:    Does not take into account the extra row.

972:    Level: advanced

974: .seealso: DSSolve(), DSSetExtraRow()
975: @*/
976: PetscErrorCode DSCond(DS ds,PetscReal *cond)
977: {
980:   DSCheckAlloc(ds,1);
982:   PetscLogEventBegin(DS_Other,ds,0,0,0);
983:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
984:   PetscUseTypeMethod(ds,cond,cond);
985:   PetscFPTrapPop();
986:   PetscLogEventEnd(DS_Other,ds,0,0,0);
987:   PetscInfo(ds,"Computed condition number = %g\n",(double)*cond);
988:   return 0;
989: }

991: /*@C
992:    DSTranslateHarmonic - Computes a translation of the dense system.

994:    Logically Collective on ds

996:    Input Parameters:
997: +  ds      - the direct solver context
998: .  tau     - the translation amount
999: .  beta    - last component of vector b
1000: -  recover - boolean flag to indicate whether to recover or not

1002:    Output Parameters:
1003: +  g       - the computed vector (optional)
1004: -  gamma   - scale factor (optional)

1006:    Notes:
1007:    This function is intended for use in the context of Krylov methods only.
1008:    It computes a translation of a Krylov decomposition in order to extract
1009:    eigenpair approximations by harmonic Rayleigh-Ritz.
1010:    The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
1011:    vector b is assumed to be beta*e_n^T.

1013:    The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
1014:    the factor by which the residual of the Krylov decomposition is scaled.

1016:    If the recover flag is activated, the computed translation undoes the
1017:    translation done previously. In that case, parameter tau is ignored.

1019:    Level: developer

1021: .seealso: DSTranslateRKS()
1022: @*/
1023: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
1024: {
1027:   DSCheckAlloc(ds,1);
1028:   if (recover) PetscInfo(ds,"Undoing the translation\n");
1029:   else PetscInfo(ds,"Computing the translation\n");
1030:   PetscLogEventBegin(DS_Other,ds,0,0,0);
1031:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1032:   PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
1033:   PetscFPTrapPop();
1034:   PetscLogEventEnd(DS_Other,ds,0,0,0);
1035:   ds->state = DS_STATE_RAW;
1036:   PetscObjectStateIncrease((PetscObject)ds);
1037:   return 0;
1038: }

1040: /*@
1041:    DSTranslateRKS - Computes a modification of the dense system corresponding
1042:    to an update of the shift in a rational Krylov method.

1044:    Logically Collective on ds

1046:    Input Parameters:
1047: +  ds    - the direct solver context
1048: -  alpha - the translation amount

1050:    Notes:
1051:    This function is intended for use in the context of Krylov methods only.
1052:    It takes the leading (k+1,k) submatrix of A, containing the truncated
1053:    Rayleigh quotient of a Krylov-Schur relation computed from a shift
1054:    sigma1 and transforms it to obtain a Krylov relation as if computed
1055:    from a different shift sigma2. The new matrix is computed as
1056:    1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
1057:    alpha = sigma1-sigma2.

1059:    Matrix Q is placed in DS_MAT_Q so that it can be used to update the
1060:    Krylov basis.

1062:    Level: developer

1064: .seealso: DSTranslateHarmonic()
1065: @*/
1066: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
1067: {
1070:   DSCheckAlloc(ds,1);
1071:   PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha));
1072:   PetscLogEventBegin(DS_Other,ds,0,0,0);
1073:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1074:   PetscUseTypeMethod(ds,transrks,alpha);
1075:   PetscFPTrapPop();
1076:   PetscLogEventEnd(DS_Other,ds,0,0,0);
1077:   ds->state   = DS_STATE_RAW;
1078:   ds->compact = PETSC_FALSE;
1079:   PetscObjectStateIncrease((PetscObject)ds);
1080:   return 0;
1081: }