Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 ARPACK package
12 : */
13 :
14 : #include <slepc/private/epsimpl.h>
15 : #include "arpack.h"
16 :
17 7 : static PetscErrorCode EPSSetUp_ARPACK(EPS eps)
18 : {
19 7 : PetscInt ncv;
20 7 : EPS_ARPACK *ar = (EPS_ARPACK*)eps->data;
21 :
22 7 : PetscFunctionBegin;
23 7 : EPSCheckDefinite(eps);
24 7 : if (eps->ncv!=PETSC_DEFAULT) {
25 3 : PetscCheck(eps->ncv>=eps->nev+2,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The value of ncv must be at least nev+2");
26 4 : } else eps->ncv = PetscMin(PetscMax(20,2*eps->nev+1),eps->n); /* set default value of ncv */
27 7 : if (eps->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
28 7 : if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(300,(PetscInt)(2*eps->n/eps->ncv));
29 7 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
30 7 : PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
31 7 : PetscCheck(eps->which!=EPS_WHICH_USER,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support user-defined ordering of eigenvalues");
32 7 : EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
33 7 : EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
34 :
35 7 : ncv = eps->ncv;
36 : #if defined(PETSC_USE_COMPLEX)
37 7 : PetscCall(PetscFree(ar->rwork));
38 7 : PetscCall(PetscMalloc1(ncv,&ar->rwork));
39 7 : ar->lworkl = 3*ncv*ncv+5*ncv;
40 7 : PetscCall(PetscFree(ar->workev));
41 7 : PetscCall(PetscMalloc1(3*ncv,&ar->workev));
42 : #else
43 : if (eps->ishermitian) {
44 : ar->lworkl = ncv*(ncv+8);
45 : } else {
46 : ar->lworkl = 3*ncv*ncv+6*ncv;
47 : PetscCall(PetscFree(ar->workev));
48 : PetscCall(PetscMalloc1(3*ncv,&ar->workev));
49 : }
50 : #endif
51 7 : PetscCall(PetscFree(ar->workl));
52 7 : PetscCall(PetscMalloc1(ar->lworkl,&ar->workl));
53 7 : PetscCall(PetscFree(ar->select));
54 7 : PetscCall(PetscMalloc1(ncv,&ar->select));
55 7 : PetscCall(PetscFree(ar->workd));
56 7 : PetscCall(PetscMalloc1(3*eps->nloc,&ar->workd));
57 :
58 7 : PetscCall(EPSAllocateSolution(eps,0));
59 7 : PetscCall(EPS_SetInnerProduct(eps));
60 7 : PetscCall(EPSSetWorkVecs(eps,2));
61 7 : PetscFunctionReturn(PETSC_SUCCESS);
62 : }
63 :
64 7 : static PetscErrorCode EPSSolve_ARPACK(EPS eps)
65 : {
66 7 : EPS_ARPACK *ar = (EPS_ARPACK*)eps->data;
67 7 : char bmat[1],howmny[] = "A";
68 7 : const char *which;
69 7 : PetscInt n,ld,iparam[11],ipntr[14],ido,info,nev,ncv,rvec;
70 : #if !defined(PETSC_HAVE_MPIUNI) && !defined(PETSC_HAVE_MSMPI)
71 7 : MPI_Fint fcomm;
72 : #endif
73 7 : PetscScalar sigmar,*pV,*resid;
74 7 : Vec x,y,w = eps->work[0];
75 7 : Mat A;
76 7 : PetscBool isSinv,isShift;
77 : #if !defined(PETSC_USE_COMPLEX)
78 : PetscScalar sigmai = 0.0;
79 : #endif
80 :
81 7 : PetscFunctionBegin;
82 7 : nev = eps->nev;
83 7 : ncv = eps->ncv;
84 : #if !defined(PETSC_HAVE_MPIUNI) && !defined(PETSC_HAVE_MSMPI)
85 7 : fcomm = MPI_Comm_c2f(PetscObjectComm((PetscObject)eps));
86 : #endif
87 7 : n = eps->nloc;
88 7 : PetscCall(EPSGetStartVector(eps,0,NULL));
89 7 : PetscCall(BVSetActiveColumns(eps->V,0,0)); /* just for deflation space */
90 7 : PetscCall(BVCopyVec(eps->V,0,eps->work[1]));
91 7 : PetscCall(BVGetLeadingDimension(eps->V,&ld));
92 7 : PetscCall(BVGetArray(eps->V,&pV));
93 7 : PetscCall(VecGetArray(eps->work[1],&resid));
94 :
95 7 : ido = 0; /* first call to reverse communication interface */
96 7 : info = 1; /* indicates an initial vector is provided */
97 7 : iparam[0] = 1; /* use exact shifts */
98 7 : iparam[2] = eps->max_it; /* max Arnoldi iterations */
99 7 : iparam[3] = 1; /* blocksize */
100 7 : iparam[4] = 0; /* number of converged Ritz values */
101 :
102 : /*
103 : Computational modes ([]=not supported):
104 : symmetric non-symmetric complex
105 : 1 1 'I' 1 'I' 1 'I'
106 : 2 3 'I' 3 'I' 3 'I'
107 : 3 2 'G' 2 'G' 2 'G'
108 : 4 3 'G' 3 'G' 3 'G'
109 : 5 [ 4 'G' ] [ 3 'G' ]
110 : 6 [ 5 'G' ] [ 4 'G' ]
111 : */
112 7 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&isSinv));
113 7 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isShift));
114 7 : PetscCall(STGetShift(eps->st,&sigmar));
115 7 : PetscCall(STGetMatrix(eps->st,0,&A));
116 7 : PetscCall(MatCreateVecsEmpty(A,&x,&y));
117 :
118 7 : if (isSinv) {
119 : /* shift-and-invert mode */
120 1 : iparam[6] = 3;
121 1 : if (eps->ispositive) bmat[0] = 'G';
122 1 : else bmat[0] = 'I';
123 6 : } else if (isShift && eps->ispositive) {
124 : /* generalized shift mode with B positive definite */
125 1 : iparam[6] = 2;
126 1 : bmat[0] = 'G';
127 : } else {
128 : /* regular mode */
129 5 : PetscCheck(!eps->ishermitian || !eps->isgeneralized,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectral transformation not supported by ARPACK hermitian solver");
130 5 : iparam[6] = 1;
131 5 : bmat[0] = 'I';
132 : }
133 :
134 : #if !defined(PETSC_USE_COMPLEX)
135 : if (eps->ishermitian) {
136 : switch (eps->which) {
137 : case EPS_TARGET_MAGNITUDE:
138 : case EPS_LARGEST_MAGNITUDE: which = "LM"; break;
139 : case EPS_SMALLEST_MAGNITUDE: which = "SM"; break;
140 : case EPS_TARGET_REAL:
141 : case EPS_LARGEST_REAL: which = "LA"; break;
142 : case EPS_SMALLEST_REAL: which = "SA"; break;
143 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Wrong value of eps->which");
144 : }
145 : } else {
146 : #endif
147 7 : switch (eps->which) {
148 : case EPS_TARGET_MAGNITUDE:
149 7 : case EPS_LARGEST_MAGNITUDE: which = "LM"; break;
150 0 : case EPS_SMALLEST_MAGNITUDE: which = "SM"; break;
151 0 : case EPS_TARGET_REAL:
152 0 : case EPS_LARGEST_REAL: which = "LR"; break;
153 2 : case EPS_SMALLEST_REAL: which = "SR"; break;
154 0 : case EPS_TARGET_IMAGINARY:
155 0 : case EPS_LARGEST_IMAGINARY: which = "LI"; break;
156 0 : case EPS_SMALLEST_IMAGINARY: which = "SI"; break;
157 0 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Wrong value of eps->which");
158 : }
159 : #if !defined(PETSC_USE_COMPLEX)
160 : }
161 : #endif
162 :
163 1373 : do {
164 :
165 : #if !defined(PETSC_USE_COMPLEX)
166 : if (eps->ishermitian) {
167 : PetscStackCallExternalVoid("ARPACKsaupd",ARPACKsaupd_(&fcomm,&ido,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&ld,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,&info));
168 : } else {
169 : PetscStackCallExternalVoid("ARPACKnaupd",ARPACKnaupd_(&fcomm,&ido,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&ld,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,&info));
170 : }
171 : #else
172 1373 : PetscStackCallExternalVoid("ARPACKnaupd",ARPACKnaupd_(&fcomm,&ido,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&ld,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,ar->rwork,&info));
173 : #endif
174 :
175 1373 : if (ido == -1 || ido == 1 || ido == 2) {
176 1366 : if (ido == 1 && iparam[6] == 3 && bmat[0] == 'G') PetscCall(VecPlaceArray(x,&ar->workd[ipntr[2]-1])); /* special case for shift-and-invert with B semi-positive definite*/
177 1366 : else PetscCall(VecPlaceArray(x,&ar->workd[ipntr[0]-1]));
178 1366 : PetscCall(VecPlaceArray(y,&ar->workd[ipntr[1]-1]));
179 :
180 1366 : if (ido == -1) {
181 : /* Y = OP * X for the initialization phase to
182 : force the starting vector into the range of OP */
183 1 : PetscCall(STApply(eps->st,x,y));
184 1365 : } else if (ido == 2) {
185 : /* Y = B * X */
186 518 : PetscCall(BVApplyMatrix(eps->V,x,y));
187 : } else { /* ido == 1 */
188 847 : if (iparam[6] == 3 && bmat[0] == 'G') {
189 : /* Y = OP * X for shift-and-invert with B semi-positive definite */
190 0 : PetscCall(STMatSolve(eps->st,x,y));
191 847 : } else if (iparam[6] == 2) {
192 : /* X=A*X Y=B^-1*X for shift with B positive definite */
193 179 : PetscCall(MatMult(A,x,y));
194 179 : if (sigmar != 0.0) {
195 0 : PetscCall(BVApplyMatrix(eps->V,x,w));
196 0 : PetscCall(VecAXPY(y,sigmar,w));
197 : }
198 179 : PetscCall(VecCopy(y,x));
199 179 : PetscCall(STMatSolve(eps->st,x,y));
200 : } else {
201 : /* Y = OP * X */
202 668 : PetscCall(STApply(eps->st,x,y));
203 : }
204 847 : PetscCall(BVOrthogonalizeVec(eps->V,y,NULL,NULL,NULL));
205 : }
206 :
207 1366 : PetscCall(VecResetArray(x));
208 1366 : PetscCall(VecResetArray(y));
209 7 : } else PetscCheck(ido==99,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Internal error in ARPACK reverse communication interface (ido=%" PetscInt_FMT ")",ido);
210 :
211 1373 : } while (ido != 99);
212 :
213 7 : eps->nconv = iparam[4];
214 7 : eps->its = iparam[2];
215 :
216 7 : PetscCheck(info!=3,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"No shift could be applied in xxAUPD. Try increasing the size of NCV relative to NEV");
217 7 : PetscCheck(info==0 || info==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by ARPACK subroutine xxAUPD (%" PetscInt_FMT ")",info);
218 :
219 7 : rvec = PETSC_TRUE;
220 :
221 7 : if (eps->nconv > 0) {
222 : #if !defined(PETSC_USE_COMPLEX)
223 : if (eps->ishermitian) {
224 : PetscStackCallExternalVoid("ARPACKseupd",ARPACKseupd_(&fcomm,&rvec,howmny,ar->select,eps->eigr,pV,&ld,&sigmar,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&ld,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,&info));
225 : } else {
226 : PetscStackCallExternalVoid("ARPACKneupd",ARPACKneupd_(&fcomm,&rvec,howmny,ar->select,eps->eigr,eps->eigi,pV,&ld,&sigmar,&sigmai,ar->workev,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&ld,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,&info));
227 : }
228 : #else
229 7 : PetscStackCallExternalVoid("ARPACKneupd",ARPACKneupd_(&fcomm,&rvec,howmny,ar->select,eps->eigr,pV,&ld,&sigmar,ar->workev,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&ld,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,ar->rwork,&info));
230 : #endif
231 7 : PetscCheck(info==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by ARPACK subroutine xxEUPD (%" PetscInt_FMT ")",info);
232 : }
233 :
234 7 : PetscCall(BVRestoreArray(eps->V,&pV));
235 7 : PetscCall(VecRestoreArray(eps->work[1],&resid));
236 7 : if (eps->nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
237 0 : else eps->reason = EPS_DIVERGED_ITS;
238 :
239 7 : PetscCall(VecDestroy(&x));
240 7 : PetscCall(VecDestroy(&y));
241 7 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 7 : static PetscErrorCode EPSBackTransform_ARPACK(EPS eps)
245 : {
246 7 : PetscBool isSinv;
247 :
248 7 : PetscFunctionBegin;
249 7 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&isSinv));
250 7 : if (!isSinv) PetscCall(EPSBackTransform_Default(eps));
251 7 : PetscFunctionReturn(PETSC_SUCCESS);
252 : }
253 :
254 6 : static PetscErrorCode EPSReset_ARPACK(EPS eps)
255 : {
256 6 : EPS_ARPACK *ar = (EPS_ARPACK*)eps->data;
257 :
258 6 : PetscFunctionBegin;
259 6 : PetscCall(PetscFree(ar->workev));
260 6 : PetscCall(PetscFree(ar->workl));
261 6 : PetscCall(PetscFree(ar->select));
262 6 : PetscCall(PetscFree(ar->workd));
263 : #if defined(PETSC_USE_COMPLEX)
264 6 : PetscCall(PetscFree(ar->rwork));
265 : #endif
266 6 : PetscFunctionReturn(PETSC_SUCCESS);
267 : }
268 :
269 6 : static PetscErrorCode EPSDestroy_ARPACK(EPS eps)
270 : {
271 6 : PetscFunctionBegin;
272 6 : PetscCall(PetscFree(eps->data));
273 6 : PetscFunctionReturn(PETSC_SUCCESS);
274 : }
275 :
276 6 : SLEPC_EXTERN PetscErrorCode EPSCreate_ARPACK(EPS eps)
277 : {
278 6 : EPS_ARPACK *ctx;
279 :
280 6 : PetscFunctionBegin;
281 6 : PetscCall(PetscNew(&ctx));
282 6 : eps->data = (void*)ctx;
283 :
284 6 : eps->ops->solve = EPSSolve_ARPACK;
285 6 : eps->ops->setup = EPSSetUp_ARPACK;
286 6 : eps->ops->setupsort = EPSSetUpSort_Basic;
287 6 : eps->ops->destroy = EPSDestroy_ARPACK;
288 6 : eps->ops->reset = EPSReset_ARPACK;
289 6 : eps->ops->backtransform = EPSBackTransform_ARPACK;
290 6 : PetscFunctionReturn(PETSC_SUCCESS);
291 : }
|