1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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: . eigr - array containing the computed eigenvalues (real part)
603: . eigi - array containing the computed eigenvalues (imaginary part)
604: . rr - (optional) array containing auxiliary values (real part)
605: - ri - (optional) array containing auxiliary values (imaginary part)
607: Input/Output Parameter:
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()
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: /*@
654: DSSynchronize - Make sure that all processes have the same data, performing
655: communication if necessary.
657: Collective on DS659: Input Parameter:
660: + ds - the direct solver context
662: Input/Output Parameters:
663: + eigr - (optional) array with the computed eigenvalues (real part)
664: - eigi - (optional) array with the computed eigenvalues (imaginary part)
666: Notes:
667: When the DS has been created with a communicator with more than one process,
668: the internal data, especially the computed matrices, may diverge in the
669: different processes. This happens when using multithreaded BLAS and may
670: cause numerical issues in some ill-conditioned problems. This function
671: performs the necessary communication among the processes so that the
672: internal data is exactly equal in all of them.
674: Depending on the parallel mode as set with DSSetParallel(), this function
675: will either do nothing or synchronize the matrices computed by DSSolve()
676: and DSSort(). The arguments eigr and eigi are typically those used in the
677: calls to DSSolve() and DSSort().
679: Level: developer
681: .seealso: DSSetParallel(), DSSolve(), DSSort()
682: @*/
683: PetscErrorCode DSSynchronize(DS ds,PetscScalar eigr[],PetscScalar eigi[])684: {
686: PetscMPIInt size;
691: DSCheckAlloc(ds,1);
692: MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
693: if (size>1 && ds->pmode==DS_PARALLEL_SYNCHRONIZED) {
694: PetscLogEventBegin(DS_Synchronize,ds,0,0,0);
695: if (ds->ops->synchronize) {
696: (*ds->ops->synchronize)(ds,eigr,eigi);
697: }
698: PetscLogEventEnd(DS_Synchronize,ds,0,0,0);
699: PetscObjectStateIncrease((PetscObject)ds);
700: }
701: return(0);
702: }
704: /*@C
705: DSVectors - Compute vectors associated to the dense system such
706: as eigenvectors.
708: Logically Collective on DS710: Input Parameters:
711: + ds - the direct solver context
712: - mat - the matrix, used to indicate which vectors are required
714: Input/Output Parameter:
715: - j - (optional) index of vector to be computed
717: Output Parameter:
718: . rnorm - (optional) computed residual norm
720: Notes:
721: Allowed values for mat are DS_MAT_X, DS_MAT_Y, DS_MAT_U and DS_MAT_VT, to
722: compute right or left eigenvectors, or left or right singular vectors,
723: respectively.
725: If NULL is passed in argument j then all vectors are computed,
726: otherwise j indicates which vector must be computed. In real non-symmetric
727: problems, on exit the index j will be incremented when a complex conjugate
728: pair is found.
730: This function can be invoked after the dense problem has been solved,
731: to get the residual norm estimate of the associated Ritz pair. In that
732: case, the relevant information is returned in rnorm.
734: For computing eigenvectors, LAPACK's _trevc is used so the matrix must
735: be in (quasi-)triangular form, or call DSSolve() first.
737: Level: intermediate
739: .seealso: DSSolve()
740: @*/
741: PetscErrorCode DSVectors(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)742: {
748: DSCheckAlloc(ds,1);
750: if (mat>=DS_NUM_MAT) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid matrix");
751: if (!ds->ops->vectors) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
752: if (rnorm && !j) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"Must give a value of j");
753: if (!ds->mat[mat]) { DSAllocateMat_Private(ds,mat); }
754: PetscLogEventBegin(DS_Vectors,ds,0,0,0);
755: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
756: (*ds->ops->vectors)(ds,mat,j,rnorm);
757: PetscFPTrapPop();
758: PetscLogEventEnd(DS_Vectors,ds,0,0,0);
759: PetscObjectStateIncrease((PetscObject)ds);
760: return(0);
761: }
763: /*@
764: DSUpdateExtraRow - Performs all necessary operations so that the extra
765: row gets up-to-date after a call to DSSolve().
767: Not Collective
769: Input Parameters:
770: . ds - the direct solver context
772: Level: advanced
774: .seealso: DSSolve(), DSSetExtraRow()
775: @*/
776: PetscErrorCode DSUpdateExtraRow(DS ds)777: {
783: DSCheckAlloc(ds,1);
784: if (!ds->ops->update) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
785: if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Should have called DSSetExtraRow");
786: PetscLogEventBegin(DS_Other,ds,0,0,0);
787: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
788: (*ds->ops->update)(ds);
789: PetscFPTrapPop();
790: PetscLogEventEnd(DS_Other,ds,0,0,0);
791: return(0);
792: }
794: /*@
795: DSCond - Compute the inf-norm condition number of the first matrix
796: as cond(A) = norm(A)*norm(inv(A)).
798: Not Collective
800: Input Parameters:
801: + ds - the direct solver context
802: - cond - the computed condition number
804: Level: advanced
806: .seealso: DSSolve()
807: @*/
808: PetscErrorCode DSCond(DS ds,PetscReal *cond)809: {
815: DSCheckAlloc(ds,1);
817: if (!ds->ops->cond) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
818: PetscLogEventBegin(DS_Other,ds,0,0,0);
819: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
820: (*ds->ops->cond)(ds,cond);
821: PetscFPTrapPop();
822: PetscLogEventEnd(DS_Other,ds,0,0,0);
823: return(0);
824: }
826: /*@C
827: DSTranslateHarmonic - Computes a translation of the dense system.
829: Logically Collective on DS831: Input Parameters:
832: + ds - the direct solver context
833: . tau - the translation amount
834: . beta - last component of vector b
835: - recover - boolean flag to indicate whether to recover or not
837: Output Parameters:
838: + g - the computed vector (optional)
839: - gamma - scale factor (optional)
841: Notes:
842: This function is intended for use in the context of Krylov methods only.
843: It computes a translation of a Krylov decomposition in order to extract
844: eigenpair approximations by harmonic Rayleigh-Ritz.
845: The matrix is updated as A + g*b' where g = (A-tau*eye(n))'\b and
846: vector b is assumed to be beta*e_n^T.
848: The gamma factor is defined as sqrt(1+g'*g) and can be interpreted as
849: the factor by which the residual of the Krylov decomposition is scaled.
851: If the recover flag is activated, the computed translation undoes the
852: translation done previously. In that case, parameter tau is ignored.
854: Level: developer
855: @*/
856: PetscErrorCode DSTranslateHarmonic(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *g,PetscReal *gamma)857: {
863: DSCheckAlloc(ds,1);
864: if (!ds->ops->transharm) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
865: PetscLogEventBegin(DS_Other,ds,0,0,0);
866: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
867: (*ds->ops->transharm)(ds,tau,beta,recover,g,gamma);
868: PetscFPTrapPop();
869: PetscLogEventEnd(DS_Other,ds,0,0,0);
870: ds->state = DS_STATE_RAW;
871: PetscObjectStateIncrease((PetscObject)ds);
872: return(0);
873: }
875: /*@
876: DSTranslateRKS - Computes a modification of the dense system corresponding
877: to an update of the shift in a rational Krylov method.
879: Logically Collective on DS881: Input Parameters:
882: + ds - the direct solver context
883: - alpha - the translation amount
885: Notes:
886: This function is intended for use in the context of Krylov methods only.
887: It takes the leading (k+1,k) submatrix of A, containing the truncated
888: Rayleigh quotient of a Krylov-Schur relation computed from a shift
889: sigma1 and transforms it to obtain a Krylov relation as if computed
890: from a different shift sigma2. The new matrix is computed as
891: 1.0/alpha*(eye(k)-Q*inv(R)), where [Q,R]=qr(eye(k)-alpha*A) and
892: alpha = sigma1-sigma2.
894: Matrix Q is placed in DS_MAT_Q so that it can be used to update the
895: Krylov basis.
897: Level: developer
898: @*/
899: PetscErrorCode DSTranslateRKS(DS ds,PetscScalar alpha)900: {
906: DSCheckAlloc(ds,1);
907: if (!ds->ops->transrks) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DS type %s",((PetscObject)ds)->type_name);
908: PetscLogEventBegin(DS_Other,ds,0,0,0);
909: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
910: (*ds->ops->transrks)(ds,alpha);
911: PetscFPTrapPop();
912: PetscLogEventEnd(DS_Other,ds,0,0,0);
913: ds->state = DS_STATE_RAW;
914: ds->compact = PETSC_FALSE;
915: PetscObjectStateIncrease((PetscObject)ds);
916: return(0);
917: }
919: /*@
920: DSCopyMat - Copies the contents of a sequential dense Mat object to
921: the indicated DS matrix, or vice versa.
923: Not Collective
925: Input Parameters:
926: + ds - the direct solver context
927: . m - the requested matrix
928: . mr - first row of m to be considered
929: . mc - first column of m to be considered
930: . A - Mat object
931: . Ar - first row of A to be considered
932: . Ac - first column of A to be considered
933: . rows - number of rows to copy
934: . cols - number of columns to copy
935: - out - whether the data is copied out of the DS937: Note:
938: If out=true, the values of the DS matrix m are copied to A, otherwise
939: the entries of A are copied to the DS.
941: Level: developer
943: .seealso: DSGetMat()
944: @*/
945: PetscErrorCode DSCopyMat(DS ds,DSMatType m,PetscInt mr,PetscInt mc,Mat A,PetscInt Ar,PetscInt Ac,PetscInt rows,PetscInt cols,PetscBool out)946: {
948: PetscInt j,mrows,mcols,arows,acols;
949: PetscScalar *pA,*M;
953: DSCheckAlloc(ds,1);
955: DSCheckValidMat(ds,m,2);
964: if (!rows || !cols) return(0);
966: DSMatGetSize(ds,m,&mrows,&mcols);
967: MatGetSize(A,&arows,&acols);
968: if (m==DS_MAT_T || m==DS_MAT_D) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Not implemented for T or D matrices");
969: if (mr<0 || mr>=mrows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in m");
970: if (mc<0 || mc>=mcols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in m");
971: if (Ar<0 || Ar>=arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial row in A");
972: if (Ac<0 || Ac>=acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial column in A");
973: if (mr+rows>mrows || Ar+rows>arows) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of rows");
974: if (mc+cols>mcols || Ac+cols>acols) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of columns");
976: M = ds->mat[m];
977: MatDenseGetArray(A,&pA);
978: for (j=0;j<cols;j++) {
979: if (out) {
980: PetscMemcpy(pA+(Ac+j)*arows+Ar,M+(mc+j)*ds->ld+mr,rows*sizeof(PetscScalar));
981: } else {
982: PetscMemcpy(M+(mc+j)*ds->ld+mr,pA+(Ac+j)*arows+Ar,rows*sizeof(PetscScalar));
983: }
984: }
985: MatDenseRestoreArray(A,&pA);
986: return(0);
987: }