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: . 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: {
132: PetscInt on,om,ol,ok;
136: DSCheckAlloc(ds,1);
141: on = ds->n; om = ds->m; 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: 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");
147: ds->n = n;
148: }
149: ds->t = ds->n; /* truncated length equal to the new dimension */
150: if (m) {
151: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
152: ds->m = ds->ld;
153: } else {
154: if (m<0 || m>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 0 and ld");
155: ds->m = m;
156: }
157: }
158: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
159: ds->l = 0;
160: } else {
161: if (l<0 || l>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
162: ds->l = l;
163: }
164: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
165: ds->k = ds->n/2;
166: } else {
167: if (k<0 || k>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
168: ds->k = k;
169: }
170: if (on!=ds->n || om!=ds->m || ol!=ds->l || ok!=ds->k) {
171: PetscInfo4(ds,"New dimensions are: n=%D, m=%D, l=%D, k=%D\n",ds->n,ds->m,ds->l,ds->k);
172: }
173: return(0);
174: }
176: /*@
177: DSGetDimensions - Returns the current dimensions.
179: Not Collective
181: Input Parameter:
182: . ds - the direct solver context
184: Output Parameter:
185: + n - the current size
186: . m - the current column size (only for DSSVD)
187: . l - number of locked (inactive) leading columns
188: . k - intermediate dimension (e.g., position of arrow)
189: - t - truncated length
191: Note:
192: The t parameter makes sense only if DSTruncate() has been called.
193: Otherwise its value equals n.
195: Level: intermediate
197: .seealso: DSSetDimensions(), DSTruncate()
198: @*/
199: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *m,PetscInt *l,PetscInt *k,PetscInt *t)200: {
203: DSCheckAlloc(ds,1);
204: if (n) *n = ds->n;
205: if (m) *m = ds->m;
206: if (l) *l = ds->l;
207: if (k) *k = ds->k;
208: if (t) *t = ds->t;
209: return(0);
210: }
212: /*@
213: DSTruncate - Truncates the system represented in the DS object.
215: Logically Collective on ds
217: Input Parameters:
218: + ds - the direct solver context
219: . n - the new size
220: - trim - a flag to indicate if the factorization must be trimmed
222: Note:
223: The new size is set to n. Note that in some cases the new size could
224: be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
225: Schur form). In cases where the extra row is meaningful, the first
226: n elements are kept as the extra row for the new system.
228: If the flag trim is turned on, it resets the locked and intermediate
229: dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
230: It also cleans the extra row if being used.
232: The typical usage of trim=true is to truncate the Schur decomposition
233: at the end of a Krylov iteration. In this case, the state must be
234: changed to RAW so that DSVectors() computes eigenvectors from scratch.
236: Level: advanced
238: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType239: @*/
240: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)241: {
243: DSStateType old;
248: DSCheckAlloc(ds,1);
251: if (!ds->ops->truncate) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
252: 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);
253: PetscLogEventBegin(DS_Other,ds,0,0,0);
254: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
255: (*ds->ops->truncate)(ds,n,trim);
256: PetscFPTrapPop();
257: PetscLogEventEnd(DS_Other,ds,0,0,0);
258: PetscInfo2(ds,"Decomposition %s to size n=%D\n",trim?"trimmed":"truncated",ds->n);
259: old = ds->state;
260: ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
261: if (old!=ds->state) {
262: PetscInfo2(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]);
263: }
264: PetscObjectStateIncrease((PetscObject)ds);
265: return(0);
266: }
268: /*@
269: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
271: Not Collective
273: Input Parameters:
274: + ds - the direct solver context
275: - t - the requested matrix
277: Output Parameters:
278: + n - the number of rows
279: - m - the number of columns
281: Note:
282: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
284: Level: developer
286: .seealso: DSSetDimensions(), DSGetMat()
287: @*/
288: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)289: {
291: PetscInt rows,cols;
296: DSCheckValidMat(ds,t,2);
297: if (ds->ops->matgetsize) {
298: (*ds->ops->matgetsize)(ds,t,&rows,&cols);
299: } else {
300: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
301: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
302: cols = ds->n;
303: }
304: if (m) *m = rows;
305: if (n) *n = cols;
306: return(0);
307: }
309: /*@
310: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
312: Not Collective
314: Input Parameters:
315: + ds - the direct solver context
316: - t - the requested matrix
318: Output Parameter:
319: . flg - the Hermitian flag
321: Note:
322: Does not check the matrix values directly. The flag is set according to the
323: problem structure. For instance, in DSHEP matrix A is Hermitian.
325: Level: developer
327: .seealso: DSGetMat()
328: @*/
329: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)330: {
336: DSCheckValidMat(ds,t,2);
338: if (ds->ops->hermitian) {
339: (*ds->ops->hermitian)(ds,t,flg);
340: } else *flg = PETSC_FALSE;
341: return(0);
342: }
344: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)345: {
346: #if !defined(PETSC_USE_COMPLEX)
347: PetscScalar *A = ds->mat[DS_MAT_A];
348: #endif
351: #if !defined(PETSC_USE_COMPLEX)
352: if (A[l+(*k)+(l+(*k)-1)*ds->ld] != 0.0) {
353: if (l+(*k)<n-1) (*k)++;
354: else (*k)--;
355: }
356: #endif
357: return(0);
358: }
360: /*@
361: DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
362: to avoid breaking a 2x2 block.
364: Not Collective
366: Input Parameters:
367: + ds - the direct solver context
368: . l - the size of the locked part (set to 0 to use ds->l)
369: . n - the total matrix size (set to 0 to use ds->n)
370: - k - the wanted truncation size
372: Output Parameter:
373: . k - the possibly modified value of the truncation size
375: Notes:
376: This should be called before DSTruncate() to make sure that the truncation
377: does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
379: The total size is n (either user-provided or ds->n if 0 is passed). The
380: size where the truncation is intended is equal to l+k (where l can be
381: equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
382: at the l+k limit, the value of k is increased (or decreased) by 1.
384: Level: advanced
386: .seealso: DSTruncate(), DSSetDimensions()
387: @*/
388: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)389: {
394: DSCheckAlloc(ds,1);
398: if (ds->ops->gettruncatesize) {
399: (*ds->ops->gettruncatesize)(ds,l?l:ds->l,n?n:ds->n,k);
400: }
401: return(0);
402: }
404: /*@
405: DSGetMat - Returns a sequential dense Mat object containing the requested
406: matrix.
408: Not Collective
410: Input Parameters:
411: + ds - the direct solver context
412: - m - the requested matrix
414: Output Parameter:
415: . A - Mat object
417: Notes:
418: The Mat is created with sizes equal to the current DS dimensions (nxm),
419: then it is filled with the values that would be obtained with DSGetArray()
420: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
421: is equal to the dimension prior to truncation, see DSTruncate().
422: The communicator is always PETSC_COMM_SELF.
424: When no longer needed, the user can either destroy the matrix or call
425: DSRestoreMat(). The latter will copy back the modified values.
427: Level: advanced
429: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
430: @*/
431: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)432: {
434: PetscInt j,rows,cols,arows,acols;
435: PetscBool create=PETSC_FALSE,flg;
436: PetscScalar *pA,*M;
440: DSCheckAlloc(ds,1);
441: DSCheckValidMat(ds,m,2);
443: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
445: DSMatGetSize(ds,m,&rows,&cols);
446: if (!ds->omat[m]) create=PETSC_TRUE;
447: else {
448: MatGetSize(ds->omat[m],&arows,&acols);
449: if (arows!=rows || acols!=cols) {
450: MatDestroy(&ds->omat[m]);
451: create=PETSC_TRUE;
452: }
453: }
454: if (create) {
455: MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);
456: }
458: /* set Hermitian flag */
459: DSMatIsHermitian(ds,m,&flg);
460: MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);
462: /* copy entries */
463: PetscObjectReference((PetscObject)ds->omat[m]);
464: *A = ds->omat[m];
465: M = ds->mat[m];
466: MatDenseGetArray(*A,&pA);
467: for (j=0;j<cols;j++) {
468: PetscArraycpy(pA+j*rows,M+j*ds->ld,rows);
469: }
470: MatDenseRestoreArray(*A,&pA);
471: return(0);
472: }
474: /*@
475: DSRestoreMat - Restores the matrix after DSGetMat() was called.
477: Not Collective
479: Input Parameters:
480: + ds - the direct solver context
481: . m - the requested matrix
482: - A - the fetched Mat object
484: Notes:
485: A call to this function must match a previous call of DSGetMat().
486: The effect is that the contents of the Mat are copied back to the
487: DS internal array, and the matrix is destroyed.
489: It is not compulsory to call this function, the matrix obtained with
490: DSGetMat() can simply be destroyed if entries need not be copied back.
492: Level: advanced
494: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
495: @*/
496: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)497: {
499: PetscInt j,rows,cols;
500: PetscScalar *pA,*M;
504: DSCheckAlloc(ds,1);
505: DSCheckValidMat(ds,m,2);
507: if (!ds->omat[m]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSRestoreMat must match a previous call to DSGetMat");
508: 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");
510: MatGetSize(*A,&rows,&cols);
511: M = ds->mat[m];
512: MatDenseGetArray(*A,&pA);
513: for (j=0;j<cols;j++) {
514: PetscArraycpy(M+j*ds->ld,pA+j*rows,rows);
515: }
516: MatDenseRestoreArray(*A,&pA);
517: MatDestroy(A);
518: return(0);
519: }
521: /*@C
522: DSGetArray - Returns a pointer to one of the internal arrays used to
523: represent matrices. You MUST call DSRestoreArray() when you no longer
524: need to access the array.
526: Not Collective
528: Input Parameters:
529: + ds - the direct solver context
530: - m - the requested matrix
532: Output Parameter:
533: . a - pointer to the values
535: Level: advanced
537: .seealso: DSRestoreArray(), DSGetArrayReal()
538: @*/
539: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])540: {
543: DSCheckAlloc(ds,1);
544: DSCheckValidMat(ds,m,2);
546: *a = ds->mat[m];
547: CHKMEMQ;
548: return(0);
549: }
551: /*@C
552: DSRestoreArray - Restores the matrix after DSGetArray() was called.
554: Not Collective
556: Input Parameters:
557: + ds - the direct solver context
558: . m - the requested matrix
559: - a - pointer to the values
561: Level: advanced
563: .seealso: DSGetArray(), DSGetArrayReal()
564: @*/
565: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])566: {
571: DSCheckAlloc(ds,1);
572: DSCheckValidMat(ds,m,2);
574: CHKMEMQ;
575: *a = 0;
576: PetscObjectStateIncrease((PetscObject)ds);
577: return(0);
578: }
580: /*@C
581: DSGetArrayReal - Returns a pointer to one of the internal arrays used to
582: represent real matrices. You MUST call DSRestoreArrayReal() when you no longer
583: need to access the array.
585: Not Collective
587: Input Parameters:
588: + ds - the direct solver context
589: - m - the requested matrix
591: Output Parameter:
592: . a - pointer to the values
594: Level: advanced
596: .seealso: DSRestoreArrayReal(), DSGetArray()
597: @*/
598: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])599: {
602: DSCheckAlloc(ds,1);
603: DSCheckValidMatReal(ds,m,2);
605: *a = ds->rmat[m];
606: CHKMEMQ;
607: return(0);
608: }
610: /*@C
611: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
613: Not Collective
615: Input Parameters:
616: + ds - the direct solver context
617: . m - the requested matrix
618: - a - pointer to the values
620: Level: advanced
622: .seealso: DSGetArrayReal(), DSGetArray()
623: @*/
624: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])625: {
630: DSCheckAlloc(ds,1);
631: DSCheckValidMatReal(ds,m,2);
633: CHKMEMQ;
634: *a = 0;
635: PetscObjectStateIncrease((PetscObject)ds);
636: return(0);
637: }
639: /*@
640: DSSolve - Solves the problem.
642: Logically Collective on ds
644: Input Parameters:
645: + ds - the direct solver context
646: . eigr - array to store the computed eigenvalues (real part)
647: - eigi - array to store the computed eigenvalues (imaginary part)
649: Note:
650: This call brings the dense system to condensed form. No ordering
651: of the eigenvalues is enforced (for this, call DSSort() afterwards).
653: Level: intermediate
655: .seealso: DSSort(), DSStateType656: @*/
657: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])658: {
664: DSCheckAlloc(ds,1);
666: if (ds->state>=DS_STATE_CONDENSED) return(0);
667: if (!ds->ops->solve[ds->method]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
668: PetscInfo3(ds,"Starting solve with problem sizes: n=%D, l=%D, k=%D\n",ds->n,ds->l,ds->k);
669: PetscLogEventBegin(DS_Solve,ds,0,0,0);
670: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
671: (*ds->ops->solve[ds->method])(ds,eigr,eigi);
672: PetscFPTrapPop();
673: PetscLogEventEnd(DS_Solve,ds,0,0,0);
674: PetscInfo1(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
675: ds->state = DS_STATE_CONDENSED;
676: PetscObjectStateIncrease((PetscObject)ds);
677: return(0);
678: }
680: /*@C
681: DSSort - Sorts the result of DSSolve() according to a given sorting
682: criterion.
684: Logically Collective on ds
686: Input Parameters:
687: + ds - the direct solver context
688: . rr - (optional) array containing auxiliary values (real part)
689: - ri - (optional) array containing auxiliary values (imaginary part)
691: Input/Output Parameters:
692: + eigr - array containing the computed eigenvalues (real part)
693: . eigi - array containing the computed eigenvalues (imaginary part)
694: - k - (optional) number of elements in the leading group
696: Notes:
697: This routine sorts the arrays provided in eigr and eigi, and also
698: sorts the dense system stored inside ds (assumed to be in condensed form).
699: The sorting criterion is specified with DSSetSlepcSC().
701: If arrays rr and ri are provided, then a (partial) reordering based on these
702: values rather than on the eigenvalues is performed. In symmetric problems
703: a total order is obtained (parameter k is ignored), but otherwise the result
704: is sorted only partially. In this latter case, it is only guaranteed that
705: all the first k elements satisfy the comparison with any of the last n-k
706: elements. The output value of parameter k is the final number of elements in
707: the first set.
709: Level: intermediate
711: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
712: @*/
713: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)714: {
716: PetscInt i;
721: DSCheckSolved(ds,1);
724: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
725: if (!ds->ops->sort) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
726: if (!ds->sc) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
727: if (k && !rr) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
729: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
730: PetscLogEventBegin(DS_Other,ds,0,0,0);
731: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
732: (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
733: PetscFPTrapPop();
734: PetscLogEventEnd(DS_Other,ds,0,0,0);
735: PetscObjectStateIncrease((PetscObject)ds);
736: PetscInfo(ds,"Finished sorting\n");
737: return(0);
738: }
740: /*@C
741: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
742: permutation.
744: Logically Collective on ds
746: Input Parameters:
747: + ds - the direct solver context
748: - perm - permutation that indicates the new ordering
750: Input/Output Parameters:
751: + eigr - array with the reordered eigenvalues (real part)
752: - eigi - array with the reordered eigenvalues (imaginary part)
754: Notes:
755: This routine reorders the arrays provided in eigr and eigi, and also the dense
756: system stored inside ds (assumed to be in condensed form). There is no sorting
757: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
759: Level: advanced
761: .seealso: DSSolve(), DSSort()
762: @*/
763: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)764: {
770: DSCheckSolved(ds,1);
773: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
774: if (!ds->ops->sortperm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
776: PetscLogEventBegin(DS_Other,ds,0,0,0);
777: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
778: (*ds->ops->sortperm)(ds,perm,eigr,eigi);
779: PetscFPTrapPop();
780: PetscLogEventEnd(DS_Other,ds,0,0,0);
781: PetscObjectStateIncrease((PetscObject)ds);
782: PetscInfo(ds,"Finished sorting\n");
783: return(0);
784: }
786: /*@
787: DSSynchronize - Make sure that all processes have the same data, performing
788: communication if necessary.
790: Collective on ds
792: Input Parameter:
793: . ds - the direct solver context
795: Input/Output Parameters:
796: + eigr - (optional) array with the computed eigenvalues (real part)
797: - eigi - (optional) array with the computed eigenvalues (imaginary part)
799: Notes:
800: When the DS has been created with a communicator with more than one process,
801: the internal data, especially the computed matrices, may diverge in the
802: different processes. This happens when using multithreaded BLAS and may
803: cause numerical issues in some ill-conditioned problems. This function
804: performs the necessary communication among the processes so that the
805: internal data is exactly equal in all of them.
807: Depending on the parallel mode as set with DSSetParallel(), this function
808: will either do nothing or synchronize the matrices computed by DSSolve()
809: and DSSort(). The arguments eigr and eigi are typically those used in the
810: calls to DSSolve() and DSSort().
812: Level: developer
814: .seealso: DSSetParallel(), DSSolve(), DSSort()
815: @*/
816: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])817: {
819: PetscMPIInt size;
824: DSCheckAlloc(ds,1);
825: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);CHKERRMPI(ierr);
826: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
827: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
828: if (ds->ops->synchronize) {
829: (*ds->ops->synchronize)(ds,eigr,eigi);
830: }
831: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
832: PetscObjectStateIncrease((PetscObject)ds);
833: PetscInfo1(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
834: }
835: return(0);
836: }
838: /*@C
839: DSVectors - Compute vectors associated to the dense system such
840: as eigenvectors.
842: Logically Collective on ds
844: Input Parameters:
845: + ds - the direct solver context
846: - mat - the matrix, used to indicate which vectors are required
848: Input/Output Parameter:
849: - j - (optional) index of vector to be computed
851: Output Parameter:
852: . rnorm - (optional) computed residual norm
854: Notes:
855: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_VT, to
856: compute right or left eigenvectors, or left or right singular vectors,
857: respectively.
859: If NULL is passed in argument j then all vectors are computed,
860: otherwise j indicates which vector must be computed. In real non-symmetric
861: problems, on exit the index j will be incremented when a complex conjugate
862: pair is found.
864: This function can be invoked after the dense problem has been solved,
865: to get the residual norm estimate of the associated Ritz pair. In that
866: case, the relevant information is returned in rnorm.
868: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
869: be in (quasi-)triangular form, or call DSSolve() first.
871: Level: intermediate
873: .seealso: DSSolve()
874: @*/
875: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)876: {
882: DSCheckAlloc(ds,1);
884: if (mat>=DS_NUM_MAT) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
885: if (!ds->ops->vectors) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
886: if (rnorm && !j) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
887: if (!ds->mat[mat]) { DSAllocateMat_Private(ds,mat); }
888: if (!j) { PetscInfo1(ds,"Computing all vectors on %s\n",DSMatName[mat]); }
889: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
890: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
891: (*ds->ops->vectors)(ds,mat,j,rnorm);
892: PetscFPTrapPop();
893: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
894: PetscObjectStateIncrease((PetscObject)ds);
895: return(0);
896: }
898: /*@
899: DSUpdateExtraRow - Performs all necessary operations so that the extra
900: row gets up-to-date after a call to DSSolve().
902: Not Collective
904: Input Parameters:
905: . ds - the direct solver context
907: Level: advanced
909: .seealso: DSSolve(), DSSetExtraRow()
910: @*/
911: PetscErrorCode DSUpdateExtraRow(DS ds)912: {
918: DSCheckAlloc(ds,1);
919: if (!ds->ops->update) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
920: if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
921: PetscInfo(ds,"Updating extra row\n");
922: PetscLogEventBegin(DS_Other,ds,0,0,0);
923: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
924: (*ds->ops->update)(ds);
925: PetscFPTrapPop();
926: PetscLogEventEnd(DS_Other,ds,0,0,0);
927: return(0);
928: }
930: /*@
931: DSCond - Compute the inf-norm condition number of the first matrix
932: as cond(A) = norm(A)*norm(inv(A)).
934: Not Collective
936: Input Parameters:
937: + ds - the direct solver context
938: - cond - the computed condition number
940: Level: advanced
942: .seealso: DSSolve()
943: @*/
944: PetscErrorCode DSCond(DS ds,PetscReal *cond)945: {
951: DSCheckAlloc(ds,1);
953: if (!ds->ops->cond) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
954: PetscLogEventBegin(DS_Other,ds,0,0,0);
955: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
956: (*ds->ops->cond)(ds,cond);
957: PetscFPTrapPop();
958: PetscLogEventEnd(DS_Other,ds,0,0,0);
959: PetscInfo1(ds,"Computed condition number = %g\n",(double)*cond);
960: return(0);
961: }
963: /*@C
964: DSTranslateHarmonic - Computes a translation of the dense system.
966: Logically Collective on ds
968: Input Parameters:
969: + ds - the direct solver context
970: . tau - the translation amount
971: . beta - last component of vector b
972: - recover - boolean flag to indicate whether to recover or not
974: Output Parameters:
975: + g - the computed vector (optional)
976: - gamma - scale factor (optional)
978: Notes:
979: This function is intended for use in the context of Krylov methods only.
980: It computes a translation of a Krylov decomposition in order to extract
981: eigenpair approximations by harmonic Rayleigh-Ritz.
982: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
983: vector b is assumed to be beta*e_n^T.
985: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
986: the factor by which the residual of the Krylov decomposition is scaled.
988: If the recover flag is activated, the computed translation undoes the
989: translation done previously. In that case, parameter tau is ignored.
991: Level: developer
992: @*/
993: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)994: {
1000: DSCheckAlloc(ds,1);
1001: if (!ds->ops->transharm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
1002: if (recover) { PetscInfo(ds,"Undoing the translation\n"); }
1003: else { PetscInfo(ds,"Computing the translation\n"); }
1004: PetscLogEventBegin(DS_Other,ds,0,0,0);
1005: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1006: (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
1007: PetscFPTrapPop();
1008: PetscLogEventEnd(DS_Other,ds,0,0,0);
1009: ds->state = DS_STATE_RAW;
1010: PetscObjectStateIncrease((PetscObject)ds);
1011: return(0);
1012: }
1014: /*@
1015: DSTranslateRKS - Computes a modification of the dense system corresponding
1016: to an update of the shift in a rational Krylov method.
1018: Logically Collective on ds
1020: Input Parameters:
1021: + ds - the direct solver context
1022: - alpha - the translation amount
1024: Notes:
1025: This function is intended for use in the context of Krylov methods only.
1026: It takes the leading (k+1,k) submatrix of A, containing the truncated
1027: Rayleigh quotient of a Krylov-Schur relation computed from a shift
1028: sigma1 and transforms it to obtain a Krylov relation as if computed
1029: from a different shift sigma2. The new matrix is computed as
1030: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
1031: alpha = sigma1-sigma2.
1033: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
1034: Krylov basis.
1036: Level: developer
1037: @*/
1038: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)1039: {
1045: DSCheckAlloc(ds,1);
1046: if (!ds->ops->transrks) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
1047: PetscInfo1(ds,"Translating with alpha=%g\n",PetscRealPart(alpha));
1048: PetscLogEventBegin(DS_Other,ds,0,0,0);
1049: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1050: (*ds->ops->transrks)(ds,alpha);
1051: PetscFPTrapPop();
1052: PetscLogEventEnd(DS_Other,ds,0,0,0);
1053: ds->state = DS_STATE_RAW;
1054: ds->compact = PETSC_FALSE;
1055: PetscObjectStateIncrease((PetscObject)ds);
1056: return(0);
1057: }
1059: /*@
1060: DSCopyMat - Copies the contents of a sequential dense Mat object to
1061: the indicated DS matrix, or vice versa.
1063: Not Collective
1065: Input Parameters:
1066: + ds - the direct solver context
1067: . m - the requested matrix
1068: . mr - first row of m to be considered
1069: . mc - first column of m to be considered
1070: . A - Mat object
1071: . Ar - first row of A to be considered
1072: . Ac - first column of A to be considered
1073: . rows - number of rows to copy
1074: . cols - number of columns to copy
1075: - out - whether the data is copied out of the DS1077: Note:
1078: If out=true, the values of the DS matrix m are copied to A, otherwise
1079: the entries of A are copied to the DS.
1081: Level: developer
1083: .seealso: DSGetMat()
1084: @*/
1085: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)1086: {
1088: PetscInt j,mrows,mcols,arows,acols;
1089: PetscScalar *pA,*M;
1093: DSCheckAlloc(ds,1);
1095: DSCheckValidMat(ds,m,2);
1104: if (!rows || !cols) return(0);
1106: DSMatGetSize(ds,m,&mrows,&mcols);
1107: MatGetSize(A,&arows,&acols);
1108: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
1109: if (mr<0 || mr>=mrows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in m");
1110: if (mc<0 || mc>=mcols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in m");
1111: if (Ar<0 || Ar>=arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in A");
1112: if (Ac<0 || Ac>=acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in A");
1113: if (mr+rows>mrows || Ar+rows>arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of rows");
1114: if (mc+cols>mcols || Ac+cols>acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of columns");
1116: M = ds->mat[m];
1117: MatDenseGetArray(A,&pA);
1118: for (j=0;j<cols;j++) {
1119: if (out) {
1120: PetscArraycpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows);
1121: } else {
1122: PetscArraycpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows);
1123: }
1124: }
1125: MatDenseRestoreArray(A,&pA);
1126: return(0);
1127: }