Actual source code: dsops.c

slepc-main 2024-11-09
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: {
 34:   PetscFunctionBegin;
 36:   PetscAssertPointer(ld,2);
 37:   *ld = ds->ld;
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

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

 44:    Logically Collective

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

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

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

 59:    Level: advanced

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

 82: /*@
 83:    DSGetState - Returns the current state.

 85:    Not Collective

 87:    Input Parameter:
 88: .  ds - the direct solver context

 90:    Output Parameter:
 91: .  state - current state

 93:    Level: advanced

 95: .seealso: DSSetState()
 96: @*/
 97: PetscErrorCode DSGetState(DS ds,DSStateType *state)
 98: {
 99:   PetscFunctionBegin;
101:   PetscAssertPointer(state,2);
102:   *state = ds->state;
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

106: /*@
107:    DSSetDimensions - Resize the matrices in the DS object.

109:    Logically Collective

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

117:    Notes:
118:    The internal arrays are not reallocated.

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

123:    Use PETSC_CURRENT to leave any of the values unchanged. Use PETSC_DETERMINE
124:    to set n to the leading dimension, l to the minimum value (0), and k to n/2.

126:    Level: intermediate

128: .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
129: @*/
130: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
131: {
132:   PetscInt       on,ol,ok;
133:   PetscBool      issvd;

135:   PetscFunctionBegin;
137:   DSCheckAlloc(ds,1);
141:   on = ds->n; ol = ds->l; ok = ds->k;
142:   if (n == PETSC_DETERMINE) {
143:     ds->n = ds->extrarow? ds->ld-1: ds->ld;
144:   } else if (n != PETSC_CURRENT) {
145:     PetscCheck(n>=0 && n<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
146:     PetscCall(PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,""));  /* SVD and GSVD have extra column instead of extra row */
147:     PetscCheck(!ds->extrarow || n<ds->ld || issvd,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"A value of n equal to ld leaves no room for extra row");
148:     ds->n = n;
149:   }
150:   ds->t = ds->n;   /* truncated length equal to the new dimension */
151:   if (l == PETSC_DETERMINE) {
152:     ds->l = 0;
153:   } else if (l != PETSC_CURRENT) {
154:     PetscCheck(l>=0 && l<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
155:     ds->l = l;
156:   }
157:   if (k == PETSC_DETERMINE) {
158:     ds->k = ds->n/2;
159:   } else if (k != PETSC_CURRENT) {
160:     PetscCheck(k>=0 || k<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
161:     ds->k = k;
162:   }
163:   if (on!=ds->n || ol!=ds->l || ok!=ds->k) PetscCall(PetscInfo(ds,"New dimensions are: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

167: /*@
168:    DSGetDimensions - Returns the current dimensions.

170:    Not Collective

172:    Input Parameter:
173: .  ds - the direct solver context

175:    Output Parameters:
176: +  n  - the current size
177: .  l  - number of locked (inactive) leading columns
178: .  k  - intermediate dimension (e.g., position of arrow)
179: -  t  - truncated length

181:    Note:
182:    The t parameter makes sense only if DSTruncate() has been called.
183:    Otherwise its value equals n.

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

188:    Level: intermediate

190: .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
191: @*/
192: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
193: {
194:   PetscFunctionBegin;
196:   DSCheckAlloc(ds,1);
197:   if (n) *n = ds->n;
198:   if (l) *l = ds->l;
199:   if (k) *k = ds->k;
200:   if (t) *t = ds->t;
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: /*@
205:    DSTruncate - Truncates the system represented in the DS object.

207:    Logically Collective

209:    Input Parameters:
210: +  ds   - the direct solver context
211: .  n    - the new size
212: -  trim - a flag to indicate if the factorization must be trimmed

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

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

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

228:    Level: advanced

230: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
231: @*/
232: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
233: {
234:   DSStateType    old;

236:   PetscFunctionBegin;
239:   DSCheckAlloc(ds,1);
242:   PetscCheck(n>=ds->l && n<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n (%" PetscInt_FMT "). Must be between l (%" PetscInt_FMT ") and n (%" PetscInt_FMT ")",n,ds->l,ds->n);
243:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
244:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
245:   PetscUseTypeMethod(ds,truncate,n,trim);
246:   PetscCall(PetscFPTrapPop());
247:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
248:   PetscCall(PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n));
249:   old = ds->state;
250:   ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
251:   if (old!=ds->state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]));
252:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

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

259:    Not Collective

261:    Input Parameters:
262: +  ds - the direct solver context
263: -  t  - the requested matrix

265:    Output Parameters:
266: +  n  - the number of rows
267: -  m  - the number of columns

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

272:    Level: developer

274: .seealso: DSSetDimensions(), DSGetMat()
275: @*/
276: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
277: {
278:   PetscInt       rows,cols;

280:   PetscFunctionBegin;
283:   DSCheckValidMat(ds,t,2);
284:   if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
285:   else {
286:     if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
287:     else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
288:     if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
289:     else if (t==DS_MAT_D) cols = 1;
290:     else cols = ds->n;
291:   }
292:   if (m) *m = rows;
293:   if (n) *n = cols;
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

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

300:    Not Collective

302:    Input Parameters:
303: +  ds - the direct solver context
304: -  t  - the requested matrix

306:    Output Parameter:
307: .  flg - the Hermitian flag

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

313:    Level: developer

315: .seealso: DSGetMat()
316: @*/
317: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
318: {
319:   PetscFunctionBegin;
322:   DSCheckValidMat(ds,t,2);
323:   PetscAssertPointer(flg,3);
324:   *flg = PETSC_FALSE;
325:   PetscTryTypeMethod(ds,hermitian,t,flg);
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

329: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
330: {
331: #if !defined(PETSC_USE_COMPLEX)
332:   PetscScalar val;
333: #endif

335:   PetscFunctionBegin;
336: #if !defined(PETSC_USE_COMPLEX)
337:   PetscCall(MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val));
338:   if (val != 0.0) {
339:     if (l+(*k)<n-1) (*k)++;
340:     else (*k)--;
341:   }
342: #endif
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }

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

350:    Not Collective

352:    Input Parameters:
353: +  ds - the direct solver context
354: .  l  - the size of the locked part (set to 0 to use ds->l)
355: -  n  - the total matrix size (set to 0 to use ds->n)

357:    Output Parameter:
358: .  k  - the wanted truncation size (possibly modified)

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

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

369:    Level: advanced

371: .seealso: DSTruncate(), DSSetDimensions()
372: @*/
373: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
374: {
375:   PetscFunctionBegin;
377:   DSCheckAlloc(ds,1);
380:   PetscAssertPointer(k,4);
381:   PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /*@
386:    DSGetMat - Returns a sequential dense Mat object containing the requested
387:    matrix.

389:    Not Collective

391:    Input Parameters:
392: +  ds - the direct solver context
393: -  m  - the requested matrix

395:    Output Parameter:
396: .  A  - Mat object

398:    Notes:
399:    The returned Mat has sizes equal to the current DS dimensions (nxm),
400:    and contains the values that would be obtained with DSGetArray()
401:    (not DSGetArrayReal()). If the DS was truncated, then the number of rows
402:    is equal to the dimension prior to truncation, see DSTruncate().
403:    The communicator is always PETSC_COMM_SELF.

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

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

413:    Level: advanced

415: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
416: @*/
417: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
418: {
419:   PetscInt  rows=0,cols=0;
420:   PetscBool flg;

422:   PetscFunctionBegin;
424:   DSCheckAlloc(ds,1);
425:   DSCheckValidMat(ds,m,2);
426:   PetscAssertPointer(A,3);

428:   PetscCall(DSMatGetSize(ds,m,&rows,&cols));
429:   PetscCheck(rows && cols,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must call DSSetDimensions() first");
430:   PetscCall(MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A));

432:   /* set Hermitian flag */
433:   PetscCall(DSMatIsHermitian(ds,m,&flg));
434:   PetscCall(MatSetOption(*A,MAT_SYMMETRIC,flg));
435:   PetscCall(MatSetOption(*A,MAT_HERMITIAN,flg));
436:   PetscFunctionReturn(PETSC_SUCCESS);
437: }

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

442:    Not Collective

444:    Input Parameters:
445: +  ds - the direct solver context
446: .  m  - the requested matrix
447: -  A  - the fetched Mat object

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

452:    Level: advanced

454: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
455: @*/
456: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
457: {
458:   PetscFunctionBegin;
460:   DSCheckAlloc(ds,1);
461:   DSCheckValidMat(ds,m,2);
462:   PetscAssertPointer(A,3);

464:   PetscCall(MatDenseRestoreSubMatrix(ds->omat[m],A));
465:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

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

473:    Not Collective

475:    Input Parameters:
476: +  ds  - the direct solver context
477: .  m   - the requested matrix
478: -  col - the requested column

480:    Output Parameters:
481: +  A   - Mat object
482: -  v   - Vec object (the column)

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

490:    Level: advanced

492: .seealso: DSRestoreMatAndColumn(), DSGetMat()
493: @*/
494: PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
495: {
496:   PetscFunctionBegin;
498:   DSCheckAlloc(ds,1);
499:   DSCheckValidMat(ds,m,2);
500:   PetscAssertPointer(A,4);
501:   PetscAssertPointer(v,5);

503:   PetscCall(DSGetMat(ds,m,A));
504:   if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
505:     const PetscScalar *as;
506:     PetscScalar       *vs;
507:     PetscReal         *ar;
508:     PetscInt          i,n,lda;
509:     PetscCall(MatCreateVecs(*A,NULL,v));
510:     PetscCall(VecGetSize(*v,&n));
511:     PetscCall(MatDenseGetLDA(*A,&lda));
512:     PetscCall(MatDenseGetArrayRead(*A,&as));
513:     PetscCall(VecGetArrayWrite(*v,&vs));
514:     ar = (PetscReal*)as;
515:     for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
516:     PetscCall(VecRestoreArrayWrite(*v,&vs));
517:     PetscCall(MatDenseRestoreArrayRead(*A,&as));
518:   } else PetscCall(MatDenseGetColumnVec(*A,col,v));
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

522: /*@
523:    DSRestoreMatAndColumn - Restores the matrix and vector after DSGetMatAndColumn()
524:    was called.

526:    Not Collective

528:    Input Parameters:
529: +  ds  - the direct solver context
530: .  m   - the requested matrix
531: .  col - the requested column
532: .  A   - the fetched Mat object
533: -  v   - the fetched Vec object

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

538:    Level: advanced

540: .seealso: DSGetMatAndColumn(), DSRestoreMat()
541: @*/
542: PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
543: {
544:   PetscFunctionBegin;
546:   DSCheckAlloc(ds,1);
547:   DSCheckValidMat(ds,m,2);
548:   PetscAssertPointer(A,4);
549:   PetscAssertPointer(v,5);

551:   if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
552:     const PetscScalar *vs;
553:     PetscScalar       *as;
554:     PetscReal         *ar;
555:     PetscInt          i,n,lda;
556:     PetscCall(VecGetSize(*v,&n));
557:     PetscCall(MatDenseGetLDA(*A,&lda));
558:     PetscCall(MatDenseGetArray(*A,&as));
559:     PetscCall(VecGetArrayRead(*v,&vs));
560:     ar = (PetscReal*)as;
561:     for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
562:     PetscCall(VecRestoreArrayRead(*v,&vs));
563:     PetscCall(MatDenseRestoreArray(*A,&as));
564:     PetscCall(VecDestroy(v));
565:   } else PetscCall(MatDenseRestoreColumnVec(*A,col,v));
566:   PetscCall(DSRestoreMat(ds,m,A));
567:   PetscFunctionReturn(PETSC_SUCCESS);
568: }

570: /*@C
571:    DSGetArray - Returns a pointer to the internal array of one of the
572:    matrices. You MUST call DSRestoreArray() when you no longer
573:    need to access the array.

575:    Not Collective

577:    Input Parameters:
578: +  ds - the direct solver context
579: -  m  - the requested matrix

581:    Output Parameter:
582: .  a  - pointer to the values

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

587:    Level: advanced

589: .seealso: DSRestoreArray(), DSGetArrayReal()
590: @*/
591: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
592: {
593:   PetscFunctionBegin;
595:   DSCheckAlloc(ds,1);
596:   DSCheckValidMat(ds,m,2);
597:   PetscAssertPointer(a,3);
598:   PetscCall(MatDenseGetArray(ds->omat[m],a));
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }

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

605:    Not Collective

607:    Input Parameters:
608: +  ds - the direct solver context
609: .  m  - the requested matrix
610: -  a  - pointer to the values

612:    Level: advanced

614: .seealso: DSGetArray(), DSGetArrayReal()
615: @*/
616: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
617: {
618:   PetscFunctionBegin;
620:   DSCheckAlloc(ds,1);
621:   DSCheckValidMat(ds,m,2);
622:   PetscAssertPointer(a,3);
623:   PetscCall(MatDenseRestoreArray(ds->omat[m],a));
624:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
625:   PetscFunctionReturn(PETSC_SUCCESS);
626: }

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

632:    Not Collective

634:    Input Parameters:
635: +  ds - the direct solver context
636: -  m  - the requested matrix

638:    Output Parameter:
639: .  a  - pointer to the values

641:    Note:
642:    This function can be used only for DS_MAT_T and DS_MAT_D. These matrices always
643:    store real values, even in complex scalars, that is why the returned pointer is
644:    PetscReal.

646:    Level: advanced

648: .seealso: DSRestoreArrayReal(), DSGetArray()
649: @*/
650: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
651: {
652: #if defined(PETSC_USE_COMPLEX)
653:   PetscScalar *as;
654: #endif

656:   PetscFunctionBegin;
658:   DSCheckAlloc(ds,1);
659:   DSCheckValidMatReal(ds,m,2);
660:   PetscAssertPointer(a,3);
661: #if defined(PETSC_USE_COMPLEX)
662:   PetscCall(MatDenseGetArray(ds->omat[m],&as));
663:   *a = (PetscReal*)as;
664: #else
665:   PetscCall(MatDenseGetArray(ds->omat[m],a));
666: #endif
667:   PetscFunctionReturn(PETSC_SUCCESS);
668: }

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

673:    Not Collective

675:    Input Parameters:
676: +  ds - the direct solver context
677: .  m  - the requested matrix
678: -  a  - pointer to the values

680:    Level: advanced

682: .seealso: DSGetArrayReal(), DSGetArray()
683: @*/
684: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
685: {
686: #if defined(PETSC_USE_COMPLEX)
687:   PetscScalar *as;
688: #endif

690:   PetscFunctionBegin;
692:   DSCheckAlloc(ds,1);
693:   DSCheckValidMatReal(ds,m,2);
694:   PetscAssertPointer(a,3);
695: #if defined(PETSC_USE_COMPLEX)
696:   PetscCall(MatDenseRestoreArray(ds->omat[m],&as));
697:   *a = NULL;
698: #else
699:   PetscCall(MatDenseRestoreArray(ds->omat[m],a));
700: #endif
701:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
702:   PetscFunctionReturn(PETSC_SUCCESS);
703: }

705: /*@
706:    DSSolve - Solves the problem.

708:    Logically Collective

710:    Input Parameters:
711: +  ds   - the direct solver context
712: .  eigr - array to store the computed eigenvalues (real part)
713: -  eigi - array to store the computed eigenvalues (imaginary part)

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

719:    Level: intermediate

721: .seealso: DSSort(), DSStateType
722: @*/
723: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
724: {
725:   PetscFunctionBegin;
728:   DSCheckAlloc(ds,1);
729:   PetscAssertPointer(eigr,2);
730:   if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(PETSC_SUCCESS);
731:   PetscCheck(ds->ops->solve[ds->method],PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
732:   PetscCall(PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
733:   PetscCall(PetscLogEventBegin(DS_Solve,ds,0,0,0));
734:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
735:   PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
736:   PetscCall(PetscFPTrapPop());
737:   PetscCall(PetscLogEventEnd(DS_Solve,ds,0,0,0));
738:   PetscCall(PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]));
739:   ds->state = DS_STATE_CONDENSED;
740:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
741:   PetscFunctionReturn(PETSC_SUCCESS);
742: }

744: /*@
745:    DSSort - Sorts the result of DSSolve() according to a given sorting
746:    criterion.

748:    Logically Collective

750:    Input Parameters:
751: +  ds   - the direct solver context
752: .  rr   - (optional) array containing auxiliary values (real part)
753: -  ri   - (optional) array containing auxiliary values (imaginary part)

755:    Input/Output Parameters:
756: +  eigr - array containing the computed eigenvalues (real part)
757: .  eigi - array containing the computed eigenvalues (imaginary part)
758: -  k    - (optional) number of elements in the leading group

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

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

773:    Level: intermediate

775: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
776: @*/
777: PetscErrorCode DSSort(DS ds,PetscScalar eigr[],PetscScalar eigi[],PetscScalar rr[],PetscScalar ri[],PetscInt *k)
778: {
779:   PetscInt       i;

781:   PetscFunctionBegin;
784:   DSCheckSolved(ds,1);
785:   PetscAssertPointer(eigr,2);
786:   if (rr) PetscAssertPointer(rr,4);
787:   PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
788:   PetscCheck(ds->sc,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
789:   PetscCheck(!k || rr,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");

791:   for (i=0;i<ds->n;i++) ds->perm[i] = i;   /* initialize to trivial permutation */
792:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
793:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
794:   PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
795:   PetscCall(PetscFPTrapPop());
796:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
797:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
798:   PetscCall(PetscInfo(ds,"Finished sorting\n"));
799:   PetscFunctionReturn(PETSC_SUCCESS);
800: }

802: /*@
803:    DSSortWithPermutation - Reorders the result of DSSolve() according to a given
804:    permutation.

806:    Logically Collective

808:    Input Parameters:
809: +  ds   - the direct solver context
810: -  perm - permutation that indicates the new ordering

812:    Input/Output Parameters:
813: +  eigr - array with the reordered eigenvalues (real part)
814: -  eigi - array with the reordered eigenvalues (imaginary part)

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

821:    Level: advanced

823: .seealso: DSSolve(), DSSort()
824: @*/
825: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt perm[],PetscScalar eigr[],PetscScalar eigi[])
826: {
827:   PetscFunctionBegin;
830:   DSCheckSolved(ds,1);
831:   PetscAssertPointer(perm,2);
832:   PetscAssertPointer(eigr,3);
833:   PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");

835:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
836:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
837:   PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
838:   PetscCall(PetscFPTrapPop());
839:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
840:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
841:   PetscCall(PetscInfo(ds,"Finished sorting\n"));
842:   PetscFunctionReturn(PETSC_SUCCESS);
843: }

845: /*@
846:    DSSynchronize - Make sure that all processes have the same data, performing
847:    communication if necessary.

849:    Collective

851:    Input Parameter:
852: .  ds   - the direct solver context

854:    Input/Output Parameters:
855: +  eigr - (optional) array with the computed eigenvalues (real part)
856: -  eigi - (optional) array with the computed eigenvalues (imaginary part)

858:    Notes:
859:    When the DS has been created with a communicator with more than one process,
860:    the internal data, especially the computed matrices, may diverge in the
861:    different processes. This happens when using multithreaded BLAS and may
862:    cause numerical issues in some ill-conditioned problems. This function
863:    performs the necessary communication among the processes so that the
864:    internal data is exactly equal in all of them.

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

871:    Level: developer

873: .seealso: DSSetParallel(), DSSolve(), DSSort()
874: @*/
875: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
876: {
877:   PetscMPIInt    size;

879:   PetscFunctionBegin;
882:   DSCheckAlloc(ds,1);
883:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
884:   if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
885:     PetscCall(PetscLogEventBegin(DS_Synchronize,ds,0,0,0));
886:     PetscUseTypeMethod(ds,synchronize,eigr,eigi);
887:     PetscCall(PetscLogEventEnd(DS_Synchronize,ds,0,0,0));
888:     PetscCall(PetscObjectStateIncrease((PetscObject)ds));
889:     PetscCall(PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]));
890:   }
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }

894: /*@
895:    DSVectors - Compute vectors associated to the dense system such
896:    as eigenvectors.

898:    Logically Collective

900:    Input Parameters:
901: +  ds  - the direct solver context
902: -  mat - the matrix, used to indicate which vectors are required

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

907:    Output Parameter:
908: .  rnorm - (optional) computed residual norm

910:    Notes:
911:    Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
912:    compute right or left eigenvectors, or left or right singular vectors,
913:    respectively.

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

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

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

927:    Level: intermediate

929: .seealso: DSSolve()
930: @*/
931: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
932: {
933:   PetscFunctionBegin;
936:   DSCheckAlloc(ds,1);
938:   PetscCheck(mat<DS_NUM_MAT,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
939:   PetscCheck(!rnorm || j,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
940:   if (!ds->omat[mat]) PetscCall(DSAllocateMat_Private(ds,mat));
941:   if (!j) PetscCall(PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]));
942:   PetscCall(PetscLogEventBegin(DS_Vectors,ds,0,0,0));
943:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
944:   PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
945:   PetscCall(PetscFPTrapPop());
946:   PetscCall(PetscLogEventEnd(DS_Vectors,ds,0,0,0));
947:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
948:   PetscFunctionReturn(PETSC_SUCCESS);
949: }

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

955:    Logically Collective

957:    Input Parameters:
958: .  ds - the direct solver context

960:    Level: advanced

962: .seealso: DSSolve(), DSSetExtraRow()
963: @*/
964: PetscErrorCode DSUpdateExtraRow(DS ds)
965: {
966:   PetscFunctionBegin;
969:   DSCheckAlloc(ds,1);
970:   PetscCheck(ds->extrarow,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
971:   PetscCall(PetscInfo(ds,"Updating extra row\n"));
972:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
973:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
974:   PetscUseTypeMethod(ds,update);
975:   PetscCall(PetscFPTrapPop());
976:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
977:   PetscFunctionReturn(PETSC_SUCCESS);
978: }

980: /*@
981:    DSCond - Compute the condition number.

983:    Logically Collective

985:    Input Parameters:
986: +  ds - the direct solver context
987: -  cond - the computed condition number

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

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

996:    Does not take into account the extra row.

998:    Level: advanced

1000: .seealso: DSSolve(), DSSetExtraRow()
1001: @*/
1002: PetscErrorCode DSCond(DS ds,PetscReal *cond)
1003: {
1004:   PetscFunctionBegin;
1007:   DSCheckAlloc(ds,1);
1008:   PetscAssertPointer(cond,2);
1009:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1010:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1011:   PetscUseTypeMethod(ds,cond,cond);
1012:   PetscCall(PetscFPTrapPop());
1013:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1014:   PetscCall(PetscInfo(ds,"Computed condition number = %g\n",(double)*cond));
1015:   PetscFunctionReturn(PETSC_SUCCESS);
1016: }

1018: /*@
1019:    DSTranslateHarmonic - Computes a translation of the dense system.

1021:    Logically Collective

1023:    Input Parameters:
1024: +  ds      - the direct solver context
1025: .  tau     - the translation amount
1026: .  beta    - last component of vector b
1027: -  recover - boolean flag to indicate whether to recover or not

1029:    Output Parameters:
1030: +  g       - the computed vector (optional)
1031: -  gamma   - scale factor (optional)

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

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

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

1046:    Level: developer

1048: .seealso: DSTranslateRKS()
1049: @*/
1050: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
1051: {
1052:   PetscFunctionBegin;
1055:   DSCheckAlloc(ds,1);
1056:   if (recover) PetscCall(PetscInfo(ds,"Undoing the translation\n"));
1057:   else PetscCall(PetscInfo(ds,"Computing the translation\n"));
1058:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1059:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1060:   PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
1061:   PetscCall(PetscFPTrapPop());
1062:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1063:   ds->state = DS_STATE_RAW;
1064:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1065:   PetscFunctionReturn(PETSC_SUCCESS);
1066: }

1068: /*@
1069:    DSTranslateRKS - Computes a modification of the dense system corresponding
1070:    to an update of the shift in a rational Krylov method.

1072:    Logically Collective

1074:    Input Parameters:
1075: +  ds    - the direct solver context
1076: -  alpha - the translation amount

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

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

1090:    Level: developer

1092: .seealso: DSTranslateHarmonic()
1093: @*/
1094: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
1095: {
1096:   PetscFunctionBegin;
1099:   DSCheckAlloc(ds,1);
1100:   PetscCall(PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha)));
1101:   PetscCall(DSSetCompact(ds,PETSC_FALSE));
1102:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1103:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1104:   PetscUseTypeMethod(ds,transrks,alpha);
1105:   PetscCall(PetscFPTrapPop());
1106:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1107:   ds->state = DS_STATE_RAW;
1108:   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1109:   PetscFunctionReturn(PETSC_SUCCESS);
1110: }