Actual source code: dsops.c

slepc-3.17.2 2022-08-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: {
 36:   *ld = ds->ld;
 37:   PetscFunctionReturn(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:   PetscFunctionReturn(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:   PetscFunctionReturn(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->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:   PetscFunctionReturn(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:   PetscFunctionReturn(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);
235:   PetscLogEventBegin(DS_Other,ds,0,0,0);
236:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
237:   (*ds->ops->truncate)(ds,n,trim);
238:   PetscFPTrapPop();
239:   PetscLogEventEnd(DS_Other,ds,0,0,0);
240:   PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n);
241:   old = ds->state;
242:   ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
243:   if (old!=ds->state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]);
244:   PetscObjectStateIncrease((PetscObject)ds);
245:   PetscFunctionReturn(0);
246: }

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

251:    Not Collective

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

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

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

264:    Level: developer

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

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

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

289:    Not Collective

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

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

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

302:    Level: developer

304: .seealso: DSGetMat()
305: @*/
306: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
307: {
310:   DSCheckValidMat(ds,t,2);
312:   if (ds->ops->hermitian) (*ds->ops->hermitian)(ds,t,flg);
313:   else *flg = PETSC_FALSE;
314:   PetscFunctionReturn(0);
315: }

317: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
318: {
319: #if !defined(PETSC_USE_COMPLEX)
320:   PetscScalar *A = ds->mat[DS_MAT_A];
321: #endif

323: #if !defined(PETSC_USE_COMPLEX)
324:   if (A[l+(*k)+(l+(*k)-1)*ds->ld] != 0.0) {
325:     if (l+(*k)<n-1) (*k)++;
326:     else (*k)--;
327:   }
328: #endif
329:   PetscFunctionReturn(0);
330: }

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

336:    Not Collective

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

344:    Output Parameter:
345: .  k  - the possibly modified value of the truncation size

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

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

356:    Level: advanced

358: .seealso: DSTruncate(), DSSetDimensions()
359: @*/
360: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
361: {
363:   DSCheckAlloc(ds,1);
367:   if (ds->ops->gettruncatesize) (*ds->ops->gettruncatesize)(ds,l?l:ds->l,n?n:ds->n,k);
368:   PetscFunctionReturn(0);
369: }

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

375:    Not Collective

377:    Input Parameters:
378: +  ds - the direct solver context
379: -  m  - the requested matrix

381:    Output Parameter:
382: .  A  - Mat object

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

391:    When no longer needed, the user can either destroy the matrix or call
392:    DSRestoreMat(). The latter will copy back the modified values.

394:    Level: advanced

396: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
397: @*/
398: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
399: {
400:   PetscInt       j,rows,cols,arows,acols;
401:   PetscBool      create=PETSC_FALSE,flg;
402:   PetscScalar    *pA,*M;

405:   DSCheckAlloc(ds,1);
406:   DSCheckValidMat(ds,m,2);

410:   DSMatGetSize(ds,m,&rows,&cols);
411:   if (!ds->omat[m]) create=PETSC_TRUE;
412:   else {
413:     MatGetSize(ds->omat[m],&arows,&acols);
414:     if (arows!=rows || acols!=cols) {
415:       MatDestroy(&ds->omat[m]);
416:       create=PETSC_TRUE;
417:     }
418:   }
419:   if (create) MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);

421:   /* set Hermitian flag */
422:   DSMatIsHermitian(ds,m,&flg);
423:   MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);

425:   /* copy entries */
426:   PetscObjectReference((PetscObject)ds->omat[m]);
427:   *A = ds->omat[m];
428:   M  = ds->mat[m];
429:   MatDenseGetArray(*A,&pA);
430:   for (j=0;j<cols;j++) PetscArraycpy(pA+j*rows,M+j*ds->ld,rows);
431:   MatDenseRestoreArray(*A,&pA);
432:   PetscFunctionReturn(0);
433: }

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

438:    Not Collective

440:    Input Parameters:
441: +  ds - the direct solver context
442: .  m  - the requested matrix
443: -  A  - the fetched Mat object

445:    Notes:
446:    A call to this function must match a previous call of DSGetMat().
447:    The effect is that the contents of the Mat are copied back to the
448:    DS internal array, and the matrix is destroyed.

450:    It is not compulsory to call this function, the matrix obtained with
451:    DSGetMat() can simply be destroyed if entries need not be copied back.

453:    Level: advanced

