Actual source code: dsops.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    DS operations: DSSolve(), DSVectors(), etc
 12: */

 14: #include <slepc/private/dsimpl.h>      /*I "slepcds.h" I*/

 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: {
 37:   *ld = ds->ld;
 38:   return(0);
 39: }

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

 44:    Logically Collective on DS

 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: {

 70:   switch (state) {
 71:     case DS_STATE_RAW:
 72:     case DS_STATE_INTERMEDIATE:
 73:     case DS_STATE_CONDENSED:
 74:     case DS_STATE_TRUNCATED:
 75:       if (ds->state<state) { PetscInfo(ds,"DS state has been increased\n"); }
 76:       ds->state = state;
 77:       break;
 78:     default:
 79:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
 80:   }
 81:   return(0);
 82: }

 84: /*@
 85:    DSGetState - Returns the current state.

 87:    Not Collective

 89:    Input Parameter:
 90: .  ds - the direct solver context

 92:    Output Parameter:
 93: .  state - current state

 95:    Level: advanced

 97: .seealso: DSSetState()
 98: @*/
 99: PetscErrorCode DSGetState(DS ds,DSStateType *state)
100: {
104:   *state = ds->state;
105:   return(0);
106: }

108: /*@
109:    DSSetDimensions - Resize the matrices in the DS object.

111:    Logically Collective on DS

113:    Input Parameters:
114: +  ds - the direct solver context
115: .  n  - the new size
116: .  m  - the new column size (only for DSSVD)
117: .  l  - number of locked (inactive) leading columns
118: -  k  - intermediate dimension (e.g., position of arrow)

120:    Notes:
121:    The internal arrays are not reallocated.

123:    The value m is not used except in the case of DSSVD, pass 0 otherwise.

125:    Level: intermediate

127: .seealso: DSGetDimensions(), DSAllocate()
128: @*/
129: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt m,PetscInt l,PetscInt k)
130: {
133:   DSCheckAlloc(ds,1);
138:   if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
139:     ds->n = ds->ld;
140:   } else {
141:     if (n<0 || n>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
142:     if (ds->extrarow && n+1>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"A value of n equal to ld leaves no room for extra row");
143:     ds->n = n;
144:   }
145:   ds->t = ds->n;   /* truncated length equal to the new dimension */
146:   if (m) {
147:     if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
148:       ds->m = ds->ld;
149:     } else {
150:       if (m<0 || m>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 0 and ld");
151:       ds->m = m;
152:     }
153:   }
154:   if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
155:     ds->l = 0;
156:   } else {
157:     if (l<0 || l>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
158:     ds->l = l;
159:   }
160:   if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
161:     ds->k = ds->n/2;
162:   } else {
163:     if (k<0 || k>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
164:     ds->k = k;
165:   }
166:   return(0);
167: }

169: /*@
170:    DSGetDimensions - Returns the current dimensions.

172:    Not Collective

174:    Input Parameter:
175: .  ds - the direct solver context

177:    Output Parameter:
178: +  n  - the current size
179: .  m  - the current column size (only for DSSVD)
180: .  l  - number of locked (inactive) leading columns
181: .  k  - intermediate dimension (e.g., position of arrow)
182: -  t  - truncated length

184:    Note:
185:    The t parameter makes sense only if DSTruncate() has been called.
186:    Otherwise its value equals n.

188:    Level: intermediate

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

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

208:    Logically Collective on DS

210:    Input Parameters:
211: +  ds - the direct solver context
212: -  n  - the new size

214:    Note:
215:    The new size is set to n. In cases where the extra row is meaningful,
216:    the first n elements are kept as the extra row for the new system.

218:    Level: advanced

220: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
221: @*/
222: PetscErrorCode DSTruncate(DS ds,PetscInt n)
223: {

229:   DSCheckAlloc(ds,1);
230:   DSCheckSolved(ds,1);
232:   if (!ds->ops->truncate) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
233:   if (n<ds->l || n>ds->n) SETERRQ3(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n (%D). Must be between l (%D) and n (%D)",n,ds->l,ds->n);
234:   PetscLogEventBegin(DS_Other,ds,0,0,0);
235:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
236:   (*ds->ops->truncate)(ds,n);
237:   PetscFPTrapPop();
238:   PetscLogEventEnd(DS_Other,ds,0,0,0);
239:   ds->state = DS_STATE_TRUNCATED;
240:   PetscObjectStateIncrease((PetscObject)ds);
241:   return(0);
242: }

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

247:    Not Collective

249:    Input Parameters:
250: +  ds - the direct solver context
251: -  t  - the requested matrix

253:    Output Parameters:
254: +  n  - the number of rows
255: -  m  - the number of columns

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

260:    Level: developer

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

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

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

288:    Not Collective

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

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

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

301:    Level: developer

303: .seealso: DSGetMat()
304: @*/
305: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
306: {

312:   DSCheckValidMat(ds,t,2);
314:   if (ds->ops->hermitian) {
315:     (*ds->ops->hermitian)(ds,t,flg);
316:   } else *flg = PETSC_FALSE;
317:   return(0);
318: }

320: /*@
321:    DSGetMat - Returns a sequential dense Mat object containing the requested
322:    matrix.

324:    Not Collective

326:    Input Parameters:
327: +  ds - the direct solver context
328: -  m  - the requested matrix

330:    Output Parameter:
331: .  A  - Mat object

333:    Notes:
334:    The Mat is created with sizes equal to the current DS dimensions (nxm),
335:    then it is filled with the values that would be obtained with DSGetArray()
336:    (not DSGetArrayReal()). If the DS was truncated, then the number of rows
337:    is equal to the dimension prior to truncation, see DSTruncate().
338:    The communicator is always PETSC_COMM_SELF.

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

343:    Level: advanced

345: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
346: @*/
347: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
348: {
350:   PetscInt       j,rows,cols,arows,acols;
351:   PetscBool      create=PETSC_FALSE,flg;
352:   PetscScalar    *pA,*M;

356:   DSCheckAlloc(ds,1);
357:   DSCheckValidMat(ds,m,2);
359:   if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");

361:   DSMatGetSize(ds,m,&rows,&cols);
362:   if (!ds->omat[m]) create=PETSC_TRUE;
363:   else {
364:     MatGetSize(ds->omat[m],&arows,&acols);
365:     if (arows!=rows || acols!=cols) {
366:       MatDestroy(&ds->omat[m]);
367:       create=PETSC_TRUE;
368:     }
369:   }
370:   if (create) {
371:     MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);
372:   }

374:   /* set Hermitian flag */
375:   DSMatIsHermitian(ds,m,&flg);
376:   MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);

378:   /* copy entries */
379:   PetscObjectReference((PetscObject)ds->omat[m]);
380:   *A = ds->omat[m];
381:   M  = ds->mat[m];
382:   MatDenseGetArray(*A,&pA);
383:   for (j=0;j<cols;j++) {
384:     PetscMemcpy(pA+j*rows,M+j*ds->ld,rows*sizeof(PetscScalar));
385:   }
386:   MatDenseRestoreArray(*A,&pA);
387:   return(0);
388: }

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

393:    Not Collective

395:    Input Parameters:
396: +  ds - the direct solver context
397: .  m  - the requested matrix
398: -  A  - the fetched Mat object

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

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

408:    Level: advanced

410: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
411: @*/
412: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
413: {
415:   PetscInt       j,rows,cols;
416:   PetscScalar    *pA,*M;

420:   DSCheckAlloc(ds,1);
421:   DSCheckValidMat(ds,m,2);
423:   if (!ds->omat[m]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSRestoreMat must match a previous call to DSGetMat");
424:   if (ds->omat[m]!=*A) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with DSGetMat");

426:   MatGetSize(*A,&rows,&cols);
427:   M  = ds->mat[m];
428:   MatDenseGetArray(*A,&pA);
429:   for (j=0;j<cols;j++) {
430:     PetscMemcpy(M+j*ds->ld,pA+j*rows,rows*sizeof(PetscScalar));
431:   }
432:   MatDenseRestoreArray(*A,&pA);
433:   MatDestroy(A);
434:   return(0);
435: }

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

442:    Not Collective

444:    Input Parameters:
445: +  ds - the direct solver context
446: -  m  - the requested matrix

448:    Output Parameter:
449: .  a  - pointer to the values

451:    Level: advanced

453: .seealso: DSRestoreArray(), DSGetArrayReal()
454: @*/
455: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
456: {
459:   DSCheckAlloc(ds,1);
460:   DSCheckValidMat(ds,m,2);
462:   *a = ds->mat[m];
463:   CHKMEMQ;
464:   return(0);
465: }

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

470:    Not Collective

472:    Input Parameters:
473: +  ds - the direct solver context
474: .  m  - the requested matrix
475: -  a  - pointer to the values

477:    Level: advanced

479: .seealso: DSGetArray(), DSGetArrayReal()
480: @*/
481: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
482: {

487:   DSCheckAlloc(ds,1);
488:   DSCheckValidMat(ds,m,2);
490:   CHKMEMQ;
491:   *a = 0;
492:   PetscObjectStateIncrease((PetscObject)ds);
493:   return(0);
494: }

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

501:    Not Collective

503:    Input Parameters:
504: +  ds - the direct solver context
505: -  m  - the requested matrix

507:    Output Parameter:
508: .  a  - pointer to the values

510:    Level: advanced

512: .seealso: DSRestoreArrayReal(), DSGetArray()
513: @*/
514: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
515: {
518:   DSCheckAlloc(ds,1);
519:   DSCheckValidMatReal(ds,m,2);
521:   *a = ds->rmat[m];
522:   CHKMEMQ;
523:   return(0);
524: }

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

529:    Not Collective

531:    Input Parameters:
532: +  ds - the direct solver context
533: .  m  - the requested matrix
534: -  a  - pointer to the values

536:    Level: advanced

538: .seealso: DSGetArrayReal(), DSGetArray()
539: @*/
540: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
541: {

546:   DSCheckAlloc(ds,1);
547:   DSCheckValidMatReal(ds,m,2);
549:   CHKMEMQ;
550:   *a = 0;
551:   PetscObjectStateIncrease((PetscObject)ds);
552:   return(0);
553: }

555: /*@
556:    DSSolve - Solves the problem.

558:    Logically Collective on DS

560:    Input Parameters:
561: +  ds   - the direct solver context
562: .  eigr - array to store the computed eigenvalues (real part)
563: -  eigi - array to store the computed eigenvalues (imaginary part)

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

569:    Level: intermediate

571: .seealso: DSSort(), DSStateType
572: @*/
573: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
574: {

580:   DSCheckAlloc(ds,1);
582:   if (ds->state>=DS_STATE_CONDENSED) return(0);
583:   if (!ds->ops->solve[ds->method]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
584:   PetscLogEventBegin(DS_Solve,ds,0,0,0);
585:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
586:   (*ds->ops->solve[ds->method])(ds,eigr,eigi);
587:   PetscFPTrapPop();
588:   PetscLogEventEnd(DS_Solve,ds,0,0,0);
589:   ds->state = DS_STATE_CONDENSED;
590:   PetscObjectStateIncrease((PetscObject)ds);
591:   return(0);
592: }

594: /*@C
595:    DSSort - Sorts the result of DSSolve() according to a given sorting
596:    criterion.

598:    Logically Collective on DS

600:    Input Parameters:
601: +  ds   - the direct solver context
602: .  rr   - (optional) array containing auxiliary values (real part)
603: -  ri   - (optional) array containing auxiliary values (imaginary part)

605:    Input/Output Parameters:
606: +  eigr - array containing the computed eigenvalues (real part)
607: .  eigi - array containing the computed eigenvalues (imaginary part)
608: -  k    - (optional) number of elements in the leading group

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

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

623:    Level: intermediate

625: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
626: @*/
627: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
628: {
630:   PetscInt       i;

635:   DSCheckSolved(ds,1);
638:   if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
639:   if (!ds->ops->sort) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
640:   if (!ds->sc) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
641:   if (k && !rr) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");

643:   for (i=0;i<ds->n;i++) ds->perm[i] = i;   /* initialize to trivial permutation */
644:   PetscLogEventBegin(DS_Other,ds,0,0,0);
645:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
646:   (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
647:   PetscFPTrapPop();
648:   PetscLogEventEnd(DS_Other,ds,0,0,0);
649:   PetscObjectStateIncrease((PetscObject)ds);
650:   return(0);
651: }

653: /*@C
654:    DSSortWithPermutation - Reorders the result of DSSolve() according to a given
655:    permutation.

657:    Logically Collective on DS

659:    Input Parameters:
660: +  ds   - the direct solver context
661: -  perm - permutation that indicates the new ordering

663:    Input/Output Parameters:
664: +  eigr - array with the reordered eigenvalues (real part)
665: -  eigi - array with the reordered eigenvalues (imaginary part)

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

672:    Level: advanced

674: .seealso: DSSolve(), DSSort()
675: @*/
676: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)
677: {

683:   DSCheckSolved(ds,1);
686:   if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
687:   if (!ds->ops->sortperm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);

689:   PetscLogEventBegin(DS_Other,ds,0,0,0);
690:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
691:   (*ds->ops->sortperm)(ds,perm,eigr,eigi);
692:   PetscFPTrapPop();
693:   PetscLogEventEnd(DS_Other,ds,0,0,0);
694:   PetscObjectStateIncrease((PetscObject)ds);
695:   return(0);
696: }

698: /*@
699:    DSSynchronize - Make sure that all processes have the same data, performing
700:    communication if necessary.

702:    Collective on DS

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

707:    Input/Output Parameters:
708: +  eigr - (optional) array with the computed eigenvalues (real part)
709: -  eigi - (optional) array with the computed eigenvalues (imaginary part)

711:    Notes:
712:    When the DS has been created with a communicator with more than one process,
713:    the internal data, especially the computed matrices, may diverge in the
714:    different processes. This happens when using multithreaded BLAS and may
715:    cause numerical issues in some ill-conditioned problems. This function
716:    performs the necessary communication among the processes so that the
717:    internal data is exactly equal in all of them.

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

724:    Level: developer

726: .seealso: DSSetParallel(), DSSolve(), DSSort()
727: @*/
728: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
729: {
731:   PetscMPIInt    size;

736:   DSCheckAlloc(ds,1);
737:   MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
738:   if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
739:     PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
740:     if (ds->ops->synchronize) {
741:       (*ds->ops->synchronize)(ds,eigr,eigi);
742:     }
743:     PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
744:     PetscObjectStateIncrease((PetscObject)ds);
745:   }
746:   return(0);
747: }

749: /*@C
750:    DSVectors - Compute vectors associated to the dense system such
751:    as eigenvectors.

753:    Logically Collective on DS

755:    Input Parameters:
756: +  ds  - the direct solver context
757: -  mat - the matrix, used to indicate which vectors are required

759:    Input/Output Parameter:
760: -  j   - (optional) index of vector to be computed

762:    Output Parameter:
763: .  rnorm - (optional) computed residual norm

765:    Notes:
766:    Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_VT, to
767:    compute right or left eigenvectors, or left or right singular vectors,
768:    respectively.

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

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

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

782:    Level: intermediate

784: .seealso: DSSolve()
785: @*/
786: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
787: {

793:   DSCheckAlloc(ds,1);
795:   if (mat>=DS_NUM_MAT) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
796:   if (!ds->ops->vectors) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
797:   if (rnorm && !j) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
798:   if (!ds->mat[mat]) { DSAllocateMat_Private(ds,mat); }
799:   PetscLogEventBegin(DS_Vectors,ds,0,0,0);
800:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
801:   (*ds->ops->vectors)(ds,mat,j,rnorm);
802:   PetscFPTrapPop();
803:   PetscLogEventEnd(DS_Vectors,ds,0,0,0);
804:   PetscObjectStateIncrease((PetscObject)ds);
805:   return(0);
806: }

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

812:    Not Collective

814:    Input Parameters:
815: .  ds - the direct solver context

817:    Level: advanced

819: .seealso: DSSolve(), DSSetExtraRow()
820: @*/
821: PetscErrorCode DSUpdateExtraRow(DS ds)
822: {

828:   DSCheckAlloc(ds,1);
829:   if (!ds->ops->update) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
830:   if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
831:   PetscLogEventBegin(DS_Other,ds,0,0,0);
832:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
833:   (*ds->ops->update)(ds);
834:   PetscFPTrapPop();
835:   PetscLogEventEnd(DS_Other,ds,0,0,0);
836:   return(0);
837: }

839: /*@
840:    DSCond - Compute the inf-norm condition number of the first matrix
841:    as cond(A) = norm(A)*norm(inv(A)).

843:    Not Collective

845:    Input Parameters:
846: +  ds - the direct solver context
847: -  cond - the computed condition number

849:    Level: advanced

851: .seealso: DSSolve()
852: @*/
853: PetscErrorCode DSCond(DS ds,PetscReal *cond)
854: {

860:   DSCheckAlloc(ds,1);
862:   if (!ds->ops->cond) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
863:   PetscLogEventBegin(DS_Other,ds,0,0,0);
864:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
865:   (*ds->ops->cond)(ds,cond);
866:   PetscFPTrapPop();
867:   PetscLogEventEnd(DS_Other,ds,0,0,0);
868:   return(0);
869: }

871: /*@C
872:    DSTranslateHarmonic - Computes a translation of the dense system.

874:    Logically Collective on DS

876:    Input Parameters:
877: +  ds      - the direct solver context
878: .  tau     - the translation amount
879: .  beta    - last component of vector b
880: -  recover - boolean flag to indicate whether to recover or not

882:    Output Parameters:
883: +  g       - the computed vector (optional)
884: -  gamma   - scale factor (optional)

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

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

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

899:    Level: developer
900: @*/
901: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
902: {

908:   DSCheckAlloc(ds,1);
909:   if (!ds->ops->transharm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
910:   PetscLogEventBegin(DS_Other,ds,0,0,0);
911:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
912:   (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
913:   PetscFPTrapPop();
914:   PetscLogEventEnd(DS_Other,ds,0,0,0);
915:   ds->state = DS_STATE_RAW;
916:   PetscObjectStateIncrease((PetscObject)ds);
917:   return(0);
918: }

920: /*@
921:    DSTranslateRKS - Computes a modification of the dense system corresponding
922:    to an update of the shift in a rational Krylov method.

924:    Logically Collective on DS

926:    Input Parameters:
927: +  ds    - the direct solver context
928: -  alpha - the translation amount

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

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

942:    Level: developer
943: @*/
944: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
945: {

951:   DSCheckAlloc(ds,1);
952:   if (!ds->ops->transrks) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
953:   PetscLogEventBegin(DS_Other,ds,0,0,0);
954:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
955:   (*ds->ops->transrks)(ds,alpha);
956:   PetscFPTrapPop();
957:   PetscLogEventEnd(DS_Other,ds,0,0,0);
958:   ds->state   = DS_STATE_RAW;
959:   ds->compact = PETSC_FALSE;
960:   PetscObjectStateIncrease((PetscObject)ds);
961:   return(0);
962: }

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

968:    Not Collective

970:    Input Parameters:
971: +  ds   - the direct solver context
972: .  m    - the requested matrix
973: .  mr   - first row of m to be considered
974: .  mc   - first column of m to be considered
975: .  A    - Mat object
976: .  Ar   - first row of A to be considered
977: .  Ac   - first column of A to be considered
978: .  rows - number of rows to copy
979: .  cols - number of columns to copy
980: -  out  - whether the data is copied out of the DS

982:    Note:
983:    If out=true, the values of the DS matrix m are copied to A, otherwise
984:    the entries of A are copied to the DS.

986:    Level: developer

988: .seealso: DSGetMat()
989: @*/
990: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)
991: {
993:   PetscInt       j,mrows,mcols,arows,acols;
994:   PetscScalar    *pA,*M;

998:   DSCheckAlloc(ds,1);
1000:   DSCheckValidMat(ds,m,2);
1009:   if (!rows || !cols) return(0);

1011:   DSMatGetSize(ds,m,&mrows,&mcols);
1012:   MatGetSize(A,&arows,&acols);
1013:   if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
1014:   if (mr<0 || mr>=mrows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in m");
1015:   if (mc<0 || mc>=mcols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in m");
1016:   if (Ar<0 || Ar>=arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in A");
1017:   if (Ac<0 || Ac>=acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in A");
1018:   if (mr+rows>mrows || Ar+rows>arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of rows");
1019:   if (mc+cols>mcols || Ac+cols>acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of columns");

1021:   M  = ds->mat[m];
1022:   MatDenseGetArray(A,&pA);
1023:   for (j=0;j<cols;j++) {
1024:     if (out) {
1025:       PetscMemcpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows*sizeof(PetscScalar));
1026:     } else {
1027:       PetscMemcpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows*sizeof(PetscScalar));
1028:     }
1029:   }
1030:   MatDenseRestoreArray(A,&pA);
1031:   return(0);
1032: }