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