Actual source code: svdprimme.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 PRIMME SVD solver
12: */
14: #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
16: #include <primme.h>
18: #if defined(PETSC_USE_COMPLEX)
19: #if defined(PETSC_USE_REAL_SINGLE)
20: #define PRIMME_DRIVER cprimme_svds
21: #else
22: #define PRIMME_DRIVER zprimme_svds
23: #endif
24: #else
25: #if defined(PETSC_USE_REAL_SINGLE)
26: #define PRIMME_DRIVER sprimme_svds
27: #else
28: #define PRIMME_DRIVER dprimme_svds
29: #endif
30: #endif
32: typedef struct {
33: primme_svds_params primme; /* param struct */
34: primme_svds_preset_method method; /* primme method */
35: SVD svd; /* reference to the solver */
36: Vec x,y; /* auxiliary vectors */
37: } SVD_PRIMME;
39: static void multMatvec_PRIMME(void*,PRIMME_INT*,void*,PRIMME_INT*,int*,int*,struct primme_svds_params*,int*);
41: static void par_GlobalSumReal(void *sendBuf,void *recvBuf,int *count,primme_svds_params *primme,int *ierr)
42: {
43: *MPI_Allreduce(sendBuf,recvBuf,*count,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)primme->commInfo));CHKERRABORT(PetscObjectComm((PetscObject)primme->commInfo),*ierr);
44: }
46: PetscErrorCode SVDSetUp_PRIMME(SVD svd)
47: {
48: PetscErrorCode ierr;
49: PetscMPIInt numProcs,procID;
50: PetscInt n,m,nloc,mloc;
51: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
52: primme_svds_params *primme = &ops->primme;
55: MPI_Comm_size(PetscObjectComm((PetscObject)svd),&numProcs);
56: MPI_Comm_rank(PetscObjectComm((PetscObject)svd),&procID);
58: /* Check some constraints and set some default values */
59: SVDMatGetSize(svd,&m,&n);
60: SVDMatGetLocalSize(svd,&mloc,&nloc);
61: SVDSetDimensions_Default(svd);
62: if (!svd->max_it) svd->max_it = PetscMax(n/svd->ncv,1000);
63: svd->leftbasis = PETSC_TRUE;
64: if (svd->stopping!=SVDStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"External packages do not support user-defined stopping test");
65: if (svd->converged != SVDConvergedAbsolute) { PetscInfo(svd,"Warning: using absolute convergence test\n"); }
67: /* Transfer SLEPc options to PRIMME options */
68: primme->m = m;
69: primme->n = n;
70: primme->mLocal = mloc;
71: primme->nLocal = nloc;
72: primme->numSvals = svd->nsv;
73: primme->matrix = ops;
74: primme->commInfo = svd;
75: primme->maxMatvecs = svd->max_it;
76: primme->eps = svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL:svd->tol;
77: primme->numProcs = numProcs;
78: primme->procID = procID;
79: primme->locking = 1;
80: primme->printLevel = 0;
82: switch (svd->which) {
83: case SVD_LARGEST:
84: primme->target = primme_svds_largest;
85: break;
86: case SVD_SMALLEST:
87: primme->target = primme_svds_smallest;
88: break;
89: default:
90: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"'which' value not supported by PRIMME");
91: break;
92: }
94: /* If user sets mpd or ncv, maxBasisSize is modified */
95: if (svd->mpd) primme->maxBasisSize = svd->mpd;
96: else if (svd->ncv) primme->maxBasisSize = svd->ncv;
98: if (primme_svds_set_method(ops->method,(primme_preset_method)EPS_PRIMME_DEFAULT_MIN_TIME,(primme_preset_method)EPS_PRIMME_DEFAULT_MIN_TIME,primme) < 0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"PRIMME method not valid");
100: svd->mpd = primme->maxBasisSize;
101: svd->ncv = svd->nsv;
103: /* Set workspace */
104: SVDAllocateSolution(svd,0);
106: /* Prepare auxiliary vectors */
107: if (!ops->x) {
108: MatCreateVecsEmpty(svd->A,&ops->x,&ops->y);
109: PetscLogObjectParent((PetscObject)svd,(PetscObject)ops->x);
110: PetscLogObjectParent((PetscObject)svd,(PetscObject)ops->y);
111: }
112: return(0);
113: }
115: PetscErrorCode SVDSolve_PRIMME(SVD svd)
116: {
118: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
119: PetscScalar *svecs, *a;
122: /* Reset some parameters left from previous runs */
123: ops->primme.aNorm = 0.0;
124: ops->primme.initSize = svd->nini;
125: ops->primme.iseed[0] = -1;
126: ops->primme.iseed[1] = -1;
127: ops->primme.iseed[2] = -1;
128: ops->primme.iseed[3] = -1;
130: /* Allocating left and right singular vectors contiguously */
131: PetscMalloc1(ops->primme.numSvals*(ops->primme.mLocal+ops->primme.nLocal),&svecs);
132: PetscLogObjectMemory((PetscObject)svd,sizeof(PetscReal)*ops->primme.numSvals*(ops->primme.mLocal+ops->primme.nLocal));
133:
134: /* Call PRIMME solver */
135: PRIMME_DRIVER(svd->sigma,svecs,svd->errest,&ops->primme);
136: if (ierr) SETERRQ1(PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"PRIMME library failed with error code=%d",ierr);
138: /* Copy left and right singular vectors into svd */
139: BVGetArray(svd->U,&a);
140: PetscMemcpy(a,svecs,sizeof(PetscScalar)*ops->primme.mLocal*ops->primme.initSize);
141: BVRestoreArray(svd->U,&a);
143: BVGetArray(svd->V,&a);
144: PetscMemcpy(a,svecs+ops->primme.mLocal*ops->primme.initSize,sizeof(PetscScalar)*ops->primme.nLocal*ops->primme.initSize);
145: BVRestoreArray(svd->V,&a);
147: PetscFree(svecs);
149: svd->nconv = ops->primme.initSize >= 0 ? ops->primme.initSize : 0;
150: svd->reason = svd->nconv >= svd->nsv ? SVD_CONVERGED_TOL: SVD_DIVERGED_ITS;
151: svd->its = ops->primme.stats.numOuterIterations;
152: return(0);
153: }
155: static void multMatvec_PRIMME(void *xa,PRIMME_INT *ldx,void *ya,PRIMME_INT *ldy,int *blockSize,int *transpose,struct primme_svds_params *primme,int *ierr)
156: {
157: PetscInt i;
158: SVD_PRIMME *ops = (SVD_PRIMME*)primme->matrix;
159: Vec x = ops->x,y = ops->y;
160: SVD svd = ops->svd;
163: for (i=0;i<*blockSize;i++) {
164: if (*transpose) {
165: *VecPlaceArray(y,(PetscScalar*)xa+(*ldx)*i);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
166: *VecPlaceArray(x,(PetscScalar*)ya+(*ldy)*i);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
167: *SVDMatMult(svd,PETSC_TRUE,y,x);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
168: } else {
169: *VecPlaceArray(x,(PetscScalar*)xa+(*ldx)*i);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
170: *VecPlaceArray(y,(PetscScalar*)ya+(*ldy)*i);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
171: *SVDMatMult(svd,PETSC_FALSE,x,y);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
172: }
173: *VecResetArray(x);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
174: *VecResetArray(y);CHKERRABORT(PetscObjectComm((PetscObject)svd),*ierr);
175: }
176: PetscFunctionReturnVoid();
177: }
179: PetscErrorCode SVDReset_PRIMME(SVD svd)
180: {
182: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
185: primme_svds_free(&ops->primme);
186: VecDestroy(&ops->x);
187: VecDestroy(&ops->y);
188: return(0);
189: }
191: PetscErrorCode SVDDestroy_PRIMME(SVD svd)
192: {
196: PetscFree(svd->data);
197: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetBlockSize_C",NULL);
198: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetBlockSize_C",NULL);
199: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetMethod_C",NULL);
200: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetMethod_C",NULL);
201: return(0);
202: }
204: PetscErrorCode SVDView_PRIMME(SVD svd,PetscViewer viewer)
205: {
206: PetscErrorCode ierr;
207: PetscBool isascii;
208: primme_svds_params *primme = &((SVD_PRIMME*)svd->data)->primme;
209: SVDPRIMMEMethod methodn;
210: PetscMPIInt rank;
213: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
214: if (isascii) {
215: PetscViewerASCIIPrintf(viewer," block size=%D\n",primme->maxBlockSize);
216: SVDPRIMMEGetMethod(svd,&methodn);
217: PetscViewerASCIIPrintf(viewer," solver method: %s\n",SVDPRIMMEMethods[methodn]);
219: /* Display PRIMME params */
220: MPI_Comm_rank(PetscObjectComm((PetscObject)svd),&rank);
221: if (!rank) primme_svds_display_params(*primme);
222: }
223: return(0);
224: }
226: PetscErrorCode SVDSetFromOptions_PRIMME(PetscOptionItems *PetscOptionsObject,SVD svd)
227: {
228: PetscErrorCode ierr;
229: SVD_PRIMME *ctx = (SVD_PRIMME*)svd->data;
230: PetscInt bs;
231: SVDPRIMMEMethod meth;
232: PetscBool flg;
235: PetscOptionsHead(PetscOptionsObject,"SVD PRIMME Options");
237: PetscOptionsInt("-svd_primme_blocksize","Maximum block size","SVDPRIMMESetBlockSize",ctx->primme.maxBlockSize,&bs,&flg);
238: if (flg) { SVDPRIMMESetBlockSize(svd,bs); }
240: PetscOptionsEnum("-svd_primme_method","Method for solving the singular value problem","SVDPRIMMESetMethod",SVDPRIMMEMethods,(PetscEnum)ctx->method,(PetscEnum*)&meth,&flg);
241: if (flg) { SVDPRIMMESetMethod(svd,meth); }
243: PetscOptionsTail();
244: return(0);
245: }
247: static PetscErrorCode SVDPRIMMESetBlockSize_PRIMME(SVD svd,PetscInt bs)
248: {
249: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
252: if (bs == PETSC_DEFAULT) ops->primme.maxBlockSize = 1;
253: else if (bs <= 0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"PRIMME: block size must be positive");
254: else ops->primme.maxBlockSize = bs;
255: return(0);
256: }
258: /*@
259: SVDPRIMMESetBlockSize - The maximum block size that PRIMME will try to use.
261: Logically Collective on SVD
263: Input Parameters:
264: + svd - the singular value solver context
265: - bs - block size
267: Options Database Key:
268: . -svd_primme_blocksize - Sets the max allowed block size value
270: Notes:
271: If the block size is not set, the value established by primme_svds_initialize
272: is used.
274: The user should set the block size based on the architecture specifics
275: of the target computer, as well as any a priori knowledge of multiplicities.
276: The code does NOT require bs > 1 to find multiple eigenvalues. For some
277: methods, keeping bs = 1 yields the best overall performance.
279: Level: advanced
281: .seealso: SVDPRIMMEGetBlockSize()
282: @*/
283: PetscErrorCode SVDPRIMMESetBlockSize(SVD svd,PetscInt bs)
284: {
290: PetscTryMethod(svd,"SVDPRIMMESetBlockSize_C",(SVD,PetscInt),(svd,bs));
291: return(0);
292: }
294: static PetscErrorCode SVDPRIMMEGetBlockSize_PRIMME(SVD svd,PetscInt *bs)
295: {
296: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
299: *bs = ops->primme.maxBlockSize;
300: return(0);
301: }
303: /*@
304: SVDPRIMMEGetBlockSize - Get the maximum block size the code will try to use.
306: Not Collective
308: Input Parameter:
309: . svd - the singular value solver context
311: Output Parameter:
312: . bs - returned block size
314: Level: advanced
316: .seealso: SVDPRIMMESetBlockSize()
317: @*/
318: PetscErrorCode SVDPRIMMEGetBlockSize(SVD svd,PetscInt *bs)
319: {
325: PetscUseMethod(svd,"SVDPRIMMEGetBlockSize_C",(SVD,PetscInt*),(svd,bs));
326: return(0);
327: }
329: static PetscErrorCode SVDPRIMMESetMethod_PRIMME(SVD svd,SVDPRIMMEMethod method)
330: {
331: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
334: ops->method = (primme_svds_preset_method)method;
335: return(0);
336: }
338: /*@
339: SVDPRIMMESetMethod - Sets the method for the PRIMME SVD solver.
341: Logically Collective on SVD
343: Input Parameters:
344: + svd - the singular value solver context
345: - method - method that will be used by PRIMME
347: Options Database Key:
348: . -svd_primme_method - Sets the method for the PRIMME SVD solver
350: Note:
351: If not set, the method defaults to SVD_PRIMME_HYBRID.
353: Level: advanced
355: .seealso: SVDPRIMMEGetMethod(), SVDPRIMMEMethod
356: @*/
357: PetscErrorCode SVDPRIMMESetMethod(SVD svd,SVDPRIMMEMethod method)
358: {
364: PetscTryMethod(svd,"SVDPRIMMESetMethod_C",(SVD,SVDPRIMMEMethod),(svd,method));
365: return(0);
366: }
368: static PetscErrorCode SVDPRIMMEGetMethod_PRIMME(SVD svd,SVDPRIMMEMethod *method)
369: {
370: SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
373: *method = (SVDPRIMMEMethod)ops->method;
374: return(0);
375: }
377: /*@
378: SVDPRIMMEGetMethod - Gets the method for the PRIMME SVD solver.
380: Not Collective
382: Input Parameter:
383: . svd - the singular value solver context
385: Output Parameter:
386: . method - method that will be used by PRIMME
388: Level: advanced
390: .seealso: SVDPRIMMESetMethod(), SVDPRIMMEMethod
391: @*/
392: PetscErrorCode SVDPRIMMEGetMethod(SVD svd,SVDPRIMMEMethod *method)
393: {
399: PetscUseMethod(svd,"SVDPRIMMEGetMethod_C",(SVD,SVDPRIMMEMethod*),(svd,method));
400: return(0);
401: }
403: SLEPC_EXTERN PetscErrorCode SVDCreate_PRIMME(SVD svd)
404: {
406: SVD_PRIMME *primme;
409: PetscNewLog(svd,&primme);
410: svd->data = (void*)primme;
412: primme_svds_initialize(&primme->primme);
413: primme->primme.matrixMatvec = multMatvec_PRIMME;
414: primme->primme.globalSumReal = par_GlobalSumReal;
415: primme->method = (primme_svds_preset_method)SVD_PRIMME_HYBRID;
416: primme->svd = svd;
418: svd->ops->solve = SVDSolve_PRIMME;
419: svd->ops->setup = SVDSetUp_PRIMME;
420: svd->ops->setfromoptions = SVDSetFromOptions_PRIMME;
421: svd->ops->destroy = SVDDestroy_PRIMME;
422: svd->ops->reset = SVDReset_PRIMME;
423: svd->ops->view = SVDView_PRIMME;
425: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetBlockSize_C",SVDPRIMMESetBlockSize_PRIMME);
426: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetBlockSize_C",SVDPRIMMEGetBlockSize_PRIMME);
427: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetMethod_C",SVDPRIMMESetMethod_PRIMME);
428: PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetMethod_C",SVDPRIMMEGetMethod_PRIMME);
429: return(0);
430: }