455: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
456: @*/
457: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
458: {
459:   PetscInt       j,rows,cols;
460:   PetscScalar    *pA,*M;

463:   DSCheckAlloc(ds,1);
464:   DSCheckValidMat(ds,m,2);

469:   MatGetSize(*A,&rows,&cols);
470:   M  = ds->mat[m];
471:   MatDenseGetArray(*A,&pA);
472:   for (j=0;j<cols;j++) PetscArraycpy(M+j*ds->ld,pA+j*rows,rows);
473:   MatDenseRestoreArray(*A,&pA);
474:   MatDestroy(A);
475:   PetscFunctionReturn(0);
476: }

478: /*@C
479:    DSGetArray - Returns a pointer to one of the internal arrays used to
480:    represent matrices. You MUST call DSRestoreArray() when you no longer
481:    need to access the array.

483:    Not Collective

485:    Input Parameters:
486: +  ds - the direct solver context
487: -  m  - the requested matrix

489:    Output Parameter:
490: .  a  - pointer to the values

492:    Level: advanced

494: .seealso: DSRestoreArray(), DSGetArrayReal()
495: @*/
496: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
497: {
499:   DSCheckAlloc(ds,1);
500:   DSCheckValidMat(ds,m,2);
502:   *a = ds->mat[m];
503:   CHKMEMQ;
504:   PetscFunctionReturn(0);
505: }

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

510:    Not Collective

512:    Input Parameters:
513: +  ds - the direct solver context
514: .  m  - the requested matrix
515: -  a  - pointer to the values

517:    Level: advanced

519: .seealso: DSGetArray(), DSGetArrayReal()
520: @*/
521: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
522: {
524:   DSCheckAlloc(ds,1);
525:   DSCheckValidMat(ds,m,2);
527:   CHKMEMQ;
528:   *a = 0;
529:   PetscObjectStateIncrease((PetscObject)ds);
530:   PetscFunctionReturn(0);
531: }

533: /*@C
534:    DSGetArrayReal - Returns a pointer to one of the internal arrays used to
535:    represent real matrices. You MUST call DSRestoreArrayReal() when you no longer
536:    need to access the array.

538:    Not Collective

540:    Input Parameters:
541: +  ds - the direct solver context
542: -  m  - the requested matrix

544:    Output Parameter:
545: .  a  - pointer to the values

547:    Level: advanced

549: .seealso: DSRestoreArrayReal(), DSGetArray()
550: @*/
551: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
552: {
554:   DSCheckAlloc(ds,1);
555:   DSCheckValidMatReal(ds,m,2);
557:   *a = ds->rmat[m];
558:   CHKMEMQ;
559:   PetscFunctionReturn(0);
560: }

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

565:    Not Collective

567:    Input Parameters:
568: +  ds - the direct solver context
569: .  m  - the requested matrix
570: -  a  - pointer to the values

572:    Level: advanced

574: .seealso: DSGetArrayReal(), DSGetArray()
575: @*/
576: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
577: {
579:   DSCheckAlloc(ds,1);
580:   DSCheckValidMatReal(ds,m,2);
582:   CHKMEMQ;
583:   *a = 0;
584:   PetscObjectStateIncrease((PetscObject)ds);
585:   PetscFunctionReturn(0);
586: }

588: /*@
589:    DSSolve - Solves the problem.

591:    Logically Collective on ds

593:    Input Parameters:
594: +  ds   - the direct solver context
595: .  eigr - array to store the computed eigenvalues (real part)
596: -  eigi - array to store the computed eigenvalues (imaginary part)

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

602:    Level: intermediate

604: .seealso: DSSort(), DSStateType
605: @*/
606: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
607: {
610:   DSCheckAlloc(ds,1);
612:   if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(0);
614:   PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
615:   PetscLogEventBegin(DS_Solve,ds,0,0,0);
616:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
617:   (*ds->ops->solve[ds->method])(ds,eigr,eigi);
618:   PetscFPTrapPop();
619:   PetscLogEventEnd(DS_Solve,ds,0,0,0);
620:   PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
621:   ds->state = DS_STATE_CONDENSED;
622:   PetscObjectStateIncrease((PetscObject)ds);
623:   PetscFunctionReturn(0);
624: }

