Actual source code: dsops.c
slepc-main 2024-12-17
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: Use PETSC_CURRENT to leave any of the values unchanged. Use PETSC_DETERMINE
124: to set n to the leading dimension, l to the minimum value (0), and k to n/2.
126: Level: intermediate
128: .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
129: @*/
130: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)
131: {
132: PetscInt on,ol,ok;
133: PetscBool issvd;
135: PetscFunctionBegin;
137: DSCheckAlloc(ds,1);
141: on = ds->n; ol = ds->l; ok = ds->k;
142: if (n == PETSC_DETERMINE) {
143: ds->n = ds->extrarow? ds->ld-1: ds->ld;
144: } else if (n != PETSC_CURRENT) {
145: PetscCheck(n>=0 && n<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
146: PetscCall(PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,"")); /* SVD and GSVD have extra column instead of extra row */
147: 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");
148: ds->n = n;
149: }
150: ds->t = ds->n; /* truncated length equal to the new dimension */
151: if (l == PETSC_DETERMINE) {
152: ds->l = 0;
153: } else if (l != PETSC_CURRENT) {
154: PetscCheck(l>=0 && l<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
155: ds->l = l;
156: }
157: if (k == PETSC_DETERMINE) {
158: ds->k = ds->n/2;
159: } else if (k != PETSC_CURRENT) {
160: PetscCheck(k>=0 || k<=ds->n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
161: ds->k = k;
162: }
163: if (on!=ds->n || ol!=ds->l || ok!=ds->k) PetscCall(PetscInfo(ds,"New dimensions are: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*@
168: DSGetDimensions - Returns the current dimensions.
170: Not Collective
172: Input Parameter:
173: . ds - the direct solver context
175: Output Parameters:
176: + n - the current size
177: . l - number of locked (inactive) leading columns
178: . k - intermediate dimension (e.g., position of arrow)
179: - t - truncated length
181: Note:
182: The t parameter makes sense only if DSTruncate() has been called.
183: Otherwise its value equals n.
185: Some DS types have additional dimensions, e.g. the number of columns
186: in DSSVD. For these, you should call a specific interface function.
188: Level: intermediate
190: .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
191: @*/
192: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)
193: {
194: PetscFunctionBegin;
196: DSCheckAlloc(ds,1);
197: if (n) *n = ds->n;
198: if (l) *l = ds->l;
199: if (k) *k = ds->k;
200: if (t) *t = ds->t;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@
205: DSTruncate - Truncates the system represented in the DS object.
207: Logically Collective
209: Input Parameters:
210: + ds - the direct solver context
211: . n - the new size
212: - trim - a flag to indicate if the factorization must be trimmed
214: Note:
215: The new size is set to n. Note that in some cases the new size could
216: be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
217: Schur form). In cases where the extra row is meaningful, the first
218: n elements are kept as the extra row for the new system.
220: If the flag trim is turned on, it resets the locked and intermediate
221: dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
222: It also cleans the extra row if being used.
224: The typical usage of trim=true is to truncate the Schur decomposition
225: at the end of a Krylov iteration. In this case, the state must be
226: changed to RAW so that DSVectors() computes eigenvectors from scratch.
228: Level: advanced
230: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType
231: @*/
232: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)
233: {
234: DSStateType old;
236: PetscFunctionBegin;
239: DSCheckAlloc(ds,1);
242: 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);
243: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
244: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
245: PetscUseTypeMethod(ds,truncate,n,trim);
246: PetscCall(PetscFPTrapPop());
247: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
248: PetscCall(PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n));
249: old = ds->state;
250: ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
251: if (old!=ds->state) PetscCall(PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]));
252: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
259: Not Collective
261: Input Parameters:
262: + ds - the direct solver context
263: - t - the requested matrix
265: Output Parameters:
266: + n - the number of rows
267: - m - the number of columns
269: Note:
270: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
272: Level: developer
274: .seealso: DSSetDimensions(), DSGetMat()
275: @*/
276: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)
277: {
278: PetscInt rows,cols;
280: PetscFunctionBegin;
283: DSCheckValidMat(ds,t,2);
284: if (ds->ops->matgetsize) PetscUseTypeMethod(ds,matgetsize,t,&rows,&cols);
285: else {
286: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
287: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
288: if (t==DS_MAT_T) cols = PetscDefined(USE_COMPLEX)? 2: 3;
289: else if (t==DS_MAT_D) cols = 1;
290: else cols = ds->n;
291: }
292: if (m) *m = rows;
293: if (n) *n = cols;
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@
298: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
300: Not Collective
302: Input Parameters:
303: + ds - the direct solver context
304: - t - the requested matrix
306: Output Parameter:
307: . flg - the Hermitian flag
309: Note:
310: Does not check the matrix values directly. The flag is set according to the
311: problem structure. For instance, in DSHEP matrix A is Hermitian.
313: Level: developer
315: .seealso: DSGetMat()
316: @*/
317: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)
318: {
319: PetscFunctionBegin;
322: DSCheckValidMat(ds,t,2);
323: PetscAssertPointer(flg,3);
324: *flg = PETSC_FALSE;
325: PetscTryTypeMethod(ds,hermitian,t,flg);
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)
330: {
331: #if !defined(PETSC_USE_COMPLEX)
332: PetscScalar val;
333: #endif
335: PetscFunctionBegin;
336: #if !defined(PETSC_USE_COMPLEX)
337: PetscCall(MatGetValue(ds->omat[DS_MAT_A],l+(*k),l+(*k)-1,&val));
338: if (val != 0.0) {
339: if (l+(*k)<n-1) (*k)++;
340: else (*k)--;
341: }
342: #endif
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@
347: DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
348: to avoid breaking a 2x2 block.
350: Not Collective
352: Input Parameters:
353: + ds - the direct solver context
354: . l - the size of the locked part (set to 0 to use ds->l)
355: - n - the total matrix size (set to 0 to use ds->n)
357: Output Parameter:
358: . k - the wanted truncation size (possibly modified)
360: Notes:
361: This should be called before DSTruncate() to make sure that the truncation
362: does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
364: The total size is n (either user-provided or ds->n if 0 is passed). The
365: size where the truncation is intended is equal to l+k (where l can be
366: equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
367: at the l+k limit, the value of k is increased (or decreased) by 1.
369: Level: advanced
371: .seealso: DSTruncate(), DSSetDimensions()
372: @*/
373: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)
374: {
375: PetscFunctionBegin;
377: DSCheckAlloc(ds,1);
380: PetscAssertPointer(k,4);
381: PetscUseTypeMethod(ds,gettruncatesize,l?l:ds->l,n?n:ds->n,k);
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: DSGetMat - Returns a sequential dense Mat object containing the requested
387: matrix.
389: Not Collective
391: Input Parameters:
392: + ds - the direct solver context
393: - m - the requested matrix
395: Output Parameter:
396: . A - Mat object
398: Notes:
399: The returned Mat has sizes equal to the current DS dimensions (nxm),
400: and contains the values that would be obtained with DSGetArray()
401: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
402: is equal to the dimension prior to truncation, see DSTruncate().
403: The communicator is always PETSC_COMM_SELF.
405: It is implemented with MatDenseGetSubMatrix(), and when no longer needed
406: the user must call DSRestoreMat() which will invoke MatDenseRestoreSubMatrix().
408: For matrices DS_MAT_T and DS_MAT_D, this function will return a Mat object
409: that cannot be used directly for computations, since it uses compact storage
410: (three and one diagonals for T and D, respectively). In complex scalars, the
411: internal array stores real values, so it is sufficient with 2 columns for T.
413: Level: advanced
415: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
416: @*/
417: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)
418: {
419: PetscInt rows=0,cols=0;
420: PetscBool flg;
422: PetscFunctionBegin;
424: DSCheckAlloc(ds,1);
425: DSCheckValidMat(ds,m,2);
426: PetscAssertPointer(A,3);
428: PetscCall(DSMatGetSize(ds,m,&rows,&cols));
429: PetscCheck(rows && cols,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must call DSSetDimensions() first");
430: PetscCall(MatDenseGetSubMatrix(ds->omat[m],0,rows,0,cols,A));
432: /* set Hermitian flag */
433: PetscCall(DSMatIsHermitian(ds,m,&flg));
434: PetscCall(MatSetOption(*A,MAT_SYMMETRIC,flg));
435: PetscCall(MatSetOption(*A,MAT_HERMITIAN,flg));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /*@
440: DSRestoreMat - Restores the matrix after DSGetMat() was called.
442: Not Collective
444: Input Parameters:
445: + ds - the direct solver context
446: . m - the requested matrix
447: - A - the fetched Mat object
449: Note:
450: A call to this function must match a previous call of DSGetMat().
452: Level: advanced
454: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
455: @*/
456: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)
457: {
458: PetscFunctionBegin;
460: DSCheckAlloc(ds,1);
461: DSCheckValidMat(ds,m,2);
462: PetscAssertPointer(A,3);
464: PetscCall(MatDenseRestoreSubMatrix(ds->omat[m],A));
465: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: DSGetMatAndColumn - Returns a sequential dense Mat object containing the requested
471: matrix and one of its columns as a Vec.
473: Not Collective
475: Input Parameters:
476: + ds - the direct solver context
477: . m - the requested matrix
478: - col - the requested column
480: Output Parameters:
481: + A - Mat object
482: - v - Vec object (the column)
484: Notes:
485: This calls DSGetMat() and then it extracts the selected column.
486: The user must call DSRestoreMatAndColumn() to recover the original state.
487: For matrices DS_MAT_T and DS_MAT_D, in complex scalars this function implies
488: copying from real values stored internally to scalar values in the Vec.
490: Level: advanced
492: .seealso: DSRestoreMatAndColumn(), DSGetMat()
493: @*/
494: PetscErrorCode DSGetMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
495: {
496: PetscFunctionBegin;
498: DSCheckAlloc(ds,1);
499: DSCheckValidMat(ds,m,2);
500: PetscAssertPointer(A,4);
501: PetscAssertPointer(v,5);
503: PetscCall(DSGetMat(ds,m,A));
504: if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
505: const PetscScalar *as;
506: PetscScalar *vs;
507: PetscReal *ar;
508: PetscInt i,n,lda;
509: PetscCall(MatCreateVecs(*A,NULL,v));
510: PetscCall(VecGetSize(*v,&n));
511: PetscCall(MatDenseGetLDA(*A,&lda));
512: PetscCall(MatDenseGetArrayRead(*A,&as));
513: PetscCall(VecGetArrayWrite(*v,&vs));
514: ar = (PetscReal*)as;
515: for (i=0;i<n;i++) vs[i] = ar[i+col*lda];
516: PetscCall(VecRestoreArrayWrite(*v,&vs));
517: PetscCall(MatDenseRestoreArrayRead(*A,&as));
518: } else PetscCall(MatDenseGetColumnVec(*A,col,v));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: DSRestoreMatAndColumn - Restores the matrix and vector after DSGetMatAndColumn()
524: was called.
526: Not Collective
528: Input Parameters:
529: + ds - the direct solver context
530: . m - the requested matrix
531: . col - the requested column
532: . A - the fetched Mat object
533: - v - the fetched Vec object
535: Note:
536: A call to this function must match a previous call of DSGetMatAndColumn().
538: Level: advanced
540: .seealso: DSGetMatAndColumn(), DSRestoreMat()
541: @*/
542: PetscErrorCode DSRestoreMatAndColumn(DS ds,DSMatType m,PetscInt col,Mat *A,Vec *v)
543: {
544: PetscFunctionBegin;
546: DSCheckAlloc(ds,1);
547: DSCheckValidMat(ds,m,2);
548: PetscAssertPointer(A,4);
549: PetscAssertPointer(v,5);
551: if (PetscDefined(USE_COMPLEX) && (m==DS_MAT_T || m==DS_MAT_D)) {
552: const PetscScalar *vs;
553: PetscScalar *as;
554: PetscReal *ar;
555: PetscInt i,n,lda;
556: PetscCall(VecGetSize(*v,&n));
557: PetscCall(MatDenseGetLDA(*A,&lda));
558: PetscCall(MatDenseGetArray(*A,&as));
559: PetscCall(VecGetArrayRead(*v,&vs));
560: ar = (PetscReal*)as;
561: for (i=0;i<n;i++) ar[i+col*lda] = PetscRealPart(vs[i]);
562: PetscCall(VecRestoreArrayRead(*v,&vs));
563: PetscCall(MatDenseRestoreArray(*A,&as));
564: PetscCall(VecDestroy(v));
565: } else PetscCall(MatDenseRestoreColumnVec(*A,col,v));
566: PetscCall(DSRestoreMat(ds,m,A));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@C
571: DSGetArray - Returns a pointer to the internal array of one of the
572: matrices. You MUST call DSRestoreArray() when you no longer
573: need to access the array.
575: Not Collective
577: Input Parameters:
578: + ds - the direct solver context
579: - m - the requested matrix
581: Output Parameter:
582: . a - pointer to the values
584: Note:
585: To get read-only access, use DSGetMat() followed by MatDenseGetArrayRead().
587: Level: advanced
589: .seealso: DSRestoreArray(), DSGetArrayReal()
590: @*/
591: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])
592: {
593: PetscFunctionBegin;
595: DSCheckAlloc(ds,1);
596: DSCheckValidMat(ds,m,2);
597: PetscAssertPointer(a,3);
598: PetscCall(MatDenseGetArray(ds->omat[m],a));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: /*@C
603: DSRestoreArray - Restores the matrix after DSGetArray() was called.
605: Not Collective
607: Input Parameters:
608: + ds - the direct solver context
609: . m - the requested matrix
610: - a - pointer to the values
612: Level: advanced
614: .seealso: DSGetArray(), DSGetArrayReal()
615: @*/
616: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])
617: {
618: PetscFunctionBegin;
620: DSCheckAlloc(ds,1);
621: DSCheckValidMat(ds,m,2);
622: PetscAssertPointer(a,3);
623: PetscCall(MatDenseRestoreArray(ds->omat[m],a));
624: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: /*@C
629: DSGetArrayReal - Returns a real pointer to the internal array of T or D.
630: You MUST call DSRestoreArrayReal() when you no longer need to access the array.
632: Not Collective
634: Input Parameters:
635: + ds - the direct solver context
636: - m - the requested matrix
638: Output Parameter:
639: . a - pointer to the values
641: Note:
642: This function can be used only for DS_MAT_T and DS_MAT_D. These matrices always
643: store real values, even in complex scalars, that is why the returned pointer is
644: PetscReal.
646: Level: advanced
648: .seealso: DSRestoreArrayReal(), DSGetArray()
649: @*/
650: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])
651: {
652: #if defined(PETSC_USE_COMPLEX)
653: PetscScalar *as;
654: #endif
656: PetscFunctionBegin;
658: DSCheckAlloc(ds,1);
659: DSCheckValidMatReal(ds,m,2);
660: PetscAssertPointer(a,3);
661: #if defined(PETSC_USE_COMPLEX)
662: PetscCall(MatDenseGetArray(ds->omat[m],&as));
663: *a = (PetscReal*)as;
664: #else
665: PetscCall(MatDenseGetArray(ds->omat[m],a));
666: #endif
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /*@C
671: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
673: Not Collective
675: Input Parameters:
676: + ds - the direct solver context
677: . m - the requested matrix
678: - a - pointer to the values
680: Level: advanced
682: .seealso: DSGetArrayReal(), DSGetArray()
683: @*/
684: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])
685: {
686: #if defined(PETSC_USE_COMPLEX)
687: PetscScalar *as;
688: #endif
690: PetscFunctionBegin;
692: DSCheckAlloc(ds,1);
693: DSCheckValidMatReal(ds,m,2);
694: PetscAssertPointer(a,3);
695: #if defined(PETSC_USE_COMPLEX)
696: PetscCall(MatDenseRestoreArray(ds->omat[m],&as));
697: *a = NULL;
698: #else
699: PetscCall(MatDenseRestoreArray(ds->omat[m],a));
700: #endif
701: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: /*@
706: DSSolve - Solves the problem.
708: Logically Collective
710: Input Parameters:
711: + ds - the direct solver context
712: . eigr - array to store the computed eigenvalues (real part)
713: - eigi - array to store the computed eigenvalues (imaginary part)
715: Note:
716: This call brings the dense system to condensed form. No ordering
717: of the eigenvalues is enforced (for this, call DSSort() afterwards).
719: Level: intermediate
721: .seealso: DSSort(), DSStateType
722: @*/
723: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])
724: {
725: PetscFunctionBegin;
728: DSCheckAlloc(ds,1);
729: PetscAssertPointer(eigr,2);
730: if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(PETSC_SUCCESS);
731: PetscCheck(ds->ops->solve[ds->method],PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
732: PetscCall(PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k));
733: PetscCall(PetscLogEventBegin(DS_Solve,ds,0,0,0));
734: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
735: PetscUseTypeMethod(ds,solve[ds->method],eigr,eigi);
736: PetscCall(PetscFPTrapPop());
737: PetscCall(PetscLogEventEnd(DS_Solve,ds,0,0,0));
738: PetscCall(PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]));
739: ds->state = DS_STATE_CONDENSED;
740: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: /*@
745: DSSort - Sorts the result of DSSolve() according to a given sorting
746: criterion.
748: Logically Collective
750: Input Parameters:
751: + ds - the direct solver context
752: . rr - (optional) array containing auxiliary values (real part)
753: - ri - (optional) array containing auxiliary values (imaginary part)
755: Input/Output Parameters:
756: + eigr - array containing the computed eigenvalues (real part)
757: . eigi - array containing the computed eigenvalues (imaginary part)
758: - k - (optional) number of elements in the leading group
760: Notes:
761: This routine sorts the arrays provided in eigr and eigi, and also
762: sorts the dense system stored inside ds (assumed to be in condensed form).
763: The sorting criterion is specified with DSSetSlepcSC().
765: If arrays rr and ri are provided, then a (partial) reordering based on these
766: values rather than on the eigenvalues is performed. In symmetric problems
767: a total order is obtained (parameter k is ignored), but otherwise the result
768: is sorted only partially. In this latter case, it is only guaranteed that
769: all the first k elements satisfy the comparison with any of the last n-k
770: elements. The output value of parameter k is the final number of elements in
771: the first set.
773: Level: intermediate
775: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
776: @*/
777: PetscErrorCode DSSort(DS ds,PetscScalar eigr[],PetscScalar eigi[],PetscScalar rr[],PetscScalar ri[],PetscInt *k)
778: {
779: PetscInt i;
781: PetscFunctionBegin;
784: DSCheckSolved(ds,1);
785: PetscAssertPointer(eigr,2);
786: if (rr) PetscAssertPointer(rr,4);
787: PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
788: PetscCheck(ds->sc,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
789: PetscCheck(!k || rr,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
791: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
792: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
793: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
794: PetscUseTypeMethod(ds,sort,eigr,eigi,rr,ri,k);
795: PetscCall(PetscFPTrapPop());
796: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
797: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
798: PetscCall(PetscInfo(ds,"Finished sorting\n"));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: /*@
803: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
804: permutation.
806: Logically Collective
808: Input Parameters:
809: + ds - the direct solver context
810: - perm - permutation that indicates the new ordering
812: Input/Output Parameters:
813: + eigr - array with the reordered eigenvalues (real part)
814: - eigi - array with the reordered eigenvalues (imaginary part)
816: Notes:
817: This routine reorders the arrays provided in eigr and eigi, and also the dense
818: system stored inside ds (assumed to be in condensed form). There is no sorting
819: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
821: Level: advanced
823: .seealso: DSSolve(), DSSort()
824: @*/
825: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt perm[],PetscScalar eigr[],PetscScalar eigi[])
826: {
827: PetscFunctionBegin;
830: DSCheckSolved(ds,1);
831: PetscAssertPointer(perm,2);
832: PetscAssertPointer(eigr,3);
833: PetscCheck(ds->state<DS_STATE_TRUNCATED,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
835: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
836: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
837: PetscUseTypeMethod(ds,sortperm,perm,eigr,eigi);
838: PetscCall(PetscFPTrapPop());
839: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
840: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
841: PetscCall(PetscInfo(ds,"Finished sorting\n"));
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }
845: /*@
846: DSSynchronize - Make sure that all processes have the same data, performing
847: communication if necessary.
849: Collective
851: Input Parameter:
852: . ds - the direct solver context
854: Input/Output Parameters:
855: + eigr - (optional) array with the computed eigenvalues (real part)
856: - eigi - (optional) array with the computed eigenvalues (imaginary part)
858: Notes:
859: When the DS has been created with a communicator with more than one process,
860: the internal data, especially the computed matrices, may diverge in the
861: different processes. This happens when using multithreaded BLAS and may
862: cause numerical issues in some ill-conditioned problems. This function
863: performs the necessary communication among the processes so that the
864: internal data is exactly equal in all of them.
866: Depending on the parallel mode as set with DSSetParallel(), this function
867: will either do nothing or synchronize the matrices computed by DSSolve()
868: and DSSort(). The arguments eigr and eigi are typically those used in the
869: calls to DSSolve() and DSSort().
871: Level: developer
873: .seealso: DSSetParallel(), DSSolve(), DSSort()
874: @*/
875: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])
876: {
877: PetscMPIInt size;
879: PetscFunctionBegin;
882: DSCheckAlloc(ds,1);
883: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
884: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
885: PetscCall(PetscLogEventBegin(DS_Synchronize,ds,0,0,0));
886: PetscUseTypeMethod(ds,synchronize,eigr,eigi);
887: PetscCall(PetscLogEventEnd(DS_Synchronize,ds,0,0,0));
888: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
889: PetscCall(PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]));
890: }
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: /*@
895: DSVectors - Compute vectors associated to the dense system such
896: as eigenvectors.
898: Logically Collective
900: Input Parameters:
901: + ds - the direct solver context
902: - mat - the matrix, used to indicate which vectors are required
904: Input/Output Parameter:
905: . j - (optional) index of vector to be computed
907: Output Parameter:
908: . rnorm - (optional) computed residual norm
910: Notes:
911: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
912: compute right or left eigenvectors, or left or right singular vectors,
913: respectively.
915: If NULL is passed in argument j then all vectors are computed,
916: otherwise j indicates which vector must be computed. In real non-symmetric
917: problems, on exit the index j will be incremented when a complex conjugate
918: pair is found.
920: This function can be invoked after the dense problem has been solved,
921: to get the residual norm estimate of the associated Ritz pair. In that
922: case, the relevant information is returned in rnorm.
924: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
925: be in (quasi-)triangular form, or call DSSolve() first.
927: Level: intermediate
929: .seealso: DSSolve()
930: @*/
931: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
932: {
933: PetscFunctionBegin;
936: DSCheckAlloc(ds,1);
938: PetscCheck(mat<DS_NUM_MAT,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
939: PetscCheck(!rnorm || j,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
940: if (!ds->omat[mat]) PetscCall(DSAllocateMat_Private(ds,mat));
941: if (!j) PetscCall(PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]));
942: PetscCall(PetscLogEventBegin(DS_Vectors,ds,0,0,0));
943: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
944: PetscUseTypeMethod(ds,vectors,mat,j,rnorm);
945: PetscCall(PetscFPTrapPop());
946: PetscCall(PetscLogEventEnd(DS_Vectors,ds,0,0,0));
947: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
948: PetscFunctionReturn(PETSC_SUCCESS);
949: }
951: /*@
952: DSUpdateExtraRow - Performs all necessary operations so that the extra
953: row gets up-to-date after a call to DSSolve().
955: Logically Collective
957: Input Parameters:
958: . ds - the direct solver context
960: Level: advanced
962: .seealso: DSSolve(), DSSetExtraRow()
963: @*/
964: PetscErrorCode DSUpdateExtraRow(DS ds)
965: {
966: PetscFunctionBegin;
969: DSCheckAlloc(ds,1);
970: PetscCheck(ds->extrarow,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
971: PetscCall(PetscInfo(ds,"Updating extra row\n"));
972: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
973: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
974: PetscUseTypeMethod(ds,update);
975: PetscCall(PetscFPTrapPop());
976: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@
981: DSCond - Compute the condition number.
983: Logically Collective
985: Input Parameters:
986: + ds - the direct solver context
987: - cond - the computed condition number
989: Notes:
990: In standard eigenvalue problems, returns the inf-norm condition number of the first
991: matrix, computed as cond(A) = norm(A)*norm(inv(A)).
993: In GSVD problems, returns the maximum of cond(A) and cond(B), where cond(.) is
994: computed as the ratio of the largest and smallest singular values.
996: Does not take into account the extra row.
998: Level: advanced
1000: .seealso: DSSolve(), DSSetExtraRow()
1001: @*/
1002: PetscErrorCode DSCond(DS ds,PetscReal *cond)
1003: {
1004: PetscFunctionBegin;
1007: DSCheckAlloc(ds,1);
1008: PetscAssertPointer(cond,2);
1009: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1010: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1011: PetscUseTypeMethod(ds,cond,cond);
1012: PetscCall(PetscFPTrapPop());
1013: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1014: PetscCall(PetscInfo(ds,"Computed condition number = %g\n",(double)*cond));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /*@
1019: DSTranslateHarmonic - Computes a translation of the dense system.
1021: Logically Collective
1023: Input Parameters:
1024: + ds - the direct solver context
1025: . tau - the translation amount
1026: . beta - last component of vector b
1027: - recover - boolean flag to indicate whether to recover or not
1029: Output Parameters:
1030: + g - the computed vector (optional)
1031: - gamma - scale factor (optional)
1033: Notes:
1034: This function is intended for use in the context of Krylov methods only.
1035: It computes a translation of a Krylov decomposition in order to extract
1036: eigenpair approximations by harmonic Rayleigh-Ritz.
1037: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
1038: vector b is assumed to be beta*e_n^T.
1040: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
1041: the factor by which the residual of the Krylov decomposition is scaled.
1043: If the recover flag is activated, the computed translation undoes the
1044: translation done previously. In that case, parameter tau is ignored.
1046: Level: developer
1048: .seealso: DSTranslateRKS()
1049: @*/
1050: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)
1051: {
1052: PetscFunctionBegin;
1055: DSCheckAlloc(ds,1);
1056: if (recover) PetscCall(PetscInfo(ds,"Undoing the translation\n"));
1057: else PetscCall(PetscInfo(ds,"Computing the translation\n"));
1058: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1059: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1060: PetscUseTypeMethod(ds,transharm,tau,beta,recover,g,gamma);
1061: PetscCall(PetscFPTrapPop());
1062: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1063: ds->state = DS_STATE_RAW;
1064: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1065: PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1068: /*@
1069: DSTranslateRKS - Computes a modification of the dense system corresponding
1070: to an update of the shift in a rational Krylov method.
1072: Logically Collective
1074: Input Parameters:
1075: + ds - the direct solver context
1076: - alpha - the translation amount
1078: Notes:
1079: This function is intended for use in the context of Krylov methods only.
1080: It takes the leading (k+1,k) submatrix of A, containing the truncated
1081: Rayleigh quotient of a Krylov-Schur relation computed from a shift
1082: sigma1 and transforms it to obtain a Krylov relation as if computed
1083: from a different shift sigma2. The new matrix is computed as
1084: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
1085: alpha = sigma1-sigma2.
1087: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
1088: Krylov basis.
1090: Level: developer
1092: .seealso: DSTranslateHarmonic()
1093: @*/
1094: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)
1095: {
1096: PetscFunctionBegin;
1099: DSCheckAlloc(ds,1);
1100: PetscCall(PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha)));
1101: PetscCall(DSSetCompact(ds,PETSC_FALSE));
1102: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
1103: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1104: PetscUseTypeMethod(ds,transrks,alpha);
1105: PetscCall(PetscFPTrapPop());
1106: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
1107: ds->state = DS_STATE_RAW;
1108: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }