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: {
36: *ld = ds->ld;
37: PetscFunctionReturn(0);
38: }
40: /*@
41: DSSetState - Change the state of the DS object.
43: Logically Collective on ds
45: Input Parameters:
46: + ds - the direct solver context
47: - state - the new state
49: Notes:
50: The state indicates that the dense system is in an initial state (raw),
51: in an intermediate state (such as tridiagonal, Hessenberg or
52: Hessenberg-triangular), in a condensed state (such as diagonal, Schur or
53: generalized Schur), or in a truncated state.
55: This function is normally used to return to the raw state when the
56: condensed structure is destroyed.
58: Level: advanced
60: .seealso: DSGetState()
61: @*/
62: PetscErrorCode DSSetState(DS ds,DSStateType state) 63: {
66: switch (state) {
67: case DS_STATE_RAW:
68: case DS_STATE_INTERMEDIATE:
69: case DS_STATE_CONDENSED:
70: case DS_STATE_TRUNCATED:
71: if (ds->state!=state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]);
72: ds->state = state;
73: break;
74: default: 75: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
76: }
77: PetscFunctionReturn(0);
78: }
80: /*@
81: DSGetState - Returns the current state.
83: Not Collective
85: Input Parameter:
86: . ds - the direct solver context
88: Output Parameter:
89: . state - current state
91: Level: advanced
93: .seealso: DSSetState()
94: @*/
95: PetscErrorCode DSGetState(DS ds,DSStateType *state) 96: {
99: *state = ds->state;
100: PetscFunctionReturn(0);
101: }
103: /*@
104: DSSetDimensions - Resize the matrices in the DS object.
106: Logically Collective on ds
108: Input Parameters:
109: + ds - the direct solver context
110: . n - the new size
111: . l - number of locked (inactive) leading columns
112: - k - intermediate dimension (e.g., position of arrow)
114: Notes:
115: The internal arrays are not reallocated.
117: Some DS types have additional dimensions, e.g. the number of columns
118: in DSSVD. For these, you should call a specific interface function.
120: Level: intermediate
122: .seealso: DSGetDimensions(), DSAllocate(), DSSVDSetDimensions()
123: @*/
124: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt l,PetscInt k)125: {
126: PetscInt on,ol,ok;
127: PetscBool issvd;
130: DSCheckAlloc(ds,1);
134: on = ds->n; ol = ds->l; ok = ds->k;
135: if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
136: ds->n = ds->ld;
137: } else {
139: PetscObjectTypeCompareAny((PetscObject)ds,&issvd,DSSVD,DSGSVD,""); /* SVD and GSVD have extra column instead of extra row */
141: ds->n = n;
142: }
143: ds->t = ds->n; /* truncated length equal to the new dimension */
144: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
145: ds->l = 0;
146: } else {
148: ds->l = l;
149: }
150: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
151: ds->k = ds->n/2;
152: } else {
154: ds->k = k;
155: }
156: if (on!=ds->n || ol!=ds->l || ok!=ds->k) PetscInfo(ds,"New dimensions are: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
157: PetscFunctionReturn(0);
158: }
160: /*@
161: DSGetDimensions - Returns the current dimensions.
163: Not Collective
165: Input Parameter:
166: . ds - the direct solver context
168: Output Parameters:
169: + n - the current size
170: . l - number of locked (inactive) leading columns
171: . k - intermediate dimension (e.g., position of arrow)
172: - t - truncated length
174: Note:
175: The t parameter makes sense only if DSTruncate() has been called.
176: Otherwise its value equals n.
178: Some DS types have additional dimensions, e.g. the number of columns
179: in DSSVD. For these, you should call a specific interface function.
181: Level: intermediate
183: .seealso: DSSetDimensions(), DSTruncate(), DSSVDGetDimensions()
184: @*/
185: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *l,PetscInt *k,PetscInt *t)186: {
188: DSCheckAlloc(ds,1);
189: if (n) *n = ds->n;
190: if (l) *l = ds->l;
191: if (k) *k = ds->k;
192: if (t) *t = ds->t;
193: PetscFunctionReturn(0);
194: }
196: /*@
197: DSTruncate - Truncates the system represented in the DS object.
199: Logically Collective on ds
201: Input Parameters:
202: + ds - the direct solver context
203: . n - the new size
204: - trim - a flag to indicate if the factorization must be trimmed
206: Note:
207: The new size is set to n. Note that in some cases the new size could
208: be n+1 or n-1 to avoid breaking a 2x2 diagonal block (e.g. in real
209: Schur form). In cases where the extra row is meaningful, the first
210: n elements are kept as the extra row for the new system.
212: If the flag trim is turned on, it resets the locked and intermediate
213: dimensions to zero, see DSSetDimensions(), and sets the state to RAW.
214: It also cleans the extra row if being used.
216: The typical usage of trim=true is to truncate the Schur decomposition
217: at the end of a Krylov iteration. In this case, the state must be
218: changed to RAW so that DSVectors() computes eigenvectors from scratch.
220: Level: advanced
222: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType223: @*/
224: PetscErrorCode DSTruncate(DS ds,PetscInt n,PetscBool trim)225: {
226: DSStateType old;
230: DSCheckAlloc(ds,1);
235: PetscLogEventBegin(DS_Other,ds,0,0,0);
236: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
237: (*ds->ops->truncate)(ds,n,trim);
238: PetscFPTrapPop();
239: PetscLogEventEnd(DS_Other,ds,0,0,0);
240: PetscInfo(ds,"Decomposition %s to size n=%" PetscInt_FMT "\n",trim?"trimmed":"truncated",ds->n);
241: old = ds->state;
242: ds->state = trim? DS_STATE_RAW: DS_STATE_TRUNCATED;
243: if (old!=ds->state) PetscInfo(ds,"State has changed from %s to %s\n",DSStateTypes[old],DSStateTypes[ds->state]);
244: PetscObjectStateIncrease((PetscObject)ds);
245: PetscFunctionReturn(0);
246: }
248: /*@
249: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
251: Not Collective
253: Input Parameters:
254: + ds - the direct solver context
255: - t - the requested matrix
257: Output Parameters:
258: + n - the number of rows
259: - m - the number of columns
261: Note:
262: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
264: Level: developer
266: .seealso: DSSetDimensions(), DSGetMat()
267: @*/
268: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)269: {
270: PetscInt rows,cols;
274: DSCheckValidMat(ds,t,2);
275: if (ds->ops->matgetsize) (*ds->ops->matgetsize)(ds,t,&rows,&cols);
276: else {
277: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
278: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
279: cols = ds->n;
280: }
281: if (m) *m = rows;
282: if (n) *n = cols;
283: PetscFunctionReturn(0);
284: }
286: /*@
287: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
289: Not Collective
291: Input Parameters:
292: + ds - the direct solver context
293: - t - the requested matrix
295: Output Parameter:
296: . flg - the Hermitian flag
298: Note:
299: Does not check the matrix values directly. The flag is set according to the
300: problem structure. For instance, in DSHEP matrix A is Hermitian.
302: Level: developer
304: .seealso: DSGetMat()
305: @*/
306: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)307: {
310: DSCheckValidMat(ds,t,2);
312: if (ds->ops->hermitian) (*ds->ops->hermitian)(ds,t,flg);
313: else *flg = PETSC_FALSE;
314: PetscFunctionReturn(0);
315: }
317: PetscErrorCode DSGetTruncateSize_Default(DS ds,PetscInt l,PetscInt n,PetscInt *k)318: {
319: #if !defined(PETSC_USE_COMPLEX)
320: PetscScalar *A = ds->mat[DS_MAT_A];
321: #endif
323: #if !defined(PETSC_USE_COMPLEX)
324: if (A[l+(*k)+(l+(*k)-1)*ds->ld] != 0.0) {
325: if (l+(*k)<n-1) (*k)++;
326: else (*k)--;
327: }
328: #endif
329: PetscFunctionReturn(0);
330: }
332: /*@
333: DSGetTruncateSize - Gets the correct size to be used in DSTruncate()
334: to avoid breaking a 2x2 block.
336: Not Collective
338: Input Parameters:
339: + ds - the direct solver context
340: . l - the size of the locked part (set to 0 to use ds->l)
341: . n - the total matrix size (set to 0 to use ds->n)
342: - k - the wanted truncation size
344: Output Parameter:
345: . k - the possibly modified value of the truncation size
347: Notes:
348: This should be called before DSTruncate() to make sure that the truncation
349: does not break a 2x2 block corresponding to a complex conjugate eigenvalue.
351: The total size is n (either user-provided or ds->n if 0 is passed). The
352: size where the truncation is intended is equal to l+k (where l can be
353: equal to the locked size ds->l if set to 0). Then if there is a 2x2 block
354: at the l+k limit, the value of k is increased (or decreased) by 1.
356: Level: advanced
358: .seealso: DSTruncate(), DSSetDimensions()
359: @*/
360: PetscErrorCode DSGetTruncateSize(DS ds,PetscInt l,PetscInt n,PetscInt *k)361: {
363: DSCheckAlloc(ds,1);
367: if (ds->ops->gettruncatesize) (*ds->ops->gettruncatesize)(ds,l?l:ds->l,n?n:ds->n,k);
368: PetscFunctionReturn(0);
369: }
371: /*@
372: DSGetMat - Returns a sequential dense Mat object containing the requested
373: matrix.
375: Not Collective
377: Input Parameters:
378: + ds - the direct solver context
379: - m - the requested matrix
381: Output Parameter:
382: . A - Mat object
384: Notes:
385: The Mat is created with sizes equal to the current DS dimensions (nxm),
386: then it is filled with the values that would be obtained with DSGetArray()
387: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
388: is equal to the dimension prior to truncation, see DSTruncate().
389: The communicator is always PETSC_COMM_SELF.
391: When no longer needed, the user can either destroy the matrix or call
392: DSRestoreMat(). The latter will copy back the modified values.
394: Level: advanced
396: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
397: @*/
398: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)399: {
400: PetscInt j,rows,cols,arows,acols;
401: PetscBool create=PETSC_FALSE,flg;
402: PetscScalar *pA,*M;
405: DSCheckAlloc(ds,1);
406: DSCheckValidMat(ds,m,2);
410: DSMatGetSize(ds,m,&rows,&cols);
411: if (!ds->omat[m]) create=PETSC_TRUE;
412: else {
413: MatGetSize(ds->omat[m],&arows,&acols);
414: if (arows!=rows || acols!=cols) {
415: MatDestroy(&ds->omat[m]);
416: create=PETSC_TRUE;
417: }
418: }
419: if (create) MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);
421: /* set Hermitian flag */
422: DSMatIsHermitian(ds,m,&flg);
423: MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);
425: /* copy entries */
426: PetscObjectReference((PetscObject)ds->omat[m]);
427: *A = ds->omat[m];
428: M = ds->mat[m];
429: MatDenseGetArray(*A,&pA);
430: for (j=0;j<cols;j++) PetscArraycpy(pA+j*rows,M+j*ds->ld,rows);
431: MatDenseRestoreArray(*A,&pA);
432: PetscFunctionReturn(0);
433: }
435: /*@
436: DSRestoreMat - Restores the matrix after DSGetMat() was called.
438: Not Collective
440: Input Parameters:
441: + ds - the direct solver context
442: . m - the requested matrix
443: - A - the fetched Mat object
445: Notes:
446: A call to this function must match a previous call of DSGetMat().
447: The effect is that the contents of the Mat are copied back to the
448: DS internal array, and the matrix is destroyed.
450: It is not compulsory to call this function, the matrix obtained with
451: DSGetMat() can simply be destroyed if entries need not be copied back.
453: Level: advanced
455: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
456: @*/
457: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)458: {
459: PetscInt j,rows,cols;
460: PetscScalar *pA,*M;
463: DSCheckAlloc(ds,1);
464: DSCheckValidMat(ds,m,2);
469: MatGetSize(*A,&rows,&cols);
470: M = ds->mat[m];
471: MatDenseGetArray(*A,&pA);
472: for (j=0;j<cols;j++) PetscArraycpy(M+j*ds->ld,pA+j*rows,rows);
473: MatDenseRestoreArray(*A,&pA);
474: MatDestroy(A);
475: PetscFunctionReturn(0);
476: }
478: /*@C
479: DSGetArray - Returns a pointer to one of the internal arrays used to
480: represent matrices. You MUST call DSRestoreArray() when you no longer
481: need to access the array.
483: Not Collective
485: Input Parameters:
486: + ds - the direct solver context
487: - m - the requested matrix
489: Output Parameter:
490: . a - pointer to the values
492: Level: advanced
494: .seealso: DSRestoreArray(), DSGetArrayReal()
495: @*/
496: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])497: {
499: DSCheckAlloc(ds,1);
500: DSCheckValidMat(ds,m,2);
502: *a = ds->mat[m];
503: CHKMEMQ;
504: PetscFunctionReturn(0);
505: }
507: /*@C
508: DSRestoreArray - Restores the matrix after DSGetArray() was called.
510: Not Collective
512: Input Parameters:
513: + ds - the direct solver context
514: . m - the requested matrix
515: - a - pointer to the values
517: Level: advanced
519: .seealso: DSGetArray(), DSGetArrayReal()
520: @*/
521: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])522: {
524: DSCheckAlloc(ds,1);
525: DSCheckValidMat(ds,m,2);
527: CHKMEMQ;
528: *a = 0;
529: PetscObjectStateIncrease((PetscObject)ds);
530: PetscFunctionReturn(0);
531: }
533: /*@C
534: DSGetArrayReal - Returns a pointer to one of the internal arrays used to
535: represent real matrices. You MUST call DSRestoreArrayReal() when you no longer
536: need to access the array.
538: Not Collective
540: Input Parameters:
541: + ds - the direct solver context
542: - m - the requested matrix
544: Output Parameter:
545: . a - pointer to the values
547: Level: advanced
549: .seealso: DSRestoreArrayReal(), DSGetArray()
550: @*/
551: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])552: {
554: DSCheckAlloc(ds,1);
555: DSCheckValidMatReal(ds,m,2);
557: *a = ds->rmat[m];
558: CHKMEMQ;
559: PetscFunctionReturn(0);
560: }
562: /*@C
563: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
565: Not Collective
567: Input Parameters:
568: + ds - the direct solver context
569: . m - the requested matrix
570: - a - pointer to the values
572: Level: advanced
574: .seealso: DSGetArrayReal(), DSGetArray()
575: @*/
576: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])577: {
579: DSCheckAlloc(ds,1);
580: DSCheckValidMatReal(ds,m,2);
582: CHKMEMQ;
583: *a = 0;
584: PetscObjectStateIncrease((PetscObject)ds);
585: PetscFunctionReturn(0);
586: }
588: /*@
589: DSSolve - Solves the problem.
591: Logically Collective on ds
593: Input Parameters:
594: + ds - the direct solver context
595: . eigr - array to store the computed eigenvalues (real part)
596: - eigi - array to store the computed eigenvalues (imaginary part)
598: Note:
599: This call brings the dense system to condensed form. No ordering
600: of the eigenvalues is enforced (for this, call DSSort() afterwards).
602: Level: intermediate
604: .seealso: DSSort(), DSStateType605: @*/
606: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])607: {
610: DSCheckAlloc(ds,1);
612: if (ds->state>=DS_STATE_CONDENSED) PetscFunctionReturn(0);
614: PetscInfo(ds,"Starting solve with problem sizes: n=%" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT "\n",ds->n,ds->l,ds->k);
615: PetscLogEventBegin(DS_Solve,ds,0,0,0);
616: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
617: (*ds->ops->solve[ds->method])(ds,eigr,eigi);
618: PetscFPTrapPop();
619: PetscLogEventEnd(DS_Solve,ds,0,0,0);
620: PetscInfo(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
621: ds->state = DS_STATE_CONDENSED;
622: PetscObjectStateIncrease((PetscObject)ds);
623: PetscFunctionReturn(0);
624: }
626: /*@C
627: DSSort - Sorts the result of DSSolve() according to a given sorting
628: criterion.
630: Logically Collective on ds
632: Input Parameters:
633: + ds - the direct solver context
634: . rr - (optional) array containing auxiliary values (real part)
635: - ri - (optional) array containing auxiliary values (imaginary part)
637: Input/Output Parameters:
638: + eigr - array containing the computed eigenvalues (real part)
639: . eigi - array containing the computed eigenvalues (imaginary part)
640: - k - (optional) number of elements in the leading group
642: Notes:
643: This routine sorts the arrays provided in eigr and eigi, and also
644: sorts the dense system stored inside ds (assumed to be in condensed form).
645: The sorting criterion is specified with DSSetSlepcSC().
647: If arrays rr and ri are provided, then a (partial) reordering based on these
648: values rather than on the eigenvalues is performed. In symmetric problems
649: a total order is obtained (parameter k is ignored), but otherwise the result
650: is sorted only partially. In this latter case, it is only guaranteed that
651: all the first k elements satisfy the comparison with any of the last n-k
652: elements. The output value of parameter k is the final number of elements in
653: the first set.
655: Level: intermediate
657: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
658: @*/
659: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)660: {
661: PetscInt i;
665: DSCheckSolved(ds,1);
673: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
674: PetscLogEventBegin(DS_Other,ds,0,0,0);
675: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
676: (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
677: PetscFPTrapPop();
678: PetscLogEventEnd(DS_Other,ds,0,0,0);
679: PetscObjectStateIncrease((PetscObject)ds);
680: PetscInfo(ds,"Finished sorting\n");
681: PetscFunctionReturn(0);
682: }
684: /*@C
685: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
686: permutation.
688: Logically Collective on ds
690: Input Parameters:
691: + ds - the direct solver context
692: - perm - permutation that indicates the new ordering
694: Input/Output Parameters:
695: + eigr - array with the reordered eigenvalues (real part)
696: - eigi - array with the reordered eigenvalues (imaginary part)
698: Notes:
699: This routine reorders the arrays provided in eigr and eigi, and also the dense
700: system stored inside ds (assumed to be in condensed form). There is no sorting
701: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
703: Level: advanced
705: .seealso: DSSolve(), DSSort()
706: @*/
707: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)708: {
711: DSCheckSolved(ds,1);
717: PetscLogEventBegin(DS_Other,ds,0,0,0);
718: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
719: (*ds->ops->sortperm)(ds,perm,eigr,eigi);
720: PetscFPTrapPop();
721: PetscLogEventEnd(DS_Other,ds,0,0,0);
722: PetscObjectStateIncrease((PetscObject)ds);
723: PetscInfo(ds,"Finished sorting\n");
724: PetscFunctionReturn(0);
725: }
727: /*@
728: DSSynchronize - Make sure that all processes have the same data, performing
729: communication if necessary.
731: Collective on ds
733: Input Parameter:
734: . ds - the direct solver context
736: Input/Output Parameters:
737: + eigr - (optional) array with the computed eigenvalues (real part)
738: - eigi - (optional) array with the computed eigenvalues (imaginary part)
740: Notes:
741: When the DS has been created with a communicator with more than one process,
742: the internal data, especially the computed matrices, may diverge in the
743: different processes. This happens when using multithreaded BLAS and may
744: cause numerical issues in some ill-conditioned problems. This function
745: performs the necessary communication among the processes so that the
746: internal data is exactly equal in all of them.
748: Depending on the parallel mode as set with DSSetParallel(), this function
749: will either do nothing or synchronize the matrices computed by DSSolve()
750: and DSSort(). The arguments eigr and eigi are typically those used in the
751: calls to DSSolve() and DSSort().
753: Level: developer
755: .seealso: DSSetParallel(), DSSolve(), DSSort()
756: @*/
757: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])758: {
759: PetscMPIInt size;
763: DSCheckAlloc(ds,1);
764: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
765: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
766: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
767: if (ds->ops->synchronize) (*ds->ops->synchronize)(ds,eigr,eigi);
768: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
769: PetscObjectStateIncrease((PetscObject)ds);
770: PetscInfo(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
771: }
772: PetscFunctionReturn(0);
773: }
775: /*@C
776: DSVectors - Compute vectors associated to the dense system such
777: as eigenvectors.
779: Logically Collective on ds
781: Input Parameters:
782: + ds - the direct solver context
783: - mat - the matrix, used to indicate which vectors are required
785: Input/Output Parameter:
786: . j - (optional) index of vector to be computed
788: Output Parameter:
789: . rnorm - (optional) computed residual norm
791: Notes:
792: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_V, to
793: compute right or left eigenvectors, or left or right singular vectors,
794: respectively.
796: If NULL is passed in argument j then all vectors are computed,
797: otherwise j indicates which vector must be computed. In real non-symmetric
798: problems, on exit the index j will be incremented when a complex conjugate
799: pair is found.
801: This function can be invoked after the dense problem has been solved,
802: to get the residual norm estimate of the associated Ritz pair. In that
803: case, the relevant information is returned in rnorm.
805: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
806: be in (quasi-)triangular form, or call DSSolve() first.
808: Level: intermediate
810: .seealso: DSSolve()
811: @*/
812: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)813: {
816: DSCheckAlloc(ds,1);
821: if (!ds->mat[mat]) DSAllocateMat_Private(ds,mat);
822: if (!j) PetscInfo(ds,"Computing all vectors on %s\n",DSMatName[mat]);
823: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
824: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
825: (*ds->ops->vectors)(ds,mat,j,rnorm);
826: PetscFPTrapPop();
827: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
828: PetscObjectStateIncrease((PetscObject)ds);
829: PetscFunctionReturn(0);
830: }
832: /*@
833: DSUpdateExtraRow - Performs all necessary operations so that the extra
834: row gets up-to-date after a call to DSSolve().
836: Not Collective
838: Input Parameters:
839: . ds - the direct solver context
841: Level: advanced
843: .seealso: DSSolve(), DSSetExtraRow()
844: @*/
845: PetscErrorCode DSUpdateExtraRow(DS ds)846: {
849: DSCheckAlloc(ds,1);
852: PetscInfo(ds,"Updating extra row\n");
853: PetscLogEventBegin(DS_Other,ds,0,0,0);
854: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
855: (*ds->ops->update)(ds);
856: PetscFPTrapPop();
857: PetscLogEventEnd(DS_Other,ds,0,0,0);
858: PetscFunctionReturn(0);
859: }
861: /*@
862: DSCond - Compute the inf-norm condition number of the first matrix
863: as cond(A) = norm(A)*norm(inv(A)).
865: Not Collective
867: Input Parameters:
868: + ds - the direct solver context
869: - cond - the computed condition number
871: Level: advanced
873: .seealso: DSSolve()
874: @*/
875: PetscErrorCode DSCond(DS ds,PetscReal *cond)876: {
879: DSCheckAlloc(ds,1);
882: PetscLogEventBegin(DS_Other,ds,0,0,0);
883: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
884: (*ds->ops->cond)(ds,cond);
885: PetscFPTrapPop();
886: PetscLogEventEnd(DS_Other,ds,0,0,0);
887: PetscInfo(ds,"Computed condition number = %g\n",(double)*cond);
888: PetscFunctionReturn(0);
889: }
891: /*@C
892: DSTranslateHarmonic - Computes a translation of the dense system.
894: Logically Collective on ds
896: Input Parameters:
897: + ds - the direct solver context
898: . tau - the translation amount
899: . beta - last component of vector b
900: - recover - boolean flag to indicate whether to recover or not
902: Output Parameters:
903: + g - the computed vector (optional)
904: - gamma - scale factor (optional)
906: Notes:
907: This function is intended for use in the context of Krylov methods only.
908: It computes a translation of a Krylov decomposition in order to extract
909: eigenpair approximations by harmonic Rayleigh-Ritz.
910: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
911: vector b is assumed to be beta*e_n^T.
913: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
914: the factor by which the residual of the Krylov decomposition is scaled.
916: If the recover flag is activated, the computed translation undoes the
917: translation done previously. In that case, parameter tau is ignored.
919: Level: developer
921: .seealso: DSTranslateRKS()
922: @*/
923: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)924: {
927: DSCheckAlloc(ds,1);
929: if (recover) PetscInfo(ds,"Undoing the translation\n");
930: else PetscInfo(ds,"Computing the translation\n");
931: PetscLogEventBegin(DS_Other,ds,0,0,0);
932: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
933: (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
934: PetscFPTrapPop();
935: PetscLogEventEnd(DS_Other,ds,0,0,0);
936: ds->state = DS_STATE_RAW;
937: PetscObjectStateIncrease((PetscObject)ds);
938: PetscFunctionReturn(0);
939: }
941: /*@
942: DSTranslateRKS - Computes a modification of the dense system corresponding
943: to an update of the shift in a rational Krylov method.
945: Logically Collective on ds
947: Input Parameters:
948: + ds - the direct solver context
949: - alpha - the translation amount
951: Notes:
952: This function is intended for use in the context of Krylov methods only.
953: It takes the leading (k+1,k) submatrix of A, containing the truncated
954: Rayleigh quotient of a Krylov-Schur relation computed from a shift
955: sigma1 and transforms it to obtain a Krylov relation as if computed
956: from a different shift sigma2. The new matrix is computed as
957: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
958: alpha = sigma1-sigma2.
960: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
961: Krylov basis.
963: Level: developer
965: .seealso: DSTranslateHarmonic()
966: @*/
967: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)968: {
971: DSCheckAlloc(ds,1);
973: PetscInfo(ds,"Translating with alpha=%g\n",(double)PetscRealPart(alpha));
974: PetscLogEventBegin(DS_Other,ds,0,0,0);
975: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
976: (*ds->ops->transrks)(ds,alpha);
977: PetscFPTrapPop();
978: PetscLogEventEnd(DS_Other,ds,0,0,0);
979: ds->state = DS_STATE_RAW;
980: ds->compact = PETSC_FALSE;
981: PetscObjectStateIncrease((PetscObject)ds);
982: PetscFunctionReturn(0);
983: }
985: /*@
986: DSCopyMat - Copies the contents of a sequential dense Mat object to
987: the indicated DS matrix, or vice versa.
989: Not Collective
991: Input Parameters:
992: + ds - the direct solver context
993: . m - the requested matrix
994: . mr - first row of m to be considered
995: . mc - first column of m to be considered
996: . A - Mat object
997: . Ar - first row of A to be considered
998: . Ac - first column of A to be considered
999: . rows - number of rows to copy
1000: . cols - number of columns to copy
1001: - out - whether the data is copied out of the DS1003: Note:
1004: If out=true, the values of the DS matrix m are copied to A, otherwise
1005: the entries of A are copied to the DS.
1007: Level: developer
1009: .seealso: DSGetMat()
1010: @*/
1011: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)1012: {
1013: PetscInt j,mrows,mcols,arows,acols;
1014: PetscScalar *pA,*M;
1017: DSCheckAlloc(ds,1);
1019: DSCheckValidMat(ds,m,2);
1028: if (!rows || !cols) PetscFunctionReturn(0);
1030: DSMatGetSize(ds,m,&mrows,&mcols);
1031: MatGetSize(A,&arows,&acols);
1040: M = ds->mat[m];
1041: MatDenseGetArray(A,&pA);
1042: for (j=0;j<cols;j++) {
1043: if (out) PetscArraycpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows);
1044: else PetscArraycpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows);
1045: }
1046: MatDenseRestoreArray(A,&pA);
1047: PetscFunctionReturn(0);
1048: }