626: /*@C
627:    DSSort - Sorts the result of DSSolve() according to a given sorting
628:    criterion.

630:    Logically Collective on ds

632:    Input Parameters:
633: +  ds   - the direct solver context
634: .  rr   - (optional) array containing auxiliary values (real part)
635: -  ri   - (optional) array containing auxiliary values (imaginary part)

637:    Input/Output Parameters:
638: +  eigr - array containing the computed eigenvalues (real part)
639: .  eigi - array containing the computed eigenvalues (imaginary part)
640: -  k    - (optional) number of elements in the leading group

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

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

655:    Level: intermediate

657: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
658: @*/
659: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
660: {
661:   PetscInt       i;

665:   DSCheckSolved(ds,1);

673:   for (i=0;i<ds->n;i++) ds->perm[i] = i;   /* initialize to trivial permutation */
674:   PetscLogEventBegin(DS_Other,ds,0,0,0);
675:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
676:   (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
677:   PetscFPTrapPop();
678:   PetscLogEventEnd(DS_Other,ds,0,0,0);
679:   PetscObjectStateIncrease((PetscObject)ds);
680:   PetscInfo(ds,"Finished sorting\n");
681:   PetscFunctionReturn(0);
682: }

684: /*@C
685:    DSSortWithPermutation - Reorders the result of DSSolve() according to a given
686:    permutation.

688:    Logically Collective on ds

690:    Input Parameters:
691: +  ds   - the direct solver context
692: -  perm - permutation that indicates the new ordering

694:    Input/Output Parameters:
695: +  eigr - array with the reordered eigenvalues (real part)
696: -  eigi - array with the reordered eigenvalues (imaginary part)

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

703:    Level: advanced

705: .seealso: DSSolve(), DSSort()
706: @*/
707: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)
708: {
711:   DSCheckSolved(ds,1);

717:   PetscLogEventBegin(DS_Other,ds,0,0,0);
718:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
719:   (*ds->ops->sortperm)(ds,perm,eigr,eigi);
720:   PetscFPTrapPop();
721:   PetscLogEventEnd(DS_Other,ds,0,0,0);
722:   PetscObjectStateIncrease((PetscObject)ds);
723:   PetscInfo(ds,"Finished sorting\n");
724:   PetscFunctionReturn(0);
725: }

727: /*@
728:    DSSynchronize - Make sure that all processes have the same data, performing
729:    communication if necessary.

731:    Collective on ds

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

736:    Input/Output Parameters:
737: +  eigr - (optional) array with the computed eigenvalues (real part)
738: -  eigi - (optional) array with the computed eigenvalues (imaginary part)

740:    Notes:
741:    When the DS has been created with a communicator with more than one process,
742:    the internal data, especially the computed matrices, may diverge in the
743:    different processes. This happens when using multithreaded BLAS and may
744:    cause numerical issues in some ill-conditioned problems. This function
745:    performs the necessary communication among the processes so that the
746:    internal data is exactly equal in all of them.

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

753:    Level: developer

755: .seealso: DSSetParallel(), DSSolve(), DSSort()
756: @*/
757: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
758: {
759:   PetscMPIInt    size;

763:   DSCheckAlloc(ds,1);
764:   MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
765:   if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
766:     PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
767:     if (ds->ops->synchronize) (*ds->ops->synchronize)(ds,eigr,eigi);
768:     PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
769:     PetscObjectStateIncrease((PetscObject)ds);
770:     PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
771:   }
772:   PetscFunctionReturn(0);
773: }

775: /*@C
776:    DSVectors - Compute vectors associated to the dense system such
777:    as eigenvectors.

779:    Logically Collective on ds

781:    Input Parameters:
782: +  ds  - the direct solver context
783: -  mat - the matrix, used to indicate which vectors are required

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

788:    Output Parameter:
789: .  rnorm - (optional) computed residual norm

791:    Notes:
792:    Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
793:    compute right or left eigenvectors, or left or right singular vectors,
794:    respectively.

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

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

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

808:    Level: intermediate

810: .seealso: DSSolve()
811: @*/
812: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
813: {
816:   DSCheckAlloc(ds,1);
821:   if (!ds->mat[mat]) DSAllocateMat_Private(ds,mat);
822:   if (!j) PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]);
823:   PetscLogEventBegin(DS_Vectors,ds,0,0,0);
824:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
825:   (*ds->ops->vectors)(ds,mat,j,rnorm);
826:   PetscFPTrapPop();
827:   PetscLogEventEnd(DS_Vectors,ds,0,0,0);
828:   PetscObjectStateIncrease((PetscObject)ds);
829:   PetscFunctionReturn(0);
830: }

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

836:    Not Collective

838:    Input Parameters:
839: .  ds - the direct solver context

841:    Level: advanced

843: .seealso: DSSolve(), DSSetExtraRow()
844: @*/
845: PetscErrorCode DSUpdateExtraRow(DS ds)
846: {
849:   DSCheckAlloc(ds,1);
852:   PetscInfo(ds,"Updating extra row\n");
853:   PetscLogEventBegin(DS_Other,ds,0,0,0);
854:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
855:   (*ds->ops->update)(ds);
856:   PetscFPTrapPop();
857:   PetscLogEventEnd(DS_Other,ds,0,0,0);
858:   PetscFunctionReturn(0);
859: }

