Actual source code: dsops.c
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: [](sec:ds), `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: The state is automatically changed in functions such as `DSSolve()` or `DSTruncate()`.
57: This function is normally used to return to the raw state when the
58: condensed structure is destroyed, or to indicate that `DSSolve()` must
59: start with a problem that already has an intermediate form.
61: Level: advanced
63: .seealso: [](sec:ds), `DSGetState()`, `DSSolve()`, `DSTruncate()`
64: @*/
65: PetscErrorCode DSSetState(DS ds,DSStateType state)
66: {
67: PetscFunctionBegin;
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) PetscCall(PetscInfo(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: PetscFunctionReturn(PETSC_SUCCESS);
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: [](sec:ds), `DSSetState()`
98: @*/
99: PetscErrorCode DSGetState(DS ds,DSStateType *state)
100: {
101: PetscFunctionBegin;
103: PetscAssertPointer(state,2);
104: *state = ds->state;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: /*@
109: DSSetDimensions - Resize the matrices in the `DS` object.
111: Logically Collective
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: Use `PETSC_CURRENT` to leave any of the values unchanged. Use `PETSC_DETERMINE`
126: to set `n` to the leading dimension, `l` to the minimum value (0), and `k` to `n/2`.
128: Level: intermediate
130: .seealso: [](sec:ds), `DSGetDimensions()`, `DSAllocate()`, `DSSVDSetDimensions()`
131: @*/
132: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
133: {
134: PetscInt on,ol,ok;
135: PetscBool issvd;
137: PetscFunctionBegin;
139: DSCheckAlloc(ds,1);
143: on = ds->n; ol = ds->l; ok = ds->k;
144: if (n == PETSC_DETERMINE) {
145: ds->n = ds->extrarow? ds->ld-1: ds->ld;
146: } else if (n != PETSC_CURRENT) {
147: PetscCheck(n>=0 && n<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
148: PetscCall(PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,"")); /* SVD and GSVD have extra column instead of extra row */
149: 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");
150: ds->n = n;
151: }
152: ds->t = ds->n; /* truncated length equal to the new dimension */
153: if (l == PETSC_DETERMINE) {
154: ds->l = 0;
155: } else if (l != PETSC_CURRENT) {
156: PetscCheck(l>=0 && l<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
157: ds->l = l;
158: }
159: if (k == PETSC_DETERMINE) {
160: ds->k = ds->n/2;
161: } else if (k != PETSC_CURRENT) {
162: PetscCheck(k>=0 || k<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
163: ds->k = k;
164: }
165: 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));
166: PetscFunctionReturn(PETSC_SUCCESS);
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: Notes:
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: [](sec:ds), `DSSetDimensions()`, `DSTruncate()`, `DSSVDGetDimensions()`
193: @*/
194: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
195: {
196: PetscFunctionBegin;
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: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*@
207: DSTruncate - Truncates the system represented in the `DS` object.
209: Logically Collective
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
224: `DS_STATE_RAW`. It also cleans the extra row if being used.
226: The typical usage of `trim=PETSC_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 `DS_STATE_RAW` so that `DSVectors()` computes eigenvectors from scratch.
230: Level: advanced
232: .seealso: [](sec:ds), `DSSetDimensions()`, `DSSetExtraRow()`, `DSStateType`
233: @*/
234: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
235: {
236: DSStateType old;
238: PetscFunctionBegin;
241: DSCheckAlloc(ds,1);
244: 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);
245: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
246: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
247: PetscUseTypeMethod(ds,truncate,n,trim);
248: PetscCall(PetscFPTrapPop());
249: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
250: PetscCall(PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n));
251: old = ds->state;
252: ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
253: if (old!=ds->state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]));
254: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*@
259: DSMatGetSize - Returns the numbers of rows and columns of one of the `DS` matrices.
261: Not Collective
263: Input Parameters:
264: + ds - the direct solver context
265: - t - the requested matrix
267: Output Parameters:
268: + n - the number of rows
269: - m - the number of columns
271: Note:
272: This is equivalent to `MatGetSize()` on a matrix obtained with `DSGetMat()`.
274: Level: developer
276: .seealso: [](sec:ds), `DSSetDimensions()`, `DSGetMat()`
277: @*/
278: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
279: {
280: PetscInt rows,cols;
282: PetscFunctionBegin;
285: DSCheckValidMat(ds,t,2);
286: if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
287: else {
288: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
289: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
290: if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
291: else if (t==DS_MAT_D) cols = 1;
292: else cols = ds->n;
293: }
294: if (m) *m = rows;
295: if (n) *n = cols;
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: /*@
300: DSMatIsHermitian - Checks if one of the `DS` matrices is known to be Hermitian.
302: Not Collective
304: Input Parameters:
305: + ds - the direct solver context
306: - t - the requested matrix
308: Output Parameter:
309: . flg - the Hermitian flag
311: Note:
312: Does not check the matrix values directly. The flag is set according to the
313: problem structure. For instance, in `DSHEP` matrix `A` is Hermitian.
315: Level: developer
317: .seealso: [](sec:ds), `DSGetMat()`
318: @*/
319: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
320: {
321: PetscFunctionBegin;
324: DSCheckValidMat(ds,t,2);
325: PetscAssertPointer(flg,3);
326: *flg = PETSC_FALSE;
327: PetscTryTypeMethod(ds,hermitian,t,flg);
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
332: {
333: #if !defined(PETSC_USE_COMPLEX)
334: PetscScalar val;
335: #endif
337: PetscFunctionBegin;
338: #if !defined(PETSC_USE_COMPLEX)
339: PetscCall(MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val));
340: if (val != 0.0) {
341: if (l+(*k)<n-1) (*k)++;
342: else (*k)--;
343: }
344: #endif
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@
349: DSGetTruncateSize - Gets the correct size to be used in `DSTruncate()`
350: to avoid breaking a 2x2 block.
352: Not Collective
354: Input Parameters:
355: + ds - the direct solver context
356: . l - the size of the locked part (set to 0 to use `ds->l`)
357: - n - the total matrix size (set to 0 to use `ds->n`)
359: Output Parameter:
360: . k - the wanted truncation size (possibly modified)
362: Notes:
363: This should be called before `DSTruncate()` to make sure that the truncation
364: does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
366: The total size is `n` (either user-provided or `ds->n` if 0 is passed). The
367: size where the truncation is intended is equal to `l+k` (where `l` can be
368: equal to the locked size `ds->l` if set to 0). Then if there is a 2x2 block
369: at the `l+k` limit, the value of `k` is increased (or decreased) by 1.
371: Level: advanced
373: .seealso: [](sec:ds), `DSTruncate()`, `DSSetDimensions()`
374: @*/
375: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
376: {
377: PetscFunctionBegin;
379: DSCheckAlloc(ds,1);
382: PetscAssertPointer(k,4);
383: PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@
388: DSGetMat - Returns a sequential dense `Mat` object containing the requested
389: matrix.
391: Not Collective
393: Input Parameters:
394: + ds - the direct solver context
395: - m - the requested matrix
397: Output Parameter:
398: . A - Mat object
400: Notes:
401: The returned `Mat` has sizes equal to the current `DS` dimensions (see `DSSetDimensions()`),
402: and contains the values that would be obtained with `DSGetArray()`
403: (not `DSGetArrayReal()`). If the `DS` was truncated, then the number of rows
404: is equal to the dimension prior to truncation, see `DSTruncate()`.
405: The communicator is always `PETSC_COMM_SELF`.
407: It is implemented with `MatDenseGetSubMatrix()`, and when no longer needed
408: the user must call `DSRestoreMat()` which will invoke `MatDenseRestoreSubMatrix()`.
410: For matrices `DS_MAT_T` and `DS_MAT_D`, this function will return a `Mat` object
411: that cannot be used directly for computations, since it uses compact storage
412: (three and one diagonals for $T$ and $D$, respectively). In complex scalars, the
413: internal array stores real values, so it is sufficient with two columns for $T$.
415: Level: advanced
417: .seealso: [](sec:ds), `DSRestoreMat()`, `DSSetDimensions()`, `DSGetArray()`, `DSGetArrayReal()`, `DSTruncate()`, `DSGetMatAndColumn()`
418: @*/
419: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
420: {
421: PetscInt rows=0,cols=0;
422: PetscBool flg;
424: PetscFunctionBegin;
426: DSCheckAlloc(ds,1);
427: DSCheckValidMat(ds,m,2);
428: PetscAssertPointer(A,3);
430: PetscCall(DSMatGetSize(ds,m,&rows,&cols));
431: PetscCheck(rows && cols,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must call DSSetDimensions() first");
432: PetscCall(MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A));
434: /* set Hermitian flag */
435: PetscCall(DSMatIsHermitian(ds,m,&flg));
436: PetscCall(MatSetOption(*A,MAT_SYMMETRIC,flg));
437: PetscCall(MatSetOption(*A,MAT_HERMITIAN,flg));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: /*@
442: DSRestoreMat - Restores the matrix after `DSGetMat()` was called.
444: Not Collective
446: Input Parameters:
447: + ds - the direct solver context
448: . m - the requested matrix
449: - A - the fetched `Mat` object
451: Note:
452: A call to this function must match a previous call of `DSGetMat()`.
454: Level: advanced
456: .seealso: [](sec:ds), `DSGetMat()`, `DSRestoreArray()`, `DSRestoreArrayReal()`
457: @*/
458: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
459: {
460: PetscFunctionBegin;
462: DSCheckAlloc(ds,1);
463: DSCheckValidMat(ds,m,2);
464: PetscAssertPointer(A,3);
466: PetscCall(MatDenseRestoreSubMatrix(ds->omat[m],A));
467: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@
472: DSGetMatAndColumn - Returns a sequential dense `Mat` object containing the requested
473: matrix and one of its columns as a `Vec`.
475: Not Collective
477: Input Parameters:
478: + ds - the direct solver context
479: . m - the requested matrix
480: - col - the index of the requested column
482: Output Parameters:
483: + A - `Mat` object
484: - v - `Vec` object (the column)
486: Notes:
487: This calls `DSGetMat()` and then it extracts the selected column.
488: The user must call `DSRestoreMatAndColumn()` to recover the original state.
489: For matrices `DS_MAT_T` and `DS_MAT_D`, in complex scalars this function implies
490: copying from real values stored internally to scalar values in the Vec.
492: Level: advanced
494: .seealso: [](sec:ds), `DSRestoreMatAndColumn()`, `DSGetMat()`
495: @*/
496: PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
497: {
498: PetscFunctionBegin;
500: DSCheckAlloc(ds,1);
501: DSCheckValidMat(ds,m,2);
502: PetscAssertPointer(A,4);
503: PetscAssertPointer(v,5);
505: PetscCall(DSGetMat(ds,m,A));
506: if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
507: const PetscScalar *as;
508: PetscScalar *vs;
509: PetscReal *ar;
510: PetscInt i,n,lda;
511: PetscCall(MatCreateVecs(*A,NULL,v));
512: PetscCall(VecGetSize(*v,&n));
513: PetscCall(MatDenseGetLDA(*A,&lda));
514: PetscCall(MatDenseGetArrayRead(*A,&as));
515: PetscCall(VecGetArrayWrite(*v,&vs));
516: ar = (PetscReal*)as;
517: for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
518: PetscCall(VecRestoreArrayWrite(*v,&vs));
519: PetscCall(MatDenseRestoreArrayRead(*A,&as));
520: } else PetscCall(MatDenseGetColumnVec(*A,col,v));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: DSRestoreMatAndColumn - Restores the matrix and vector after `DSGetMatAndColumn()`
526: was called.
528: Not Collective
530: Input Parameters:
531: + ds - the direct solver context
532: . m - the requested matrix
533: . col - the requested column
534: . A - the fetched `Mat` object
535: - v - the fetched `Vec` object
537: Note:
538: A call to this function must match a previous call of `DSGetMatAndColumn()`.
540: Level: advanced
542: .seealso: [](sec:ds), `DSGetMatAndColumn()`, `DSRestoreMat()`
543: @*/
544: PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
545: {
546: PetscFunctionBegin;
548: DSCheckAlloc(ds,1);
549: DSCheckValidMat(ds,m,2);
550: PetscAssertPointer(A,4);
551: PetscAssertPointer(v,5);
553: if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
554: const PetscScalar *vs;
555: PetscScalar *as;
556: PetscReal *ar;
557: PetscInt i,n,lda;
558: PetscCall(VecGetSize(*v,&n));
559: PetscCall(MatDenseGetLDA(*A,&lda));
560: PetscCall(MatDenseGetArray(*A,&as));
561: PetscCall(VecGetArrayRead(*v,&vs));
562: ar = (PetscReal*)as;
563: for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
564: PetscCall(VecRestoreArrayRead(*v,&vs));
565: PetscCall(MatDenseRestoreArray(*A,&as));
566: PetscCall(VecDestroy(v));
567: } else PetscCall(MatDenseRestoreColumnVec(*A,col,v));
568: PetscCall(DSRestoreMat(ds,m,A));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@C
573: DSGetArray - Returns a pointer to the internal array of one of the
574: matrices. You MUST call `DSRestoreArray()` when you no longer
575: need to access the array.
577: Not Collective
579: Input Parameters:
580: + ds - the direct solver context
581: - m - the requested matrix
583: Output Parameter:
584: . a - pointer to the values
586: Note:
587: To get read-only access, use `DSGetMat()` followed by `MatDenseGetArrayRead()`.
589: Level: advanced
591: .seealso: [](sec:ds), `DSRestoreArray()`, `DSGetArrayReal()`
592: @*/
593: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
594: {
595: PetscFunctionBegin;
597: DSCheckAlloc(ds,1);
598: DSCheckValidMat(ds,m,2);
599: PetscAssertPointer(a,3);
600: PetscCall(MatDenseGetArray(ds->omat[m],a));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@C
605: DSRestoreArray - Restores the matrix after `DSGetArray()` 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: [](sec:ds), `DSGetArray()`, `DSGetArrayReal()`
617: @*/
618: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
619: {
620: PetscFunctionBegin;
622: DSCheckAlloc(ds,1);
623: DSCheckValidMat(ds,m,2);
624: PetscAssertPointer(a,3);
625: PetscCall(MatDenseRestoreArray(ds->omat[m],a));
626: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: /*@C
631: DSGetArrayReal - Returns a real pointer to the internal array of $T$ or $D$.
632: You MUST call `DSRestoreArrayReal()` when you no longer need to access the array.
634: Not Collective
636: Input Parameters:
637: + ds - the direct solver context
638: - m - the requested matrix
640: Output Parameter:
641: . a - pointer to the values
643: Note:
644: This function can be used only for `DS_MAT_T` and `DS_MAT_D`. These matrices always
645: store real values, even in complex scalars, that is why the returned pointer is
646: `PetscReal`.
648: Level: advanced
650: .seealso: [](sec:ds), `DSRestoreArrayReal()`, `DSGetArray()`
651: @*/
652: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
653: {
654: #if defined(PETSC_USE_COMPLEX)
655: PetscScalar *as;
656: #endif
658: PetscFunctionBegin;
660: DSCheckAlloc(ds,1);
661: DSCheckValidMatReal(ds,m,2);
662: PetscAssertPointer(a,3);
663: #if defined(PETSC_USE_COMPLEX)
664: PetscCall(MatDenseGetArray(ds->omat[m],&as));
665: *a = (PetscReal*)as;
666: #else
667: PetscCall(MatDenseGetArray(ds->omat[m],a));
668: #endif
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /*@C
673: DSRestoreArrayReal - Restores the matrix after `DSGetArrayReal()` was called.
675: Not Collective
677: Input Parameters:
678: + ds - the direct solver context
679: . m - the requested matrix
680: - a - pointer to the values
682: Level: advanced
684: .seealso: [](sec:ds), `DSGetArrayReal()`, `DSGetArray()`
685: @*/
686: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
687: {
688: #if defined(PETSC_USE_COMPLEX)
689: PetscScalar *as;
690: #endif
692: PetscFunctionBegin;
694: DSCheckAlloc(ds,1);
695: DSCheckValidMatReal(ds,m,2);
696: PetscAssertPointer(a,3);
697: #if defined(PETSC_USE_COMPLEX)
698: PetscCall(MatDenseRestoreArray(ds->omat[m],&as));
699: *a = NULL;
700: #else
701: PetscCall(MatDenseRestoreArray(ds->omat[m],a));
702: #endif
703: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@
708: DSSolve - Solve the problem.
710: Logically Collective
712: Input Parameters:
713: + ds - the direct solver context
714: . eigr - array to store the computed eigenvalues (real part)
715: - eigi - array to store the computed eigenvalues (imaginary part)
717: Notes:
718: This call brings the dense system to condensed form. No ordering
719: of the eigenvalues is enforced (for this, call `DSSort()` afterwards).
721: In some `DS` types, the arguments will hold singular values,
722: instead of eigenvalues. Singular values are always real.
724: Level: intermediate
726: .seealso: [](sec:ds), `DSSort()`, `DSStateType`
727: @*/
728: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
729: {
730: PetscFunctionBegin;
733: DSCheckAlloc(ds,1);
734: PetscAssertPointer(eigr,2);
735: if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(PETSC_SUCCESS);
736: PetscCheck(ds->ops->solve[ds->method],PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
737: PetscCall(PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
738: PetscCall(PetscLogEventBegin(DS_Solve,ds,0,0,0));
739: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
740: PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
741: PetscCall(PetscFPTrapPop());
742: PetscCall(PetscLogEventEnd(DS_Solve,ds,0,0,0));
743: PetscCall(PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]));
744: ds->state = DS_STATE_CONDENSED;
745: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /*@
750: DSSort - Sort the result of `DSSolve()` according to a given sorting
751: criterion.
753: Logically Collective
755: Input Parameters:
756: + ds - the direct solver context
757: . rr - (optional) array containing auxiliary values (real part)
758: - ri - (optional) array containing auxiliary values (imaginary part)
760: Input/Output Parameters:
761: + eigr - array containing the computed eigenvalues (real part)
762: . eigi - array containing the computed eigenvalues (imaginary part)
763: - k - (optional) number of elements in the leading group
765: Notes:
766: This routine sorts the arrays provided in `eigr` and `eigi`, and also
767: sorts the dense system stored inside `ds` (assumed to be in condensed form).
768: The sorting criterion is specified with `DSSetSlepcSC()`.
770: If arrays `rr` and `ri` are provided, then a (partial) reordering based on these
771: values rather than on the eigenvalues is performed. In symmetric problems
772: a total order is obtained (parameter `k` is ignored), but otherwise the result
773: is sorted only partially. In this latter case, it is only guaranteed that
774: all the first `k` elements satisfy the comparison with any of the last `n-k`
775: elements. The output value of parameter `k` is the final number of elements in
776: the first set.
778: Level: intermediate
780: .seealso: [](sec:ds), `DSSolve()`, `DSSetSlepcSC()`, `DSSortWithPermutation()`
781: @*/
782: PetscErrorCode DSSort(DS ds,PetscScalar eigr[],PetscScalar eigi[],PetscScalar rr[],PetscScalar ri[],PetscInt *k)
783: {
784: PetscInt i;
786: PetscFunctionBegin;
789: DSCheckSolved(ds,1);
790: PetscAssertPointer(eigr,2);
791: if (rr) PetscAssertPointer(rr,4);
792: PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
793: PetscCheck(ds->sc,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
794: PetscCheck(!k || rr,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
796: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
797: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
798: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
799: PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
800: PetscCall(PetscFPTrapPop());
801: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
802: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
803: PetscCall(PetscInfo(ds,"Finished sorting\n"));
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: /*@
808: DSSortWithPermutation - Reorders the result of `DSSolve()` according to a given
809: permutation.
811: Logically Collective
813: Input Parameters:
814: + ds - the direct solver context
815: - perm - permutation that indicates the new ordering
817: Input/Output Parameters:
818: + eigr - array with the reordered eigenvalues (real part)
819: - eigi - array with the reordered eigenvalues (imaginary part)
821: Notes:
822: This routine reorders the arrays provided in `eigr` and `eigi`, and also the dense
823: system stored inside `ds` (assumed to be in condensed form). There is no sorting
824: criterion, as opposed to `DSSort()`. Instead, the new ordering is given in argument
825: `perm`.
827: Level: advanced
829: .seealso: [](sec:ds), `DSSolve()`, `DSSort()`
830: @*/
831: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt perm[],PetscScalar eigr[],PetscScalar eigi[])
832: {
833: PetscFunctionBegin;
836: DSCheckSolved(ds,1);
837: PetscAssertPointer(perm,2);
838: PetscAssertPointer(eigr,3);
839: PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
841: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
842: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
843: PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
844: PetscCall(PetscFPTrapPop());
845: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
846: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
847: PetscCall(PetscInfo(ds,"Finished sorting\n"));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: /*@
852: DSSynchronize - Make sure that all processes have the same data, performing
853: communication if necessary.
855: Collective
857: Input Parameter:
858: . ds - the direct solver context
860: Input/Output Parameters:
861: + eigr - (optional) array with the computed eigenvalues (real part)
862: - eigi - (optional) array with the computed eigenvalues (imaginary part)
864: Notes:
865: When the `DS` has been created with a communicator with more than one process,
866: the internal data, especially the computed matrices, may diverge in the
867: different processes. This happens when using multithreaded BLAS and may
868: cause numerical issues in some ill-conditioned problems. This function
869: performs the necessary communication among the processes so that the
870: internal data is exactly equal in all of them.
872: Depending on the parallel mode as set with `DSSetParallel()`, this function
873: will either do nothing or synchronize the matrices computed by `DSSolve()`
874: and `DSSort()`. The arguments `eigr` and `eigi` are typically those used in the
875: calls to `DSSolve()` and `DSSort()`.
877: Level: developer
879: .seealso: [](sec:ds), `DSSetParallel()`, `DSSolve()`, `DSSort()`
880: @*/
881: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
882: {
883: PetscMPIInt size;
885: PetscFunctionBegin;
888: DSCheckAlloc(ds,1);
889: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
890: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
891: PetscCall(PetscLogEventBegin(DS_Synchronize,ds,0,0,0));
892: PetscUseTypeMethod(ds,synchronize,eigr,eigi);
893: PetscCall(PetscLogEventEnd(DS_Synchronize,ds,0,0,0));
894: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
895: PetscCall(PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]));
896: }
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*@
901: DSVectors - Compute vectors associated to the dense system such
902: as eigenvectors.
904: Logically Collective
906: Input Parameters:
907: + ds - the direct solver context
908: - mat - the matrix, used to indicate which vectors are required
910: Input/Output Parameter:
911: . j - (optional) index of vector to be computed
913: Output Parameter:
914: . rnorm - (optional) computed residual norm
916: Notes:
917: Allowed values for mat are `DS_MAT_X`, `DS_MAT_Y`, `DS_MAT_U` and `DS_MAT_V`, to
918: compute right or left eigenvectors, or left or right singular vectors,
919: respectively.
921: If `NULL` is passed in argument `j` then all vectors are computed,
922: otherwise `j` indicates which vector must be computed. In real non-symmetric
923: problems, on exit the index `j` will be incremented when a complex conjugate
924: pair is found.
926: This function can be invoked after the dense problem has been solved,
927: to get the residual norm estimate of the associated Ritz pair. In that
928: case, the relevant information is returned in `rnorm`.
930: For computing eigenvectors, LAPACK's `_trevc` is used so the matrix must
931: be in (quasi-)triangular form, or call `DSSolve()` first.
933: Level: intermediate
935: .seealso: [](sec:ds), `DSSolve()`
936: @*/
937: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
938: {
939: PetscFunctionBegin;
942: DSCheckAlloc(ds,1);
944: PetscCheck(mat<DS_NUM_MAT,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
945: PetscCheck(!rnorm || j,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
946: if (!ds->omat[mat]) PetscCall(DSAllocateMat_Private(ds,mat));
947: if (!j) PetscCall(PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]));
948: PetscCall(PetscLogEventBegin(DS_Vectors,ds,0,0,0));
949: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
950: PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
951: PetscCall(PetscFPTrapPop());
952: PetscCall(PetscLogEventEnd(DS_Vectors,ds,0,0,0));
953: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: /*@
958: DSUpdateExtraRow - Performs all necessary operations so that the extra
959: row gets up-to-date after a call to `DSSolve()`.
961: Logically Collective
963: Input Parameter:
964: . ds - the direct solver context
966: Level: advanced
968: .seealso: [](sec:ds), `DSSolve()`, `DSSetExtraRow()`
969: @*/
970: PetscErrorCode DSUpdateExtraRow(DS ds)
971: {
972: PetscFunctionBegin;
975: DSCheckAlloc(ds,1);
976: PetscCheck(ds->extrarow,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
977: PetscCall(PetscInfo(ds,"Updating extra row\n"));
978: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
979: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
980: PetscUseTypeMethod(ds,update);
981: PetscCall(PetscFPTrapPop());
982: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
983: PetscFunctionReturn(PETSC_SUCCESS);
984: }
986: /*@
987: DSCond - Compute the condition number.
989: Logically Collective
991: Input Parameters:
992: + ds - the direct solver context
993: - cond - the computed condition number
995: Notes:
996: In standard eigenvalue problems, returns the $\infty$-norm condition number of the first
997: matrix, computed as $\kappa_\infty(A) = \|A\|_\infty\|A^{-1}\|_\infty$.
999: In GSVD problems, returns the maximum of $\kappa_2(A)$ and $\kappa_2(B)$, where $\kappa_2(\cdot)$ is
1000: computed as the ratio of the largest and smallest singular values.
1002: Does not take into account the extra row.
1004: Level: advanced
1006: .seealso: [](sec:ds), `DSSolve()`, `DSSetExtraRow()`
1007: @*/
1008: PetscErrorCode DSCond(DS ds,PetscReal *cond)
1009: {
1010: PetscFunctionBegin;
1013: DSCheckAlloc(ds,1);
1014: PetscAssertPointer(cond,2);
1015: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1016: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1017: PetscUseTypeMethod(ds,cond,cond);
1018: PetscCall(PetscFPTrapPop());
1019: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1020: PetscCall(PetscInfo(ds,"Computed condition number = %g\n",(double)*cond));
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: /*@
1025: DSTranslateHarmonic - Computes a translation of the dense system.
1027: Logically Collective
1029: Input Parameters:
1030: + ds - the direct solver context
1031: . tau - the translation amount
1032: . beta - last component of vector $b$
1033: - recover - boolean flag to indicate whether to recover or not
1035: Output Parameters:
1036: + g - the computed vector (optional)
1037: - gamma - scale factor (optional)
1039: Notes:
1040: This function is intended for use in the context of Krylov methods only.
1041: It computes a translation of a Krylov decomposition in order to extract
1042: eigenpair approximations by harmonic Rayleigh-Ritz.
1043: The matrix is updated as $A + gb^*$ where $g = (A-\tau I)^{-*}b$ and
1044: vector $b$ is assumed to be $\beta e_n^*$.
1046: The $\gamma$ factor is defined as $\sqrt{1+g^*g}$ and can be interpreted as
1047: the factor by which the residual of the Krylov decomposition is scaled.
1049: If the `recover` flag is activated, the computed translation undoes the
1050: translation done previously. In that case, parameter `tau` is ignored.
1052: Level: developer
1054: .seealso: [](sec:ds), `DSTranslateRKS()`
1055: @*/
1056: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar g[],PetscReal *gamma)
1057: {
1058: PetscFunctionBegin;
1061: DSCheckAlloc(ds,1);
1062: if (recover) PetscCall(PetscInfo(ds,"Undoing the translation\n"));
1063: else PetscCall(PetscInfo(ds,"Computing the translation\n"));
1064: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1065: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1066: PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
1067: PetscCall(PetscFPTrapPop());
1068: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1069: ds->state = DS_STATE_RAW;
1070: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*@
1075: DSTranslateRKS - Computes a modification of the dense system corresponding
1076: to an update of the shift in a rational Krylov method.
1078: Logically Collective
1080: Input Parameters:
1081: + ds - the direct solver context
1082: - alpha - the translation amount
1084: Notes:
1085: This function is intended for use in the context of Krylov methods only.
1086: It takes the leading $(k+1,k)$ submatrix of $A$, containing the truncated
1087: Rayleigh quotient of a Krylov-Schur relation computed from a shift
1088: $\sigma_1$ and transforms it to obtain a Krylov relation as if computed
1089: from a different shift $\sigma_2$. The new matrix is computed as
1090: $\alpha^{-1}(I-QR^{-1})$, where $QR=I-\alpha A$ and
1091: $\alpha = \sigma_1-\sigma_2$.
1093: Matrix $Q$ is placed in `DS_MAT_Q` so that it can be used to update the
1094: Krylov basis.
1096: Level: developer
1098: .seealso: [](sec:ds), `DSTranslateHarmonic()`
1099: @*/
1100: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
1101: {
1102: PetscFunctionBegin;
1105: DSCheckAlloc(ds,1);
1106: PetscCall(PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha)));
1107: PetscCall(DSSetCompact(ds,PETSC_FALSE));
1108: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1109: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1110: PetscUseTypeMethod(ds,transrks,alpha);
1111: PetscCall(PetscFPTrapPop());
1112: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1113: ds->state = DS_STATE_RAW;
1114: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }