Actual source code: dsops.c
slepc-3.18.2 2023-01-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: DS operations: DSSolve(), DSVectors(), etc
12: */
14: #include <slepc/private/dsimpl.h>
16: /*@
17: DSGetLeadingDimension - Returns the leading dimension of the allocated
18: matrices.
20: Not Collective
22: Input Parameter:
23: . ds - the direct solver context
25: Output Parameter:
26: . ld - leading dimension (maximum allowed dimension for the matrices)
28: Level: advanced
30: .seealso: DSAllocate(), DSSetDimensions()
31: @*/
32: PetscErrorCode DSGetLeadingDimension(DS ds,PetscInt *ld)
33: {
36: *ld = ds->ld;
37: return 0;
38: }
40: /*@
41: DSSetState - Change the state of the DS object.
43: Logically Collective on ds
45: Input Parameters:
46: + ds - the direct solver context
47: - state - the new state
49: Notes:
50: The state indicates that the dense system is in an initial state (raw),
51: in an intermediate state (such as tridiagonal, Hessenberg or
52: Hessenberg-triangular), in a condensed state (such as diagonal, Schur or
53: generalized Schur), or in a truncated state.
55: This function is normally used to return to the raw state when the
56: condensed structure is destroyed.
58: Level: advanced
60: .seealso: DSGetState()
61: @*/
62: PetscErrorCode DSSetState(DS ds,DSStateType state)
63: {
66: switch (state) {
67: case DS_STATE_RAW:
68: case DS_STATE_INTERMEDIATE:
69: case DS_STATE_CONDENSED:
70: case DS_STATE_TRUNCATED:
71: if (ds->state!=state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]);
72: ds->state = state;
73: break;
74: default:
75: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
76: }
77: return 0;
78: }
80: /*@
81: DSGetState - Returns the current state.
83: Not Collective
85: Input Parameter:
86: . ds - the direct solver context
88: Output Parameter:
89: . state - current state
91: Level: advanced
93: .seealso: DSSetState()
94: @*/
95: PetscErrorCode DSGetState(DS ds,DSStateType *state)
96: {
99: *state = ds->state;
100: return 0;
101: }
103: /*@
104: DSSetDimensions - Resize the matrices in the DS object.
106: Logically Collective on ds
108: Input Parameters:
109: + ds - the direct solver context
110: . n - the new size
111: . l - number of locked (inactive) leading columns
112: - k - intermediate dimension (e.g., position of arrow)
114: Notes:
115: The internal arrays are not reallocated.
117: Some DS types have additional dimensions, e.g. the number of columns
118: in DSSVD. For these, you should call a specific interface function.
120: Level: intermediate
122: .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
123: @*/
124: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
125: {
126: PetscInt on,ol,ok;
127: PetscBool issvd;
130: DSCheckAlloc(ds,1);
134: on = ds->n; ol = ds->l; ok = ds->k;
135: if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
136: ds->n = ds->extrarow? ds->ld-1: ds->ld;
137: } else {
139: PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,""); /* SVD and GSVD have extra column instead of extra row */
141: ds->n = n;
142: }
143: ds->t = ds->n; /* truncated length equal to the new dimension */
144: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
145: ds->l = 0;
146: } else {
148: ds->l = l;
149: }
150: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
151: ds->k = ds->n/2;
152: } else {
154: ds->k = k;
155: }
156: if (on!=ds->n || ol!=ds->l || ok!=ds->k) PetscInfo(ds,"New dimensions are: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
157: return 0;
158: }
160: /*@
161: DSGetDimensions - Returns the current dimensions.
163: Not Collective
165: Input Parameter:
166: . ds - the direct solver context
168: Output Parameters:
169: + n - the current size
170: . l - number of locked (inactive) leading columns
171: . k - intermediate dimension (e.g., position of arrow)
172: - t - truncated length
174: Note:
175: The t parameter makes sense only if DSTruncate() has been called.
176: Otherwise its value equals n.
178: Some DS types have additional dimensions, e.g. the number of columns
179: in DSSVD. For these, you should call a specific interface function.
181: Level: intermediate
183: .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
184: @*/
185: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
186: {
188: DSCheckAlloc(ds,1);
189: if (n) *n = ds->n;
190: if (l) *l = ds->l;
191: if (k) *k = ds->k;
192: if (t) *t = ds->t;
193: return 0;
194: }
196: /*@
197: DSTruncate - Truncates the system represented in the DS object.
199: Logically Collective on ds
201: Input Parameters:
202: + ds - the direct solver context
203: . n - the new size
204: - trim - a flag to indicate if the factorization must be trimmed
206: Note:
207: The new size is set to n. Note that in some cases the new size could
208: be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
209: Schur form). In cases where the extra row is meaningful, the first
210: n elements are kept as the extra row for the new system.
212: If the flag trim is turned on, it resets the locked and intermediate
213: dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
214: It also cleans the extra row if being used.
216: The typical usage of trim=true is to truncate the Schur decomposition
217: at the end of a Krylov iteration. In this case, the state must be
218: changed to RAW so that DSVectors() computes eigenvectors from scratch.
220: Level: advanced
222: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
223: @*/
224: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
225: {
226: DSStateType old;
230: DSCheckAlloc(ds,1);
234: PetscLogEventBegin(DS_Other,ds,0,0,0);
235: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
236: PetscUseTypeMethod(ds,truncate,n,trim);
237: PetscFPTrapPop();
238: PetscLogEventEnd(DS_Other,ds,0,0,0);
239: PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n);
240: old = ds->state;
241: ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
242: if (old!=ds->state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]);
243: PetscObjectStateIncrease((PetscObject)ds);
244: return 0;
245: }
247: /*@
248: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
250: Not Collective
252: Input Parameters:
253: + ds - the direct solver context
254: - t - the requested matrix
256: Output Parameters:
257: + n - the number of rows
258: - m - the number of columns
260: Note:
261: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
263: Level: developer
265: .seealso: DSSetDimensions(), DSGetMat()
266: @*/
267: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
268: {
269: PetscInt rows,cols;
273: DSCheckValidMat(ds,t,2);
274: if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
275: else {
276: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
277: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
278: if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
279: else if (t==DS_MAT_D) cols = 1;
280: else cols = ds->n;
281: }
282: if (m) *m = rows;
283: if (n) *n = cols;
284: return 0;
285: }
287: /*@
288: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
290: Not Collective
292: Input Parameters:
293: + ds - the direct solver context
294: - t - the requested matrix
296: Output Parameter:
297: . flg - the Hermitian flag
299: Note:
300: Does not check the matrix values directly. The flag is set according to the
301: problem structure. For instance, in DSHEP matrix A is Hermitian.
303: Level: developer
305: .seealso: DSGetMat()
306: @*/
307: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
308: {
311: DSCheckValidMat(ds,t,2);
313: *flg = PETSC_FALSE;
314: PetscTryTypeMethod(ds,hermitian,t,flg);
315: return 0;
316: }
318: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
319: {
320: #if !defined(PETSC_USE_COMPLEX)
321: PetscScalar val;
322: #endif
324: #if !defined(PETSC_USE_COMPLEX)
325: MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val);
326: if (val != 0.0) {
327: if (l+(*k)<n-1) (*k)++;
328: else (*k)--;
329: }
330: #endif
331: return 0;
332: }
334: /*@
335: DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
336: to avoid breaking a 2x2 block.
338: Not Collective
340: Input Parameters:
341: + ds - the direct solver context
342: . l - the size of the locked part (set to 0 to use ds->l)
343: . n - the total matrix size (set to 0 to use ds->n)
344: - k - the wanted truncation size
346: Output Parameter:
347: . k - the possibly modified value of the truncation size
349: Notes:
350: This should be called before DSTruncate() to make sure that the truncation
351: does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
353: The total size is n (either user-provided or ds->n if 0 is passed). The
354: size where the truncation is intended is equal to l+k (where l can be
355: equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
356: at the l+k limit, the value of k is increased (or decreased) by 1.
358: Level: advanced
360: .seealso: DSTruncate(), DSSetDimensions()
361: @*/
362: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
363: {
365: DSCheckAlloc(ds,1);
369: PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
370: return 0;
371: }
373: /*@
374: DSGetMat - Returns a sequential dense Mat object containing the requested
375: matrix.
377: Not Collective
379: Input Parameters:
380: + ds - the direct solver context
381: - m - the requested matrix
383: Output Parameter:
384: . A - Mat object
386: Notes:
387: The returned Mat has sizes equal to the current DS dimensions (nxm),
388: and contains the values that would be obtained with DSGetArray()
389: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
390: is equal to the dimension prior to truncation, see DSTruncate().
391: The communicator is always PETSC_COMM_SELF.
393: It is implemented with MatDenseGetSubMatrix(), and when no longer needed
394: the user must call DSRestoreMat() which will invoke MatDenseRestoreSubMatrix().
396: For matrices DS_MAT_T and DS_MAT_D, this function will return a Mat object
397: that cannot be used directly for computations, since it uses compact storage
398: (three and one diagonals for T and D, respectively). In complex scalars, the
399: internal array stores real values, so it is sufficient with 2 columns for T.
401: Level: advanced
403: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
404: @*/
405: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
406: {
407: PetscInt rows,cols;
408: PetscBool flg;
411: DSCheckAlloc(ds,1);
412: DSCheckValidMat(ds,m,2);
415: DSMatGetSize(ds,m,&rows,&cols);
417: MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A);
419: /* set Hermitian flag */
420: DSMatIsHermitian(ds,m,&flg);
421: MatSetOption(*A,MAT_SYMMETRIC,flg);
422: MatSetOption(*A,MAT_HERMITIAN,flg);
423: return 0;
424: }
426: /*@
427: DSRestoreMat - Restores the matrix after DSGetMat() was called.
429: Not Collective
431: Input Parameters:
432: + ds - the direct solver context
433: . m - the requested matrix
434: - A - the fetched Mat object
436: Note:
437: A call to this function must match a previous call of DSGetMat().
439: Level: advanced
441: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
442: @*/
443: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
444: {
446: DSCheckAlloc(ds,1);
447: DSCheckValidMat(ds,m,2);
450: MatDenseRestoreSubMatrix(ds->omat[m],A);
451: PetscObjectStateIncrease((PetscObject)ds);
452: return 0;
453: }
455: /*@
456: DSGetMatAndColumn - Returns a sequential dense Mat object containing the requested
457: matrix and one of its columns as a Vec.
459: Not Collective
461: Input Parameters:
462: + ds - the direct solver context
463: . m - the requested matrix
464: - col - the requested column
466: Output Parameters:
467: + A - Mat object
468: - v - Vec object (the column)
470: Notes:
471: This calls DSGetMat() and then it extracts the selected column.
472: The user must call DSRestoreMatAndColumn() to recover the original state.
473: For matrices DS_MAT_T and DS_MAT_D, in complex scalars this function implies
474: copying from real values stored internally to scalar values in the Vec.
476: Level: advanced
478: .seealso: DSRestoreMatAndColumn(), DSGetMat()
479: @*/
480: PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
481: {
483: DSCheckAlloc(ds,1);
484: DSCheckValidMat(ds,m,2);
488: DSGetMat(ds,m,A);
489: if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
490: const PetscScalar *as;
491: PetscScalar *vs;
492: PetscReal *ar;
493: PetscInt i,n,lda;
494: MatCreateVecs(*A,NULL,v);
495: VecGetSize(*v,&n);
496: MatDenseGetLDA(*A,&lda);
497: MatDenseGetArrayRead(*A,&as);
498: VecGetArrayWrite(*v,&vs);
499: ar = (PetscReal*)as;
500: for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
501: VecRestoreArrayWrite(*v,&vs);
502: MatDenseRestoreArrayRead(*A,&as);
503: } else MatDenseGetColumnVec(*A,col,v);
504: return 0;
505: }
507: /*@
508: DSRestoreMatAndColumn - Restores the matrix and vector after DSGetMatAndColumn()
509: was called.
511: Not Collective
513: Input Parameters:
514: + ds - the direct solver context
515: . m - the requested matrix
516: . col - the requested column
517: . A - the fetched Mat object
518: - v - the fetched Vec object
520: Note:
521: A call to this function must match a previous call of DSGetMatAndColumn().
523: Level: advanced
525: .seealso: DSGetMatAndColumn(), DSRestoreMat()
526: @*/
527: PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
528: {
530: DSCheckAlloc(ds,1);
531: DSCheckValidMat(ds,m,2);
535: if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
536: const PetscScalar *vs;
537: PetscScalar *as;
538: PetscReal *ar;
539: PetscInt i,n,lda;
540: VecGetSize(*v,&n);
541: MatDenseGetLDA(*A,&lda);
542: MatDenseGetArray(*A,&as);
543: VecGetArrayRead(*v,&vs);
544: ar = (PetscReal*)as;
545: for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
546: VecRestoreArrayRead(*v,&vs);
547: MatDenseRestoreArray(*A,&as);
548: VecDestroy(v);
549: } else MatDenseRestoreColumnVec(*A,col,v);
550: DSRestoreMat(ds,m,A);
551: return 0;
552: }
554: /*@C
555: DSGetArray - Returns a pointer to the internal array of one of the
556: matrices. You MUST call DSRestoreArray() when you no longer
557: need to access the array.
559: Not Collective
561: Input Parameters:
562: + ds - the direct solver context
563: - m - the requested matrix
565: Output Parameter:
566: . a - pointer to the values
568: Note:
569: To get read-only access, use DSGetMat() followed by MatDenseGetArrayRead().
571: Level: advanced
573: .seealso: DSRestoreArray(), DSGetArrayReal()
574: @*/
575: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
576: {
578: DSCheckAlloc(ds,1);
579: DSCheckValidMat(ds,m,2);
581: MatDenseGetArray(ds->omat[m],a);
582: return 0;
583: }
585: /*@C
586: DSRestoreArray - Restores the matrix after DSGetArray() was called.
588: Not Collective
590: Input Parameters:
591: + ds - the direct solver context
592: . m - the requested matrix
593: - a - pointer to the values
595: Level: advanced
597: .seealso: DSGetArray(), DSGetArrayReal()
598: @*/
599: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
600: {
602: DSCheckAlloc(ds,1);
603: DSCheckValidMat(ds,m,2);
605: MatDenseRestoreArray(ds->omat[m],a);
606: PetscObjectStateIncrease((PetscObject)ds);
607: return 0;
608: }
610: /*@C
611: DSGetArrayReal - Returns a real pointer to the internal array of T or D.
612: You MUST call DSRestoreArrayReal() when you no longer need to access the array.
614: Not Collective
616: Input Parameters:
617: + ds - the direct solver context
618: - m - the requested matrix
620: Output Parameter:
621: . a - pointer to the values
623: Note:
624: This function can be used only for DS_MAT_T and DS_MAT_D. These matrices always
625: store real values, even in complex scalars, that is why the returned pointer is
626: PetscReal.
628: Level: advanced
630: .seealso: DSRestoreArrayReal(), DSGetArray()
631: @*/
632: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
633: {
634: #if defined(PETSC_USE_COMPLEX)
635: PetscScalar *as;
636: #endif
639: DSCheckAlloc(ds,1);
640: DSCheckValidMatReal(ds,m,2);
642: #if defined(PETSC_USE_COMPLEX)
643: MatDenseGetArray(ds->omat[m],&as);
644: *a = (PetscReal*)as;
645: #else
646: MatDenseGetArray(ds->omat[m],a);
647: #endif
648: return 0;
649: }
651: /*@C
652: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
654: Not Collective
656: Input Parameters:
657: + ds - the direct solver context
658: . m - the requested matrix
659: - a - pointer to the values
661: Level: advanced
663: .seealso: DSGetArrayReal(), DSGetArray()
664: @*/
665: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
666: {
667: #if defined(PETSC_USE_COMPLEX)
668: PetscScalar *as;
669: #endif
672: DSCheckAlloc(ds,1);
673: DSCheckValidMatReal(ds,m,2);
675: #if defined(PETSC_USE_COMPLEX)
676: MatDenseRestoreArray(ds->omat[m],&as);
677: *a = NULL;
678: #else
679: MatDenseRestoreArray(ds->omat[m],a);
680: #endif
681: PetscObjectStateIncrease((PetscObject)ds);
682: return 0;
683: }
685: /*@
686: DSSolve - Solves the problem.
688: Logically Collective on ds
690: Input Parameters:
691: + ds - the direct solver context
692: . eigr - array to store the computed eigenvalues (real part)
693: - eigi - array to store the computed eigenvalues (imaginary part)
695: Note:
696: This call brings the dense system to condensed form. No ordering
697: of the eigenvalues is enforced (for this, call DSSort() afterwards).
699: Level: intermediate
701: .seealso: DSSort(), DSStateType
702: @*/
703: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
704: {
707: DSCheckAlloc(ds,1);
709: if (ds->state>=DS_STATE_CONDENSED) return 0;
711: PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
712: PetscLogEventBegin(DS_Solve,ds,0,0,0);
713: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
714: PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
715: PetscFPTrapPop();
716: PetscLogEventEnd(DS_Solve,ds,0,0,0);
717: PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
718: ds->state = DS_STATE_CONDENSED;
719: PetscObjectStateIncrease((PetscObject)ds);
720: return 0;
721: }
723: /*@C
724: DSSort - Sorts the result of DSSolve() according to a given sorting
725: criterion.
727: Logically Collective on ds
729: Input Parameters:
730: + ds - the direct solver context
731: . rr - (optional) array containing auxiliary values (real part)
732: - ri - (optional) array containing auxiliary values (imaginary part)
734: Input/Output Parameters:
735: + eigr - array containing the computed eigenvalues (real part)
736: . eigi - array containing the computed eigenvalues (imaginary part)
737: - k - (optional) number of elements in the leading group
739: Notes:
740: This routine sorts the arrays provided in eigr and eigi, and also
741: sorts the dense system stored inside ds (assumed to be in condensed form).
742: The sorting criterion is specified with DSSetSlepcSC().
744: If arrays rr and ri are provided, then a (partial) reordering based on these
745: values rather than on the eigenvalues is performed. In symmetric problems
746: a total order is obtained (parameter k is ignored), but otherwise the result
747: is sorted only partially. In this latter case, it is only guaranteed that
748: all the first k elements satisfy the comparison with any of the last n-k
749: elements. The output value of parameter k is the final number of elements in
750: the first set.
752: Level: intermediate
754: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
755: @*/
756: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
757: {
758: PetscInt i;
762: DSCheckSolved(ds,1);
769: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
770: PetscLogEventBegin(DS_Other,ds,0,0,0);
771: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
772: PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
773: PetscFPTrapPop();
774: PetscLogEventEnd(DS_Other,ds,0,0,0);
775: PetscObjectStateIncrease((PetscObject)ds);
776: PetscInfo(ds,"Finished sorting\n");
777: return 0;
778: }
780: /*@C
781: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
782: permutation.
784: Logically Collective on ds
786: Input Parameters:
787: + ds - the direct solver context
788: - perm - permutation that indicates the new ordering
790: Input/Output Parameters:
791: + eigr - array with the reordered eigenvalues (real part)
792: - eigi - array with the reordered eigenvalues (imaginary part)
794: Notes:
795: This routine reorders the arrays provided in eigr and eigi, and also the dense
796: system stored inside ds (assumed to be in condensed form). There is no sorting
797: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
799: Level: advanced
801: .seealso: DSSolve(), DSSort()
802: @*/
803: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)
804: {
807: DSCheckSolved(ds,1);
812: PetscLogEventBegin(DS_Other,ds,0,0,0);
813: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
814: PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
815: PetscFPTrapPop();
816: PetscLogEventEnd(DS_Other,ds,0,0,0);
817: PetscObjectStateIncrease((PetscObject)ds);
818: PetscInfo(ds,"Finished sorting\n");
819: return 0;
820: }
822: /*@
823: DSSynchronize - Make sure that all processes have the same data, performing
824: communication if necessary.
826: Collective on ds
828: Input Parameter:
829: . ds - the direct solver context
831: Input/Output Parameters:
832: + eigr - (optional) array with the computed eigenvalues (real part)
833: - eigi - (optional) array with the computed eigenvalues (imaginary part)
835: Notes:
836: When the DS has been created with a communicator with more than one process,
837: the internal data, especially the computed matrices, may diverge in the
838: different processes. This happens when using multithreaded BLAS and may
839: cause numerical issues in some ill-conditioned problems. This function
840: performs the necessary communication among the processes so that the
841: internal data is exactly equal in all of them.
843: Depending on the parallel mode as set with DSSetParallel(), this function
844: will either do nothing or synchronize the matrices computed by DSSolve()
845: and DSSort(). The arguments eigr and eigi are typically those used in the
846: calls to DSSolve() and DSSort().
848: Level: developer
850: .seealso: DSSetParallel(), DSSolve(), DSSort()
851: @*/
852: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
853: {
854: PetscMPIInt size;
858: DSCheckAlloc(ds,1);
859: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
860: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
861: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
862: PetscUseTypeMethod(ds,synchronize,eigr,eigi);
863: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
864: PetscObjectStateIncrease((PetscObject)ds);
865: PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
866: }
867: return 0;
868: }
870: /*@C
871: DSVectors - Compute vectors associated to the dense system such
872: as eigenvectors.
874: Logically Collective on ds
876: Input Parameters:
877: + ds - the direct solver context
878: - mat - the matrix, used to indicate which vectors are required
880: Input/Output Parameter:
881: . j - (optional) index of vector to be computed
883: Output Parameter:
884: . rnorm - (optional) computed residual norm
886: Notes:
887: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
888: compute right or left eigenvectors, or left or right singular vectors,
889: respectively.
891: If NULL is passed in argument j then all vectors are computed,
892: otherwise j indicates which vector must be computed. In real non-symmetric
893: problems, on exit the index j will be incremented when a complex conjugate
894: pair is found.
896: This function can be invoked after the dense problem has been solved,
897: to get the residual norm estimate of the associated Ritz pair. In that
898: case, the relevant information is returned in rnorm.
900: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
901: be in (quasi-)triangular form, or call DSSolve() first.
903: Level: intermediate
905: .seealso: DSSolve()
906: @*/
907: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
908: {
911: DSCheckAlloc(ds,1);
915: if (!ds->omat[mat]) DSAllocateMat_Private(ds,mat);
916: if (!j) PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]);
917: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
918: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
919: PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
920: PetscFPTrapPop();
921: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
922: PetscObjectStateIncrease((PetscObject)ds);
923: return 0;
924: }
926: /*@
927: DSUpdateExtraRow - Performs all necessary operations so that the extra
928: row gets up-to-date after a call to DSSolve().
930: Not Collective
932: Input Parameters:
933: . ds - the direct solver context
935: Level: advanced
937: .seealso: DSSolve(), DSSetExtraRow()
938: @*/
939: PetscErrorCode DSUpdateExtraRow(DS ds)
940: {
943: DSCheckAlloc(ds,1);
945: PetscInfo(ds,"Updating extra row\n");
946: PetscLogEventBegin(DS_Other,ds,0,0,0);
947: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
948: PetscUseTypeMethod(ds,update);
949: PetscFPTrapPop();
950: PetscLogEventEnd(DS_Other,ds,0,0,0);
951: return 0;
952: }
954: /*@
955: DSCond - Compute the condition number.
957: Not Collective
959: Input Parameters:
960: + ds - the direct solver context
961: - cond - the computed condition number
963: Notes:
964: In standard eigenvalue problems, returns the inf-norm condition number of the first
965: matrix, computed as cond(A) = norm(A)*norm(inv(A)).
967: In GSVD problems, returns the maximum of cond(A) and cond(B), where cond(.) is
968: computed as the ratio of the largest and smallest singular values.
970: Does not take into account the extra row.
972: Level: advanced
974: .seealso: DSSolve(), DSSetExtraRow()
975: @*/
976: PetscErrorCode DSCond(DS ds,PetscReal *cond)
977: {
980: DSCheckAlloc(ds,1);
982: PetscLogEventBegin(DS_Other,ds,0,0,0);
983: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
984: PetscUseTypeMethod(ds,cond,cond);
985: PetscFPTrapPop();
986: PetscLogEventEnd(DS_Other,ds,0,0,0);
987: PetscInfo(ds,"Computed condition number = %g\n",(double)*cond);
988: return 0;
989: }
991: /*@C
992: DSTranslateHarmonic - Computes a translation of the dense system.
994: Logically Collective on ds
996: Input Parameters:
997: + ds - the direct solver context
998: . tau - the translation amount
999: . beta - last component of vector b
1000: - recover - boolean flag to indicate whether to recover or not
1002: Output Parameters:
1003: + g - the computed vector (optional)
1004: - gamma - scale factor (optional)
1006: Notes:
1007: This function is intended for use in the context of Krylov methods only.
1008: It computes a translation of a Krylov decomposition in order to extract
1009: eigenpair approximations by harmonic Rayleigh-Ritz.
1010: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
1011: vector b is assumed to be beta*e_n^T.
1013: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
1014: the factor by which the residual of the Krylov decomposition is scaled.
1016: If the recover flag is activated, the computed translation undoes the
1017: translation done previously. In that case, parameter tau is ignored.
1019: Level: developer
1021: .seealso: DSTranslateRKS()
1022: @*/
1023: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
1024: {
1027: DSCheckAlloc(ds,1);
1028: if (recover) PetscInfo(ds,"Undoing the translation\n");
1029: else PetscInfo(ds,"Computing the translation\n");
1030: PetscLogEventBegin(DS_Other,ds,0,0,0);
1031: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1032: PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
1033: PetscFPTrapPop();
1034: PetscLogEventEnd(DS_Other,ds,0,0,0);
1035: ds->state = DS_STATE_RAW;
1036: PetscObjectStateIncrease((PetscObject)ds);
1037: return 0;
1038: }
1040: /*@
1041: DSTranslateRKS - Computes a modification of the dense system corresponding
1042: to an update of the shift in a rational Krylov method.
1044: Logically Collective on ds
1046: Input Parameters:
1047: + ds - the direct solver context
1048: - alpha - the translation amount
1050: Notes:
1051: This function is intended for use in the context of Krylov methods only.
1052: It takes the leading (k+1,k) submatrix of A, containing the truncated
1053: Rayleigh quotient of a Krylov-Schur relation computed from a shift
1054: sigma1 and transforms it to obtain a Krylov relation as if computed
1055: from a different shift sigma2. The new matrix is computed as
1056: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
1057: alpha = sigma1-sigma2.
1059: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
1060: Krylov basis.
1062: Level: developer
1064: .seealso: DSTranslateHarmonic()
1065: @*/
1066: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
1067: {
1070: DSCheckAlloc(ds,1);
1071: PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha));
1072: PetscLogEventBegin(DS_Other,ds,0,0,0);
1073: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1074: PetscUseTypeMethod(ds,transrks,alpha);
1075: PetscFPTrapPop();
1076: PetscLogEventEnd(DS_Other,ds,0,0,0);
1077: ds->state = DS_STATE_RAW;
1078: ds->compact = PETSC_FALSE;
1079: PetscObjectStateIncrease((PetscObject)ds);
1080: return 0;
1081: }