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 PRIMME package
12 : */
13 :
14 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 :
16 : #include <primme.h>
17 :
18 : #if defined(PETSC_USE_COMPLEX)
19 : #if defined(PETSC_USE_REAL_SINGLE)
20 : #define PRIMME_DRIVER cprimme
21 : #else
22 : #define PRIMME_DRIVER zprimme
23 : #endif
24 : #else
25 : #if defined(PETSC_USE_REAL_SINGLE)
26 : #define PRIMME_DRIVER sprimme
27 : #else
28 : #define PRIMME_DRIVER dprimme
29 : #endif
30 : #endif
31 :
32 : #if defined(PRIMME_VERSION_MAJOR) && PRIMME_VERSION_MAJOR*100+PRIMME_VERSION_MINOR >= 202
33 : #define SLEPC_HAVE_PRIMME2p2
34 : #endif
35 :
36 : typedef struct {
37 : primme_params primme; /* param struct */
38 : PetscInt bs; /* block size */
39 : primme_preset_method method; /* primme method */
40 : Mat A,B; /* problem matrices */
41 : KSP ksp; /* linear solver and preconditioner */
42 : Vec x,y; /* auxiliary vectors */
43 : double target; /* a copy of eps's target */
44 : } EPS_PRIMME;
45 :
46 852 : static void par_GlobalSumReal(void *sendBuf,void *recvBuf,int *count,primme_params *primme,int *ierr)
47 : {
48 852 : if (sendBuf == recvBuf) {
49 1704 : *ierr = MPI_Allreduce(MPI_IN_PLACE,recvBuf,*count,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)primme->commInfo));
50 : } else {
51 0 : *ierr = MPI_Allreduce(sendBuf,recvBuf,*count,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)primme->commInfo));
52 : }
53 852 : }
54 :
55 : #if defined(SLEPC_HAVE_PRIMME3)
56 2252 : static void par_broadcastReal(void *buf,int *count,primme_params *primme,int *ierr)
57 : {
58 2252 : *ierr = MPI_Bcast(buf,*count,MPIU_REAL,0/*root*/,PetscObjectComm((PetscObject)primme->commInfo));
59 2252 : }
60 : #endif
61 :
62 : #if defined(SLEPC_HAVE_PRIMME2p2)
63 7413 : static void convTestFun(double *eval,void *evec,double *resNorm,int *isconv,primme_params *primme,int *err)
64 : {
65 7413 : PetscErrorCode ierr;
66 7413 : EPS eps = (EPS)primme->commInfo;
67 7413 : PetscScalar eigvr = eval?*eval:0.0;
68 7413 : PetscReal r = resNorm?*resNorm:HUGE_VAL,errest;
69 :
70 7413 : ierr = (*eps->converged)(eps,eigvr,0.0,r,&errest,eps->convergedctx);
71 7413 : if (ierr) *err = 1;
72 : else {
73 7413 : *isconv = (errest<=eps->tol?1:0);
74 7413 : *err = 0;
75 : }
76 7413 : }
77 :
78 6591 : static void monitorFun(void *basisEvals,int *basisSize,int *basisFlags,int *iblock,int *blockSize,void *basisNorms,int *numConverged,void *lockedEvals,int *numLocked,int *lockedFlags,void *lockedNorms,int *inner_its,void *LSRes,
79 : #if defined(SLEPC_HAVE_PRIMME3)
80 : const char *msg,double *time,
81 : #endif
82 : primme_event *event,struct primme_params *primme,int *err)
83 : {
84 6591 : PetscErrorCode ierr = PETSC_SUCCESS;
85 6591 : EPS eps = (EPS)primme->commInfo;
86 6591 : PetscInt i,k,nerrest;
87 :
88 6591 : switch (*event) {
89 795 : case primme_event_outer_iteration:
90 : /* Update EPS */
91 795 : PetscCallVoid(PetscIntCast(primme->stats.numOuterIterations,&eps->its));
92 795 : eps->nconv = primme->initSize;
93 795 : k=0;
94 1411 : if (lockedEvals && numLocked) for (i=0; i<*numLocked && k<eps->ncv; i++) eps->eigr[k++] = ((PetscReal*)lockedEvals)[i];
95 795 : nerrest = k;
96 795 : if (iblock && blockSize) {
97 2112 : for (i=0; i<*blockSize && k+iblock[i]<eps->ncv; i++) eps->errest[k+iblock[i]] = ((PetscReal*)basisNorms)[i];
98 795 : nerrest = k+(*blockSize>0?1+iblock[*blockSize-1]:0);
99 : }
100 16867 : if (basisEvals && basisSize) for (i=0; i<*basisSize && k<eps->ncv; i++) eps->eigr[k++] = ((PetscReal*)basisEvals)[i];
101 : /* Show progress */
102 795 : ierr = EPSMonitor(eps,eps->its,numConverged?*numConverged:0,eps->eigr,eps->eigi,eps->errest,nerrest);
103 795 : break;
104 : #if defined(SLEPC_HAVE_PRIMME3)
105 0 : case primme_event_message:
106 : /* Print PRIMME information messages */
107 0 : ierr = PetscInfo(eps,"%s\n",msg);
108 0 : break;
109 : #endif
110 : default:
111 : break;
112 : }
113 6591 : *err = (ierr!=0)? 1: 0;
114 : }
115 : #endif /* SLEPC_HAVE_PRIMME2p2 */
116 :
117 5940 : static void matrixMatvec_PRIMME(void *xa,PRIMME_INT *ldx,void *ya,PRIMME_INT *ldy,int *blockSize,struct primme_params *primme,int *ierr)
118 : {
119 5940 : PetscInt i;
120 5940 : EPS_PRIMME *ops = (EPS_PRIMME*)primme->matrix;
121 5940 : Vec x = ops->x,y = ops->y;
122 5940 : Mat A = ops->A;
123 :
124 5940 : PetscFunctionBegin;
125 13752 : for (i=0;i<*blockSize;i++) {
126 7812 : PetscCallAbort(PetscObjectComm((PetscObject)A),VecPlaceArray(x,(PetscScalar*)xa+(*ldx)*i));
127 7812 : PetscCallAbort(PetscObjectComm((PetscObject)A),VecPlaceArray(y,(PetscScalar*)ya+(*ldy)*i));
128 7812 : PetscCallAbort(PetscObjectComm((PetscObject)A),MatMult(A,x,y));
129 7812 : PetscCallAbort(PetscObjectComm((PetscObject)A),VecResetArray(x));
130 7812 : PetscCallAbort(PetscObjectComm((PetscObject)A),VecResetArray(y));
131 : }
132 5940 : PetscFunctionReturnVoid();
133 : }
134 :
135 : #if defined(SLEPC_HAVE_PRIMME3)
136 5506 : static void massMatrixMatvec_PRIMME(void *xa,PRIMME_INT *ldx,void *ya,PRIMME_INT *ldy,int *blockSize,struct primme_params *primme,int *ierr)
137 : {
138 5506 : PetscInt i;
139 5506 : EPS_PRIMME *ops = (EPS_PRIMME*)primme->massMatrix;
140 5506 : Vec x = ops->x,y = ops->y;
141 5506 : Mat B = ops->B;
142 :
143 5506 : PetscFunctionBegin;
144 11012 : for (i=0;i<*blockSize;i++) {
145 5506 : PetscCallAbort(PetscObjectComm((PetscObject)B),VecPlaceArray(x,(PetscScalar*)xa+(*ldx)*i));
146 5506 : PetscCallAbort(PetscObjectComm((PetscObject)B),VecPlaceArray(y,(PetscScalar*)ya+(*ldy)*i));
147 5506 : PetscCallAbort(PetscObjectComm((PetscObject)B),MatMult(B,x,y));
148 5506 : PetscCallAbort(PetscObjectComm((PetscObject)B),VecResetArray(x));
149 5506 : PetscCallAbort(PetscObjectComm((PetscObject)B),VecResetArray(y));
150 : }
151 5506 : PetscFunctionReturnVoid();
152 : }
153 : #endif
154 :
155 5336 : static void applyPreconditioner_PRIMME(void *xa,PRIMME_INT *ldx,void *ya,PRIMME_INT *ldy,int *blockSize,struct primme_params *primme,int *ierr)
156 : {
157 5336 : PetscInt i;
158 5336 : EPS_PRIMME *ops = (EPS_PRIMME*)primme->matrix;
159 5336 : Vec x = ops->x,y = ops->y;
160 :
161 5336 : PetscFunctionBegin;
162 12375 : for (i=0;i<*blockSize;i++) {
163 7039 : PetscCallAbort(PetscObjectComm((PetscObject)ops->ksp),VecPlaceArray(x,(PetscScalar*)xa+(*ldx)*i));
164 7039 : PetscCallAbort(PetscObjectComm((PetscObject)ops->ksp),VecPlaceArray(y,(PetscScalar*)ya+(*ldy)*i));
165 7039 : PetscCallAbort(PetscObjectComm((PetscObject)ops->ksp),KSPSolve(ops->ksp,x,y));
166 7039 : PetscCallAbort(PetscObjectComm((PetscObject)ops->ksp),VecResetArray(x));
167 7039 : PetscCallAbort(PetscObjectComm((PetscObject)ops->ksp),VecResetArray(y));
168 : }
169 5336 : PetscFunctionReturnVoid();
170 : }
171 :
172 10 : static PetscErrorCode EPSSetUp_PRIMME(EPS eps)
173 : {
174 10 : PetscMPIInt numProcs,procID;
175 10 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
176 10 : primme_params *primme = &ops->primme;
177 10 : PetscBool flg;
178 :
179 10 : PetscFunctionBegin;
180 10 : EPSCheckHermitianDefinite(eps);
181 10 : EPSCheckNotStructured(eps);
182 10 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&numProcs));
183 10 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&procID));
184 :
185 : /* Check some constraints and set some default values */
186 10 : if (eps->max_it==PETSC_DETERMINE) eps->max_it = PETSC_INT_MAX;
187 10 : PetscCall(STGetMatrix(eps->st,0,&ops->A));
188 10 : if (eps->isgeneralized) {
189 : #if defined(SLEPC_HAVE_PRIMME3)
190 1 : PetscCall(STGetMatrix(eps->st,1,&ops->B));
191 : #else
192 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This version of PRIMME is not available for generalized problems");
193 : #endif
194 : }
195 10 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
196 10 : EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
197 10 : if (!eps->which) eps->which = EPS_LARGEST_REAL;
198 : #if !defined(SLEPC_HAVE_PRIMME2p2)
199 : if (eps->converged != EPSConvergedAbsolute) PetscCall(PetscInfo(eps,"Warning: using absolute convergence test\n"));
200 : #else
201 10 : EPSCheckIgnored(eps,EPS_FEATURE_CONVERGENCE);
202 : #endif
203 :
204 : /* Transfer SLEPc options to PRIMME options */
205 10 : primme_free(primme);
206 10 : primme_initialize(primme);
207 10 : primme->n = (PRIMME_INT)eps->n;
208 10 : primme->nLocal = (PRIMME_INT)eps->nloc;
209 10 : primme->numEvals = (int)eps->nev;
210 10 : primme->matrix = ops;
211 10 : primme->matrixMatvec = matrixMatvec_PRIMME;
212 : #if defined(SLEPC_HAVE_PRIMME3)
213 10 : if (eps->isgeneralized) {
214 1 : primme->massMatrix = ops;
215 1 : primme->massMatrixMatvec = massMatrixMatvec_PRIMME;
216 : }
217 : #endif
218 10 : primme->commInfo = eps;
219 10 : primme->maxOuterIterations = (PRIMME_INT)eps->max_it;
220 : #if !defined(SLEPC_HAVE_PRIMME2p2)
221 : primme->eps = SlepcDefaultTol(eps->tol);
222 : #endif
223 10 : primme->numProcs = numProcs;
224 10 : primme->procID = procID;
225 10 : primme->printLevel = 1;
226 10 : primme->correctionParams.precondition = 1;
227 10 : primme->globalSumReal = par_GlobalSumReal;
228 : #if defined(SLEPC_HAVE_PRIMME3)
229 10 : primme->broadcastReal = par_broadcastReal;
230 : #endif
231 : #if defined(SLEPC_HAVE_PRIMME2p2)
232 10 : primme->convTestFun = convTestFun;
233 10 : primme->monitorFun = monitorFun;
234 : #endif
235 10 : if (ops->bs > 0) primme->maxBlockSize = (int)ops->bs;
236 :
237 10 : switch (eps->which) {
238 3 : case EPS_LARGEST_REAL:
239 3 : primme->target = primme_largest;
240 3 : break;
241 2 : case EPS_SMALLEST_REAL:
242 2 : primme->target = primme_smallest;
243 2 : break;
244 1 : case EPS_LARGEST_MAGNITUDE:
245 1 : primme->target = primme_largest_abs;
246 1 : ops->target = 0.0;
247 1 : primme->numTargetShifts = 1;
248 1 : primme->targetShifts = &ops->target;
249 1 : break;
250 1 : case EPS_SMALLEST_MAGNITUDE:
251 1 : primme->target = primme_closest_abs;
252 1 : ops->target = 0.0;
253 1 : primme->numTargetShifts = 1;
254 1 : primme->targetShifts = &ops->target;
255 1 : break;
256 3 : case EPS_TARGET_MAGNITUDE:
257 : case EPS_TARGET_REAL:
258 3 : primme->target = primme_closest_abs;
259 3 : primme->numTargetShifts = 1;
260 3 : ops->target = PetscRealPart(eps->target);
261 3 : primme->targetShifts = &ops->target;
262 3 : break;
263 0 : default:
264 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"'which' value not supported by PRIMME");
265 : }
266 :
267 10 : switch (eps->extraction) {
268 7 : case EPS_RITZ:
269 7 : primme->projectionParams.projection = primme_proj_RR;
270 7 : break;
271 1 : case EPS_HARMONIC:
272 1 : primme->projectionParams.projection = primme_proj_harmonic;
273 1 : break;
274 2 : case EPS_REFINED:
275 2 : primme->projectionParams.projection = primme_proj_refined;
276 2 : break;
277 0 : default:
278 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"'extraction' value not supported by PRIMME");
279 : }
280 :
281 : /* If user sets mpd or ncv, maxBasisSize is modified */
282 10 : if (eps->mpd!=PETSC_DETERMINE) {
283 1 : primme->maxBasisSize = (int)eps->mpd;
284 1 : if (eps->ncv!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: 'ncv' is ignored by PRIMME\n"));
285 9 : } else if (eps->ncv!=PETSC_DETERMINE) primme->maxBasisSize = (int)eps->ncv;
286 :
287 10 : PetscCheck(primme_set_method(ops->method,primme)>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"PRIMME method not valid");
288 :
289 10 : eps->mpd = (PetscInt)primme->maxBasisSize;
290 10 : eps->ncv = (PetscInt)(primme->locking?eps->nev:0)+primme->maxBasisSize;
291 10 : ops->bs = (PetscInt)primme->maxBlockSize;
292 :
293 : /* Set workspace */
294 10 : PetscCall(EPSAllocateSolution(eps,0));
295 :
296 : /* Setup the preconditioner */
297 10 : if (primme->correctionParams.precondition) {
298 10 : PetscCall(STGetKSP(eps->st,&ops->ksp));
299 10 : PetscCall(PetscObjectTypeCompare((PetscObject)ops->ksp,KSPPREONLY,&flg));
300 10 : if (!flg) PetscCall(PetscInfo(eps,"Warning: ignoring KSP, should use KSPPREONLY\n"));
301 10 : primme->preconditioner = NULL;
302 10 : primme->applyPreconditioner = applyPreconditioner_PRIMME;
303 : }
304 :
305 : /* Prepare auxiliary vectors */
306 10 : if (!ops->x) PetscCall(MatCreateVecsEmpty(ops->A,&ops->x,&ops->y));
307 10 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 :
310 10 : static PetscErrorCode EPSSolve_PRIMME(EPS eps)
311 : {
312 10 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
313 10 : PetscScalar *a;
314 10 : PetscInt i,ld,ierrprimme;
315 10 : PetscReal *evals,*rnorms;
316 :
317 10 : PetscFunctionBegin;
318 : /* Reset some parameters left from previous runs */
319 : #if defined(SLEPC_HAVE_PRIMME2p2)
320 10 : ops->primme.aNorm = 0.0;
321 : #else
322 : /* Force PRIMME to stop by absolute error */
323 : ops->primme.aNorm = 1.0;
324 : #endif
325 10 : ops->primme.initSize = (int)eps->nini;
326 10 : ops->primme.iseed[0] = -1;
327 10 : ops->primme.iseed[1] = -1;
328 10 : ops->primme.iseed[2] = -1;
329 10 : ops->primme.iseed[3] = -1;
330 10 : PetscCall(BVGetLeadingDimension(eps->V,&ld));
331 10 : ops->primme.ldevecs = (PRIMME_INT)ld;
332 :
333 : /* Call PRIMME solver */
334 10 : PetscCall(BVGetArray(eps->V,&a));
335 10 : PetscCall(PetscMalloc2(eps->ncv,&evals,eps->ncv,&rnorms));
336 10 : ierrprimme = PRIMME_DRIVER(evals,a,rnorms,&ops->primme);
337 267 : for (i=0;i<eps->ncv;i++) eps->eigr[i] = evals[i];
338 257 : for (i=0;i<eps->ncv;i++) eps->errest[i] = rnorms[i];
339 10 : PetscCall(PetscFree2(evals,rnorms));
340 10 : PetscCall(BVRestoreArray(eps->V,&a));
341 :
342 10 : eps->nconv = ops->primme.initSize >= 0 ? (PetscInt)ops->primme.initSize : 0;
343 10 : eps->reason = eps->nconv >= eps->nev ? EPS_CONVERGED_TOL: EPS_DIVERGED_ITS;
344 10 : PetscCall(PetscIntCast(ops->primme.stats.numOuterIterations,&eps->its));
345 :
346 : /* Process PRIMME error code */
347 10 : switch (ierrprimme) {
348 : case 0: /* no error */
349 : break;
350 0 : case -1:
351 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": unexpected error",ierrprimme);
352 0 : case -2:
353 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": allocation error",ierrprimme);
354 : case -3: /* stop due to maximum number of iterations or matvecs */
355 : break;
356 0 : default:
357 0 : PetscCheck(ierrprimme<-39,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": configuration error; check PRIMME's manual",ierrprimme);
358 0 : PetscCheck(ierrprimme>=-39,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": runtime error; check PRIMME's manual",ierrprimme);
359 : }
360 10 : PetscFunctionReturn(PETSC_SUCCESS);
361 : }
362 :
363 9 : static PetscErrorCode EPSReset_PRIMME(EPS eps)
364 : {
365 9 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
366 :
367 9 : PetscFunctionBegin;
368 9 : primme_free(&ops->primme);
369 9 : PetscCall(VecDestroy(&ops->x));
370 9 : PetscCall(VecDestroy(&ops->y));
371 9 : PetscFunctionReturn(PETSC_SUCCESS);
372 : }
373 :
374 9 : static PetscErrorCode EPSDestroy_PRIMME(EPS eps)
375 : {
376 9 : PetscFunctionBegin;
377 9 : PetscCall(PetscFree(eps->data));
378 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMESetBlockSize_C",NULL));
379 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMESetMethod_C",NULL));
380 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMEGetBlockSize_C",NULL));
381 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMEGetMethod_C",NULL));
382 9 : PetscFunctionReturn(PETSC_SUCCESS);
383 : }
384 :
385 0 : static PetscErrorCode EPSView_PRIMME(EPS eps,PetscViewer viewer)
386 : {
387 0 : PetscBool isascii;
388 0 : EPS_PRIMME *ctx = (EPS_PRIMME*)eps->data;
389 0 : PetscMPIInt rank;
390 :
391 0 : PetscFunctionBegin;
392 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
393 0 : if (isascii) {
394 0 : PetscCall(PetscViewerASCIIPrintf(viewer," block size=%" PetscInt_FMT "\n",ctx->bs));
395 0 : PetscCall(PetscViewerASCIIPrintf(viewer," solver method: %s\n",EPSPRIMMEMethods[(EPSPRIMMEMethod)ctx->method]));
396 :
397 : /* Display PRIMME params */
398 0 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)eps),&rank));
399 0 : if (!rank) primme_display_params(ctx->primme);
400 : }
401 0 : PetscFunctionReturn(PETSC_SUCCESS);
402 : }
403 :
404 8 : static PetscErrorCode EPSSetFromOptions_PRIMME(EPS eps,PetscOptionItems *PetscOptionsObject)
405 : {
406 8 : EPS_PRIMME *ctx = (EPS_PRIMME*)eps->data;
407 8 : PetscInt bs;
408 8 : EPSPRIMMEMethod meth;
409 8 : PetscBool flg;
410 :
411 8 : PetscFunctionBegin;
412 8 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS PRIMME Options");
413 :
414 8 : PetscCall(PetscOptionsInt("-eps_primme_blocksize","Maximum block size","EPSPRIMMESetBlockSize",ctx->bs,&bs,&flg));
415 8 : if (flg) PetscCall(EPSPRIMMESetBlockSize(eps,bs));
416 :
417 8 : PetscCall(PetscOptionsEnum("-eps_primme_method","Method for solving the eigenproblem","EPSPRIMMESetMethod",EPSPRIMMEMethods,(PetscEnum)ctx->method,(PetscEnum*)&meth,&flg));
418 8 : if (flg) PetscCall(EPSPRIMMESetMethod(eps,meth));
419 :
420 8 : PetscOptionsHeadEnd();
421 8 : PetscFunctionReturn(PETSC_SUCCESS);
422 : }
423 :
424 5 : static PetscErrorCode EPSPRIMMESetBlockSize_PRIMME(EPS eps,PetscInt bs)
425 : {
426 5 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
427 :
428 5 : PetscFunctionBegin;
429 5 : if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) ops->bs = 0;
430 : else {
431 5 : PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"PRIMME: block size must be positive");
432 5 : ops->bs = bs;
433 : }
434 5 : PetscFunctionReturn(PETSC_SUCCESS);
435 : }
436 :
437 : /*@
438 : EPSPRIMMESetBlockSize - The maximum block size that PRIMME will try to use.
439 :
440 : Logically Collective
441 :
442 : Input Parameters:
443 : + eps - the eigenproblem solver context
444 : - bs - block size
445 :
446 : Options Database Key:
447 : . -eps_primme_blocksize - Sets the max allowed block size value
448 :
449 : Notes:
450 : If the block size is not set, the value established by primme_initialize
451 : is used.
452 :
453 : The user should set the block size based on the architecture specifics
454 : of the target computer, as well as any a priori knowledge of multiplicities.
455 : The code does NOT require bs > 1 to find multiple eigenvalues. For some
456 : methods, keeping bs = 1 yields the best overall performance.
457 :
458 : Level: advanced
459 :
460 : .seealso: EPSPRIMMEGetBlockSize()
461 : @*/
462 5 : PetscErrorCode EPSPRIMMESetBlockSize(EPS eps,PetscInt bs)
463 : {
464 5 : PetscFunctionBegin;
465 5 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
466 15 : PetscValidLogicalCollectiveInt(eps,bs,2);
467 5 : PetscTryMethod(eps,"EPSPRIMMESetBlockSize_C",(EPS,PetscInt),(eps,bs));
468 5 : PetscFunctionReturn(PETSC_SUCCESS);
469 : }
470 :
471 4 : static PetscErrorCode EPSPRIMMEGetBlockSize_PRIMME(EPS eps,PetscInt *bs)
472 : {
473 4 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
474 :
475 4 : PetscFunctionBegin;
476 4 : *bs = ops->bs;
477 4 : PetscFunctionReturn(PETSC_SUCCESS);
478 : }
479 :
480 : /*@
481 : EPSPRIMMEGetBlockSize - Get the maximum block size the code will try to use.
482 :
483 : Not Collective
484 :
485 : Input Parameter:
486 : . eps - the eigenproblem solver context
487 :
488 : Output Parameter:
489 : . bs - returned block size
490 :
491 : Level: advanced
492 :
493 : .seealso: EPSPRIMMESetBlockSize()
494 : @*/
495 4 : PetscErrorCode EPSPRIMMEGetBlockSize(EPS eps,PetscInt *bs)
496 : {
497 4 : PetscFunctionBegin;
498 4 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
499 4 : PetscAssertPointer(bs,2);
500 4 : PetscUseMethod(eps,"EPSPRIMMEGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
501 4 : PetscFunctionReturn(PETSC_SUCCESS);
502 : }
503 :
504 5 : static PetscErrorCode EPSPRIMMESetMethod_PRIMME(EPS eps,EPSPRIMMEMethod method)
505 : {
506 5 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
507 :
508 5 : PetscFunctionBegin;
509 5 : ops->method = (primme_preset_method)method;
510 5 : PetscFunctionReturn(PETSC_SUCCESS);
511 : }
512 :
513 : /*@
514 : EPSPRIMMESetMethod - Sets the method for the PRIMME library.
515 :
516 : Logically Collective
517 :
518 : Input Parameters:
519 : + eps - the eigenproblem solver context
520 : - method - method that will be used by PRIMME
521 :
522 : Options Database Key:
523 : . -eps_primme_method - Sets the method for the PRIMME library
524 :
525 : Note:
526 : If not set, the method defaults to EPS_PRIMME_DEFAULT_MIN_TIME.
527 :
528 : Level: advanced
529 :
530 : .seealso: EPSPRIMMEGetMethod(), EPSPRIMMEMethod
531 : @*/
532 5 : PetscErrorCode EPSPRIMMESetMethod(EPS eps,EPSPRIMMEMethod method)
533 : {
534 5 : PetscFunctionBegin;
535 5 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
536 15 : PetscValidLogicalCollectiveEnum(eps,method,2);
537 5 : PetscTryMethod(eps,"EPSPRIMMESetMethod_C",(EPS,EPSPRIMMEMethod),(eps,method));
538 5 : PetscFunctionReturn(PETSC_SUCCESS);
539 : }
540 :
541 4 : static PetscErrorCode EPSPRIMMEGetMethod_PRIMME(EPS eps,EPSPRIMMEMethod *method)
542 : {
543 4 : EPS_PRIMME *ops = (EPS_PRIMME*)eps->data;
544 :
545 4 : PetscFunctionBegin;
546 4 : *method = (EPSPRIMMEMethod)ops->method;
547 4 : PetscFunctionReturn(PETSC_SUCCESS);
548 : }
549 :
550 : /*@
551 : EPSPRIMMEGetMethod - Gets the method for the PRIMME library.
552 :
553 : Not Collective
554 :
555 : Input Parameter:
556 : . eps - the eigenproblem solver context
557 :
558 : Output Parameter:
559 : . method - method that will be used by PRIMME
560 :
561 : Level: advanced
562 :
563 : .seealso: EPSPRIMMESetMethod(), EPSPRIMMEMethod
564 : @*/
565 4 : PetscErrorCode EPSPRIMMEGetMethod(EPS eps,EPSPRIMMEMethod *method)
566 : {
567 4 : PetscFunctionBegin;
568 4 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
569 4 : PetscAssertPointer(method,2);
570 4 : PetscUseMethod(eps,"EPSPRIMMEGetMethod_C",(EPS,EPSPRIMMEMethod*),(eps,method));
571 4 : PetscFunctionReturn(PETSC_SUCCESS);
572 : }
573 :
574 9 : SLEPC_EXTERN PetscErrorCode EPSCreate_PRIMME(EPS eps)
575 : {
576 9 : EPS_PRIMME *primme;
577 :
578 9 : PetscFunctionBegin;
579 9 : PetscCall(PetscNew(&primme));
580 9 : eps->data = (void*)primme;
581 :
582 9 : primme_initialize(&primme->primme);
583 9 : primme->primme.globalSumReal = par_GlobalSumReal;
584 : #if defined(SLEPC_HAVE_PRIMME3)
585 9 : primme->primme.broadcastReal = par_broadcastReal;
586 : #endif
587 : #if defined(SLEPC_HAVE_PRIMME2p2)
588 9 : primme->primme.convTestFun = convTestFun;
589 9 : primme->primme.monitorFun = monitorFun;
590 : #endif
591 9 : primme->method = (primme_preset_method)EPS_PRIMME_DEFAULT_MIN_TIME;
592 :
593 9 : eps->categ = EPS_CATEGORY_PRECOND;
594 :
595 9 : eps->ops->solve = EPSSolve_PRIMME;
596 9 : eps->ops->setup = EPSSetUp_PRIMME;
597 9 : eps->ops->setupsort = EPSSetUpSort_Basic;
598 9 : eps->ops->setfromoptions = EPSSetFromOptions_PRIMME;
599 9 : eps->ops->destroy = EPSDestroy_PRIMME;
600 9 : eps->ops->reset = EPSReset_PRIMME;
601 9 : eps->ops->view = EPSView_PRIMME;
602 9 : eps->ops->backtransform = EPSBackTransform_Default;
603 9 : eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
604 :
605 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMESetBlockSize_C",EPSPRIMMESetBlockSize_PRIMME));
606 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMESetMethod_C",EPSPRIMMESetMethod_PRIMME));
607 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMEGetBlockSize_C",EPSPRIMMEGetBlockSize_PRIMME));
608 9 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSPRIMMEGetMethod_C",EPSPRIMMEGetMethod_PRIMME));
609 9 : PetscFunctionReturn(PETSC_SUCCESS);
610 : }
|
