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 8 : static PetscErrorCode EPSSetUp_ARPACK(EPS eps)
18 : {
19 8 : PetscInt ncv;
20 8 : EPS_ARPACK *ar = (EPS_ARPACK*)eps->data;
21 :
22 8 : PetscFunctionBegin;
23 8 : EPSCheckDefinite(eps);
24 8 : 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 5 : } else eps->ncv = PetscMin(PetscMax(20,2*eps->nev+1),eps->n); /* set default value of ncv */
27 8 : if (eps->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
28 8 : if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(300,(PetscInt)(2*eps->n/eps->ncv));
29 8 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
30 8 : PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
31 8 : PetscCheck(eps->which!=EPS_WHICH_USER,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support user-defined ordering of eigenvalues");
32 8 : EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING | EPS_FEATURE_TWOSIDED);
33 8 : EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION);
34 :
35 8 : ncv = eps->ncv;
36 : #if defined(PETSC_USE_COMPLEX)
37 : PetscCall(PetscFree(ar->rwork));
38 : PetscCall(PetscMalloc1(ncv,&ar->rwork));
39 : ar->lworkl = 3*ncv*ncv+5*ncv;
40 : PetscCall(PetscFree(ar->workev));
41 : PetscCall(PetscMalloc1(3*ncv,&ar->workev));
42 : #else
43 8 : if (eps->ishermitian) {
44 6 : ar->lworkl = ncv*(ncv+8);
45 : } else {
46 2 : ar->lworkl = 3*ncv*ncv+6*ncv;
47 2 : PetscCall(PetscFree(ar->workev));
48 2 : PetscCall(PetscMalloc1(3*ncv,&ar->workev));
49 : }
50 : #endif
51 8 : PetscCall(PetscFree(ar->workl));
52 8 : PetscCall(PetscMalloc1(ar->lworkl,&ar->workl));
53 8 : PetscCall(PetscFree(ar->select));
54 8 : PetscCall(PetscMalloc1(ncv,&ar->select));
55 8 : PetscCall(PetscFree(ar->workd));
56 8 : PetscCall(PetscMalloc1(3*eps->nloc,&ar->workd));
57 :
58 8 : PetscCall(EPSAllocateSolution(eps,0));
59 8 : PetscCall(EPS_SetInnerProduct(eps));
60 8 : PetscCall(EPSSetWorkVecs(eps,2));
61 8 : PetscFunctionReturn(PETSC_SUCCESS);
62 : }
63 :
64 8 : static PetscErrorCode EPSSolve_ARPACK(EPS eps)
65 : {
66 8 : EPS_ARPACK *ar = (EPS_ARPACK*)eps->data;
67 8 : char bmat[1],howmny[] = "A";
68 8 : const char *which;
69 8 : PetscInt n,ld,iparam[11],ipntr[14],ido,info,nev,ncv,rvec;
70 : #if !defined(PETSC_HAVE_MPIUNI) && !defined(PETSC_HAVE_MSMPI)
71 8 : MPI_Fint fcomm;
72 : #endif
73 8 : PetscScalar sigmar,*pV,*resid;
74 8 : Vec x,y,w = eps->work[0];
75 8 : Mat A;
76 8 : PetscBool isSinv,isShift;
77 : #if !defined(PETSC_USE_COMPLEX)
78 8 : PetscScalar sigmai = 0.0;
79 : #endif
80 :
81 8 : PetscFunctionBegin;
82 8 : nev = eps->nev;
83 8 : ncv = eps->ncv;
84 : #if !defined(PETSC_HAVE_MPIUNI) && !defined(PETSC_HAVE_MSMPI)
85 8 : fcomm = MPI_Comm_c2f(PetscObjectComm((PetscObject)eps));
86 : #endif
87 8 : n = eps->nloc;
88 8 : PetscCall(EPSGetStartVector(eps,0,NULL));
89 8 : PetscCall(BVSetActiveColumns(eps->V,0,0)); /* just for deflation space */
90 8 : PetscCall(BVCopyVec(eps->V,0,eps->work[1]));
91 8 : PetscCall(BVGetLeadingDimension(eps->V,&ld));
92 8 : PetscCall(BVGetArray(eps->V,&pV));
93 8 : PetscCall(VecGetArray(eps->work[1],&resid));
94 :
95 8 : ido = 0; /* first call to reverse communication interface */
96 8 : info = 1; /* indicates an initial vector is provided */
97 8 : iparam[0] = 1; /* use exact shifts */
98 8 : iparam[2] = eps->max_it; /* max Arnoldi iterations */
99 8 : iparam[3] = 1; /* blocksize */
100 8 : 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 8 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&isSinv));
113 8 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isShift));
114 8 : PetscCall(STGetShift(eps->st,&sigmar));
115 8 : PetscCall(STGetMatrix(eps->st,0,&A));
116 8 : PetscCall(MatCreateVecsEmpty(A,&x,&y));
117 :
118 8 : 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 7 : } 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 6 : PetscCheck(!eps->ishermitian || !eps->isgeneralized,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectral transformation not supported by ARPACK hermitian solver");
130 6 : iparam[6] = 1;
131 6 : bmat[0] = 'I';
132 : }
133 :
134 : #if !defined(PETSC_USE_COMPLEX)
135 8 : if (eps->ishermitian) {
136 6 : switch (eps->which) {
137 : case EPS_TARGET_MAGNITUDE:
138 8 : case EPS_LARGEST_MAGNITUDE: which = "LM"; break;
139 0 : case EPS_SMALLEST_MAGNITUDE: which = "SM"; break;
140 0 : case EPS_TARGET_REAL:
141 0 : case EPS_LARGEST_REAL: which = "LA"; break;
142 2 : case EPS_SMALLEST_REAL: which = "SA"; break;
143 0 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Wrong value of eps->which");
144 : }
145 : } else {
146 : #endif
147 2 : switch (eps->which) {
148 : case EPS_TARGET_MAGNITUDE:
149 : case EPS_LARGEST_MAGNITUDE: which = "LM"; break;
150 0 : case EPS_SMALLEST_MAGNITUDE: which = "SM"; break;
151 1 : case EPS_TARGET_REAL:
152 1 : case EPS_LARGEST_REAL: which = "LR"; break;
153 0 : 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 1449 : }
161 : #endif
162 :
163 1457 : do {
164 :
165 : #if !defined(PETSC_USE_COMPLEX)
166 1457 : if (eps->ishermitian) {
167 664 : 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 793 : 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 : 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 1457 : if (ido == -1 || ido == 1 || ido == 2) {
176 1449 : 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 1449 : else PetscCall(VecPlaceArray(x,&ar->workd[ipntr[0]-1]));
178 1449 : PetscCall(VecPlaceArray(y,&ar->workd[ipntr[1]-1]));
179 :
180 1449 : 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 1448 : } else if (ido == 2) {
185 : /* Y = B * X */
186 579 : PetscCall(BVApplyMatrix(eps->V,x,y));
187 : } else { /* ido == 1 */
188 869 : 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 869 : } else if (iparam[6] == 2) {
192 : /* X=A*X Y=B^-1*X for shift with B positive definite */
193 189 : PetscCall(MatMult(A,x,y));
194 189 : if (sigmar != 0.0) {
195 0 : PetscCall(BVApplyMatrix(eps->V,x,w));
196 0 : PetscCall(VecAXPY(y,sigmar,w));
197 : }
198 189 : PetscCall(VecCopy(y,x));
199 189 : PetscCall(STMatSolve(eps->st,x,y));
200 : } else {
201 : /* Y = OP * X */
202 680 : PetscCall(STApply(eps->st,x,y));
203 : }
204 869 : PetscCall(BVOrthogonalizeVec(eps->V,y,NULL,NULL,NULL));
205 : }
206 :
207 1449 : PetscCall(VecResetArray(x));
208 1449 : PetscCall(VecResetArray(y));
209 8 : } else PetscCheck(ido==99,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Internal error in ARPACK reverse communication interface (ido=%" PetscInt_FMT ")",ido);
210 :
211 1457 : } while (ido != 99);
212 :
213 8 : eps->nconv = iparam[4];
214 8 : eps->its = iparam[2];
215 :
216 8 : 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 8 : PetscCheck(info==0 || info==1,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by ARPACK subroutine xxAUPD (%" PetscInt_FMT ")",info);
218 :
219 8 : rvec = PETSC_TRUE;
220 :
221 8 : if (eps->nconv > 0) {
222 : #if !defined(PETSC_USE_COMPLEX)
223 8 : if (eps->ishermitian) {
224 6 : 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 2 : 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 : 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 8 : PetscCheck(info==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error reported by ARPACK subroutine xxEUPD (%" PetscInt_FMT ")",info);
232 : }
233 :
234 8 : PetscCall(BVRestoreArray(eps->V,&pV));
235 8 : PetscCall(VecRestoreArray(eps->work[1],&resid));
236 8 : if (eps->nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
237 0 : else eps->reason = EPS_DIVERGED_ITS;
238 :
239 8 : PetscCall(VecDestroy(&x));
240 8 : PetscCall(VecDestroy(&y));
241 8 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 8 : static PetscErrorCode EPSBackTransform_ARPACK(EPS eps)
245 : {
246 8 : PetscBool isSinv;
247 :
248 8 : PetscFunctionBegin;
249 8 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&isSinv));
250 8 : if (!isSinv) PetscCall(EPSBackTransform_Default(eps));
251 8 : PetscFunctionReturn(PETSC_SUCCESS);
252 : }
253 :
254 7 : static PetscErrorCode EPSReset_ARPACK(EPS eps)
255 : {
256 7 : EPS_ARPACK *ar = (EPS_ARPACK*)eps->data;
257 :
258 7 : PetscFunctionBegin;
259 7 : PetscCall(PetscFree(ar->workev));
260 7 : PetscCall(PetscFree(ar->workl));
261 7 : PetscCall(PetscFree(ar->select));
262 7 : PetscCall(PetscFree(ar->workd));
263 : #if defined(PETSC_USE_COMPLEX)
264 : PetscCall(PetscFree(ar->rwork));
265 : #endif
266 7 : PetscFunctionReturn(PETSC_SUCCESS);
267 : }
268 :
269 7 : static PetscErrorCode EPSDestroy_ARPACK(EPS eps)
270 : {
271 7 : PetscFunctionBegin;
272 7 : PetscCall(PetscFree(eps->data));
273 7 : PetscFunctionReturn(PETSC_SUCCESS);
274 : }
275 :
276 7 : SLEPC_EXTERN PetscErrorCode EPSCreate_ARPACK(EPS eps)
277 : {
278 7 : EPS_ARPACK *ctx;
279 :
280 7 : PetscFunctionBegin;
281 7 : PetscCall(PetscNew(&ctx));
282 7 : eps->data = (void*)ctx;
283 :
284 7 : eps->ops->solve = EPSSolve_ARPACK;
285 7 : eps->ops->setup = EPSSetUp_ARPACK;
286 7 : eps->ops->setupsort = EPSSetUpSort_Basic;
287 7 : eps->ops->destroy = EPSDestroy_ARPACK;
288 7 : eps->ops->reset = EPSReset_ARPACK;
289 7 : eps->ops->backtransform = EPSBackTransform_ARPACK;
290 7 : PetscFunctionReturn(PETSC_SUCCESS);
291 : }
|