Actual source code: arpack.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  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:   PetscLogObjectMemory((PetscObject)eps,ncv*sizeof(PetscReal));
 39:   ar->lworkl = 3*ncv*ncv+5*ncv;
 40:   PetscFree(ar->workev);
 41:   PetscMalloc1(3*ncv,&ar->workev);
 42:   PetscLogObjectMemory((PetscObject)eps,3*ncv*sizeof(PetscScalar));
 43: #else
 44:   if (eps->ishermitian) {
 45:     ar->lworkl = ncv*(ncv+8);
 46:   } else {
 47:     ar->lworkl = 3*ncv*ncv+6*ncv;
 48:     PetscFree(ar->workev);
 49:     PetscMalloc1(3*ncv,&ar->workev);
 50:     PetscLogObjectMemory((PetscObject)eps,3*ncv*sizeof(PetscScalar));
 51:   }
 52: #endif
 53:   PetscFree(ar->workl);
 54:   PetscMalloc1(ar->lworkl,&ar->workl);
 55:   PetscLogObjectMemory((PetscObject)eps,ar->lworkl*sizeof(PetscScalar));
 56:   PetscFree(ar->select);
 57:   PetscMalloc1(ncv,&ar->select);
 58:   PetscLogObjectMemory((PetscObject)eps,ncv*sizeof(PetscInt));
 59:   PetscFree(ar->workd);
 60:   PetscMalloc1(3*eps->nloc,&ar->workd);
 61:   PetscLogObjectMemory((PetscObject)eps,3*eps->nloc*sizeof(PetscScalar));

 63:   EPSAllocateSolution(eps,0);
 64:   EPS_SetInnerProduct(eps);
 65:   EPSSetWorkVecs(eps,2);
 66:   PetscFunctionReturn(0);
 67: }

 69: PetscErrorCode EPSSolve_ARPACK(EPS eps)
 70: {
 71:   EPS_ARPACK     *ar = (EPS_ARPACK*)eps->data;
 72:   char           bmat[1],howmny[] = "A";
 73:   const char     *which;
 74:   PetscInt       n,iparam[11],ipntr[14],ido,info,nev,ncv,rvec;
 75: #if !defined(PETSC_HAVE_MPIUNI) && !defined(PETSC_HAVE_MSMPI)
 76:   MPI_Fint       fcomm;
 77: #endif
 78:   PetscScalar    sigmar,*pV,*resid;
 79:   Vec            x,y,w = eps->work[0];
 80:   Mat            A;
 81:   PetscBool      isSinv,isShift;
 82: #if !defined(PETSC_USE_COMPLEX)
 83:   PetscScalar    sigmai = 0.0;
 84: #endif

 86:   nev = eps->nev;
 87:   ncv = eps->ncv;
 88: #if !defined(PETSC_HAVE_MPIUNI) && !defined(PETSC_HAVE_MSMPI)
 89:   fcomm = MPI_Comm_c2f(PetscObjectComm((PetscObject)eps));
 90: #endif
 91:   n = eps->nloc;
 92:   EPSGetStartVector(eps,0,NULL);
 93:   BVSetActiveColumns(eps->V,0,0);  /* just for deflation space */
 94:   BVCopyVec(eps->V,0,eps->work[1]);
 95:   BVGetArray(eps->V,&pV);
 96:   VecGetArray(eps->work[1],&resid);

 98:   ido  = 0;            /* first call to reverse communication interface */
 99:   info = 1;            /* indicates an initial vector is provided */
100:   iparam[0] = 1;       /* use exact shifts */
101:   iparam[2] = eps->max_it;  /* max Arnoldi iterations */
102:   iparam[3] = 1;       /* blocksize */
103:   iparam[4] = 0;       /* number of converged Ritz values */

105:   /*
106:      Computational modes ([]=not supported):
107:             symmetric    non-symmetric    complex
108:         1     1  'I'        1  'I'         1  'I'
109:         2     3  'I'        3  'I'         3  'I'
110:         3     2  'G'        2  'G'         2  'G'
111:         4     3  'G'        3  'G'         3  'G'
112:         5   [ 4  'G' ]    [ 3  'G' ]
113:         6   [ 5  'G' ]    [ 4  'G' ]
114:    */
115:   PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&isSinv);
116:   PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isShift);
117:   STGetShift(eps->st,&sigmar);
118:   STGetMatrix(eps->st,0,&A);
119:   MatCreateVecsEmpty(A,&x,&y);

121:   if (isSinv) {
122:     /* shift-and-invert mode */
123:     iparam[6] = 3;
124:     if (eps->ispositive) bmat[0] = 'G';
125:     else bmat[0] = 'I';
126:   } else if (isShift && eps->ispositive) {
127:     /* generalized shift mode with B positive definite */
128:     iparam[6] = 2;
129:     bmat[0] = 'G';
130:   } else {
131:     /* regular mode */
133:     iparam[6] = 1;
134:     bmat[0] = 'I';
135:   }

