Actual source code: blzpack.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepceps.h" I*/
15: #include <../src/eps/impls/external/blzpack/blzpackp.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,istrivial,flg;
61: if (eps->ncv) {
62: 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)");
63: } else eps->ncv = PetscMin(eps->nev+10,eps->nev*2);
64: if (eps->mpd) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
65: if (!eps->max_it) eps->max_it = PetscMax(1000,eps->n);
67: if (!blz->block_size) blz->block_size = 3;
68: if (!eps->ishermitian) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");
69: if (eps->which==EPS_ALL) {
70: if (eps->inta==0.0 && eps->intb==0.0) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Must define a computational interval when using EPS_ALL");
71: blz->slice = 1;
72: }
73: PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&issinv);
74: if (blz->slice || eps->isgeneralized) {
75: if (!issinv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Shift-and-invert ST is needed for generalized problems or spectrum slicing");
76: }
77: if (blz->slice) {
78: if (eps->intb >= PETSC_MAX_REAL) { /* right-open interval */
79: if (eps->inta <= PETSC_MIN_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The defined computational interval should have at least one of their sides bounded");
80: STSetDefaultShift(eps->st,eps->inta);
81: } else {
82: STSetDefaultShift(eps->st,eps->intb);
83: }
84: }
85: if (!eps->which) {
86: if (issinv) eps->which = EPS_TARGET_REAL;
87: else eps->which = EPS_SMALLEST_REAL;
88: }
89: 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),1,"Wrong value of eps->which");
90: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
92: k1 = PetscMin(eps->n,180);
93: k2 = blz->block_size;
94: k4 = PetscMin(eps->ncv,eps->n);
95: k3 = 484+k1*(13+k1*2+k2+PetscMax(18,k2+2))+k2*k2*3+k4*2;
97: listor = 123+k1*12;
98: PetscFree(blz->istor);
99: PetscMalloc1((17+listor),&blz->istor);
100: PetscLogObjectMemory((PetscObject)eps,(17+listor)*sizeof(PetscBLASInt));
101: PetscBLASIntCast(listor,&blz->istor[14]);
103: if (blz->slice) lrstor = eps->nloc*(k2*4+k1*2+k4)+k3;
104: else lrstor = eps->nloc*(k2*4+k1)+k3;
105: lrstor*=10;
106: PetscFree(blz->rstor);
107: PetscMalloc1((4+lrstor),&blz->rstor);
108: PetscLogObjectMemory((PetscObject)eps,(4+lrstor)*sizeof(PetscReal));
109: blz->rstor[3] = lrstor;
111: ncuv = PetscMax(3,blz->block_size);
112: PetscFree(blz->u);
113: PetscMalloc1(ncuv*eps->nloc,&blz->u);
114: PetscLogObjectMemory((PetscObject)eps,ncuv*eps->nloc*sizeof(PetscScalar));
115: PetscFree(blz->v);
116: PetscMalloc1(ncuv*eps->nloc,&blz->v);
117: PetscLogObjectMemory((PetscObject)eps,ncuv*eps->nloc*sizeof(PetscScalar));
119: PetscFree(blz->eig);
120: PetscMalloc1(2*eps->ncv,&blz->eig);
121: PetscLogObjectMemory((PetscObject)eps,2*eps->ncv*sizeof(PetscReal));
123: if (eps->extraction) { PetscInfo(eps,"Warning: extraction type ignored\n"); }
125: EPSAllocateSolution(eps,0);
126: EPS_SetInnerProduct(eps);
127: PetscObjectTypeCompare((PetscObject)eps->V,BVVECS,&flg);
128: if (flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires a BV with contiguous storage");
129: RGIsTrivial(eps->rg,&istrivial);
130: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support region filtering");
131: if (eps->stopping!=EPSStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"External packages do not support user-defined stopping test");
132: return(0);
133: }
135: PetscErrorCode EPSSolve_BLZPACK(EPS eps)
136: {
138: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
139: PetscInt nn;
140: PetscBLASInt i,nneig,lflag,nvopu;
141: Vec x,y,v0;
142: PetscScalar sigma,*pV;
143: Mat A,M;
144: KSP ksp;
145: PC pc;
148: STGetMatrix(eps->st,0,&A);
149: MatCreateVecsEmpty(A,&x,&y);
150: EPSGetStartVector(eps,0,NULL);
151: BVSetActiveColumns(eps->V,0,0); /* just for deflation space */
152: BVGetColumn(eps->V,0,&v0);
153: VecGetArray(v0,&pV);
155: if (eps->isgeneralized && !blz->slice) {
156: STGetShift(eps->st,&sigma); /* shift of origin */
157: blz->rstor[0] = sigma; /* lower limit of eigenvalue interval */
158: blz->rstor[1] = sigma; /* upper limit of eigenvalue interval */
159: } else {
160: sigma = 0.0;
161: blz->rstor[0] = eps->inta; /* lower limit of eigenvalue interval */
162: blz->rstor[1] = eps->intb; /* upper limit of eigenvalue interval */
163: }
164: nneig = 0; /* no. of eigs less than sigma */
166: PetscBLASIntCast(eps->nloc,&blz->istor[0]); /* no. of rows of U, V, X */
167: PetscBLASIntCast(eps->nloc,&blz->istor[1]); /* leading dim of U, V, X */
168: PetscBLASIntCast(eps->nev,&blz->istor[2]); /* required eigenpairs */
169: PetscBLASIntCast(eps->ncv,&blz->istor[3]); /* working eigenpairs */
170: blz->istor[4] = blz->block_size; /* number of vectors in a block */
171: blz->istor[5] = blz->nsteps; /* maximun number of steps per run */
172: blz->istor[6] = 1; /* number of starting vectors as input */
173: blz->istor[7] = 0; /* number of eigenpairs given as input */
174: blz->istor[8] = (blz->slice || eps->isgeneralized) ? 1 : 0; /* problem type */
175: blz->istor[9] = blz->slice; /* spectrum slicing */
176: blz->istor[10] = eps->isgeneralized ? 1 : 0; /* solutions refinement (purify) */
177: blz->istor[11] = 0; /* level of printing */
178: blz->istor[12] = 6; /* file unit for output */
179: PetscBLASIntCast(MPI_Comm_c2f(PetscObjectComm((PetscObject)eps)),&blz->istor[13]);
181: blz->rstor[2] = eps->tol; /* threshold for convergence */
183: lflag = 0; /* reverse communication interface flag */
185: do {
186: BLZpack_(blz->istor,blz->rstor,&sigma,&nneig,blz->u,blz->v,&lflag,&nvopu,blz->eig,pV);
188: switch (lflag) {
189: case 1:
190: /* compute v = OP u */
191: for (i=0;i<nvopu;i++) {
192: VecPlaceArray(x,blz->u+i*eps->nloc);
193: VecPlaceArray(y,blz->v+i*eps->nloc);
194: if (blz->slice || eps->isgeneralized) {
195: STMatSolve(eps->st,x,y);
196: } else {
197: STApply(eps->st,x,y);
198: }
199: BVOrthogonalizeVec(eps->V,y,NULL,NULL,NULL);
200: VecResetArray(x);
201: VecResetArray(y);
202: }
203: /* monitor */
204: eps->nconv = BLZistorr_(blz->istor,"NTEIG",5);
205: EPSMonitor(eps,eps->its,eps->nconv,
206: blz->rstor+BLZistorr_(blz->istor,"IRITZ",5),
207: eps->eigi,
208: blz->rstor+BLZistorr_(blz->istor,"IRITZ",5)+BLZistorr_(blz->istor,"JT",2),
209: BLZistorr_(blz->istor,"NRITZ",5));
210: eps->its = eps->its + 1;
211: if (eps->its >= eps->max_it || eps->nconv >= eps->nev) lflag = 5;
212: break;
213: case 2:
214: /* compute v = B u */
215: for (i=0;i<nvopu;i++) {
216: VecPlaceArray(x,blz->u+i*eps->nloc);
217: VecPlaceArray(y,blz->v+i*eps->nloc);
218: BVApplyMatrix(eps->V,x,y);
219: VecResetArray(x);
220: VecResetArray(y);
221: }
222: break;
223: case 3:
224: /* update shift */
225: PetscInfo1(eps,"Factorization update (sigma=%g)\n",(double)sigma);
226: STSetShift(eps->st,sigma);
227: STGetKSP(eps->st,&ksp);
228: KSPGetPC(ksp,&pc);
229: PCFactorGetMatrix(pc,&M);
230: MatGetInertia(M,&nn,NULL,NULL);
231: PetscBLASIntCast(nn,&nneig);
232: break;
233: case 4:
234: /* copy the initial vector */
235: VecPlaceArray(x,blz->v);
236: VecCopy(v0,x);
237: VecResetArray(x);
238: break;
239: }
241: } while (lflag > 0);
243: VecRestoreArray(v0,&pV);
244: BVRestoreColumn(eps->V,0,&v0);
246: eps->nconv = BLZistorr_(blz->istor,"NTEIG",5);
247: eps->reason = EPS_CONVERGED_TOL;
249: for (i=0;i<eps->nconv;i++) {
250: eps->eigr[i]=blz->eig[i];
251: }
253: if (lflag!=0) {
254: char msg[2048] = "";
255: for (i = 0; i < 33; i++) {
256: if (blz->istor[15] & (1 << i)) PetscStrcat(msg,blzpack_error[i]);
257: }
258: SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error in BLZPACK (code=%d): '%s'",blz->istor[15],msg);
259: }
260: VecDestroy(&x);
261: VecDestroy(&y);
262: return(0);
263: }
265: PetscErrorCode EPSBackTransform_BLZPACK(EPS eps)
266: {
268: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
271: if (!blz->slice && !eps->isgeneralized) {
272: EPSBackTransform_Default(eps);
273: }
274: return(0);
275: }
277: PetscErrorCode EPSReset_BLZPACK(EPS eps)
278: {
280: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
283: PetscFree(blz->istor);
284: PetscFree(blz->rstor);
285: PetscFree(blz->u);
286: PetscFree(blz->v);
287: PetscFree(blz->eig);
288: return(0);
289: }
291: PetscErrorCode EPSDestroy_BLZPACK(EPS eps)
292: {
296: PetscFree(eps->data);
297: PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetBlockSize_C",NULL);
298: PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetNSteps_C",NULL);
299: return(0);
300: }
302: PetscErrorCode EPSView_BLZPACK(EPS eps,PetscViewer viewer)
303: {
305: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
306: PetscBool isascii;
309: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
310: if (isascii) {
311: PetscViewerASCIIPrintf(viewer," block size=%d\n",blz->block_size);
312: PetscViewerASCIIPrintf(viewer," maximum number of steps per run=%d\n",blz->nsteps);
313: if (blz->slice) {
314: PetscViewerASCIIPrintf(viewer," computational interval [%f,%f]\n",eps->inta,eps->intb);
315: }
316: }
317: return(0);
318: }
320: PetscErrorCode EPSSetFromOptions_BLZPACK(PetscOptionItems *PetscOptionsObject,EPS eps)
321: {
323: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
324: PetscInt bs,n;
325: PetscBool flg;
328: PetscOptionsHead(PetscOptionsObject,"EPS BLZPACK Options");
330: bs = blz->block_size;
331: PetscOptionsInt("-eps_blzpack_blocksize","Block size","EPSBlzpackSetBlockSize",bs,&bs,&flg);
332: if (flg) { EPSBlzpackSetBlockSize(eps,bs); }
334: n = blz->nsteps;
335: PetscOptionsInt("-eps_blzpack_nsteps","Number of steps","EPSBlzpackSetNSteps",n,&n,&flg);
336: if (flg) { EPSBlzpackSetNSteps(eps,n); }
338: PetscOptionsTail();
339: return(0);
340: }
342: static PetscErrorCode EPSBlzpackSetBlockSize_BLZPACK(EPS eps,PetscInt bs)
343: {
345: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
348: if (bs == PETSC_DEFAULT) blz->block_size = 3;
349: else if (bs <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Block size must be positive");
350: else {
351: PetscBLASIntCast(bs,&blz->block_size);
352: }
353: return(0);
354: }
356: /*@
357: EPSBlzpackSetBlockSize - Sets the block size for the BLZPACK package.
359: Collective on EPS
361: Input Parameters:
362: + eps - the eigenproblem solver context
363: - bs - block size
365: Options Database Key:
366: . -eps_blzpack_blocksize - Sets the value of the block size
368: Level: advanced
369: @*/
370: PetscErrorCode EPSBlzpackSetBlockSize(EPS eps,PetscInt bs)
371: {
377: PetscTryMethod(eps,"EPSBlzpackSetBlockSize_C",(EPS,PetscInt),(eps,bs));
378: return(0);
379: }
381: static PetscErrorCode EPSBlzpackSetNSteps_BLZPACK(EPS eps,PetscInt nsteps)
382: {
384: EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;
387: if (nsteps == PETSC_DEFAULT) blz->nsteps = 0;
388: else {
389: PetscBLASIntCast(nsteps,&blz->nsteps);
390: }
391: return(0);
392: }
394: /*@
395: EPSBlzpackSetNSteps - Sets the maximum number of steps per run for the BLZPACK
396: package.
398: Collective on EPS
400: Input Parameters:
401: + eps - the eigenproblem solver context
402: - nsteps - maximum number of steps
404: Options Database Key:
405: . -eps_blzpack_nsteps - Sets the maximum number of steps per run
407: Level: advanced
409: @*/
410: PetscErrorCode EPSBlzpackSetNSteps(EPS eps,PetscInt nsteps)
411: {
417: PetscTryMethod(eps,"EPSBlzpackSetNSteps_C",(EPS,PetscInt),(eps,nsteps));
418: return(0);
419: }
421: SLEPC_EXTERN PetscErrorCode EPSCreate_BLZPACK(EPS eps)
422: {
424: EPS_BLZPACK *blzpack;
427: PetscNewLog(eps,&blzpack);
428: eps->data = (void*)blzpack;
430: eps->ops->solve = EPSSolve_BLZPACK;
431: eps->ops->setup = EPSSetUp_BLZPACK;
432: eps->ops->setfromoptions = EPSSetFromOptions_BLZPACK;
433: eps->ops->destroy = EPSDestroy_BLZPACK;
434: eps->ops->reset = EPSReset_BLZPACK;
435: eps->ops->view = EPSView_BLZPACK;
436: eps->ops->backtransform = EPSBackTransform_BLZPACK;
438: PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetBlockSize_C",EPSBlzpackSetBlockSize_BLZPACK);
439: PetscObjectComposeFunction((PetscObject)eps,"EPSBlzpackSetNSteps_C",EPSBlzpackSetNSteps_BLZPACK);
440: return(0);
441: }