Actual source code: blzpack.c

slepc-3.15.2 2021-09-20
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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 BLZPACK package
 12: */

 14: #include <slepc/private/epsimpl.h>
 15: #include "blzpack.h"

 17: const char* blzpack_error[33] = {
 18:   "",
 19:   "illegal data, LFLAG ",
 20:   "illegal data, dimension of (U), (V), (X) ",
 21:   "illegal data, leading dimension of (U), (V), (X) ",
 22:   "illegal data, leading dimension of (EIG) ",
 23:   "illegal data, number of required eigenpairs ",
 24:   "illegal data, Lanczos algorithm block size ",
 25:   "illegal data, maximum number of steps ",
 26:   "illegal data, number of starting vectors ",
 27:   "illegal data, number of eigenpairs provided ",
 28:   "illegal data, problem type flag ",
 29:   "illegal data, spectrum slicing flag ",
 30:   "illegal data, eigenvectors purification flag ",
 31:   "illegal data, level of output ",
 32:   "illegal data, output file unit ",
 33:   "illegal data, LCOMM (MPI or PVM) ",
 34:   "illegal data, dimension of ISTOR ",
 35:   "illegal data, convergence threshold ",
 36:   "illegal data, dimension of RSTOR ",
 37:   "illegal data on at least one PE ",
 38:   "ISTOR(3:14) must be equal on all PEs ",
 39:   "RSTOR(1:3) must be equal on all PEs ",
 40:   "not enough space in ISTOR to start eigensolution ",
 41:   "not enough space in RSTOR to start eigensolution ",
 42:   "illegal data, number of negative eigenvalues ",
 43:   "illegal data, entries of V ",
 44:   "illegal data, entries of X ",
 45:   "failure in computational subinterval ",
 46:   "file I/O error, blzpack.__.BQ ",
 47:   "file I/O error, blzpack.__.BX ",
 48:   "file I/O error, blzpack.__.Q ",
 49:   "file I/O error, blzpack.__.X ",
 50:   "parallel interface error "
 51: };

 53: PetscErrorCode EPSSetUp_BLZPACK(EPS eps)
 54: {
 56:   PetscInt       listor,lrstor,ncuv,k1,k2,k3,k4;
 57:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;
 58:   PetscBool      issinv;

 61:   EPSCheckHermitianDefinite(eps);
 62:   if (eps->ncv!=PETSC_DEFAULT) {
 63:     if (eps->ncv < PetscMin(eps->nev+10,eps->nev*2)) SETERRQ(PetscObjectComm((PetscObject)eps),0,"Warning: BLZpack recommends that ncv be larger than min(nev+10,nev*2)");
 64:   } else eps->ncv = PetscMin(eps->nev+10,eps->nev*2);
 65:   if (eps->mpd!=PETSC_DEFAULT) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
 66:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(1000,eps->n);

 68:   PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&issinv);
 69:   if (!eps->which) {
 70:     if (issinv) eps->which = EPS_TARGET_REAL;
 71:     else eps->which = EPS_SMALLEST_REAL;
 72:   }
 73:   if ((issinv && eps->which!=EPS_TARGET_REAL && eps->which!=EPS_TARGET_MAGNITUDE && eps->which!=EPS_ALL) || (!issinv && eps->which!=EPS_SMALLEST_REAL)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Wrong value of eps->which");
 74:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING);
 75:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);

 77:   if (!blz->block_size) blz->block_size = 3;
 78:   if (eps->which==EPS_ALL) {
 79:     if (eps->inta==eps->intb) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Must define a computational interval when using EPS_ALL");
 80:     blz->slice = 1;
 81:   }
 82:   if (blz->slice || eps->isgeneralized) {
 83:     if (!issinv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Shift-and-invert ST is needed for generalized problems or spectrum slicing");
 84:   }
 85:   if (blz->slice) {
 86:     if (eps->intb >= PETSC_MAX_REAL) { /* right-open interval */
 87:       if (eps->inta <= PETSC_MIN_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The defined computational interval should have at least one of their sides bounded");
 88:       STSetDefaultShift(eps->st,eps->inta);
 89:     } else {
 90:       STSetDefaultShift(eps->st,eps->intb);
 91:     }
 92:   }

 94:   k1 = PetscMin(eps->n,180);
 95:   k2 = blz->block_size;
 96:   k4 = PetscMin(eps->ncv,eps->n);
 97:   k3 = 484+k1*(13+k1*2+k2+PetscMax(18,k2+2))+k2*k2*3+k4*2;

 99:   listor = 123+k1*12;
100:   PetscFree(blz->istor);
101:   PetscMalloc1((17+listor),&blz->istor);
102:   PetscLogObjectMemory((PetscObject)eps,(17+listor)*sizeof(PetscBLASInt));
103:   PetscBLASIntCast(listor,&blz->istor[14]);

105:   if (blz->slice) lrstor = eps->nloc*(k2*4+k1*2+k4)+k3;
106:   else lrstor = eps->nloc*(k2*4+k1)+k3;
107: lrstor*=10;
108:   PetscFree(blz->rstor);
109:   PetscMalloc1((4+lrstor),&blz->rstor);
110:   PetscLogObjectMemory((PetscObject)eps,(4+lrstor)*sizeof(PetscReal));
111:   blz->rstor[3] = lrstor;

113:   ncuv = PetscMax(3,blz->block_size);
114:   PetscFree(blz->u);
115:   PetscMalloc1(ncuv*eps->nloc,&blz->u);
116:   PetscLogObjectMemory((PetscObject)eps,ncuv*eps->nloc*sizeof(PetscScalar));
117:   PetscFree(blz->v);
118:   PetscMalloc1(ncuv*eps->nloc,&blz->v);
119:   PetscLogObjectMemory((PetscObject)eps,ncuv*eps->nloc*sizeof(PetscScalar));

121:   PetscFree(blz->eig);
122:   PetscMalloc1(2*eps->ncv,&blz->eig);
123:   PetscLogObjectMemory((PetscObject)eps,2*eps->ncv*sizeof(PetscReal));

125:   EPSAllocateSolution(eps,0);
126:   EPS_SetInnerProduct(eps);
127:   return(0);
128: }

130: PetscErrorCode EPSSolve_BLZPACK(EPS eps)
131: {
133:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;
134:   PetscInt       nn;
135:   PetscBLASInt   i,nneig,lflag,nvopu;
136:   Vec            x,y;
137:   PetscScalar    sigma,*pV;
138:   Mat            A,M;
139:   KSP            ksp;
140:   PC             pc;
141:   MPI_Fint       fcomm;

144:   STGetMatrix(eps->st,0,&A);
145:   MatCreateVecsEmpty(A,&x,&y);
146:   EPSGetStartVector(eps,0,NULL);
147:   BVSetActiveColumns(eps->V,0,0);  /* just for deflation space */
148:   BVGetArray(eps->V,&pV);

150:   if (eps->isgeneralized && !blz->slice) {
151:     STGetShift(eps->st,&sigma); /* shift of origin */
152:     blz->rstor[0]  = sigma;        /* lower limit of eigenvalue interval */
153:     blz->rstor[1]  = sigma;        /* upper limit of eigenvalue interval */
154:   } else {
155:     sigma = 0.0;
156:     blz->rstor[0]  = eps->inta;    /* lower limit of eigenvalue interval */
157:     blz->rstor[1]  = eps->intb;    /* upper limit of eigenvalue interval */
158:   }
159:   nneig = 0;                       /* no. of eigs less than sigma */

161:   PetscBLASIntCast(eps->nloc,&blz->istor[0]); /* no. of rows of U, V, X */
162:   PetscBLASIntCast(eps->nloc,&blz->istor[1]); /* leading dim of U, V, X */
163:   PetscBLASIntCast(eps->nev,&blz->istor[2]);  /* required eigenpairs */
164:   PetscBLASIntCast(eps->ncv,&blz->istor[3]);  /* working eigenpairs */
165:   blz->istor[4]  = blz->block_size;    /* number of vectors in a block */
166:   blz->istor[5]  = blz->nsteps;        /* maximun number of steps per run */
167:   blz->istor[6]  = 1;                  /* number of starting vectors as input */
168:   blz->istor[7]  = 0;                  /* number of eigenpairs given as input */
169:   blz->istor[8]  = (blz->slice || eps->isgeneralized) ? 1 : 0;   /* problem type */
170:   blz->istor[9]  = blz->slice;         /* spectrum slicing */
171:   blz->istor[10] = eps->isgeneralized ? 1 : 0;   /* solutions refinement (purify) */
172:   blz->istor[11] = 0;                  /* level of printing */
173:   blz->istor[12] = 6;                  /* file unit for output */
174:   fcomm = MPI_Comm_c2f(PetscObjectComm((PetscObject)eps));
175:   blz->istor[13] = fcomm;

177:   blz->rstor[2]  = eps->tol;           /* threshold for convergence */

179:   lflag = 0;           /* reverse communication interface flag */

181:   do {
182:     BLZpack_(blz->istor,blz->rstor,&sigma,&nneig,blz->u,blz->v,&lflag,&nvopu,blz->eig,pV);

184:     switch (lflag) {
185:     case 1:
186:       /* compute v = OP u */
187:       for (i=0;i<nvopu;i++) {
188:         VecPlaceArray(x,blz->u+i*eps->nloc);
189:         VecPlaceArray(y,blz->v+i*eps->nloc);
190:         if (blz->slice || eps->isgeneralized) {
191:           STMatSolve(eps->st,x,y);
192:         } else {
193:           STApply(eps->st,x,y);
194:         }
195:         BVOrthogonalizeVec(eps->V,y,NULL,NULL,NULL);
196:         VecResetArray(x);
197:         VecResetArray(y);
198:       }
199:       /* monitor */
200:       eps->nconv  = BLZistorr_(blz->istor,"NTEIG",5);
201:       EPSMonitor(eps,eps->its,eps->nconv,blz->rstor+BLZistorr_(blz->istor,"IRITZ",5),eps->eigi,blz->rstor+BLZistorr_(blz->istor,"IRITZ",5)+BLZistorr_(blz->istor,"JT",2),BLZistorr_(blz->istor,"NRITZ",5));
202:       eps->its = eps->its + 1;
203:       if (eps->its >= eps->max_it || eps->nconv >= eps->nev) lflag = 5;
204:       break;
205:     case 2:
206:       /* compute v = B u */
207:       for (i=0;i<nvopu;i++) {
208:         VecPlaceArray(x,blz->u+i*eps->nloc);
209:         VecPlaceArray(y,blz->v+i*eps->nloc);
210:         BVApplyMatrix(eps->V,x,y);
211:         VecResetArray(x);
212:         VecResetArray(y);
213:       }
214:       break;
215:     case 3:
216:       /* update shift */
217:       PetscInfo1(eps,"Factorization update (sigma=%g)\n",(double)sigma);
218:       STSetShift(eps->st,sigma);
219:       STGetKSP(eps->st,&ksp);
220:       KSPGetPC(ksp,&pc);
221:       PCFactorGetMatrix(pc,&M);
222:       MatGetInertia(M,&nn,NULL,NULL);
223:       PetscBLASIntCast(nn,&nneig);
224:       break;
225:     case 4:
226:       /* copy the initial vector */
227:       VecPlaceArray(x,blz->v);
228:       BVCopyVec(eps->V,0,x);
229:       VecResetArray(x);
230:       break;
231:     }

233:   } while (lflag > 0);

235:   BVRestoreArray(eps->V,&pV);

237:   eps->nconv  = BLZistorr_(blz->istor,"NTEIG",5);
238:   eps->reason = EPS_CONVERGED_TOL;
239:   if (blz->slice) eps->nev = eps->nconv;
240:   for (i=0;i<eps->nconv;i++) eps->eigr[i]=blz->eig[i];

242:   if (lflag!=0) {
243:     char msg[2048] = "";
244:     for (i = 0; i < 33; i++) {
245:       if (blz->istor[15] & (1 << i)) { PetscStrlcat(msg,blzpack_error[i],sizeof(msg)); }
246:     }
247:     SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error in BLZPACK (code=%d): '%s'",blz->istor[15],msg);
248:   }
249:   VecDestroy(&x);
250:   VecDestroy(&y);
251:   return(0);
252: }

254: PetscErrorCode EPSBackTransform_BLZPACK(EPS eps)
255: {
257:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;

260:   if (!blz->slice && !eps->isgeneralized) {
261:     EPSBackTransform_Default(eps);
262:   }
263:   return(0);
264: }

266: PetscErrorCode EPSReset_BLZPACK(EPS eps)
267: {
269:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;

272:   PetscFree(blz->istor);
273:   PetscFree(blz->rstor);
274:   PetscFree(blz->u);
275:   PetscFree(blz->v);
276:   PetscFree(blz->eig);
277:   return(0);
278: }

280: PetscErrorCode EPSDestroy_BLZPACK(EPS eps)
281: {

285:   PetscFree(eps->data);
286:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetBlockSize_C",NULL);
287:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackGetBlockSize_C",NULL);
288:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetNSteps_C",NULL);
289:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackGetNSteps_C",NULL);
290:   return(0);
291: }

293: PetscErrorCode EPSView_BLZPACK(EPS eps,PetscViewer viewer)
294: {
296:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;
297:   PetscBool      isascii;

300:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
301:   if (isascii) {
302:     PetscViewerASCIIPrintf(viewer,"  block size=%d\n",blz->block_size);
303:     PetscViewerASCIIPrintf(viewer,"  maximum number of steps per run=%d\n",blz->nsteps);
304:     if (blz->slice) {
305:       PetscViewerASCIIPrintf(viewer,"  computational interval [%f,%f]\n",eps->inta,eps->intb);
306:     }
307:   }
308:   return(0);
309: }

311: PetscErrorCode EPSSetFromOptions_BLZPACK(PetscOptionItems *PetscOptionsObject,EPS eps)
312: {
314:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;
315:   PetscInt       bs,n;
316:   PetscBool      flg;

319:   PetscOptionsHead(PetscOptionsObject,"EPS BLZPACK Options");

321:     bs = blz->block_size;
322:     PetscOptionsInt("-eps_blzpack_blocksize","Block size","EPSBlzpackSetBlockSize",bs,&bs,&flg);
323:     if (flg) { EPSBlzpackSetBlockSize(eps,bs); }

325:     n = blz->nsteps;
326:     PetscOptionsInt("-eps_blzpack_nsteps","Number of steps","EPSBlzpackSetNSteps",n,&n,&flg);
327:     if (flg) { EPSBlzpackSetNSteps(eps,n); }

329:   PetscOptionsTail();
330:   return(0);
331: }

333: static PetscErrorCode EPSBlzpackSetBlockSize_BLZPACK(EPS eps,PetscInt bs)
334: {
336:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;

339:   if (bs == PETSC_DEFAULT) blz->block_size = 3;
340:   else if (bs <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Block size must be positive");
341:   else {
342:     PetscBLASIntCast(bs,&blz->block_size);
343:   }
344:   return(0);
345: }

347: /*@
348:    EPSBlzpackSetBlockSize - Sets the block size for the BLZPACK package.

350:    Collective on eps

352:    Input Parameters:
353: +  eps - the eigenproblem solver context
354: -  bs - block size

356:    Options Database Key:
357: .  -eps_blzpack_blocksize - Sets the value of the block size

359:    Level: advanced
360: @*/
361: PetscErrorCode EPSBlzpackSetBlockSize(EPS eps,PetscInt bs)
362: {

368:   PetscTryMethod(eps,"EPSBlzpackSetBlockSize_C",(EPS,PetscInt),(eps,bs));
369:   return(0);
370: }

372: static PetscErrorCode EPSBlzpackGetBlockSize_BLZPACK(EPS eps,PetscInt *bs)
373: {
374:   EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;

377:   *bs = blz->block_size;
378:   return(0);
379: }

381: /*@
382:    EPSBlzpackGetBlockSize - Gets the block size used in the BLZPACK solver.

384:    Not Collective

386:    Input Parameter:
387: .  eps - the eigenproblem solver context

389:    Output Parameter:
390: .  bs - the block size

392:    Level: advanced

394: .seealso: EPSBlzpackSetBlockSize()
395: @*/
396: PetscErrorCode EPSBlzpackGetBlockSize(EPS eps,PetscInt *bs)
397: {

403:   PetscUseMethod(eps,"EPSBlzpackGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
404:   return(0);
405: }

407: static PetscErrorCode EPSBlzpackSetNSteps_BLZPACK(EPS eps,PetscInt nsteps)
408: {
410:   EPS_BLZPACK    *blz = (EPS_BLZPACK*)eps->data;

413:   if (nsteps == PETSC_DEFAULT) blz->nsteps = 0;
414:   else {
415:     PetscBLASIntCast(nsteps,&blz->nsteps);
416:   }
417:   return(0);
418: }

420: /*@
421:    EPSBlzpackSetNSteps - Sets the maximum number of steps per run for the BLZPACK
422:    package.

424:    Collective on eps

426:    Input Parameters:
427: +  eps     - the eigenproblem solver context
428: -  nsteps  - maximum number of steps

430:    Options Database Key:
431: .  -eps_blzpack_nsteps - Sets the maximum number of steps per run

433:    Level: advanced

435: @*/
436: PetscErrorCode EPSBlzpackSetNSteps(EPS eps,PetscInt nsteps)
437: {

443:   PetscTryMethod(eps,"EPSBlzpackSetNSteps_C",(EPS,PetscInt),(eps,nsteps));
444:   return(0);
445: }

447: static PetscErrorCode EPSBlzpackGetNSteps_BLZPACK(EPS eps,PetscInt *nsteps)
448: {
449:   EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;

452:   *nsteps = blz->nsteps;
453:   return(0);
454: }

456: /*@
457:    EPSBlzpackGetNSteps - Gets the maximum number of steps per run in the BLZPACK solver.

459:    Not Collective

461:    Input Parameter:
462: .  eps - the eigenproblem solver context

464:    Output Parameter:
465: .  nsteps - the maximum number of steps

467:    Level: advanced

469: .seealso: EPSBlzpackSetNSteps()
470: @*/
471: PetscErrorCode EPSBlzpackGetNSteps(EPS eps,PetscInt *nsteps)
472: {

478:   PetscUseMethod(eps,"EPSBlzpackGetNSteps_C",(EPS,PetscInt*),(eps,nsteps));
479:   return(0);
480: }

482: SLEPC_EXTERN PetscErrorCode EPSCreate_BLZPACK(EPS eps)
483: {
485:   EPS_BLZPACK    *blzpack;

488:   PetscNewLog(eps,&blzpack);
489:   eps->data = (void*)blzpack;

491:   eps->ops->solve          = EPSSolve_BLZPACK;
492:   eps->ops->setup          = EPSSetUp_BLZPACK;
493:   eps->ops->setupsort      = EPSSetUpSort_Basic;
494:   eps->ops->setfromoptions = EPSSetFromOptions_BLZPACK;
495:   eps->ops->destroy        = EPSDestroy_BLZPACK;
496:   eps->ops->reset          = EPSReset_BLZPACK;
497:   eps->ops->view           = EPSView_BLZPACK;
498:   eps->ops->backtransform  = EPSBackTransform_BLZPACK;

500:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetBlockSize_C",EPSBlzpackSetBlockSize_BLZPACK);
501:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackGetBlockSize_C",EPSBlzpackGetBlockSize_BLZPACK);
502:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetNSteps_C",EPSBlzpackSetNSteps_BLZPACK);
503:   PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackGetNSteps_C",EPSBlzpackGetNSteps_BLZPACK);
504:   return(0);
505: }