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