LCOV - code coverage report
Current view: top level - sys/classes/ds/interface - dsops.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 339 342 99.1 %
Date: 2024-03-29 00:55:19 Functions: 27 27 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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             : */
      13             : 
      14             : #include <slepc/private/dsimpl.h>      /*I "slepcds.h" I*/
      15             : 
      16             : /*@
      17             :    DSGetLeadingDimension - Returns the leading dimension of the allocated
      18             :    matrices.
      19             : 
      20             :    Not Collective
      21             : 
      22             :    Input Parameter:
      23             : .  ds - the direct solver context
      24             : 
      25             :    Output Parameter:
      26             : .  ld - leading dimension (maximum allowed dimension for the matrices)
      27             : 
      28             :    Level: advanced
      29             : 
      30             : .seealso: DSAllocate(), DSSetDimensions()
      31             : @*/
      32       37175 : PetscErrorCode DSGetLeadingDimension(DS ds,PetscInt *ld)
      33             : {
      34       37175 :   PetscFunctionBegin;
      35       37175 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
      36       37175 :   PetscAssertPointer(ld,2);
      37       37175 :   *ld = ds->ld;
      38       37175 :   PetscFunctionReturn(PETSC_SUCCESS);
      39             : }
      40             : 
      41             : /*@
      42             :    DSSetState - Change the state of the DS object.
      43             : 
      44             :    Logically Collective
      45             : 
      46             :    Input Parameters:
      47             : +  ds    - the direct solver context
      48             : -  state - the new state
      49             : 
      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.
      55             : 
      56             :    This function is normally used to return to the raw state when the
      57             :    condensed structure is destroyed.
      58             : 
      59             :    Level: advanced
      60             : 
      61             : .seealso: DSGetState()
      62             : @*/
      63       22434 : PetscErrorCode DSSetState(DS ds,DSStateType state)
      64             : {
      65       22434 :   PetscFunctionBegin;
      66       22434 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
      67       89736 :   PetscValidLogicalCollectiveEnum(ds,state,2);
      68       22434 :   switch (state) {
      69       22434 :     case DS_STATE_RAW:
      70             :     case DS_STATE_INTERMEDIATE:
      71             :     case DS_STATE_CONDENSED:
      72             :     case DS_STATE_TRUNCATED:
      73       22434 :       if (ds->state!=state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]));
      74       22434 :       ds->state = state;
      75       22434 :       break;
      76           0 :     default:
      77           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
      78             :   }
      79       22434 :   PetscFunctionReturn(PETSC_SUCCESS);
      80             : }
      81             : 
      82             : /*@
      83             :    DSGetState - Returns the current state.
      84             : 
      85             :    Not Collective
      86             : 
      87             :    Input Parameter:
      88             : .  ds - the direct solver context
      89             : 
      90             :    Output Parameter:
      91             : .  state - current state
      92             : 
      93             :    Level: advanced
      94             : 
      95             : .seealso: DSSetState()
      96             : @*/
      97           2 : PetscErrorCode DSGetState(DS ds,DSStateType *state)
      98             : {
      99           2 :   PetscFunctionBegin;
     100           2 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     101           2 :   PetscAssertPointer(state,2);
     102           2 :   *state = ds->state;
     103           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     104             : }
     105             : 
     106             : /*@
     107             :    DSSetDimensions - Resize the matrices in the DS object.
     108             : 
     109             :    Logically Collective
     110             : 
     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)
     116             : 
     117             :    Notes:
     118             :    The internal arrays are not reallocated.
     119             : 
     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.
     122             : 
     123             :    Level: intermediate
     124             : 
     125             : .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
     126             : @*/
     127       31234 : PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
     128             : {
     129       31234 :   PetscInt       on,ol,ok;
     130       31234 :   PetscBool      issvd;
     131             : 
     132       31234 :   PetscFunctionBegin;
     133       31234 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     134       31234 :   DSCheckAlloc(ds,1);
     135      124936 :   PetscValidLogicalCollectiveInt(ds,n,2);
     136      124936 :   PetscValidLogicalCollectiveInt(ds,l,3);
     137      124936 :   PetscValidLogicalCollectiveInt(ds,k,4);
     138       31234 :   on = ds->n; ol = ds->l; ok = ds->k;
     139       31234 :   if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
     140         131 :     ds->n = ds->extrarow? ds->ld-1: ds->ld;
     141             :   } else {
     142       31103 :     PetscCheck(n>=0 && n<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
     143       31103 :     PetscCall(PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,""));  /* SVD and GSVD have extra column instead of extra row */
     144       31103 :     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");
     145       31103 :     ds->n = n;
     146             :   }
     147       31234 :   ds->t = ds->n;   /* truncated length equal to the new dimension */
     148       31234 :   if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
     149         487 :     ds->l = 0;
     150             :   } else {
     151       30747 :     PetscCheck(l>=0 && l<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
     152       30747 :     ds->l = l;
     153             :   }
     154       31234 :   if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
     155        1419 :     ds->k = ds->n/2;
     156             :   } else {
     157       29815 :     PetscCheck(k>=0 || k<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
     158       29815 :     ds->k = k;
     159             :   }
     160       31234 :   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));
     161       31234 :   PetscFunctionReturn(PETSC_SUCCESS);
     162             : }
     163             : 
     164             : /*@
     165             :    DSGetDimensions - Returns the current dimensions.
     166             : 
     167             :    Not Collective
     168             : 
     169             :    Input Parameter:
     170             : .  ds - the direct solver context
     171             : 
     172             :    Output Parameters:
     173             : +  n  - the current size
     174             : .  l  - number of locked (inactive) leading columns
     175             : .  k  - intermediate dimension (e.g., position of arrow)
     176             : -  t  - truncated length
     177             : 
     178             :    Note:
     179             :    The t parameter makes sense only if DSTruncate() has been called.
     180             :    Otherwise its value equals n.
     181             : 
     182             :    Some DS types have additional dimensions, e.g. the number of columns
     183             :    in DSSVD. For these, you should call a specific interface function.
     184             : 
     185             :    Level: intermediate
     186             : 
     187             : .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
     188             : @*/
     189       22251 : PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
     190             : {
     191       22251 :   PetscFunctionBegin;
     192       22251 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     193       22251 :   DSCheckAlloc(ds,1);
     194       22251 :   if (n) *n = ds->n;
     195       22251 :   if (l) *l = ds->l;
     196       22251 :   if (k) *k = ds->k;
     197       22251 :   if (t) *t = ds->t;
     198       22251 :   PetscFunctionReturn(PETSC_SUCCESS);
     199             : }
     200             : 
     201             : /*@
     202             :    DSTruncate - Truncates the system represented in the DS object.
     203             : 
     204             :    Logically Collective
     205             : 
     206             :    Input Parameters:
     207             : +  ds   - the direct solver context
     208             : .  n    - the new size
     209             : -  trim - a flag to indicate if the factorization must be trimmed
     210             : 
     211             :    Note:
     212             :    The new size is set to n. Note that in some cases the new size could
     213             :    be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
     214             :    Schur form). In cases where the extra row is meaningful, the first
     215             :    n elements are kept as the extra row for the new system.
     216             : 
     217             :    If the flag trim is turned on, it resets the locked and intermediate
     218             :    dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
     219             :    It also cleans the extra row if being used.
     220             : 
     221             :    The typical usage of trim=true is to truncate the Schur decomposition
     222             :    at the end of a Krylov iteration. In this case, the state must be
     223             :    changed to RAW so that DSVectors() computes eigenvectors from scratch.
     224             : 
     225             :    Level: advanced
     226             : 
     227             : .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
     228             : @*/
     229        9026 : PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
     230             : {
     231        9026 :   DSStateType    old;
     232             : 
     233        9026 :   PetscFunctionBegin;
     234        9026 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     235        9026 :   PetscValidType(ds,1);
     236        9026 :   DSCheckAlloc(ds,1);
     237       36104 :   PetscValidLogicalCollectiveInt(ds,n,2);
     238       36104 :   PetscValidLogicalCollectiveBool(ds,trim,3);
     239        9026 :   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);
     240        9026 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
     241        9026 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     242        9026 :   PetscUseTypeMethod(ds,truncate,n,trim);
     243        9026 :   PetscCall(PetscFPTrapPop());
     244        9026 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
     245       17363 :   PetscCall(PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n));
     246        9026 :   old = ds->state;
     247        9026 :   ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
     248        9026 :   if (old!=ds->state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]));
     249        9026 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     250        9026 :   PetscFunctionReturn(PETSC_SUCCESS);
     251             : }
     252             : 
     253             : /*@
     254             :    DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
     255             : 
     256             :    Not Collective
     257             : 
     258             :    Input Parameters:
     259             : +  ds - the direct solver context
     260             : -  t  - the requested matrix
     261             : 
     262             :    Output Parameters:
     263             : +  n  - the number of rows
     264             : -  m  - the number of columns
     265             : 
     266             :    Note:
     267             :    This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
     268             : 
     269             :    Level: developer
     270             : 
     271             : .seealso: DSSetDimensions(), DSGetMat()
     272             : @*/
     273       58634 : PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
     274             : {
     275       58634 :   PetscInt       rows,cols;
     276             : 
     277       58634 :   PetscFunctionBegin;
     278       58634 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     279       58634 :   PetscValidType(ds,1);
     280       58634 :   DSCheckValidMat(ds,t,2);
     281       58634 :   if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
     282             :   else {
     283       46391 :     if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
     284       39680 :     else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
     285       46391 :     if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
     286       44026 :     else if (t==DS_MAT_D) cols = 1;
     287       40700 :     else cols = ds->n;
     288             :   }
     289       58634 :   if (m) *m = rows;
     290       58634 :   if (n) *n = cols;
     291       58634 :   PetscFunctionReturn(PETSC_SUCCESS);
     292             : }
     293             : 
     294             : /*@
     295             :    DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
     296             : 
     297             :    Not Collective
     298             : 
     299             :    Input Parameters:
     300             : +  ds - the direct solver context
     301             : -  t  - the requested matrix
     302             : 
     303             :    Output Parameter:
     304             : .  flg - the Hermitian flag
     305             : 
     306             :    Note:
     307             :    Does not check the matrix values directly. The flag is set according to the
     308             :    problem structure. For instance, in DSHEP matrix A is Hermitian.
     309             : 
     310             :    Level: developer
     311             : 
     312             : .seealso: DSGetMat()
     313             : @*/
     314       58606 : PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
     315             : {
     316       58606 :   PetscFunctionBegin;
     317       58606 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     318       58606 :   PetscValidType(ds,1);
     319       58606 :   DSCheckValidMat(ds,t,2);
     320       58606 :   PetscAssertPointer(flg,3);
     321       58606 :   *flg = PETSC_FALSE;
     322       58606 :   PetscTryTypeMethod(ds,hermitian,t,flg);
     323       58606 :   PetscFunctionReturn(PETSC_SUCCESS);
     324             : }
     325             : 
     326        1632 : PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
     327             : {
     328             : #if !defined(PETSC_USE_COMPLEX)
     329        1632 :   PetscScalar val;
     330             : #endif
     331             : 
     332        1632 :   PetscFunctionBegin;
     333             : #if !defined(PETSC_USE_COMPLEX)
     334        1632 :   PetscCall(MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val));
     335        1632 :   if (val != 0.0) {
     336         433 :     if (l+(*k)<n-1) (*k)++;
     337           0 :     else (*k)--;
     338             :   }
     339             : #endif
     340        1632 :   PetscFunctionReturn(PETSC_SUCCESS);
     341             : }
     342             : 
     343             : /*@
     344             :    DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
     345             :    to avoid breaking a 2x2 block.
     346             : 
     347             :    Not Collective
     348             : 
     349             :    Input Parameters:
     350             : +  ds - the direct solver context
     351             : .  l  - the size of the locked part (set to 0 to use ds->l)
     352             : -  n  - the total matrix size (set to 0 to use ds->n)
     353             : 
     354             :    Output Parameter:
     355             : .  k  - the wanted truncation size (possibly modified)
     356             : 
     357             :    Notes:
     358             :    This should be called before DSTruncate() to make sure that the truncation
     359             :    does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
     360             : 
     361             :    The total size is n (either user-provided or ds->n if 0 is passed). The
     362             :    size where the truncation is intended is equal to l+k (where l can be
     363             :    equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
     364             :    at the l+k limit, the value of k is increased (or decreased) by 1.
     365             : 
     366             :    Level: advanced
     367             : 
     368             : .seealso: DSTruncate(), DSSetDimensions()
     369             : @*/
     370        4939 : PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
     371             : {
     372        4939 :   PetscFunctionBegin;
     373        4939 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     374        4939 :   DSCheckAlloc(ds,1);
     375       19756 :   PetscValidLogicalCollectiveInt(ds,l,2);
     376       19756 :   PetscValidLogicalCollectiveInt(ds,n,3);
     377        4939 :   PetscAssertPointer(k,4);
     378        4939 :   PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
     379        4939 :   PetscFunctionReturn(PETSC_SUCCESS);
     380             : }
     381             : 
     382             : /*@
     383             :    DSGetMat - Returns a sequential dense Mat object containing the requested
     384             :    matrix.
     385             : 
     386             :    Not Collective
     387             : 
     388             :    Input Parameters:
     389             : +  ds - the direct solver context
     390             : -  m  - the requested matrix
     391             : 
     392             :    Output Parameter:
     393             : .  A  - Mat object
     394             : 
     395             :    Notes:
     396             :    The returned Mat has sizes equal to the current DS dimensions (nxm),
     397             :    and contains the values that would be obtained with DSGetArray()
     398             :    (not DSGetArrayReal()). If the DS was truncated, then the number of rows
     399             :    is equal to the dimension prior to truncation, see DSTruncate().
     400             :    The communicator is always PETSC_COMM_SELF.
     401             : 
     402             :    It is implemented with MatDenseGetSubMatrix(), and when no longer needed
     403             :    the user must call DSRestoreMat() which will invoke MatDenseRestoreSubMatrix().
     404             : 
     405             :    For matrices DS_MAT_T and DS_MAT_D, this function will return a Mat object
     406             :    that cannot be used directly for computations, since it uses compact storage
     407             :    (three and one diagonals for T and D, respectively). In complex scalars, the
     408             :    internal array stores real values, so it is sufficient with 2 columns for T.
     409             : 
     410             :    Level: advanced
     411             : 
     412             : .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
     413             : @*/
     414       58606 : PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
     415             : {
     416       58606 :   PetscInt  rows,cols;
     417       58606 :   PetscBool flg;
     418             : 
     419       58606 :   PetscFunctionBegin;
     420       58606 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     421       58606 :   DSCheckAlloc(ds,1);
     422       58606 :   DSCheckValidMat(ds,m,2);
     423       58606 :   PetscAssertPointer(A,3);
     424             : 
     425       58606 :   PetscCall(DSMatGetSize(ds,m,&rows,&cols));
     426       58606 :   PetscCheck(rows && cols,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must call DSSetDimensions() first");
     427       58606 :   PetscCall(MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A));
     428             : 
     429             :   /* set Hermitian flag */
     430       58606 :   PetscCall(DSMatIsHermitian(ds,m,&flg));
     431       58606 :   PetscCall(MatSetOption(*A,MAT_SYMMETRIC,flg));
     432       58606 :   PetscCall(MatSetOption(*A,MAT_HERMITIAN,flg));
     433       58606 :   PetscFunctionReturn(PETSC_SUCCESS);
     434             : }
     435             : 
     436             : /*@
     437             :    DSRestoreMat - Restores the matrix after DSGetMat() was called.
     438             : 
     439             :    Not Collective
     440             : 
     441             :    Input Parameters:
     442             : +  ds - the direct solver context
     443             : .  m  - the requested matrix
     444             : -  A  - the fetched Mat object
     445             : 
     446             :    Note:
     447             :    A call to this function must match a previous call of DSGetMat().
     448             : 
     449             :    Level: advanced
     450             : 
     451             : .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
     452             : @*/
     453       58606 : PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
     454             : {
     455       58606 :   PetscFunctionBegin;
     456       58606 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     457       58606 :   DSCheckAlloc(ds,1);
     458       58606 :   DSCheckValidMat(ds,m,2);
     459       58606 :   PetscAssertPointer(A,3);
     460             : 
     461       58606 :   PetscCall(MatDenseRestoreSubMatrix(ds->omat[m],A));
     462       58606 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     463       58606 :   PetscFunctionReturn(PETSC_SUCCESS);
     464             : }
     465             : 
     466             : /*@
     467             :    DSGetMatAndColumn - Returns a sequential dense Mat object containing the requested
     468             :    matrix and one of its columns as a Vec.
     469             : 
     470             :    Not Collective
     471             : 
     472             :    Input Parameters:
     473             : +  ds  - the direct solver context
     474             : .  m   - the requested matrix
     475             : -  col - the requested column
     476             : 
     477             :    Output Parameters:
     478             : +  A   - Mat object
     479             : -  v   - Vec object (the column)
     480             : 
     481             :    Notes:
     482             :    This calls DSGetMat() and then it extracts the selected column.
     483             :    The user must call DSRestoreMatAndColumn() to recover the original state.
     484             :    For matrices DS_MAT_T and DS_MAT_D, in complex scalars this function implies
     485             :    copying from real values stored internally to scalar values in the Vec.
     486             : 
     487             :    Level: advanced
     488             : 
     489             : .seealso: DSRestoreMatAndColumn(), DSGetMat()
     490             : @*/
     491        3690 : PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
     492             : {
     493        3690 :   PetscFunctionBegin;
     494        3690 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     495        3690 :   DSCheckAlloc(ds,1);
     496        3690 :   DSCheckValidMat(ds,m,2);
     497        3690 :   PetscAssertPointer(A,4);
     498        3690 :   PetscAssertPointer(v,5);
     499             : 
     500        3690 :   PetscCall(DSGetMat(ds,m,A));
     501        3690 :   if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
     502             :     const PetscScalar *as;
     503             :     PetscScalar       *vs;
     504             :     PetscReal         *ar;
     505             :     PetscInt          i,n,lda;
     506             :     PetscCall(MatCreateVecs(*A,NULL,v));
     507             :     PetscCall(VecGetSize(*v,&n));
     508             :     PetscCall(MatDenseGetLDA(*A,&lda));
     509             :     PetscCall(MatDenseGetArrayRead(*A,&as));
     510             :     PetscCall(VecGetArrayWrite(*v,&vs));
     511             :     ar = (PetscReal*)as;
     512             :     for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
     513             :     PetscCall(VecRestoreArrayWrite(*v,&vs));
     514             :     PetscCall(MatDenseRestoreArrayRead(*A,&as));
     515        3690 :   } else PetscCall(MatDenseGetColumnVec(*A,col,v));
     516        3690 :   PetscFunctionReturn(PETSC_SUCCESS);
     517             : }
     518             : 
     519             : /*@
     520             :    DSRestoreMatAndColumn - Restores the matrix and vector after DSGetMatAndColumn()
     521             :    was called.
     522             : 
     523             :    Not Collective
     524             : 
     525             :    Input Parameters:
     526             : +  ds  - the direct solver context
     527             : .  m   - the requested matrix
     528             : .  col - the requested column
     529             : .  A   - the fetched Mat object
     530             : -  v   - the fetched Vec object
     531             : 
     532             :    Note:
     533             :    A call to this function must match a previous call of DSGetMatAndColumn().
     534             : 
     535             :    Level: advanced
     536             : 
     537             : .seealso: DSGetMatAndColumn(), DSRestoreMat()
     538             : @*/
     539        3690 : PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
     540             : {
     541        3690 :   PetscFunctionBegin;
     542        3690 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     543        3690 :   DSCheckAlloc(ds,1);
     544        3690 :   DSCheckValidMat(ds,m,2);
     545        3690 :   PetscAssertPointer(A,4);
     546        3690 :   PetscAssertPointer(v,5);
     547             : 
     548        3690 :   if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
     549             :     const PetscScalar *vs;
     550             :     PetscScalar       *as;
     551             :     PetscReal         *ar;
     552             :     PetscInt          i,n,lda;
     553             :     PetscCall(VecGetSize(*v,&n));
     554             :     PetscCall(MatDenseGetLDA(*A,&lda));
     555             :     PetscCall(MatDenseGetArray(*A,&as));
     556             :     PetscCall(VecGetArrayRead(*v,&vs));
     557             :     ar = (PetscReal*)as;
     558             :     for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
     559             :     PetscCall(VecRestoreArrayRead(*v,&vs));
     560             :     PetscCall(MatDenseRestoreArray(*A,&as));
     561             :     PetscCall(VecDestroy(v));
     562        3690 :   } else PetscCall(MatDenseRestoreColumnVec(*A,col,v));
     563        3690 :   PetscCall(DSRestoreMat(ds,m,A));
     564        3690 :   PetscFunctionReturn(PETSC_SUCCESS);
     565             : }
     566             : 
     567             : /*@C
     568             :    DSGetArray - Returns a pointer to the internal array of one of the
     569             :    matrices. You MUST call DSRestoreArray() when you no longer
     570             :    need to access the array.
     571             : 
     572             :    Not Collective
     573             : 
     574             :    Input Parameters:
     575             : +  ds - the direct solver context
     576             : -  m  - the requested matrix
     577             : 
     578             :    Output Parameter:
     579             : .  a  - pointer to the values
     580             : 
     581             :    Note:
     582             :    To get read-only access, use DSGetMat() followed by MatDenseGetArrayRead().
     583             : 
     584             :    Level: advanced
     585             : 
     586             : .seealso: DSRestoreArray(), DSGetArrayReal()
     587             : @*/
     588       62860 : PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
     589             : {
     590       62860 :   PetscFunctionBegin;
     591       62860 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     592       62860 :   DSCheckAlloc(ds,1);
     593       62860 :   DSCheckValidMat(ds,m,2);
     594       62860 :   PetscAssertPointer(a,3);
     595       62860 :   PetscCall(MatDenseGetArray(ds->omat[m],a));
     596       62860 :   PetscFunctionReturn(PETSC_SUCCESS);
     597             : }
     598             : 
     599             : /*@C
     600             :    DSRestoreArray - Restores the matrix after DSGetArray() was called.
     601             : 
     602             :    Not Collective
     603             : 
     604             :    Input Parameters:
     605             : +  ds - the direct solver context
     606             : .  m  - the requested matrix
     607             : -  a  - pointer to the values
     608             : 
     609             :    Level: advanced
     610             : 
     611             : .seealso: DSGetArray(), DSGetArrayReal()
     612             : @*/
     613       62759 : PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
     614             : {
     615       62759 :   PetscFunctionBegin;
     616       62759 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     617       62759 :   DSCheckAlloc(ds,1);
     618       62759 :   DSCheckValidMat(ds,m,2);
     619       62759 :   PetscAssertPointer(a,3);
     620       62759 :   PetscCall(MatDenseRestoreArray(ds->omat[m],a));
     621       62759 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     622       62759 :   PetscFunctionReturn(PETSC_SUCCESS);
     623             : }
     624             : 
     625             : /*@C
     626             :    DSGetArrayReal - Returns a real pointer to the internal array of T or D.
     627             :    You MUST call DSRestoreArrayReal() when you no longer need to access the array.
     628             : 
     629             :    Not Collective
     630             : 
     631             :    Input Parameters:
     632             : +  ds - the direct solver context
     633             : -  m  - the requested matrix
     634             : 
     635             :    Output Parameter:
     636             : .  a  - pointer to the values
     637             : 
     638             :    Note:
     639             :    This function can be used only for DS_MAT_T and DS_MAT_D. These matrices always
     640             :    store real values, even in complex scalars, that is why the returned pointer is
     641             :    PetscReal.
     642             : 
     643             :    Level: advanced
     644             : 
     645             : .seealso: DSRestoreArrayReal(), DSGetArray()
     646             : @*/
     647      137795 : PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
     648             : {
     649             : #if defined(PETSC_USE_COMPLEX)
     650             :   PetscScalar *as;
     651             : #endif
     652             : 
     653      137795 :   PetscFunctionBegin;
     654      137795 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     655      137795 :   DSCheckAlloc(ds,1);
     656      137795 :   DSCheckValidMatReal(ds,m,2);
     657      137795 :   PetscAssertPointer(a,3);
     658             : #if defined(PETSC_USE_COMPLEX)
     659             :   PetscCall(MatDenseGetArray(ds->omat[m],&as));
     660             :   *a = (PetscReal*)as;
     661             : #else
     662      137795 :   PetscCall(MatDenseGetArray(ds->omat[m],a));
     663             : #endif
     664      137795 :   PetscFunctionReturn(PETSC_SUCCESS);
     665             : }
     666             : 
     667             : /*@C
     668             :    DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
     669             : 
     670             :    Not Collective
     671             : 
     672             :    Input Parameters:
     673             : +  ds - the direct solver context
     674             : .  m  - the requested matrix
     675             : -  a  - pointer to the values
     676             : 
     677             :    Level: advanced
     678             : 
     679             : .seealso: DSGetArrayReal(), DSGetArray()
     680             : @*/
     681      137795 : PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
     682             : {
     683             : #if defined(PETSC_USE_COMPLEX)
     684             :   PetscScalar *as;
     685             : #endif
     686             : 
     687      137795 :   PetscFunctionBegin;
     688      137795 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     689      137795 :   DSCheckAlloc(ds,1);
     690      137795 :   DSCheckValidMatReal(ds,m,2);
     691      137795 :   PetscAssertPointer(a,3);
     692             : #if defined(PETSC_USE_COMPLEX)
     693             :   PetscCall(MatDenseRestoreArray(ds->omat[m],&as));
     694             :   *a = NULL;
     695             : #else
     696      137795 :   PetscCall(MatDenseRestoreArray(ds->omat[m],a));
     697             : #endif
     698      137795 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     699      137795 :   PetscFunctionReturn(PETSC_SUCCESS);
     700             : }
     701             : 
     702             : /*@
     703             :    DSSolve - Solves the problem.
     704             : 
     705             :    Logically Collective
     706             : 
     707             :    Input Parameters:
     708             : +  ds   - the direct solver context
     709             : .  eigr - array to store the computed eigenvalues (real part)
     710             : -  eigi - array to store the computed eigenvalues (imaginary part)
     711             : 
     712             :    Note:
     713             :    This call brings the dense system to condensed form. No ordering
     714             :    of the eigenvalues is enforced (for this, call DSSort() afterwards).
     715             : 
     716             :    Level: intermediate
     717             : 
     718             : .seealso: DSSort(), DSStateType
     719             : @*/
     720       22423 : PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     721             : {
     722       22423 :   PetscFunctionBegin;
     723       22423 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     724       22423 :   PetscValidType(ds,1);
     725       22423 :   DSCheckAlloc(ds,1);
     726       22423 :   PetscAssertPointer(eigr,2);
     727       22423 :   if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(PETSC_SUCCESS);
     728       22398 :   PetscCheck(ds->ops->solve[ds->method],PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
     729       22398 :   PetscCall(PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
     730       22398 :   PetscCall(PetscLogEventBegin(DS_Solve,ds,0,0,0));
     731       22398 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     732       22398 :   PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
     733       22398 :   PetscCall(PetscFPTrapPop());
     734       22398 :   PetscCall(PetscLogEventEnd(DS_Solve,ds,0,0,0));
     735       22398 :   PetscCall(PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]));
     736       22398 :   ds->state = DS_STATE_CONDENSED;
     737       22398 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     738       22398 :   PetscFunctionReturn(PETSC_SUCCESS);
     739             : }
     740             : 
     741             : /*@C
     742             :    DSSort - Sorts the result of DSSolve() according to a given sorting
     743             :    criterion.
     744             : 
     745             :    Logically Collective
     746             : 
     747             :    Input Parameters:
     748             : +  ds   - the direct solver context
     749             : .  rr   - (optional) array containing auxiliary values (real part)
     750             : -  ri   - (optional) array containing auxiliary values (imaginary part)
     751             : 
     752             :    Input/Output Parameters:
     753             : +  eigr - array containing the computed eigenvalues (real part)
     754             : .  eigi - array containing the computed eigenvalues (imaginary part)
     755             : -  k    - (optional) number of elements in the leading group
     756             : 
     757             :    Notes:
     758             :    This routine sorts the arrays provided in eigr and eigi, and also
     759             :    sorts the dense system stored inside ds (assumed to be in condensed form).
     760             :    The sorting criterion is specified with DSSetSlepcSC().
     761             : 
     762             :    If arrays rr and ri are provided, then a (partial) reordering based on these
     763             :    values rather than on the eigenvalues is performed. In symmetric problems
     764             :    a total order is obtained (parameter k is ignored), but otherwise the result
     765             :    is sorted only partially. In this latter case, it is only guaranteed that
     766             :    all the first k elements satisfy the comparison with any of the last n-k
     767             :    elements. The output value of parameter k is the final number of elements in
     768             :    the first set.
     769             : 
     770             :    Level: intermediate
     771             : 
     772             : .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
     773             : @*/
     774       29819 : PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     775             : {
     776       29819 :   PetscInt       i;
     777             : 
     778       29819 :   PetscFunctionBegin;
     779       29819 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     780       29819 :   PetscValidType(ds,1);
     781       29819 :   DSCheckSolved(ds,1);
     782       29819 :   PetscAssertPointer(eigr,2);
     783       29819 :   if (rr) PetscAssertPointer(rr,4);
     784       29819 :   PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
     785       29819 :   PetscCheck(ds->sc,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
     786       29819 :   PetscCheck(!k || rr,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
     787             : 
     788      476754 :   for (i=0;i<ds->n;i++) ds->perm[i] = i;   /* initialize to trivial permutation */
     789       29819 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
     790       29819 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     791       29819 :   PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
     792       29819 :   PetscCall(PetscFPTrapPop());
     793       29819 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
     794       29819 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     795       29819 :   PetscCall(PetscInfo(ds,"Finished sorting\n"));
     796       29819 :   PetscFunctionReturn(PETSC_SUCCESS);
     797             : }
     798             : 
     799             : /*@C
     800             :    DSSortWithPermutation - Reorders the result of DSSolve() according to a given
     801             :    permutation.
     802             : 
     803             :    Logically Collective
     804             : 
     805             :    Input Parameters:
     806             : +  ds   - the direct solver context
     807             : -  perm - permutation that indicates the new ordering
     808             : 
     809             :    Input/Output Parameters:
     810             : +  eigr - array with the reordered eigenvalues (real part)
     811             : -  eigi - array with the reordered eigenvalues (imaginary part)
     812             : 
     813             :    Notes:
     814             :    This routine reorders the arrays provided in eigr and eigi, and also the dense
     815             :    system stored inside ds (assumed to be in condensed form). There is no sorting
     816             :    criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
     817             : 
     818             :    Level: advanced
     819             : 
     820             : .seealso: DSSolve(), DSSort()
     821             : @*/
     822           1 : PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)
     823             : {
     824           1 :   PetscFunctionBegin;
     825           1 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     826           1 :   PetscValidType(ds,1);
     827           1 :   DSCheckSolved(ds,1);
     828           1 :   PetscAssertPointer(perm,2);
     829           1 :   PetscAssertPointer(eigr,3);
     830           1 :   PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
     831             : 
     832           1 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
     833           1 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     834           1 :   PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
     835           1 :   PetscCall(PetscFPTrapPop());
     836           1 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
     837           1 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     838           1 :   PetscCall(PetscInfo(ds,"Finished sorting\n"));
     839           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     840             : }
     841             : 
     842             : /*@
     843             :    DSSynchronize - Make sure that all processes have the same data, performing
     844             :    communication if necessary.
     845             : 
     846             :    Collective
     847             : 
     848             :    Input Parameter:
     849             : .  ds   - the direct solver context
     850             : 
     851             :    Input/Output Parameters:
     852             : +  eigr - (optional) array with the computed eigenvalues (real part)
     853             : -  eigi - (optional) array with the computed eigenvalues (imaginary part)
     854             : 
     855             :    Notes:
     856             :    When the DS has been created with a communicator with more than one process,
     857             :    the internal data, especially the computed matrices, may diverge in the
     858             :    different processes. This happens when using multithreaded BLAS and may
     859             :    cause numerical issues in some ill-conditioned problems. This function
     860             :    performs the necessary communication among the processes so that the
     861             :    internal data is exactly equal in all of them.
     862             : 
     863             :    Depending on the parallel mode as set with DSSetParallel(), this function
     864             :    will either do nothing or synchronize the matrices computed by DSSolve()
     865             :    and DSSort(). The arguments eigr and eigi are typically those used in the
     866             :    calls to DSSolve() and DSSort().
     867             : 
     868             :    Level: developer
     869             : 
     870             : .seealso: DSSetParallel(), DSSolve(), DSSort()
     871             : @*/
     872       22171 : PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     873             : {
     874       22171 :   PetscMPIInt    size;
     875             : 
     876       22171 :   PetscFunctionBegin;
     877       22171 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     878       22171 :   PetscValidType(ds,1);
     879       22171 :   DSCheckAlloc(ds,1);
     880       22171 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
     881       22171 :   if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
     882         340 :     PetscCall(PetscLogEventBegin(DS_Synchronize,ds,0,0,0));
     883         340 :     PetscUseTypeMethod(ds,synchronize,eigr,eigi);
     884         340 :     PetscCall(PetscLogEventEnd(DS_Synchronize,ds,0,0,0));
     885         340 :     PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     886         340 :     PetscCall(PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]));
     887             :   }
     888       22171 :   PetscFunctionReturn(PETSC_SUCCESS);
     889             : }
     890             : 
     891             : /*@C
     892             :    DSVectors - Compute vectors associated to the dense system such
     893             :    as eigenvectors.
     894             : 
     895             :    Logically Collective
     896             : 
     897             :    Input Parameters:
     898             : +  ds  - the direct solver context
     899             : -  mat - the matrix, used to indicate which vectors are required
     900             : 
     901             :    Input/Output Parameter:
     902             : .  j   - (optional) index of vector to be computed
     903             : 
     904             :    Output Parameter:
     905             : .  rnorm - (optional) computed residual norm
     906             : 
     907             :    Notes:
     908             :    Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
     909             :    compute right or left eigenvectors, or left or right singular vectors,
     910             :    respectively.
     911             : 
     912             :    If NULL is passed in argument j then all vectors are computed,
     913             :    otherwise j indicates which vector must be computed. In real non-symmetric
     914             :    problems, on exit the index j will be incremented when a complex conjugate
     915             :    pair is found.
     916             : 
     917             :    This function can be invoked after the dense problem has been solved,
     918             :    to get the residual norm estimate of the associated Ritz pair. In that
     919             :    case, the relevant information is returned in rnorm.
     920             : 
     921             :    For computing eigenvectors, LAPACK's _trevc is used so the matrix must
     922             :    be in (quasi-)triangular form, or call DSSolve() first.
     923             : 
     924             :    Level: intermediate
     925             : 
     926             : .seealso: DSSolve()
     927             : @*/
     928       46428 : PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
     929             : {
     930       46428 :   PetscFunctionBegin;
     931       46428 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     932       46428 :   PetscValidType(ds,1);
     933       46428 :   DSCheckAlloc(ds,1);
     934      185712 :   PetscValidLogicalCollectiveEnum(ds,mat,2);
     935       46428 :   PetscCheck(mat<DS_NUM_MAT,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
     936       46428 :   PetscCheck(!rnorm || j,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
     937       46428 :   if (!ds->omat[mat]) PetscCall(DSAllocateMat_Private(ds,mat));
     938       46428 :   if (!j) PetscCall(PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]));
     939       46428 :   PetscCall(PetscLogEventBegin(DS_Vectors,ds,0,0,0));
     940       46428 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     941       46428 :   PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
     942       46428 :   PetscCall(PetscFPTrapPop());
     943       46428 :   PetscCall(PetscLogEventEnd(DS_Vectors,ds,0,0,0));
     944       46428 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
     945       46428 :   PetscFunctionReturn(PETSC_SUCCESS);
     946             : }
     947             : 
     948             : /*@
     949             :    DSUpdateExtraRow - Performs all necessary operations so that the extra
     950             :    row gets up-to-date after a call to DSSolve().
     951             : 
     952             :    Logically Collective
     953             : 
     954             :    Input Parameters:
     955             : .  ds - the direct solver context
     956             : 
     957             :    Level: advanced
     958             : 
     959             : .seealso: DSSolve(), DSSetExtraRow()
     960             : @*/
     961       10488 : PetscErrorCode DSUpdateExtraRow(DS ds)
     962             : {
     963       10488 :   PetscFunctionBegin;
     964       10488 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     965       10488 :   PetscValidType(ds,1);
     966       10488 :   DSCheckAlloc(ds,1);
     967       10488 :   PetscCheck(ds->extrarow,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
     968       10488 :   PetscCall(PetscInfo(ds,"Updating extra row\n"));
     969       10488 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
     970       10488 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     971       10488 :   PetscUseTypeMethod(ds,update);
     972       10488 :   PetscCall(PetscFPTrapPop());
     973       10488 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
     974       10488 :   PetscFunctionReturn(PETSC_SUCCESS);
     975             : }
     976             : 
     977             : /*@
     978             :    DSCond - Compute the condition number.
     979             : 
     980             :    Logically Collective
     981             : 
     982             :    Input Parameters:
     983             : +  ds - the direct solver context
     984             : -  cond - the computed condition number
     985             : 
     986             :    Notes:
     987             :    In standard eigenvalue problems, returns the inf-norm condition number of the first
     988             :    matrix, computed as cond(A) = norm(A)*norm(inv(A)).
     989             : 
     990             :    In GSVD problems, returns the maximum of cond(A) and cond(B), where cond(.) is
     991             :    computed as the ratio of the largest and smallest singular values.
     992             : 
     993             :    Does not take into account the extra row.
     994             : 
     995             :    Level: advanced
     996             : 
     997             : .seealso: DSSolve(), DSSetExtraRow()
     998             : @*/
     999         164 : PetscErrorCode DSCond(DS ds,PetscReal *cond)
    1000             : {
    1001         164 :   PetscFunctionBegin;
    1002         164 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
    1003         164 :   PetscValidType(ds,1);
    1004         164 :   DSCheckAlloc(ds,1);
    1005         164 :   PetscAssertPointer(cond,2);
    1006         164 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
    1007         164 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1008         164 :   PetscUseTypeMethod(ds,cond,cond);
    1009         164 :   PetscCall(PetscFPTrapPop());
    1010         164 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
    1011         164 :   PetscCall(PetscInfo(ds,"Computed condition number = %g\n",(double)*cond));
    1012         164 :   PetscFunctionReturn(PETSC_SUCCESS);
    1013             : }
    1014             : 
    1015             : /*@C
    1016             :    DSTranslateHarmonic - Computes a translation of the dense system.
    1017             : 
    1018             :    Logically Collective
    1019             : 
    1020             :    Input Parameters:
    1021             : +  ds      - the direct solver context
    1022             : .  tau     - the translation amount
    1023             : .  beta    - last component of vector b
    1024             : -  recover - boolean flag to indicate whether to recover or not
    1025             : 
    1026             :    Output Parameters:
    1027             : +  g       - the computed vector (optional)
    1028             : -  gamma   - scale factor (optional)
    1029             : 
    1030             :    Notes:
    1031             :    This function is intended for use in the context of Krylov methods only.
    1032             :    It computes a translation of a Krylov decomposition in order to extract
    1033             :    eigenpair approximations by harmonic Rayleigh-Ritz.
    1034             :    The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
    1035             :    vector b is assumed to be beta*e_n^T.
    1036             : 
    1037             :    The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
    1038             :    the factor by which the residual of the Krylov decomposition is scaled.
    1039             : 
    1040             :    If the recover flag is activated, the computed translation undoes the
    1041             :    translation done previously. In that case, parameter tau is ignored.
    1042             : 
    1043             :    Level: developer
    1044             : 
    1045             : .seealso: DSTranslateRKS()
    1046             : @*/
    1047          81 : PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
    1048             : {
    1049          81 :   PetscFunctionBegin;
    1050          81 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
    1051          81 :   PetscValidType(ds,1);
    1052          81 :   DSCheckAlloc(ds,1);
    1053          81 :   if (recover) PetscCall(PetscInfo(ds,"Undoing the translation\n"));
    1054          57 :   else PetscCall(PetscInfo(ds,"Computing the translation\n"));
    1055          81 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
    1056          81 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1057          81 :   PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
    1058          81 :   PetscCall(PetscFPTrapPop());
    1059          81 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
    1060          81 :   ds->state = DS_STATE_RAW;
    1061          81 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
    1062          81 :   PetscFunctionReturn(PETSC_SUCCESS);
    1063             : }
    1064             : 
    1065             : /*@
    1066             :    DSTranslateRKS - Computes a modification of the dense system corresponding
    1067             :    to an update of the shift in a rational Krylov method.
    1068             : 
    1069             :    Logically Collective
    1070             : 
    1071             :    Input Parameters:
    1072             : +  ds    - the direct solver context
    1073             : -  alpha - the translation amount
    1074             : 
    1075             :    Notes:
    1076             :    This function is intended for use in the context of Krylov methods only.
    1077             :    It takes the leading (k+1,k) submatrix of A, containing the truncated
    1078             :    Rayleigh quotient of a Krylov-Schur relation computed from a shift
    1079             :    sigma1 and transforms it to obtain a Krylov relation as if computed
    1080             :    from a different shift sigma2. The new matrix is computed as
    1081             :    1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
    1082             :    alpha = sigma1-sigma2.
    1083             : 
    1084             :    Matrix Q is placed in DS_MAT_Q so that it can be used to update the
    1085             :    Krylov basis.
    1086             : 
    1087             :    Level: developer
    1088             : 
    1089             : .seealso: DSTranslateHarmonic()
    1090             : @*/
    1091          12 : PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
    1092             : {
    1093          12 :   PetscFunctionBegin;
    1094          12 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
    1095          12 :   PetscValidType(ds,1);
    1096          12 :   DSCheckAlloc(ds,1);
    1097          12 :   PetscCall(PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha)));
    1098          12 :   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
    1099          12 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1100          12 :   PetscUseTypeMethod(ds,transrks,alpha);
    1101          12 :   PetscCall(PetscFPTrapPop());
    1102          12 :   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
    1103          12 :   ds->state   = DS_STATE_RAW;
    1104          12 :   ds->compact = PETSC_FALSE;
    1105          12 :   PetscCall(PetscObjectStateIncrease((PetscObject)ds));
    1106          12 :   PetscFunctionReturn(PETSC_SUCCESS);
    1107             : }

Generated by: LCOV version 1.14