861: /*@
862:    DSCond - Compute the inf-norm condition number of the first matrix
863:    as cond(A) = norm(A)*norm(inv(A)).

865:    Not Collective

867:    Input Parameters:
868: +  ds - the direct solver context
869: -  cond - the computed condition number

871:    Level: advanced

873: .seealso: DSSolve()
874: @*/
875: PetscErrorCode DSCond(DS ds,PetscReal *cond)
876: {
879:   DSCheckAlloc(ds,1);
882:   PetscLogEventBegin(DS_Other,ds,0,0,0);
883:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
884:   (*ds->ops->cond)(ds,cond);
885:   PetscFPTrapPop();
886:   PetscLogEventEnd(DS_Other,ds,0,0,0);
887:   PetscInfo(ds,"Computed condition number = %g\n",(double)*cond);
888:   PetscFunctionReturn(0);
889: }

891: /*@C
892:    DSTranslateHarmonic - Computes a translation of the dense system.

894:    Logically Collective on ds

896:    Input Parameters:
897: +  ds      - the direct solver context
898: .  tau     - the translation amount
899: .  beta    - last component of vector b
900: -  recover - boolean flag to indicate whether to recover or not

902:    Output Parameters:
903: +  g       - the computed vector (optional)
904: -  gamma   - scale factor (optional)

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

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

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

919:    Level: developer

921: .seealso: DSTranslateRKS()
922: @*/
923: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
924: {
927:   DSCheckAlloc(ds,1);
929:   if (recover) PetscInfo(ds,"Undoing the translation\n");
930:   else PetscInfo(ds,"Computing the translation\n");
931:   PetscLogEventBegin(DS_Other,ds,0,0,0);
932:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
933:   (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
934:   PetscFPTrapPop();
935:   PetscLogEventEnd(DS_Other,ds,0,0,0);
936:   ds->state = DS_STATE_RAW;
937:   PetscObjectStateIncrease((PetscObject)ds);
938:   PetscFunctionReturn(0);
939: }

941: /*@
942:    DSTranslateRKS - Computes a modification of the dense system corresponding
943:    to an update of the shift in a rational Krylov method.

945:    Logically Collective on ds

947:    Input Parameters:
948: +  ds    - the direct solver context
949: -  alpha - the translation amount

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

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

963:    Level: developer

965: .seealso: DSTranslateHarmonic()
966: @*/
967: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
968: {
971:   DSCheckAlloc(ds,1);
973:   PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha));
974:   PetscLogEventBegin(DS_Other,ds,0,0,0);
975:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
976:   (*ds->ops->transrks)(ds,alpha);
977:   PetscFPTrapPop();
978:   PetscLogEventEnd(DS_Other,ds,0,0,0);
979:   ds->state   = DS_STATE_RAW;
980:   ds->compact = PETSC_FALSE;
981:   PetscObjectStateIncrease((PetscObject)ds);
982:   PetscFunctionReturn(0);
983: }

985: /*@
986:    DSCopyMat - Copies the contents of a sequential dense Mat object to
987:    the indicated DS matrix, or vice versa.

989:    Not Collective

991:    Input Parameters:
992: +  ds   - the direct solver context
993: .  m    - the requested matrix
994: .  mr   - first row of m to be considered
995: .  mc   - first column of m to be considered
996: .  A    - Mat object
997: .  Ar   - first row of A to be considered
998: .  Ac   - first column of A to be considered
999: .  rows - number of rows to copy
1000: .  cols - number of columns to copy
1001: -  out  - whether the data is copied out of the DS

1003:    Note:
1004:    If out=true, the values of the DS matrix m are copied to A, otherwise
1005:    the entries of A are copied to the DS.

1007:    Level: developer

1009: .seealso: DSGetMat()
1010: @*/
1011: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)
1012: {
1013:   PetscInt       j,mrows,mcols,arows,acols;
1014:   PetscScalar    *pA,*M;

1017:   DSCheckAlloc(ds,1);
1019:   DSCheckValidMat(ds,m,2);
1028:   if (!rows || !cols) PetscFunctionReturn(0);

1030:   DSMatGetSize(ds,m,&mrows,&mcols);
1031:   MatGetSize(A,&arows,&acols);

1040:   M  = ds->mat[m];
1041:   MatDenseGetArray(A,&pA);
1042:   for (j=0;j<cols;j++) {
1043:     if (out) PetscArraycpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows);
1044:     else PetscArraycpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows);
1045:   }
1046:   MatDenseRestoreArray(A,&pA);
1047:   PetscFunctionReturn(0);
1048: }