1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: {
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) { PetscInfo2(ds,"State has changed from %s to %s\n",DSStateTypes[ds->state],DSStateTypes[state]); }
76: ds->state = state;
77: break;
78: default: 79: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Wrong state");
80: }
81: 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 ds
113: 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: {
132: PetscInt on,om,ol,ok;
136: DSCheckAlloc(ds,1);
141: on = ds->n; om = ds->m; ol = ds->l; ok = ds->k;
142: if (n==PETSC_DECIDE || n==PETSC_DEFAULT) {
143: ds->n = ds->ld;
144: } else {
145: if (n<0 || n>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of n. Must be between 0 and ld");
146: 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");
147: ds->n = n;
148: }
149: ds->t = ds->n; /* truncated length equal to the new dimension */
150: if (m) {
151: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
152: ds->m = ds->ld;
153: } else {
154: if (m<0 || m>ds->ld) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 0 and ld");
155: ds->m = m;
156: }
157: }
158: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
159: ds->l = 0;
160: } else {
161: if (l<0 || l>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and n");
162: ds->l = l;
163: }
164: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
165: ds->k = ds->n/2;
166: } else {
167: if (k<0 || k>ds->n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and n");
168: ds->k = k;
169: }
170: if (on!=ds->n || om!=ds->m || ol!=ds->l || ok!=ds->k) {
171: PetscInfo4(ds,"New dimensions are: n=%D, m=%D, l=%D, k=%D\n",ds->n,ds->m,ds->l,ds->k);
172: }
173: return(0);
174: }
176: /*@
177: DSGetDimensions - Returns the current dimensions.
179: Not Collective
181: Input Parameter:
182: . ds - the direct solver context
184: Output Parameter:
185: + n - the current size
186: . m - the current column size (only for DSSVD)
187: . l - number of locked (inactive) leading columns
188: . k - intermediate dimension (e.g., position of arrow)
189: - t - truncated length
191: Note:
192: The t parameter makes sense only if DSTruncate() has been called.
193: Otherwise its value equals n.
195: Level: intermediate
197: .seealso: DSSetDimensions(), DSTruncate()
198: @*/
199: PetscErrorCode DSGetDimensions(DS ds,PetscInt *n,PetscInt *m,PetscInt *l,PetscInt *k,PetscInt *t)200: {
203: DSCheckAlloc(ds,1);
204: if (n) *n = ds->n;
205: if (m) *m = ds->m;
206: if (l) *l = ds->l;
207: if (k) *k = ds->k;
208: if (t) *t = ds->t;
209: return(0);
210: }
212: /*@
213: DSTruncate - Truncates the system represented in the DS object.
215: Logically Collective on ds
217: Input Parameters:
218: + ds - the direct solver context
219: - n - the new size
221: Note:
222: The new size is set to n. In cases where the extra row is meaningful,
223: the first n elements are kept as the extra row for the new system.
225: Level: advanced
227: .seealso: DSSetDimensions(), DSSetExtraRow(), DSStateType228: @*/
229: PetscErrorCode DSTruncate(DS ds,PetscInt n)230: {
236: DSCheckAlloc(ds,1);
237: DSCheckSolved(ds,1);
239: if (!ds->ops->truncate) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
240: 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);
241: PetscLogEventBegin(DS_Other,ds,0,0,0);
242: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
243: (*ds->ops->truncate)(ds,n);
244: PetscFPTrapPop();
245: PetscLogEventEnd(DS_Other,ds,0,0,0);
246: PetscInfo1(ds,"State has changed from %s to TRUNCATED\n",DSStateTypes[ds->state]);
247: ds->state = DS_STATE_TRUNCATED;
248: PetscObjectStateIncrease((PetscObject)ds);
249: return(0);
250: }
252: /*@
253: DSMatGetSize - Returns the numbers of rows and columns of one of the DS matrices.
255: Not Collective
257: Input Parameters:
258: + ds - the direct solver context
259: - t - the requested matrix
261: Output Parameters:
262: + n - the number of rows
263: - m - the number of columns
265: Note:
266: This is equivalent to MatGetSize() on a matrix obtained with DSGetMat().
268: Level: developer
270: .seealso: DSSetDimensions(), DSGetMat()
271: @*/
272: PetscErrorCode DSMatGetSize(DS ds,DSMatType t,PetscInt *m,PetscInt *n)273: {
275: PetscInt rows,cols;
280: DSCheckValidMat(ds,t,2);
281: if (ds->ops->matgetsize) {
282: (*ds->ops->matgetsize)(ds,t,&rows,&cols);
283: } else {
284: if (ds->state==DS_STATE_TRUNCATED && t>=DS_MAT_Q) rows = ds->t;
285: else rows = (t==DS_MAT_A && ds->extrarow)? ds->n+1: ds->n;
286: cols = ds->n;
287: }
288: if (m) *m = rows;
289: if (n) *n = cols;
290: return(0);
291: }
293: /*@
294: DSMatIsHermitian - Checks if one of the DS matrices is known to be Hermitian.
296: Not Collective
298: Input Parameters:
299: + ds - the direct solver context
300: - t - the requested matrix
302: Output Parameter:
303: . flg - the Hermitian flag
305: Note:
306: Does not check the matrix values directly. The flag is set according to the
307: problem structure. For instance, in DSHEP matrix A is Hermitian.
309: Level: developer
311: .seealso: DSGetMat()
312: @*/
313: PetscErrorCode DSMatIsHermitian(DS ds,DSMatType t,PetscBool *flg)314: {
320: DSCheckValidMat(ds,t,2);
322: if (ds->ops->hermitian) {
323: (*ds->ops->hermitian)(ds,t,flg);
324: } else *flg = PETSC_FALSE;
325: return(0);
326: }
328: /*@
329: DSGetMat - Returns a sequential dense Mat object containing the requested
330: matrix.
332: Not Collective
334: Input Parameters:
335: + ds - the direct solver context
336: - m - the requested matrix
338: Output Parameter:
339: . A - Mat object
341: Notes:
342: The Mat is created with sizes equal to the current DS dimensions (nxm),
343: then it is filled with the values that would be obtained with DSGetArray()
344: (not DSGetArrayReal()). If the DS was truncated, then the number of rows
345: is equal to the dimension prior to truncation, see DSTruncate().
346: The communicator is always PETSC_COMM_SELF.
348: When no longer needed, the user can either destroy the matrix or call
349: DSRestoreMat(). The latter will copy back the modified values.
351: Level: advanced
353: .seealso: DSRestoreMat(), DSSetDimensions(), DSGetArray(), DSGetArrayReal(), DSTruncate()
354: @*/
355: PetscErrorCode DSGetMat(DS ds,DSMatType m,Mat *A)356: {
358: PetscInt j,rows,cols,arows,acols;
359: PetscBool create=PETSC_FALSE,flg;
360: PetscScalar *pA,*M;
364: DSCheckAlloc(ds,1);
365: DSCheckValidMat(ds,m,2);
367: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
369: DSMatGetSize(ds,m,&rows,&cols);
370: if (!ds->omat[m]) create=PETSC_TRUE;
371: else {
372: MatGetSize(ds->omat[m],&arows,&acols);
373: if (arows!=rows || acols!=cols) {
374: MatDestroy(&ds->omat[m]);
375: create=PETSC_TRUE;
376: }
377: }
378: if (create) {
379: MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]);
380: }
382: /* set Hermitian flag */
383: DSMatIsHermitian(ds,m,&flg);
384: MatSetOption(ds->omat[m],MAT_HERMITIAN,flg);
386: /* copy entries */
387: PetscObjectReference((PetscObject)ds->omat[m]);
388: *A = ds->omat[m];
389: M = ds->mat[m];
390: MatDenseGetArray(*A,&pA);
391: for (j=0;j<cols;j++) {
392: PetscArraycpy(pA+j*rows,M+j*ds->ld,rows);
393: }
394: MatDenseRestoreArray(*A,&pA);
395: return(0);
396: }
398: /*@
399: DSRestoreMat - Restores the matrix after DSGetMat() was called.
401: Not Collective
403: Input Parameters:
404: + ds - the direct solver context
405: . m - the requested matrix
406: - A - the fetched Mat object
408: Notes:
409: A call to this function must match a previous call of DSGetMat().
410: The effect is that the contents of the Mat are copied back to the
411: DS internal array, and the matrix is destroyed.
413: It is not compulsory to call this function, the matrix obtained with
414: DSGetMat() can simply be destroyed if entries need not be copied back.
416: Level: advanced
418: .seealso: DSGetMat(), DSRestoreArray(), DSRestoreArrayReal()
419: @*/
420: PetscErrorCode DSRestoreMat(DS ds,DSMatType m,Mat *A)421: {
423: PetscInt j,rows,cols;
424: PetscScalar *pA,*M;
428: DSCheckAlloc(ds,1);
429: DSCheckValidMat(ds,m,2);
431: if (!ds->omat[m]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSRestoreMat must match a previous call to DSGetMat");
432: 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");
434: MatGetSize(*A,&rows,&cols);
435: M = ds->mat[m];
436: MatDenseGetArray(*A,&pA);
437: for (j=0;j<cols;j++) {
438: PetscArraycpy(M+j*ds->ld,pA+j*rows,rows);
439: }
440: MatDenseRestoreArray(*A,&pA);
441: MatDestroy(A);
442: return(0);
443: }
445: /*@C
446: DSGetArray - Returns a pointer to one of the internal arrays used to
447: represent matrices. You MUST call DSRestoreArray() when you no longer
448: need to access the array.
450: Not Collective
452: Input Parameters:
453: + ds - the direct solver context
454: - m - the requested matrix
456: Output Parameter:
457: . a - pointer to the values
459: Level: advanced
461: .seealso: DSRestoreArray(), DSGetArrayReal()
462: @*/
463: PetscErrorCode DSGetArray(DS ds,DSMatType m,PetscScalar *a[])464: {
467: DSCheckAlloc(ds,1);
468: DSCheckValidMat(ds,m,2);
470: *a = ds->mat[m];
471: CHKMEMQ;
472: return(0);
473: }
475: /*@C
476: DSRestoreArray - Restores the matrix after DSGetArray() was called.
478: Not Collective
480: Input Parameters:
481: + ds - the direct solver context
482: . m - the requested matrix
483: - a - pointer to the values
485: Level: advanced
487: .seealso: DSGetArray(), DSGetArrayReal()
488: @*/
489: PetscErrorCode DSRestoreArray(DS ds,DSMatType m,PetscScalar *a[])490: {
495: DSCheckAlloc(ds,1);
496: DSCheckValidMat(ds,m,2);
498: CHKMEMQ;
499: *a = 0;
500: PetscObjectStateIncrease((PetscObject)ds);
501: return(0);
502: }
504: /*@C
505: DSGetArrayReal - Returns a pointer to one of the internal arrays used to
506: represent real matrices. You MUST call DSRestoreArrayReal() when you no longer
507: need to access the array.
509: Not Collective
511: Input Parameters:
512: + ds - the direct solver context
513: - m - the requested matrix
515: Output Parameter:
516: . a - pointer to the values
518: Level: advanced
520: .seealso: DSRestoreArrayReal(), DSGetArray()
521: @*/
522: PetscErrorCode DSGetArrayReal(DS ds,DSMatType m,PetscReal *a[])523: {
526: DSCheckAlloc(ds,1);
527: DSCheckValidMatReal(ds,m,2);
529: *a = ds->rmat[m];
530: CHKMEMQ;
531: return(0);
532: }
534: /*@C
535: DSRestoreArrayReal - Restores the matrix after DSGetArrayReal() was called.
537: Not Collective
539: Input Parameters:
540: + ds - the direct solver context
541: . m - the requested matrix
542: - a - pointer to the values
544: Level: advanced
546: .seealso: DSGetArrayReal(), DSGetArray()
547: @*/
548: PetscErrorCode DSRestoreArrayReal(DS ds,DSMatType m,PetscReal *a[])549: {
554: DSCheckAlloc(ds,1);
555: DSCheckValidMatReal(ds,m,2);
557: CHKMEMQ;
558: *a = 0;
559: PetscObjectStateIncrease((PetscObject)ds);
560: return(0);
561: }
563: /*@
564: DSSolve - Solves the problem.
566: Logically Collective on ds
568: Input Parameters:
569: + ds - the direct solver context
570: . eigr - array to store the computed eigenvalues (real part)
571: - eigi - array to store the computed eigenvalues (imaginary part)
573: Note:
574: This call brings the dense system to condensed form. No ordering
575: of the eigenvalues is enforced (for this, call DSSort() afterwards).
577: Level: intermediate
579: .seealso: DSSort(), DSStateType580: @*/
581: PetscErrorCode DSSolve(DS ds,PetscScalar eigr[],PetscScalar eigi[])582: {
588: DSCheckAlloc(ds,1);
590: if (ds->state>=DS_STATE_CONDENSED) return(0);
591: if (!ds->ops->solve[ds->method]) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The specified method number does not exist for this DS");
592: PetscInfo3(ds,"Starting solve with problem sizes: n=%D, l=%D, k=%D\n",ds->n,ds->l,ds->k);
593: PetscLogEventBegin(DS_Solve,ds,0,0,0);
594: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
595: (*ds->ops->solve[ds->method])(ds,eigr,eigi);
596: PetscFPTrapPop();
597: PetscLogEventEnd(DS_Solve,ds,0,0,0);
598: PetscInfo1(ds,"State has changed from %s to CONDENSED\n",DSStateTypes[ds->state]);
599: ds->state = DS_STATE_CONDENSED;
600: PetscObjectStateIncrease((PetscObject)ds);
601: return(0);
602: }
604: /*@C
605: DSSort - Sorts the result of DSSolve() according to a given sorting
606: criterion.
608: Logically Collective on ds
610: Input Parameters:
611: + ds - the direct solver context
612: . rr - (optional) array containing auxiliary values (real part)
613: - ri - (optional) array containing auxiliary values (imaginary part)
615: Input/Output Parameters:
616: + eigr - array containing the computed eigenvalues (real part)
617: . eigi - array containing the computed eigenvalues (imaginary part)
618: - k - (optional) number of elements in the leading group
620: Notes:
621: This routine sorts the arrays provided in eigr and eigi, and also
622: sorts the dense system stored inside ds (assumed to be in condensed form).
623: The sorting criterion is specified with DSSetSlepcSC().
625: If arrays rr and ri are provided, then a (partial) reordering based on these
626: values rather than on the eigenvalues is performed. In symmetric problems
627: a total order is obtained (parameter k is ignored), but otherwise the result
628: is sorted only partially. In this latter case, it is only guaranteed that
629: all the first k elements satisfy the comparison with any of the last n-k
630: elements. The output value of parameter k is the final number of elements in
631: the first set.
633: Level: intermediate
635: .seealso: DSSolve(), DSSetSlepcSC(), DSSortWithPermutation()
636: @*/
637: PetscErrorCode DSSort(DS ds,PetscScalar *eigr,PetscScalar *eigi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)638: {
640: PetscInt i;
645: DSCheckSolved(ds,1);
648: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
649: if (!ds->ops->sort) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
650: if (!ds->sc) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must provide a sorting criterion first");
651: if (k && !rr) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Argument k can only be used together with rr");
653: for (i=0;i<ds->n;i++) ds->perm[i] = i; /* initialize to trivial permutation */
654: PetscLogEventBegin(DS_Other,ds,0,0,0);
655: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
656: (*ds->ops->sort)(ds,eigr,eigi,rr,ri,k);
657: PetscFPTrapPop();
658: PetscLogEventEnd(DS_Other,ds,0,0,0);
659: PetscObjectStateIncrease((PetscObject)ds);
660: PetscInfo(ds,"Finished sorting\n");
661: return(0);
662: }
664: /*@C
665: DSSortWithPermutation - Reorders the result of DSSolve() according to a given
666: permutation.
668: Logically Collective on ds
670: Input Parameters:
671: + ds - the direct solver context
672: - perm - permutation that indicates the new ordering
674: Input/Output Parameters:
675: + eigr - array with the reordered eigenvalues (real part)
676: - eigi - array with the reordered eigenvalues (imaginary part)
678: Notes:
679: This routine reorders the arrays provided in eigr and eigi, and also the dense
680: system stored inside ds (assumed to be in condensed form). There is no sorting
681: criterion, as opposed to DSSort(). Instead, the new ordering is given in argument perm.
683: Level: advanced
685: .seealso: DSSolve(), DSSort()
686: @*/
687: PetscErrorCode DSSortWithPermutation(DS ds,PetscInt *perm,PetscScalar *eigr,PetscScalar *eigi)688: {
694: DSCheckSolved(ds,1);
697: if (ds->state==DS_STATE_TRUNCATED) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Cannot sort a truncated DS");
698: if (!ds->ops->sortperm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
700: PetscLogEventBegin(DS_Other,ds,0,0,0);
701: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
702: (*ds->ops->sortperm)(ds,perm,eigr,eigi);
703: PetscFPTrapPop();
704: PetscLogEventEnd(DS_Other,ds,0,0,0);
705: PetscObjectStateIncrease((PetscObject)ds);
706: PetscInfo(ds,"Finished sorting\n");
707: return(0);
708: }
710: /*@
711: DSSynchronize - Make sure that all processes have the same data, performing
712: communication if necessary.
714: Collective on ds
716: Input Parameter:
717: . ds - the direct solver context
719: Input/Output Parameters:
720: + eigr - (optional) array with the computed eigenvalues (real part)
721: - eigi - (optional) array with the computed eigenvalues (imaginary part)
723: Notes:
724: When the DS has been created with a communicator with more than one process,
725: the internal data, especially the computed matrices, may diverge in the
726: different processes. This happens when using multithreaded BLAS and may
727: cause numerical issues in some ill-conditioned problems. This function
728: performs the necessary communication among the processes so that the
729: internal data is exactly equal in all of them.
731: Depending on the parallel mode as set with DSSetParallel(), this function
732: will either do nothing or synchronize the matrices computed by DSSolve()
733: and DSSort(). The arguments eigr and eigi are typically those used in the
734: calls to DSSolve() and DSSort().
736: Level: developer
738: .seealso: DSSetParallel(), DSSolve(), DSSort()
739: @*/
740: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])741: {
743: PetscMPIInt size;
748: DSCheckAlloc(ds,1);
749: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
750: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
751: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
752: if (ds->ops->synchronize) {
753: (*ds->ops->synchronize)(ds,eigr,eigi);
754: }
755: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
756: PetscObjectStateIncrease((PetscObject)ds);
757: PetscInfo1(ds,"Synchronization completed (%s)\n",DSParallelTypes[ds->pmode]);
758: }
759: return(0);
760: }
762: /*@C
763: DSVectors - Compute vectors associated to the dense system such
764: as eigenvectors.
766: Logically Collective on ds
768: Input Parameters:
769: + ds - the direct solver context
770: - mat - the matrix, used to indicate which vectors are required
772: Input/Output Parameter:
773: - j - (optional) index of vector to be computed
775: Output Parameter:
776: . rnorm - (optional) computed residual norm
778: Notes:
779: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_VT, to
780: compute right or left eigenvectors, or left or right singular vectors,
781: respectively.
783: If NULL is passed in argument j then all vectors are computed,
784: otherwise j indicates which vector must be computed. In real non-symmetric
785: problems, on exit the index j will be incremented when a complex conjugate
786: pair is found.
788: This function can be invoked after the dense problem has been solved,
789: to get the residual norm estimate of the associated Ritz pair. In that
790: case, the relevant information is returned in rnorm.
792: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
793: be in (quasi-)triangular form, or call DSSolve() first.
795: Level: intermediate
797: .seealso: DSSolve()
798: @*/
799: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)800: {
806: DSCheckAlloc(ds,1);
808: if (mat>=DS_NUM_MAT) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
809: if (!ds->ops->vectors) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
810: if (rnorm && !j) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
811: if (!ds->mat[mat]) { DSAllocateMat_Private(ds,mat); }
812: if (!j) { PetscInfo1(ds,"Computing all vectors on %s\n",DSMatName[mat]); }
813: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
814: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
815: (*ds->ops->vectors)(ds,mat,j,rnorm);
816: PetscFPTrapPop();
817: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
818: PetscObjectStateIncrease((PetscObject)ds);
819: return(0);
820: }
822: /*@
823: DSUpdateExtraRow - Performs all necessary operations so that the extra
824: row gets up-to-date after a call to DSSolve().
826: Not Collective
828: Input Parameters:
829: . ds - the direct solver context
831: Level: advanced
833: .seealso: DSSolve(), DSSetExtraRow()
834: @*/
835: PetscErrorCode DSUpdateExtraRow(DS ds)836: {
842: DSCheckAlloc(ds,1);
843: if (!ds->ops->update) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
844: if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
845: PetscInfo(ds,"Updating extra row\n");
846: PetscLogEventBegin(DS_Other,ds,0,0,0);
847: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
848: (*ds->ops->update)(ds);
849: PetscFPTrapPop();
850: PetscLogEventEnd(DS_Other,ds,0,0,0);
851: return(0);
852: }
854: /*@
855: DSCond - Compute the inf-norm condition number of the first matrix
856: as cond(A) = norm(A)*norm(inv(A)).
858: Not Collective
860: Input Parameters:
861: + ds - the direct solver context
862: - cond - the computed condition number
864: Level: advanced
866: .seealso: DSSolve()
867: @*/
868: PetscErrorCode DSCond(DS ds,PetscReal *cond)869: {
875: DSCheckAlloc(ds,1);
877: if (!ds->ops->cond) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
878: PetscLogEventBegin(DS_Other,ds,0,0,0);
879: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
880: (*ds->ops->cond)(ds,cond);
881: PetscFPTrapPop();
882: PetscLogEventEnd(DS_Other,ds,0,0,0);
883: PetscInfo1(ds,"Computed condition number = %g\n",(double)*cond);
884: return(0);
885: }
887: /*@C
888: DSTranslateHarmonic - Computes a translation of the dense system.
890: Logically Collective on ds
892: Input Parameters:
893: + ds - the direct solver context
894: . tau - the translation amount
895: . beta - last component of vector b
896: - recover - boolean flag to indicate whether to recover or not
898: Output Parameters:
899: + g - the computed vector (optional)
900: - gamma - scale factor (optional)
902: Notes:
903: This function is intended for use in the context of Krylov methods only.
904: It computes a translation of a Krylov decomposition in order to extract
905: eigenpair approximations by harmonic Rayleigh-Ritz.
906: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
907: vector b is assumed to be beta*e_n^T.
909: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
910: the factor by which the residual of the Krylov decomposition is scaled.
912: If the recover flag is activated, the computed translation undoes the
913: translation done previously. In that case, parameter tau is ignored.
915: Level: developer
916: @*/
917: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)918: {
924: DSCheckAlloc(ds,1);
925: if (!ds->ops->transharm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
926: if (recover) { PetscInfo(ds,"Undoing the translation\n"); }
927: else { PetscInfo(ds,"Computing the translation\n"); }
928: PetscLogEventBegin(DS_Other,ds,0,0,0);
929: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
930: (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
931: PetscFPTrapPop();
932: PetscLogEventEnd(DS_Other,ds,0,0,0);
933: ds->state = DS_STATE_RAW;
934: PetscObjectStateIncrease((PetscObject)ds);
935: return(0);
936: }
938: /*@
939: DSTranslateRKS - Computes a modification of the dense system corresponding
940: to an update of the shift in a rational Krylov method.
942: Logically Collective on ds
944: Input Parameters:
945: + ds - the direct solver context
946: - alpha - the translation amount
948: Notes:
949: This function is intended for use in the context of Krylov methods only.
950: It takes the leading (k+1,k) submatrix of A, containing the truncated
951: Rayleigh quotient of a Krylov-Schur relation computed from a shift
952: sigma1 and transforms it to obtain a Krylov relation as if computed
953: from a different shift sigma2. The new matrix is computed as
954: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
955: alpha = sigma1-sigma2.
957: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
958: Krylov basis.
960: Level: developer
961: @*/
962: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)963: {
969: DSCheckAlloc(ds,1);
970: if (!ds->ops->transrks) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
971: PetscInfo1(ds,"Translating with alpha=%g\n",PetscRealPart(alpha));
972: PetscLogEventBegin(DS_Other,ds,0,0,0);
973: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
974: (*ds->ops->transrks)(ds,alpha);
975: PetscFPTrapPop();
976: PetscLogEventEnd(DS_Other,ds,0,0,0);
977: ds->state = DS_STATE_RAW;
978: ds->compact = PETSC_FALSE;
979: PetscObjectStateIncrease((PetscObject)ds);
980: return(0);
981: }
983: /*@
984: DSCopyMat - Copies the contents of a sequential dense Mat object to
985: the indicated DS matrix, or vice versa.
987: Not Collective
989: Input Parameters:
990: + ds - the direct solver context
991: . m - the requested matrix
992: . mr - first row of m to be considered
993: . mc - first column of m to be considered
994: . A - Mat object
995: . Ar - first row of A to be considered
996: . Ac - first column of A to be considered
997: . rows - number of rows to copy
998: . cols - number of columns to copy
999: - out - whether the data is copied out of the DS1001: Note:
1002: If out=true, the values of the DS matrix m are copied to A, otherwise
1003: the entries of A are copied to the DS.
1005: Level: developer
1007: .seealso: DSGetMat()
1008: @*/
1009: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)1010: {
1012: PetscInt j,mrows,mcols,arows,acols;
1013: PetscScalar *pA,*M;
1017: DSCheckAlloc(ds,1);
1019: DSCheckValidMat(ds,m,2);
1028: if (!rows || !cols) return(0);
1030: DSMatGetSize(ds,m,&mrows,&mcols);
1031: MatGetSize(A,&arows,&acols);
1032: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
1033: if (mr<0 || mr>=mrows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in m");
1034: if (mc<0 || mc>=mcols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in m");
1035: if (Ar<0 || Ar>=arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in A");
1036: if (Ac<0 || Ac>=acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in A");
1037: if (mr+rows>mrows || Ar+rows>arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of rows");
1038: if (mc+cols>mcols || Ac+cols>acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of columns");
1040: M = ds->mat[m];
1041: MatDenseGetArray(A,&pA);
1042: for (j=0;j<cols;j++) {
1043: if (out) {
1044: PetscArraycpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows);
1045: } else {
1046: PetscArraycpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows);
1047: }
1048: }
1049: MatDenseRestoreArray(A,&pA);
1050: return(0);
1051: }