137: #if !defined(PETSC_USE_COMPLEX)
138:   if (eps->ishermitian) {
139:     switch (eps->which) {
140:       case EPS_TARGET_MAGNITUDE:
141:       case EPS_LARGEST_MAGNITUDE:  which = "LM"; break;
142:       case EPS_SMALLEST_MAGNITUDE: which = "SM"; break;
143:       case EPS_TARGET_REAL:
144:       case EPS_LARGEST_REAL:       which = "LA"; break;
145:       case EPS_SMALLEST_REAL:      which = "SA"; break;
146:       default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Wrong value of eps->which");
147:     }
148:   } else {
149: #endif
150:     switch (eps->which) {
151:       case EPS_TARGET_MAGNITUDE:
152:       case EPS_LARGEST_MAGNITUDE:  which = "LM"; break;
153:       case EPS_SMALLEST_MAGNITUDE: which = "SM"; break;
154:       case EPS_TARGET_REAL:
155:       case EPS_LARGEST_REAL:       which = "LR"; break;
156:       case EPS_SMALLEST_REAL:      which = "SR"; break;
157:       case EPS_TARGET_IMAGINARY:
158:       case EPS_LARGEST_IMAGINARY:  which = "LI"; break;
159:       case EPS_SMALLEST_IMAGINARY: which = "SI"; break;
160:       default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Wrong value of eps->which");
161:     }
162: #if !defined(PETSC_USE_COMPLEX)
163:   }
164: #endif

166:   do {

168: #if !defined(PETSC_USE_COMPLEX)
169:     if (eps->ishermitian) {
170:       PetscStackCall("ARPACKsaupd",ARPACKsaupd_(&fcomm,&ido,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&n,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,&info));
171:     } else {
172:       PetscStackCall("ARPACKnaupd",ARPACKnaupd_(&fcomm,&ido,bmat,&n,which,&nev,&eps->tol,resid,&ncv,pV,&n,iparam,ipntr,ar->workd,ar->workl,&ar->lworkl,&info));
173:     }
174: #else
175:     PetscStackCall("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));
176: #endif

178:     if (ido == -1 || ido == 1 || ido == 2) {
179:       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*/
180:       else VecPlaceArray(x,&ar->workd[ipntr[0]-1]);
181:       VecPlaceArray(y,&ar->workd[ipntr[1]-1]);

183:       if (ido == -1) {
184:         /* Y = OP * X for for the initialization phase to
185:            force the starting vector into the range of OP */
186:         STApply(eps->st,x,y);
187:       } else if (ido == 2) {
188:         /* Y = B * X */
189:         BVApplyMatrix(eps->V,x,y);
190:       } else { /* ido == 1 */
191:         if (iparam[6] == 3 && bmat[0] == 'G') {
192:           /* Y = OP * X for shift-and-invert with B semi-positive definite */
193:           STMatSolve(eps->st,x,y);
194:         } else if (iparam[6] == 2) {
195:           /* X=A*X Y=B^-1*X for shift with B positive definite */
196:           MatMult(A,x,y);
197:           if (sigmar != 0.0) {
198:             BVApplyMatrix(eps->V,x,w);
199:             VecAXPY(y,sigmar,w);
200:           }
201:           VecCopy(y,x);
202:           STMatSolve(eps->st,x,y);
203:         } else {
204:           /* Y = OP * X */
205:           STApply(eps->st,x,y);
206:         }
207:         BVOrthogonalizeVec(eps->V,y,NULL,NULL,NULL);
208:       }

210:       VecResetArray(x);
211:       VecResetArray(y);

214:   } while (ido != 99);

216:   eps->nconv = iparam[4];
217:   eps->its = iparam[2];


222:   rvec = PETSC_TRUE;

224:   if (eps->nconv > 0) {
225: #if !defined(PETSC_USE_COMPLEX)
226:     if (eps->ishermitian) {
227:       PetscStackCall("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));
228:     } else {
229:       PetscStackCall("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));
230:     }
231: #else
232:     PetscStackCall("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));
233: #endif
235:   }

237:   BVRestoreArray(eps->V,&pV);
238:   VecRestoreArray(eps->work[1],&resid);
239:   if (eps->nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
240:   else eps->reason = EPS_DIVERGED_ITS;

242:   VecDestroy(&x);
243:   VecDestroy(&y);
244:   PetscFunctionReturn(0);
245: }

247: PetscErrorCode EPSBackTransform_ARPACK(EPS eps)
248: {
249:   PetscBool      isSinv;

251:   PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&isSinv);
252:   if (!isSinv) EPSBackTransform_Default(eps);
253:   PetscFunctionReturn(0);
254: }

256: PetscErrorCode EPSReset_ARPACK(EPS eps)
257: {
258:   EPS_ARPACK     *ar = (EPS_ARPACK*)eps->data;

260:   PetscFree(ar->workev);
261:   PetscFree(ar->workl);
262:   PetscFree(ar->select);
263:   PetscFree(ar->workd);
264: #if defined(PETSC_USE_COMPLEX)
265:   PetscFree(ar->rwork);
266: #endif
267:   PetscFunctionReturn(0);
268: }

270: PetscErrorCode EPSDestroy_ARPACK(EPS eps)
271: {
272:   PetscFree(eps->data);
273:   PetscFunctionReturn(0);
274: }

276: SLEPC_EXTERN PetscErrorCode EPSCreate_ARPACK(EPS eps)
277: {
278:   EPS_ARPACK     *ctx;

280:   PetscNewLog(eps,&ctx);
281:   eps->data = (void*)ctx;

283:   eps->ops->solve          = EPSSolve_ARPACK;
284:   eps->ops->setup          = EPSSetUp_ARPACK;
285:   eps->ops->setupsort      = EPSSetUpSort_Basic;
286:   eps->ops->destroy        = EPSDestroy_ARPACK;
287:   eps->ops->reset          = EPSReset_ARPACK;
288:   eps->ops->backtransform  = EPSBackTransform_ARPACK;
289:   PetscFunctionReturn(0);
290: }