Actual source code: lapack.c
slepc-3.13.4 2020-09-02
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: This file implements a wrapper to the LAPACK eigenvalue subroutines.
12: Generalized problems are transformed to standard ones only if necessary.
13: */
15: #include <slepc/private/epsimpl.h>
17: PetscErrorCode EPSSetUp_LAPACK(EPS eps)
18: {
19: PetscErrorCode ierr,ierra,ierrb;
20: PetscBool isshift,flg,denseok=PETSC_FALSE;
21: Mat A,B,OP,shell,Ar,Br,Adense=NULL,Bdense=NULL;
22: PetscScalar shift,*Ap,*Bp;
23: PetscInt i,ld,nmat;
24: KSP ksp;
25: PC pc;
26: Vec v;
29: eps->ncv = eps->n;
30: if (eps->mpd!=PETSC_DEFAULT) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
31: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
32: if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
33: if (eps->which==EPS_ALL && eps->inta!=eps->intb) SETERRQ(PetscObjectComm((PetscObject)eps),1,"This solver does not support interval computation");
34: if (eps->balance!=EPS_BALANCE_NONE) { PetscInfo(eps,"Warning: balancing ignored\n"); }
35: if (eps->stopping!=EPSStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"User-defined stopping test not supported");
36: if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
37: EPSAllocateSolution(eps,0);
39: /* attempt to get dense representations of A and B separately */
40: PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift);
41: if (isshift) {
42: STGetNumMatrices(eps->st,&nmat);
43: STGetMatrix(eps->st,0,&A);
44: MatHasOperation(A,MATOP_CREATE_SUBMATRICES,&flg);
45: if (flg) {
46: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
47: ierra = MatCreateRedundantMatrix(A,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar);
48: if (!ierra) { ierra |= MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&Adense); }
49: ierra |= MatDestroy(&Ar);
50: PetscPopErrorHandler();
51: } else ierra = 1;
52: if (nmat>1) {
53: STGetMatrix(eps->st,1,&B);
54: MatHasOperation(B,MATOP_CREATE_SUBMATRICES,&flg);
55: if (flg) {
56: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
57: ierrb = MatCreateRedundantMatrix(B,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Br);
58: if (!ierrb) { ierrb |= MatConvert(Br,MATSEQDENSE,MAT_INITIAL_MATRIX,&Bdense); }
59: ierrb |= MatDestroy(&Br);
60: PetscPopErrorHandler();
61: } else ierrb = 1;
62: } else ierrb = 0;
63: denseok = PetscNot(ierra || ierrb);
64: }
66: /* setup DS */
67: if (denseok) {
68: if (eps->isgeneralized) {
69: if (eps->ishermitian) {
70: if (eps->ispositive) {
71: DSSetType(eps->ds,DSGHEP);
72: } else {
73: DSSetType(eps->ds,DSGNHEP); /* TODO: should be DSGHIEP */
74: }
75: } else {
76: DSSetType(eps->ds,DSGNHEP);
77: }
78: } else {
79: if (eps->ishermitian) {
80: DSSetType(eps->ds,DSHEP);
81: } else {
82: DSSetType(eps->ds,DSNHEP);
83: }
84: }
85: } else {
86: DSSetType(eps->ds,DSNHEP);
87: }
88: DSAllocate(eps->ds,eps->ncv);
89: DSGetLeadingDimension(eps->ds,&ld);
90: DSSetDimensions(eps->ds,eps->ncv,0,0,0);
92: if (denseok) {
93: STGetShift(eps->st,&shift);
94: if (shift != 0.0) {
95: MatShift(Adense,shift);
96: }
97: /* use dummy pc and ksp to avoid problems when B is not positive definite */
98: STGetKSP(eps->st,&ksp);
99: KSPSetType(ksp,KSPPREONLY);
100: KSPGetPC(ksp,&pc);
101: PCSetType(pc,PCNONE);
102: } else {
103: PetscInfo(eps,"Using slow explicit operator\n");
104: STGetOperator(eps->st,&shell);
105: MatComputeOperator(shell,MATDENSE,&OP);
106: STRestoreOperator(eps->st,&shell);
107: MatDestroy(&Adense);
108: MatCreateRedundantMatrix(OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Adense);
109: MatDestroy(&OP);
110: }
112: /* fill DS matrices */
113: VecCreateSeqWithArray(PETSC_COMM_SELF,1,ld,NULL,&v);
114: DSGetArray(eps->ds,DS_MAT_A,&Ap);
115: for (i=0;i<ld;i++) {
116: VecPlaceArray(v,Ap+i*ld);
117: MatGetColumnVector(Adense,v,i);
118: VecResetArray(v);
119: }
120: DSRestoreArray(eps->ds,DS_MAT_A,&Ap);
121: if (denseok && eps->isgeneralized) {
122: DSGetArray(eps->ds,DS_MAT_B,&Bp);
123: for (i=0;i<ld;i++) {
124: VecPlaceArray(v,Bp+i*ld);
125: MatGetColumnVector(Bdense,v,i);
126: VecResetArray(v);
127: }
128: DSRestoreArray(eps->ds,DS_MAT_B,&Bp);
129: }
130: VecDestroy(&v);
131: DSSetState(eps->ds,DS_STATE_RAW);
132: MatDestroy(&Adense);
133: MatDestroy(&Bdense);
134: return(0);
135: }
137: PetscErrorCode EPSSolve_LAPACK(EPS eps)
138: {
140: PetscInt n=eps->n,i,low,high;
141: PetscScalar *array,*pX,*pY;
142: Vec v,w;
145: DSSolve(eps->ds,eps->eigr,eps->eigi);
146: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
147: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
149: /* right eigenvectors */
150: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
151: DSGetArray(eps->ds,DS_MAT_X,&pX);
152: for (i=0;i<eps->ncv;i++) {
153: BVGetColumn(eps->V,i,&v);
154: VecGetOwnershipRange(v,&low,&high);
155: VecGetArray(v,&array);
156: PetscArraycpy(array,pX+i*n+low,high-low);
157: VecRestoreArray(v,&array);
158: BVRestoreColumn(eps->V,i,&v);
159: }
160: DSRestoreArray(eps->ds,DS_MAT_X,&pX);
162: /* left eigenvectors */
163: if (eps->twosided) {
164: DSVectors(eps->ds,DS_MAT_Y,NULL,NULL);
165: DSGetArray(eps->ds,DS_MAT_Y,&pY);
166: for (i=0;i<eps->ncv;i++) {
167: BVGetColumn(eps->W,i,&w);
168: VecGetOwnershipRange(w,&low,&high);
169: VecGetArray(w,&array);
170: PetscArraycpy(array,pY+i*n+low,high-low);
171: VecRestoreArray(w,&array);
172: BVRestoreColumn(eps->W,i,&w);
173: }
174: DSRestoreArray(eps->ds,DS_MAT_Y,&pY);
175: }
177: eps->nconv = eps->ncv;
178: eps->its = 1;
179: eps->reason = EPS_CONVERGED_TOL;
180: return(0);
181: }
183: SLEPC_EXTERN PetscErrorCode EPSCreate_LAPACK(EPS eps)
184: {
186: eps->useds = PETSC_TRUE;
187: eps->hasts = PETSC_TRUE;
188: eps->categ = EPS_CATEGORY_OTHER;
190: eps->ops->solve = EPSSolve_LAPACK;
191: eps->ops->setup = EPSSetUp_LAPACK;
192: eps->ops->backtransform = EPSBackTransform_Default;
193: return(0);
194: }