1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: {
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) { PetscInfo2(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]); }
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: . l - number of locked (inactive) leading columns
117: - k - intermediate dimension (e.g., position of arrow)
119: Notes:
120: The internal arrays are not reallocated.
122: Some DS types have additional dimensions, e.g. the number of columns
123: in DSSVD. For these, you should call a specific interface function.
125: Level: intermediate
127: .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
128: @*/
129: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)130: {
132: PetscInt on,ol,ok;
133: PetscBool issvd;
137: DSCheckAlloc(ds,1);
141: on = ds->n; ol = ds->l; ok = ds->k;
142: if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
143: ds->n = ds->ld;
144: } else {
145: if (n<0 || n>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
146: PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,""); /* SVD and GSVD have extra column instead of extra row */
147: if (ds->extrarow && n+1>ds->ld && !issvd) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"A value of n equal to ld leaves no room for extra row");
148: ds->n = n;
149: }
150: ds->t = ds->n; /* truncated length equal to the new dimension */
151: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
152: ds->l = 0;
153: } else {
154: if (l<0 || l>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
155: ds->l = l;
156: }
157: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
158: ds->k = ds->n/2;
159: } else {
160: if (k<0 || k>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
161: ds->k = k;
162: }
163: if (on!=ds->n || ol!=ds->l || ok!=ds->k) {
164: PetscInfo3(ds,"New dimensions are: n=%D, l=%D, k=%D\n",ds->n,ds->l,ds->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 Parameters:
178: + n - the current size
179: . l - number of locked (inactive) leading columns
180: . k - intermediate dimension (e.g., position of arrow)
181: - t - truncated length
183: Note:
184: The t parameter makes sense only if DSTruncate() has been called.
185: Otherwise its value equals n.
187: Some DS types have additional dimensions, e.g. the number of columns
188: in DSSVD. For these, you should call a specific interface function.
190: Level: intermediate
192: .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
193: @*/
194: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)195: {
198: DSCheckAlloc(ds,1);
199: if (n) *n = ds->n;
200: if (l) *l = ds->l;
201: if (k) *k = ds->k;
202: if (t) *t = ds->t;
203: return(0);
204: }
206: /*@
207: DSTruncate - Truncates the system represented in the DS object.
209: Logically Collective on ds
211: Input Parameters:
212: + ds - the direct solver context
213: . n - the new size
214: - trim - a flag to indicate if the factorization must be trimmed
216: Note:
217: The new size is set to n. Note that in some cases the new size could
218: be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
219: Schur form). In cases where the extra row is meaningful, the first
220: n elements are kept as the extra row for the new system.
222: If the flag trim is turned on, it resets the locked and intermediate
223: dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
224: It also cleans the extra row if being used.
226: The typical usage of trim=true is to truncate the Schur decomposition
227: at the end of a Krylov iteration. In this case, the state must be
228: changed to RAW so that DSVectors() computes eigenvectors from scratch.
230: Level: advanced
232: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType233: @*/
234: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)235: {
237: DSStateType old;
242: DSCheckAlloc(ds,1);
245: if (!ds->ops->truncate) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
246: 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);
247: PetscLogEventBegin(DS_Other,ds,0,0,0);
248: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
249: (*ds->ops->truncate)(ds,n,trim);
250: PetscFPTrapPop();
251: PetscLogEventEnd(DS_Other,ds,0,0,0);
252: PetscInfo2(ds,"Decomposition %s to size n=%D\n",trim?"trimmed":"truncated",ds->n);
253: old = ds->state;
254: ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
255: if (old!=ds->state) {
256: PetscInfo2(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]);
257: }
258: PetscObjectStateIncrease((PetscObject)ds);
259: return(0);
260: }
262: /*@
263: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
265: Not Collective
267: Input Parameters:
268: + ds - the direct solver context
269: - t - the requested matrix
271: Output Parameters:
272: + n - the number of rows
273: - m - the number of columns
275: Note:
276: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
278: Level: developer
280: .seealso: DSSetDimensions(), DSGetMat()
281: @*/
282: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)283: {
285: PetscInt rows,cols;
290: DSCheckValidMat(ds,t,2);
291: if (ds->ops->matgetsize) {
292: (*ds->ops->matgetsize)(ds,t,&rows,&cols);
293: } else {
294: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
295: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
296: cols = ds->n;
297: }
298: if (m) *m = rows;
299: if (n) *n = cols;
300: return(0);
301: }
303: /*@
304: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
306: Not Collective
308: Input Parameters:
309: + ds - the direct solver context
310: - t - the requested matrix
312: Output Parameter:
313: . flg - the Hermitian flag
315: Note:
316: Does not check the matrix values directly. The flag is set according to the
317: problem structure. For instance, in DSHEP matrix A is Hermitian.
319: Level: developer
321: .seealso: DSGetMat()
322: @*/
323: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)324: {
330: DSCheckValidMat(ds,t,2);
332: if (ds->ops->hermitian) {
333: (*ds->ops->hermitian)(ds,t,flg);
334: } else *flg = PETSC_FALSE;
335: return(0);
336: }
338: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)339: {
340: #if !defined(PETSC_USE_COMPLEX)
341: PetscScalar *A = ds->mat[DS_MAT_A];
342: #endif
345: #if !defined(PETSC_USE_COMPLEX)
346: if (A[l+(*k)+(l+(*k)-1)*ds->ld] != 0.0) {
347: if (l+(*k)<n-1) (*k)++;
348: else (*k)--;
349: }
350: #endif
351: return(0);
352: }
354: /*@
355: DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
356: to avoid breaking a 2x2 block.
358: Not Collective
360: Input Parameters:
361: + ds - the direct solver context
362: . l - the size of the locked part (set to 0 to use ds->l)
363: . n - the total matrix size (set to 0 to use ds->n)
364: - k - the wanted truncation size
366: Output Parameter:
367: . k - the possibly modified value of the truncation size
369: Notes:
370: This should be called before DSTruncate() to make sure that the truncation
371: does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
373: The total size is n (either user-provided or ds->n if 0 is passed). The
374: size where the truncation is intended is equal to l+k (where l can be
375: equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
376: at the l+k limit, the value of k is increased (or decreased) by 1.
378: Level: advanced
380: .seealso: DSTruncate(), DSSetDimensions()
381: @*/
382: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)383: {
388: DSCheckAlloc(ds,1);
392: if (ds->ops->gettruncatesize) {
393: (*ds->ops->gettruncatesize)(ds,l?l:ds->l,n?n:ds->n,k);
394: }
395: return(0);
396: }
398: /*@
399: DSGetMat - Returns a sequential dense Mat object containing the requested
400: matrix.
402: Not Collective
404: Input Parameters:
405: + ds - the direct solver context
406: - m - the requested matrix
408: Output Parameter:
409: . A - Mat object
411: Notes:
412: The Mat is created with sizes equal to the current DS dimensions (nxm),
413: then it is filled with the values that would be obtained with DSGetArray()
414: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
415: is equal to the dimension prior to truncation, see DSTruncate().
416: The communicator is always PETSC_COMM_SELF.
418: When no longer needed, the user can either destroy the matrix or call
419: DSRestoreMat(). The latter will copy back the modified values.
421: Level: advanced
423: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
424: @*/
425: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)426: {
428: PetscInt j,rows,cols,arows,acols;
429: PetscBool create=PETSC_FALSE,flg;
430: PetscScalar *pA,*M;
434: DSCheckAlloc(ds,1);
435: DSCheckValidMat(ds,m,2);
437: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
439: DSMatGetSize(ds,m,&rows,&cols);
440: if (!ds->omat[m]) create=PETSC_TRUE;
441: else {
442: MatGetSize(ds->omat[m],&arows,&acols);
443: if (arows!=rows || acols!=cols) {
444: MatDestroy(&ds->omat[m]);
445: create=PETSC_TRUE;
446: }
447: }
448: if (create) {
449: MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);
450: }
452: /* set Hermitian flag */
453: DSMatIsHermitian(ds,m,&flg);
454: MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);
456: /* copy entries */
457: PetscObjectReference((PetscObject)ds->omat[m]);
458: *A = ds->omat[m];
459: M = ds->mat[m];
460: MatDenseGetArray(*A,&pA);
461: for (j=0;j<cols;j++) {
462: PetscArraycpy(pA+j*rows,M+j*ds->ld,rows);
463: }
464: MatDenseRestoreArray(*A,&pA);
465: return(0);
466: }
468: /*@
469: DSRestoreMat - Restores the matrix after DSGetMat() was called.
471: Not Collective
473: Input Parameters:
474: + ds - the direct solver context
475: . m - the requested matrix
476: - A - the fetched Mat object
478: Notes:
479: A call to this function must match a previous call of DSGetMat().
480: The effect is that the contents of the Mat are copied back to the
481: DS internal array, and the matrix is destroyed.
483: It is not compulsory to call this function, the matrix obtained with
484: DSGetMat() can simply be destroyed if entries need not be copied back.
486: Level: advanced
488: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
489: @*/
490: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)491: {
493: PetscInt j,rows,cols;
494: PetscScalar *pA,*M;
498: DSCheckAlloc(ds,1);
499: DSCheckValidMat(ds,m,2);
501: if (!ds->omat[m]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSRestoreMat must match a previous call to DSGetMat");
502: 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");
504: MatGetSize(*A,&rows,&cols);
505: M = ds->mat[m];
506: MatDenseGetArray(*A,&pA);
507: for (j=0;j<cols;j++) {
508: PetscArraycpy(M+j*ds->ld,pA+j*rows,rows);
509: }
510: MatDenseRestoreArray(*A,&pA);
511: MatDestroy(A);
512: return(0);
513: }
515: /*@C
516: DSGetArray - Returns a pointer to one of the internal arrays used to
517: represent matrices. You MUST call DSRestoreArray() when you no longer
518: need to access the array.
520: Not Collective
522: Input Parameters:
523: + ds - the direct solver context
524: - m - the requested matrix
526: Output Parameter:
527: . a - pointer to the values
529: Level: advanced
531: .seealso: DSRestoreArray(), DSGetArrayReal()
532: @*/
533: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])534: {
537: DSCheckAlloc(ds,1);
538: DSCheckValidMat(ds,m,2);
540: *a = ds->mat[m];
541: CHKMEMQ;
542: return(0);
543: }
545: /*@C
546: DSRestoreArray - Restores the matrix after DSGetArray() was called.
548: Not Collective
550: Input Parameters:
551: + ds - the direct solver context
552: . m - the requested matrix
553: - a - pointer to the values
555: Level: advanced
557: .seealso: DSGetArray(), DSGetArrayReal()
558: @*/
559: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])560: {
565: DSCheckAlloc(ds,1);
566: DSCheckValidMat(ds,m,2);
568: CHKMEMQ;
569: *a = 0;
570: PetscObjectStateIncrease((PetscObject)ds);
571: return(0);
572: }
574: /*@C
575: DSGetArrayReal - Returns a pointer to one of the internal arrays used to
576: represent real matrices. You MUST call DSRestoreArrayReal() when you no longer
577: need to access the array.
579: Not Collective
581: Input Parameters:
582: + ds - the direct solver context
583: - m - the requested matrix
585: Output Parameter:
586: . a - pointer to the values
588: Level: advanced
590: .seealso: DSRestoreArrayReal(), DSGetArray()
591: @*/
592: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])593: {
596: DSCheckAlloc(ds,1);
597: DSCheckValidMatReal(ds,m,2);
599: *a = ds->rmat[m];
600: CHKMEMQ;
601: return(0);
602: }
604: /*@C
605: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
607: Not Collective
609: Input Parameters:
610: + ds - the direct solver context
611: . m - the requested matrix
612: - a - pointer to the values
614: Level: advanced
616: .seealso: DSGetArrayReal(), DSGetArray()
617: @*/
618: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])619: {
624: DSCheckAlloc(ds,1);
625: DSCheckValidMatReal(ds,m,2);
627: CHKMEMQ;
628: *a = 0;
629: PetscObjectStateIncrease((PetscObject)ds);
630: return(0);
631: }
633: /*@
634: DSSolve - Solves the problem.
636: Logically Collective on ds
638: Input Parameters:
639: + ds - the direct solver context
640: . eigr - array to store the computed eigenvalues (real part)
641: - eigi - array to store the computed eigenvalues (imaginary part)
643: Note:
644: This call brings the dense system to condensed form. No ordering
645: of the eigenvalues is enforced (for this, call DSSort() afterwards).
647: Level: intermediate
649: .seealso: DSSort(), DSStateType650: @*/
651: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])652: {
658: DSCheckAlloc(ds,1);
660: if (ds->state>=DS_STATE_CONDENSED) return(0);
661: if (!ds->ops->solve[ds->method]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
662: PetscInfo3(ds,"Starting solve with problem sizes: n=%D, l=%D, k=%D\n",ds->n,ds->l,ds->k);
663: PetscLogEventBegin(DS_Solve,ds,0,0,0);
664: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
665: (*ds->ops->solve[ds->method])(ds,eigr,eigi);
666: PetscFPTrapPop();
667: PetscLogEventEnd(DS_Solve,ds,0,0,0);
668: PetscInfo1(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
669: ds->state = DS_STATE_CONDENSED;
670: PetscObjectStateIncrease((PetscObject)ds);
671: return(0);
672: }
674: /*@C
675: DSSort - Sorts the result of DSSolve() according to a given sorting
676: criterion.
678: Logically Collective on ds
680: Input Parameters:
681: + ds - the direct solver context
682: . rr - (optional) array containing auxiliary values (real part)
683: - ri - (optional) array containing auxiliary values (imaginary part)
685: Input/Output Parameters:
686: + eigr - array containing the computed eigenvalues (real part)
687: . eigi - array containing the computed eigenvalues (imaginary part)
688: - k - (optional) number of elements in the leading group
690: Notes:
691: This routine sorts the arrays provided in eigr and eigi, and also
692: sorts the dense system stored inside ds (assumed to be in condensed form).
693: The sorting criterion is specified with DSSetSlepcSC().
695: If arrays rr and ri are provided, then a (partial) reordering based on these
696: values rather than on the eigenvalues is performed. In symmetric problems
697: a total order is obtained (parameter k is ignored), but otherwise the result
698: is sorted only partially. In this latter case, it is only guaranteed that
699: all the first k elements satisfy the comparison with any of the last n-k
700: elements. The output value of parameter k is the final number of elements in
701: the first set.
703: Level: intermediate
705: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
706: @*/
707: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)708: {
710: PetscInt i;
715: DSCheckSolved(ds,1);
718: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
719: if (!ds->ops->sort) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
720: if (!ds->sc) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
721: if (k && !rr) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
723: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
724: PetscLogEventBegin(DS_Other,ds,0,0,0);
725: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
726: (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
727: PetscFPTrapPop();
728: PetscLogEventEnd(DS_Other,ds,0,0,0);
729: PetscObjectStateIncrease((PetscObject)ds);
730: PetscInfo(ds,"Finished sorting\n");
731: return(0);
732: }
734: /*@C
735: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
736: permutation.
738: Logically Collective on ds
740: Input Parameters:
741: + ds - the direct solver context
742: - perm - permutation that indicates the new ordering
744: Input/Output Parameters:
745: + eigr - array with the reordered eigenvalues (real part)
746: - eigi - array with the reordered eigenvalues (imaginary part)
748: Notes:
749: This routine reorders the arrays provided in eigr and eigi, and also the dense
750: system stored inside ds (assumed to be in condensed form). There is no sorting
751: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
753: Level: advanced
755: .seealso: DSSolve(), DSSort()
756: @*/
757: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)758: {
764: DSCheckSolved(ds,1);
767: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
768: if (!ds->ops->sortperm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
770: PetscLogEventBegin(DS_Other,ds,0,0,0);
771: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
772: (*ds->ops->sortperm)(ds,perm,eigr,eigi);
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: /*@
781: DSSynchronize - Make sure that all processes have the same data, performing
782: communication if necessary.
784: Collective on ds
786: Input Parameter:
787: . ds - the direct solver context
789: Input/Output Parameters:
790: + eigr - (optional) array with the computed eigenvalues (real part)
791: - eigi - (optional) array with the computed eigenvalues (imaginary part)
793: Notes:
794: When the DS has been created with a communicator with more than one process,
795: the internal data, especially the computed matrices, may diverge in the
796: different processes. This happens when using multithreaded BLAS and may
797: cause numerical issues in some ill-conditioned problems. This function
798: performs the necessary communication among the processes so that the
799: internal data is exactly equal in all of them.
801: Depending on the parallel mode as set with DSSetParallel(), this function
802: will either do nothing or synchronize the matrices computed by DSSolve()
803: and DSSort(). The arguments eigr and eigi are typically those used in the
804: calls to DSSolve() and DSSort().
806: Level: developer
808: .seealso: DSSetParallel(), DSSolve(), DSSort()
809: @*/
810: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])811: {
813: PetscMPIInt size;
818: DSCheckAlloc(ds,1);
819: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
820: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
821: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
822: if (ds->ops->synchronize) {
823: (*ds->ops->synchronize)(ds,eigr,eigi);
824: }
825: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
826: PetscObjectStateIncrease((PetscObject)ds);
827: PetscInfo1(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
828: }
829: return(0);
830: }
832: /*@C
833: DSVectors - Compute vectors associated to the dense system such
834: as eigenvectors.
836: Logically Collective on ds
838: Input Parameters:
839: + ds - the direct solver context
840: - mat - the matrix, used to indicate which vectors are required
842: Input/Output Parameter:
843: - j - (optional) index of vector to be computed
845: Output Parameter:
846: . rnorm - (optional) computed residual norm
848: Notes:
849: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
850: compute right or left eigenvectors, or left or right singular vectors,
851: respectively.
853: If NULL is passed in argument j then all vectors are computed,
854: otherwise j indicates which vector must be computed. In real non-symmetric
855: problems, on exit the index j will be incremented when a complex conjugate
856: pair is found.
858: This function can be invoked after the dense problem has been solved,
859: to get the residual norm estimate of the associated Ritz pair. In that
860: case, the relevant information is returned in rnorm.
862: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
863: be in (quasi-)triangular form, or call DSSolve() first.
865: Level: intermediate
867: .seealso: DSSolve()
868: @*/
869: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)870: {
876: DSCheckAlloc(ds,1);
878: if (mat>=DS_NUM_MAT) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
879: if (!ds->ops->vectors) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
880: if (rnorm && !j) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
881: if (!ds->mat[mat]) { DSAllocateMat_Private(ds,mat); }
882: if (!j) { PetscInfo1(ds,"Computing all vectors on %s\n",DSMatName[mat]); }
883: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
884: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
885: (*ds->ops->vectors)(ds,mat,j,rnorm);
886: PetscFPTrapPop();
887: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
888: PetscObjectStateIncrease((PetscObject)ds);
889: return(0);
890: }
892: /*@
893: DSUpdateExtraRow - Performs all necessary operations so that the extra
894: row gets up-to-date after a call to DSSolve().
896: Not Collective
898: Input Parameters:
899: . ds - the direct solver context
901: Level: advanced
903: .seealso: DSSolve(), DSSetExtraRow()
904: @*/
905: PetscErrorCode DSUpdateExtraRow(DS ds)906: {
912: DSCheckAlloc(ds,1);
913: if (!ds->ops->update) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
914: if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
915: PetscInfo(ds,"Updating extra row\n");
916: PetscLogEventBegin(DS_Other,ds,0,0,0);
917: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
918: (*ds->ops->update)(ds);
919: PetscFPTrapPop();
920: PetscLogEventEnd(DS_Other,ds,0,0,0);
921: return(0);
922: }
924: /*@
925: DSCond - Compute the inf-norm condition number of the first matrix
926: as cond(A) = norm(A)*norm(inv(A)).
928: Not Collective
930: Input Parameters:
931: + ds - the direct solver context
932: - cond - the computed condition number
934: Level: advanced
936: .seealso: DSSolve()
937: @*/
938: PetscErrorCode DSCond(DS ds,PetscReal *cond)939: {
945: DSCheckAlloc(ds,1);
947: if (!ds->ops->cond) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
948: PetscLogEventBegin(DS_Other,ds,0,0,0);
949: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
950: (*ds->ops->cond)(ds,cond);
951: PetscFPTrapPop();
952: PetscLogEventEnd(DS_Other,ds,0,0,0);
953: PetscInfo1(ds,"Computed condition number = %g\n",(double)*cond);
954: return(0);
955: }
957: /*@C
958: DSTranslateHarmonic - Computes a translation of the dense system.
960: Logically Collective on ds
962: Input Parameters:
963: + ds - the direct solver context
964: . tau - the translation amount
965: . beta - last component of vector b
966: - recover - boolean flag to indicate whether to recover or not
968: Output Parameters:
969: + g - the computed vector (optional)
970: - gamma - scale factor (optional)
972: Notes:
973: This function is intended for use in the context of Krylov methods only.
974: It computes a translation of a Krylov decomposition in order to extract
975: eigenpair approximations by harmonic Rayleigh-Ritz.
976: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
977: vector b is assumed to be beta*e_n^T.
979: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
980: the factor by which the residual of the Krylov decomposition is scaled.
982: If the recover flag is activated, the computed translation undoes the
983: translation done previously. In that case, parameter tau is ignored.
985: Level: developer
986: @*/
987: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)988: {
994: DSCheckAlloc(ds,1);
995: if (!ds->ops->transharm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
996: if (recover) { PetscInfo(ds,"Undoing the translation\n"); }
997: else { PetscInfo(ds,"Computing the translation\n"); }
998: PetscLogEventBegin(DS_Other,ds,0,0,0);
999: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1000: (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
1001: PetscFPTrapPop();
1002: PetscLogEventEnd(DS_Other,ds,0,0,0);
1003: ds->state = DS_STATE_RAW;
1004: PetscObjectStateIncrease((PetscObject)ds);
1005: return(0);
1006: }
1008: /*@
1009: DSTranslateRKS - Computes a modification of the dense system corresponding
1010: to an update of the shift in a rational Krylov method.
1012: Logically Collective on ds
1014: Input Parameters:
1015: + ds - the direct solver context
1016: - alpha - the translation amount
1018: Notes:
1019: This function is intended for use in the context of Krylov methods only.
1020: It takes the leading (k+1,k) submatrix of A, containing the truncated
1021: Rayleigh quotient of a Krylov-Schur relation computed from a shift
1022: sigma1 and transforms it to obtain a Krylov relation as if computed
1023: from a different shift sigma2. The new matrix is computed as
1024: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
1025: alpha = sigma1-sigma2.
1027: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
1028: Krylov basis.
1030: Level: developer
1031: @*/
1032: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)1033: {
1039: DSCheckAlloc(ds,1);
1040: if (!ds->ops->transrks) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
1041: PetscInfo1(ds,"Translating with alpha=%g\n",PetscRealPart(alpha));
1042: PetscLogEventBegin(DS_Other,ds,0,0,0);
1043: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1044: (*ds->ops->transrks)(ds,alpha);
1045: PetscFPTrapPop();
1046: PetscLogEventEnd(DS_Other,ds,0,0,0);
1047: ds->state = DS_STATE_RAW;
1048: ds->compact = PETSC_FALSE;
1049: PetscObjectStateIncrease((PetscObject)ds);
1050: return(0);
1051: }
1053: /*@
1054: DSCopyMat - Copies the contents of a sequential dense Mat object to
1055: the indicated DS matrix, or vice versa.
1057: Not Collective
1059: Input Parameters:
1060: + ds - the direct solver context
1061: . m - the requested matrix
1062: . mr - first row of m to be considered
1063: . mc - first column of m to be considered
1064: . A - Mat object
1065: . Ar - first row of A to be considered
1066: . Ac - first column of A to be considered
1067: . rows - number of rows to copy
1068: . cols - number of columns to copy
1069: - out - whether the data is copied out of the DS1071: Note:
1072: If out=true, the values of the DS matrix m are copied to A, otherwise
1073: the entries of A are copied to the DS.
1075: Level: developer
1077: .seealso: DSGetMat()
1078: @*/
1079: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)1080: {
1082: PetscInt j,mrows,mcols,arows,acols;
1083: PetscScalar *pA,*M;
1087: DSCheckAlloc(ds,1);
1089: DSCheckValidMat(ds,m,2);
1098: if (!rows || !cols) return(0);
1100: DSMatGetSize(ds,m,&mrows,&mcols);
1101: MatGetSize(A,&arows,&acols);
1102: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
1103: if (mr<0 || mr>=mrows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in m");
1104: if (mc<0 || mc>=mcols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in m");
1105: if (Ar<0 || Ar>=arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in A");
1106: if (Ac<0 || Ac>=acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in A");
1107: if (mr+rows>mrows || Ar+rows>arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of rows");
1108: if (mc+cols>mcols || Ac+cols>acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of columns");
1110: M = ds->mat[m];
1111: MatDenseGetArray(A,&pA);
1112: for (j=0;j<cols;j++) {
1113: if (out) {
1114: PetscArraycpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows);
1115: } else {
1116: PetscArraycpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows);
1117: }
1118: }
1119: MatDenseRestoreArray(A,&pA);
1120: return(0);
1121: }