1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcds.h" I*/
16: /*@
17: DSGetLeadingDimension - Returns the leading dimension of the allocated
18: matrices.
20: Not Collective
22: Input Parameter:
23: . ds - the direct solver context
25: Output Parameter:
26: . ld - leading dimension (maximum allowed dimension for the matrices)
28: Level: advanced
30: .seealso: DSAllocate(), DSSetDimensions()
31: @*/
32: PetscErrorCode DSGetLeadingDimension(DS ds,PetscInt *ld) 33: {
37: *ld = ds->ld;
38: return(0);
39: }
41: /*@
42: DSSetState - Change the state of the DS object.
44: Logically Collective on DS 46: Input Parameters:
47: + ds - the direct solver context
48: - state - the new state
50: Notes:
51: The state indicates that the dense system is in an initial state (raw),
52: in an intermediate state (such as tridiagonal, Hessenberg or
53: Hessenberg-triangular), in a condensed state (such as diagonal, Schur or
54: generalized Schur), or in a truncated state.
56: This function is normally used to return to the raw state when the
57: condensed structure is destroyed.
59: Level: advanced
61: .seealso: DSGetState()
62: @*/
63: PetscErrorCode DSSetState(DS ds,DSStateType state) 64: {
70: switch (state) {
71: case DS_STATE_RAW:
72: case DS_STATE_INTERMEDIATE:
73: case DS_STATE_CONDENSED:
74: case DS_STATE_TRUNCATED:
75: if (ds->state<state) { PetscInfo(ds,"DS state has been increased\n"); }
76: ds->state = state;
77: break;
78: default: 79: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
80: }
81: return(0);
82: }
84: /*@
85: DSGetState - Returns the current state.
87: Not Collective
89: Input Parameter:
90: . ds - the direct solver context
92: Output Parameter:
93: . state - current state
95: Level: advanced
97: .seealso: DSSetState()
98: @*/
99: PetscErrorCode DSGetState(DS ds,DSStateType *state)100: {
104: *state = ds->state;
105: return(0);
106: }
108: /*@
109: DSSetDimensions - Resize the matrices in the DS object.
111: Logically Collective on DS113: Input Parameters:
114: + ds - the direct solver context
115: . n - the new size
116: . m - the new column size (only for DSSVD)
117: . l - number of locked (inactive) leading columns
118: - k - intermediate dimension (e.g., position of arrow)
120: Notes:
121: The internal arrays are not reallocated.
123: The value m is not used except in the case of DSSVD, pass 0 otherwise.
125: Level: intermediate
127: .seealso: DSGetDimensions(), DSAllocate()
128: @*/
129: PetscErrorCode DSSetDimensions(DS ds,PetscInt n,PetscInt m,PetscInt l,PetscInt k)130: {
133: DSCheckAlloc(ds,1);
138: if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
139: ds->n = ds->ld;
140: } else {
141: if (n<0 || n>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
142: if (ds->extrarow && n+1>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"A value of n equal to ld leaves no room for extra row");
143: ds->n = n;
144: }
145: ds->t = ds->n; /* truncated length equal to the new dimension */
146: if (m) {
147: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
148: ds->m = ds->ld;
149: } else {
150: if (m<0 || m>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 0 and ld");
151: ds->m = m;
152: }
153: }
154: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
155: ds->l = 0;
156: } else {
157: if (l<0 || l>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
158: ds->l = l;
159: }
160: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
161: ds->k = ds->n/2;
162: } else {
163: if (k<0 || k>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
164: ds->k = k;
165: }
166: return(0);
167: }
169: /*@
170: DSGetDimensions - Returns the current dimensions.
172: Not Collective
174: Input Parameter:
175: . ds - the direct solver context
177: Output Parameter:
178: + n - the current size
179: . m - the current column size (only for DSSVD)
180: . l - number of locked (inactive) leading columns
181: . k - intermediate dimension (e.g., position of arrow)
182: - t - truncated length
184: Note:
185: The t parameter makes sense only if DSTruncate() has been called.
186: Otherwise its value equals n.
188: Level: intermediate
190: .seealso: DSSetDimensions(), DSTruncate()
191: @*/
192: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *m,PetscInt *l,PetscInt *k,PetscInt *t)193: {
196: DSCheckAlloc(ds,1);
197: if (n) *n = ds->n;
198: if (m) *m = ds->m;
199: if (l) *l = ds->l;
200: if (k) *k = ds->k;
201: if (t) *t = ds->t;
202: return(0);
203: }
205: /*@
206: DSTruncate - Truncates the system represented in the DS object.
208: Logically Collective on DS210: Input Parameters:
211: + ds - the direct solver context
212: - n - the new size
214: Note:
215: The new size is set to n. In cases where the extra row is meaningful,
216: the first n elements are kept as the extra row for the new system.
218: Level: advanced
220: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType221: @*/
222: PetscErrorCode DSTruncate(DS ds,PetscInt n)223: {
229: DSCheckAlloc(ds,1);
230: DSCheckSolved(ds,1);
232: if (!ds->ops->truncate) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
233: if (n<ds->l || n>ds->n) SETERRQ3(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n (%D). Must be between l (%D) and n (%D)",n,ds->l,ds->n);
234: PetscLogEventBegin(DS_Other,ds,0,0,0);
235: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
236: (*ds->ops->truncate)(ds,n);
237: PetscFPTrapPop();
238: PetscLogEventEnd(DS_Other,ds,0,0,0);
239: ds->state = DS_STATE_TRUNCATED;
240: PetscObjectStateIncrease((PetscObject)ds);
241: return(0);
242: }
244: /*@
245: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
247: Not Collective
249: Input Parameters:
250: + ds - the direct solver context
251: - t - the requested matrix
253: Output Parameters:
254: + n - the number of rows
255: - m - the number of columns
257: Note:
258: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
260: Level: developer
262: .seealso: DSSetDimensions(), DSGetMat()
263: @*/
264: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)265: {
267: PetscInt rows,cols;
272: DSCheckValidMat(ds,t,2);
273: if (ds->ops->matgetsize) {
274: (*ds->ops->matgetsize)(ds,t,&rows,&cols);
275: } else {
276: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
277: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
278: cols = ds->n;
279: }
280: if (m) *m = rows;
281: if (n) *n = cols;
282: return(0);
283: }
285: /*@
286: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
288: Not Collective
290: Input Parameters:
291: + ds - the direct solver context
292: - t - the requested matrix
294: Output Parameter:
295: . flg - the Hermitian flag
297: Note:
298: Does not check the matrix values directly. The flag is set according to the
299: problem structure. For instance, in DSHEP matrix A is Hermitian.
301: Level: developer
303: .seealso: DSGetMat()
304: @*/
305: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)306: {
312: DSCheckValidMat(ds,t,2);
314: if (ds->ops->hermitian) {
315: (*ds->ops->hermitian)(ds,t,flg);
316: } else *flg = PETSC_FALSE;
317: return(0);
318: }
320: /*@
321: DSGetMat - Returns a sequential dense Mat object containing the requested
322: matrix.
324: Not Collective
326: Input Parameters:
327: + ds - the direct solver context
328: - m - the requested matrix
330: Output Parameter:
331: . A - Mat object
333: Notes:
334: The Mat is created with sizes equal to the current DS dimensions (nxm),
335: then it is filled with the values that would be obtained with DSGetArray()
336: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
337: is equal to the dimension prior to truncation, see DSTruncate().
338: The communicator is always PETSC_COMM_SELF.
340: When no longer needed, the user can either destroy the matrix or call
341: DSRestoreMat(). The latter will copy back the modified values.
343: Level: advanced
345: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
346: @*/
347: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)348: {
350: PetscInt j,rows,cols,arows,acols;
351: PetscBool create=PETSC_FALSE,flg;
352: PetscScalar *pA,*M;
356: DSCheckAlloc(ds,1);
357: DSCheckValidMat(ds,m,2);
359: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
361: DSMatGetSize(ds,m,&rows,&cols);
362: if (!ds->omat[m]) create=PETSC_TRUE;
363: else {
364: MatGetSize(ds->omat[m],&arows,&acols);
365: if (arows!=rows || acols!=cols) {
366: MatDestroy(&ds->omat[m]);
367: create=PETSC_TRUE;
368: }
369: }
370: if (create) {
371: MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);
372: }
374: /* set Hermitian flag */
375: DSMatIsHermitian(ds,m,&flg);
376: MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);
378: /* copy entries */
379: PetscObjectReference((PetscObject)ds->omat[m]);
380: *A = ds->omat[m];
381: M = ds->mat[m];
382: MatDenseGetArray(*A,&pA);
383: for (j=0;j<cols;j++) {
384: PetscMemcpy(pA+j*rows,M+j*ds->ld,rows*sizeof(PetscScalar));
385: }
386: MatDenseRestoreArray(*A,&pA);
387: return(0);
388: }
390: /*@
391: DSRestoreMat - Restores the matrix after DSGetMat() was called.
393: Not Collective
395: Input Parameters:
396: + ds - the direct solver context
397: . m - the requested matrix
398: - A - the fetched Mat object
400: Notes:
401: A call to this function must match a previous call of DSGetMat().
402: The effect is that the contents of the Mat are copied back to the
403: DS internal array, and the matrix is destroyed.
405: It is not compulsory to call this function, the matrix obtained with
406: DSGetMat() can simply be destroyed if entries need not be copied back.
408: Level: advanced
410: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
411: @*/
412: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)413: {
415: PetscInt j,rows,cols;
416: PetscScalar *pA,*M;
420: DSCheckAlloc(ds,1);
421: DSCheckValidMat(ds,m,2);
423: if (!ds->omat[m]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSRestoreMat must match a previous call to DSGetMat");
424: if (ds->omat[m]!=*A) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with DSGetMat");
426: MatGetSize(*A,&rows,&cols);
427: M = ds->mat[m];
428: MatDenseGetArray(*A,&pA);
429: for (j=0;j<cols;j++) {
430: PetscMemcpy(M+j*ds->ld,pA+j*rows,rows*sizeof(PetscScalar));
431: }
432: MatDenseRestoreArray(*A,&pA);
433: MatDestroy(A);
434: return(0);
435: }
437: /*@C
438: DSGetArray - Returns a pointer to one of the internal arrays used to
439: represent matrices. You MUST call DSRestoreArray() when you no longer
440: need to access the array.
442: Not Collective
444: Input Parameters:
445: + ds - the direct solver context
446: - m - the requested matrix
448: Output Parameter:
449: . a - pointer to the values
451: Level: advanced
453: .seealso: DSRestoreArray(), DSGetArrayReal()
454: @*/
455: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])456: {
459: DSCheckAlloc(ds,1);
460: DSCheckValidMat(ds,m,2);
462: *a = ds->mat[m];
463: CHKMEMQ;
464: return(0);
465: }
467: /*@C
468: DSRestoreArray - Restores the matrix after DSGetArray() was called.
470: Not Collective
472: Input Parameters:
473: + ds - the direct solver context
474: . m - the requested matrix
475: - a - pointer to the values
477: Level: advanced
479: .seealso: DSGetArray(), DSGetArrayReal()
480: @*/
481: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])482: {
487: DSCheckAlloc(ds,1);
488: DSCheckValidMat(ds,m,2);
490: CHKMEMQ;
491: *a = 0;
492: PetscObjectStateIncrease((PetscObject)ds);
493: return(0);
494: }
496: /*@C
497: DSGetArrayReal - Returns a pointer to one of the internal arrays used to
498: represent real matrices. You MUST call DSRestoreArrayReal() when you no longer
499: need to access the array.
501: Not Collective
503: Input Parameters:
504: + ds - the direct solver context
505: - m - the requested matrix
507: Output Parameter:
508: . a - pointer to the values
510: Level: advanced
512: .seealso: DSRestoreArrayReal(), DSGetArray()
513: @*/
514: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])515: {
518: DSCheckAlloc(ds,1);
519: DSCheckValidMatReal(ds,m,2);
521: *a = ds->rmat[m];
522: CHKMEMQ;
523: return(0);
524: }
526: /*@C
527: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
529: Not Collective
531: Input Parameters:
532: + ds - the direct solver context
533: . m - the requested matrix
534: - a - pointer to the values
536: Level: advanced
538: .seealso: DSGetArrayReal(), DSGetArray()
539: @*/
540: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])541: {
546: DSCheckAlloc(ds,1);
547: DSCheckValidMatReal(ds,m,2);
549: CHKMEMQ;
550: *a = 0;
551: PetscObjectStateIncrease((PetscObject)ds);
552: return(0);
553: }
555: /*@
556: DSSolve - Solves the problem.
558: Logically Collective on DS560: Input Parameters:
561: + ds - the direct solver context
562: . eigr - array to store the computed eigenvalues (real part)
563: - eigi - array to store the computed eigenvalues (imaginary part)
565: Note:
566: This call brings the dense system to condensed form. No ordering
567: of the eigenvalues is enforced (for this, call DSSort() afterwards).
569: Level: intermediate
571: .seealso: DSSort(), DSStateType572: @*/
573: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])574: {
580: DSCheckAlloc(ds,1);
582: if (ds->state>=DS_STATE_CONDENSED) return(0);
583: if (!ds->ops->solve[ds->method]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
584: PetscLogEventBegin(DS_Solve,ds,0,0,0);
585: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
586: (*ds->ops->solve[ds->method])(ds,eigr,eigi);
587: PetscFPTrapPop();
588: PetscLogEventEnd(DS_Solve,ds,0,0,0);
589: ds->state = DS_STATE_CONDENSED;
590: PetscObjectStateIncrease((PetscObject)ds);
591: return(0);
592: }
594: /*@C
595: DSSort - Sorts the result of DSSolve() according to a given sorting
596: criterion.
598: Logically Collective on DS600: Input Parameters:
601: + ds - the direct solver context
602: . rr - (optional) array containing auxiliary values (real part)
603: - ri - (optional) array containing auxiliary values (imaginary part)
605: Input/Output Parameters:
606: + eigr - array containing the computed eigenvalues (real part)
607: . eigi - array containing the computed eigenvalues (imaginary part)
608: - k - (optional) number of elements in the leading group
610: Notes:
611: This routine sorts the arrays provided in eigr and eigi, and also
612: sorts the dense system stored inside ds (assumed to be in condensed form).
613: The sorting criterion is specified with DSSetSlepcSC().
615: If arrays rr and ri are provided, then a (partial) reordering based on these
616: values rather than on the eigenvalues is performed. In symmetric problems
617: a total order is obtained (parameter k is ignored), but otherwise the result
618: is sorted only partially. In this latter case, it is only guaranteed that
619: all the first k elements satisfy the comparison with any of the last n-k
620: elements. The output value of parameter k is the final number of elements in
621: the first set.
623: Level: intermediate
625: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
626: @*/
627: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)628: {
630: PetscInt i;
635: DSCheckSolved(ds,1);
638: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
639: if (!ds->ops->sort) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
640: if (!ds->sc) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
641: if (k && !rr) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
643: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
644: PetscLogEventBegin(DS_Other,ds,0,0,0);
645: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
646: (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
647: PetscFPTrapPop();
648: PetscLogEventEnd(DS_Other,ds,0,0,0);
649: PetscObjectStateIncrease((PetscObject)ds);
650: return(0);
651: }
653: /*@C
654: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
655: permutation.
657: Logically Collective on DS659: Input Parameters:
660: + ds - the direct solver context
661: - perm - permutation that indicates the new ordering
663: Input/Output Parameters:
664: + eigr - array with the reordered eigenvalues (real part)
665: - eigi - array with the reordered eigenvalues (imaginary part)
667: Notes:
668: This routine reorders the arrays provided in eigr and eigi, and also the dense
669: system stored inside ds (assumed to be in condensed form). There is no sorting
670: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
672: Level: advanced
674: .seealso: DSSolve(), DSSort()
675: @*/
676: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)677: {
683: DSCheckSolved(ds,1);
686: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
687: if (!ds->ops->sortperm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
689: PetscLogEventBegin(DS_Other,ds,0,0,0);
690: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
691: (*ds->ops->sortperm)(ds,perm,eigr,eigi);
692: PetscFPTrapPop();
693: PetscLogEventEnd(DS_Other,ds,0,0,0);
694: PetscObjectStateIncrease((PetscObject)ds);
695: return(0);
696: }
698: /*@
699: DSSynchronize - Make sure that all processes have the same data, performing
700: communication if necessary.
702: Collective on DS704: Input Parameter:
705: + ds - the direct solver context
707: Input/Output Parameters:
708: + eigr - (optional) array with the computed eigenvalues (real part)
709: - eigi - (optional) array with the computed eigenvalues (imaginary part)
711: Notes:
712: When the DS has been created with a communicator with more than one process,
713: the internal data, especially the computed matrices, may diverge in the
714: different processes. This happens when using multithreaded BLAS and may
715: cause numerical issues in some ill-conditioned problems. This function
716: performs the necessary communication among the processes so that the
717: internal data is exactly equal in all of them.
719: Depending on the parallel mode as set with DSSetParallel(), this function
720: will either do nothing or synchronize the matrices computed by DSSolve()
721: and DSSort(). The arguments eigr and eigi are typically those used in the
722: calls to DSSolve() and DSSort().
724: Level: developer
726: .seealso: DSSetParallel(), DSSolve(), DSSort()
727: @*/
728: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])729: {
731: PetscMPIInt size;
736: DSCheckAlloc(ds,1);
737: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
738: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
739: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
740: if (ds->ops->synchronize) {
741: (*ds->ops->synchronize)(ds,eigr,eigi);
742: }
743: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
744: PetscObjectStateIncrease((PetscObject)ds);
745: }
746: return(0);
747: }
749: /*@C
750: DSVectors - Compute vectors associated to the dense system such
751: as eigenvectors.
753: Logically Collective on DS755: Input Parameters:
756: + ds - the direct solver context
757: - mat - the matrix, used to indicate which vectors are required
759: Input/Output Parameter:
760: - j - (optional) index of vector to be computed
762: Output Parameter:
763: . rnorm - (optional) computed residual norm
765: Notes:
766: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_VT, to
767: compute right or left eigenvectors, or left or right singular vectors,
768: respectively.
770: If NULL is passed in argument j then all vectors are computed,
771: otherwise j indicates which vector must be computed. In real non-symmetric
772: problems, on exit the index j will be incremented when a complex conjugate
773: pair is found.
775: This function can be invoked after the dense problem has been solved,
776: to get the residual norm estimate of the associated Ritz pair. In that
777: case, the relevant information is returned in rnorm.
779: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
780: be in (quasi-)triangular form, or call DSSolve() first.
782: Level: intermediate
784: .seealso: DSSolve()
785: @*/
786: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)787: {
793: DSCheckAlloc(ds,1);
795: if (mat>=DS_NUM_MAT) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
796: if (!ds->ops->vectors) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
797: if (rnorm && !j) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
798: if (!ds->mat[mat]) { DSAllocateMat_Private(ds,mat); }
799: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
800: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
801: (*ds->ops->vectors)(ds,mat,j,rnorm);
802: PetscFPTrapPop();
803: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
804: PetscObjectStateIncrease((PetscObject)ds);
805: return(0);
806: }
808: /*@
809: DSUpdateExtraRow - Performs all necessary operations so that the extra
810: row gets up-to-date after a call to DSSolve().
812: Not Collective
814: Input Parameters:
815: . ds - the direct solver context
817: Level: advanced
819: .seealso: DSSolve(), DSSetExtraRow()
820: @*/
821: PetscErrorCode DSUpdateExtraRow(DS ds)822: {
828: DSCheckAlloc(ds,1);
829: if (!ds->ops->update) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
830: if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
831: PetscLogEventBegin(DS_Other,ds,0,0,0);
832: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
833: (*ds->ops->update)(ds);
834: PetscFPTrapPop();
835: PetscLogEventEnd(DS_Other,ds,0,0,0);
836: return(0);
837: }
839: /*@
840: DSCond - Compute the inf-norm condition number of the first matrix
841: as cond(A) = norm(A)*norm(inv(A)).
843: Not Collective
845: Input Parameters:
846: + ds - the direct solver context
847: - cond - the computed condition number
849: Level: advanced
851: .seealso: DSSolve()
852: @*/
853: PetscErrorCode DSCond(DS ds,PetscReal *cond)854: {
860: DSCheckAlloc(ds,1);
862: if (!ds->ops->cond) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
863: PetscLogEventBegin(DS_Other,ds,0,0,0);
864: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
865: (*ds->ops->cond)(ds,cond);
866: PetscFPTrapPop();
867: PetscLogEventEnd(DS_Other,ds,0,0,0);
868: return(0);
869: }
871: /*@C
872: DSTranslateHarmonic - Computes a translation of the dense system.
874: Logically Collective on DS876: Input Parameters:
877: + ds - the direct solver context
878: . tau - the translation amount
879: . beta - last component of vector b
880: - recover - boolean flag to indicate whether to recover or not
882: Output Parameters:
883: + g - the computed vector (optional)
884: - gamma - scale factor (optional)
886: Notes:
887: This function is intended for use in the context of Krylov methods only.
888: It computes a translation of a Krylov decomposition in order to extract
889: eigenpair approximations by harmonic Rayleigh-Ritz.
890: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
891: vector b is assumed to be beta*e_n^T.
893: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
894: the factor by which the residual of the Krylov decomposition is scaled.
896: If the recover flag is activated, the computed translation undoes the
897: translation done previously. In that case, parameter tau is ignored.
899: Level: developer
900: @*/
901: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)902: {
908: DSCheckAlloc(ds,1);
909: if (!ds->ops->transharm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
910: PetscLogEventBegin(DS_Other,ds,0,0,0);
911: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
912: (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
913: PetscFPTrapPop();
914: PetscLogEventEnd(DS_Other,ds,0,0,0);
915: ds->state = DS_STATE_RAW;
916: PetscObjectStateIncrease((PetscObject)ds);
917: return(0);
918: }
920: /*@
921: DSTranslateRKS - Computes a modification of the dense system corresponding
922: to an update of the shift in a rational Krylov method.
924: Logically Collective on DS926: Input Parameters:
927: + ds - the direct solver context
928: - alpha - the translation amount
930: Notes:
931: This function is intended for use in the context of Krylov methods only.
932: It takes the leading (k+1,k) submatrix of A, containing the truncated
933: Rayleigh quotient of a Krylov-Schur relation computed from a shift
934: sigma1 and transforms it to obtain a Krylov relation as if computed
935: from a different shift sigma2. The new matrix is computed as
936: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
937: alpha = sigma1-sigma2.
939: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
940: Krylov basis.
942: Level: developer
943: @*/
944: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)945: {
951: DSCheckAlloc(ds,1);
952: if (!ds->ops->transrks) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
953: PetscLogEventBegin(DS_Other,ds,0,0,0);
954: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
955: (*ds->ops->transrks)(ds,alpha);
956: PetscFPTrapPop();
957: PetscLogEventEnd(DS_Other,ds,0,0,0);
958: ds->state = DS_STATE_RAW;
959: ds->compact = PETSC_FALSE;
960: PetscObjectStateIncrease((PetscObject)ds);
961: return(0);
962: }
964: /*@
965: DSCopyMat - Copies the contents of a sequential dense Mat object to
966: the indicated DS matrix, or vice versa.
968: Not Collective
970: Input Parameters:
971: + ds - the direct solver context
972: . m - the requested matrix
973: . mr - first row of m to be considered
974: . mc - first column of m to be considered
975: . A - Mat object
976: . Ar - first row of A to be considered
977: . Ac - first column of A to be considered
978: . rows - number of rows to copy
979: . cols - number of columns to copy
980: - out - whether the data is copied out of the DS982: Note:
983: If out=true, the values of the DS matrix m are copied to A, otherwise
984: the entries of A are copied to the DS.
986: Level: developer
988: .seealso: DSGetMat()
989: @*/
990: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)991: {
993: PetscInt j,mrows,mcols,arows,acols;
994: PetscScalar *pA,*M;
998: DSCheckAlloc(ds,1);
1000: DSCheckValidMat(ds,m,2);
1009: if (!rows || !cols) return(0);
1011: DSMatGetSize(ds,m,&mrows,&mcols);
1012: MatGetSize(A,&arows,&acols);
1013: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
1014: if (mr<0 || mr>=mrows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in m");
1015: if (mc<0 || mc>=mcols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in m");
1016: if (Ar<0 || Ar>=arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in A");
1017: if (Ac<0 || Ac>=acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in A");
1018: if (mr+rows>mrows || Ar+rows>arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of rows");
1019: if (mc+cols>mcols || Ac+cols>acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of columns");
1021: M = ds->mat[m];
1022: MatDenseGetArray(A,&pA);
1023: for (j=0;j<cols;j++) {
1024: if (out) {
1025: PetscMemcpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows*sizeof(PetscScalar));
1026: } else {
1027: PetscMemcpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows*sizeof(PetscScalar));
1028: }
1029: }
1030: MatDenseRestoreArray(A,&pA);
1031: return(0);
1032: }