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 SVD solver
      12             : */
      13             : 
      14             : #include <slepc/private/svdimpl.h>    /*I "slepcsvd.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_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
      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_svds_params        primme;   /* param struct */
      38             :   PetscInt                  bs;       /* block size */
      39             :   primme_svds_preset_method method;   /* primme method */
      40             :   SVD                       svd;      /* reference to the solver */
      41             :   Vec                       x,y;      /* auxiliary vectors */
      42             : } SVD_PRIMME;
      43             : 
      44             : static void multMatvec_PRIMME(void*,PRIMME_INT*,void*,PRIMME_INT*,int*,int*,struct primme_svds_params*,int*);
      45             : 
      46        2313 : static void par_GlobalSumReal(void *sendBuf,void *recvBuf,int *count,primme_svds_params *primme,int *ierr)
      47             : {
      48        2313 :   if (sendBuf == recvBuf) {
      49        4626 :     *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        2313 : }
      54             : 
      55             : #if defined(SLEPC_HAVE_PRIMME3)
      56        1238 : static void par_broadcastReal(void *buf,int *count,primme_svds_params *primme,int *ierr)
      57             : {
      58        1238 :   *ierr = MPI_Bcast(buf,*count,MPIU_REAL,0/*root*/,PetscObjectComm((PetscObject)primme->commInfo));
      59        1238 : }
      60             : #endif
      61             : 
      62             : #if defined(SLEPC_HAVE_PRIMME2p2)
      63         964 : static void convTestFun(double *sval,void *leftsvec,void *rightsvec,double *resNorm,
      64             : #if defined(SLEPC_HAVE_PRIMME3)
      65             :                         int *method,
      66             : #endif
      67             :                         int *isconv,struct primme_svds_params *primme,int *err)
      68             : {
      69         964 :   PetscErrorCode ierr;
      70         964 :   SVD            svd = (SVD)primme->commInfo;
      71         964 :   PetscReal      sigma = sval?*sval:0.0;
      72         964 :   PetscReal      r = resNorm?*resNorm:HUGE_VAL,errest;
      73             : 
      74         964 :   ierr = (*svd->converged)(svd,sigma,r,&errest,svd->convergedctx);
      75         964 :   if (ierr) *err = 1;
      76             :   else {
      77         964 :     *isconv = (errest<=svd->tol?1:0);
      78         964 :     *err = 0;
      79             :   }
      80         964 : }
      81             : 
      82         783 : static void monitorFun(void *basisSvals,int *basisSize,int *basisFlags,int *iblock,int *blockSize,void *basisNorms,int *numConverged,void *lockedSvals,int *numLocked,int *lockedFlags,void *lockedNorms,int *inner_its,void *LSRes,
      83             : #if defined(SLEPC_HAVE_PRIMME3)
      84             :                        const char *msg,double *time,
      85             : #endif
      86             :                        primme_event *event,int *stage,struct primme_svds_params *primme,int *err)
      87             : {
      88             : 
      89         783 :   PetscErrorCode ierr = PETSC_SUCCESS;
      90         783 :   SVD            svd = (SVD)primme->commInfo;
      91         783 :   PetscInt       i,k,nerrest;
      92             : 
      93         783 :   *err = 1;
      94         783 :   switch (*event) {
      95         143 :     case primme_event_outer_iteration:
      96             :       /* Update SVD */
      97         143 :       PetscCallVoid(PetscIntCast(primme->stats.numOuterIterations,&svd->its));
      98         143 :       if (numConverged) svd->nconv = *numConverged;
      99         143 :       k = 0;
     100         143 :       if (lockedSvals && numLocked) for (i=0; i<*numLocked && k<svd->ncv; i++) svd->sigma[k++] = ((PetscReal*)lockedSvals)[i];
     101         143 :       nerrest = k;
     102         143 :       if (iblock && blockSize) {
     103         388 :         for (i=0; i<*blockSize && k+iblock[i]<svd->ncv; i++) svd->errest[k+iblock[i]] = ((PetscReal*)basisNorms)[i];
     104         143 :         nerrest = k+(*blockSize>0?1+iblock[*blockSize-1]:0);
     105             :       }
     106        1065 :       if (basisSvals && basisSize) for (i=0; i<*basisSize && k<svd->ncv; i++) svd->sigma[k++] = ((PetscReal*)basisSvals)[i];
     107             :       /* Show progress */
     108         143 :       ierr = SVDMonitor(svd,svd->its,numConverged?*numConverged:0,svd->sigma,svd->errest,nerrest);
     109         143 :       break;
     110             : #if defined(SLEPC_HAVE_PRIMME3)
     111           0 :     case primme_event_message:
     112             :       /* Print PRIMME information messages */
     113           0 :       ierr = PetscInfo(svd,"%s\n",msg);
     114           0 :       break;
     115             : #endif
     116             :     default:
     117             :       break;
     118             :   }
     119         783 :   *err = (ierr!=0)? 1: 0;
     120             : }
     121             : #endif /* SLEPC_HAVE_PRIMME2p2 */
     122             : 
     123        1437 : 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)
     124             : {
     125        1437 :   PetscInt   i;
     126        1437 :   SVD_PRIMME *ops = (SVD_PRIMME*)primme->matrix;
     127        1437 :   Vec        x = ops->x,y = ops->y;
     128        1437 :   SVD        svd = ops->svd;
     129             : 
     130        1437 :   PetscFunctionBegin;
     131        3649 :   for (i=0;i<*blockSize;i++) {
     132        2212 :     if (*transpose) {
     133        1092 :       PetscCallAbort(PetscObjectComm((PetscObject)svd),VecPlaceArray(y,(PetscScalar*)xa+(*ldx)*i));
     134        1092 :       PetscCallAbort(PetscObjectComm((PetscObject)svd),VecPlaceArray(x,(PetscScalar*)ya+(*ldy)*i));
     135        1092 :       PetscCallAbort(PetscObjectComm((PetscObject)svd),MatMult(svd->AT,y,x));
     136             :     } else {
     137        1120 :       PetscCallAbort(PetscObjectComm((PetscObject)svd),VecPlaceArray(x,(PetscScalar*)xa+(*ldx)*i));
     138        1120 :       PetscCallAbort(PetscObjectComm((PetscObject)svd),VecPlaceArray(y,(PetscScalar*)ya+(*ldy)*i));
     139        1120 :       PetscCallAbort(PetscObjectComm((PetscObject)svd),MatMult(svd->A,x,y));
     140             :     }
     141        2212 :     PetscCallAbort(PetscObjectComm((PetscObject)svd),VecResetArray(x));
     142        2212 :     PetscCallAbort(PetscObjectComm((PetscObject)svd),VecResetArray(y));
     143             :   }
     144        1437 :   PetscFunctionReturnVoid();
     145             : }
     146             : 
     147           7 : static PetscErrorCode SVDSetUp_PRIMME(SVD svd)
     148             : {
     149           7 :   PetscMPIInt        numProcs,procID;
     150           7 :   PetscInt           n,m,nloc,mloc;
     151           7 :   SVD_PRIMME         *ops = (SVD_PRIMME*)svd->data;
     152           7 :   primme_svds_params *primme = &ops->primme;
     153             : 
     154           7 :   PetscFunctionBegin;
     155           7 :   SVDCheckStandard(svd);
     156           7 :   SVDCheckDefinite(svd);
     157           7 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)svd),&numProcs));
     158           7 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)svd),&procID));
     159             : 
     160             :   /* Check some constraints and set some default values */
     161           7 :   PetscCall(MatGetSize(svd->A,&m,&n));
     162           7 :   PetscCall(MatGetLocalSize(svd->A,&mloc,&nloc));
     163           7 :   PetscCall(SVDSetDimensions_Default(svd));
     164           7 :   if (svd->max_it==PETSC_DETERMINE) svd->max_it = PETSC_INT_MAX;
     165           7 :   svd->leftbasis = PETSC_TRUE;
     166           7 :   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
     167             : #if !defined(SLEPC_HAVE_PRIMME2p2)
     168             :   if (svd->converged != SVDConvergedAbsolute) PetscCall(PetscInfo(svd,"Warning: using absolute convergence test\n"));
     169             : #endif
     170             : 
     171             :   /* Transfer SLEPc options to PRIMME options */
     172           7 :   primme_svds_free(primme);
     173           7 :   primme_svds_initialize(primme);
     174           7 :   primme->m             = (PRIMME_INT)m;
     175           7 :   primme->n             = (PRIMME_INT)n;
     176           7 :   primme->mLocal        = (PRIMME_INT)mloc;
     177           7 :   primme->nLocal        = (PRIMME_INT)nloc;
     178           7 :   primme->numSvals      = (int)svd->nsv;
     179           7 :   primme->matrix        = ops;
     180           7 :   primme->commInfo      = svd;
     181           7 :   primme->maxMatvecs    = (PRIMME_INT)svd->max_it;
     182             : #if !defined(SLEPC_HAVE_PRIMME2p2)
     183             :   primme->eps           = SlepcDefaultTol(svd->tol);
     184             : #endif
     185           7 :   primme->numProcs      = numProcs;
     186           7 :   primme->procID        = procID;
     187           7 :   primme->printLevel    = 1;
     188           7 :   primme->matrixMatvec  = multMatvec_PRIMME;
     189           7 :   primme->globalSumReal = par_GlobalSumReal;
     190             : #if defined(SLEPC_HAVE_PRIMME3)
     191           7 :   primme->broadcastReal = par_broadcastReal;
     192             : #endif
     193             : #if defined(SLEPC_HAVE_PRIMME2p2)
     194           7 :   primme->convTestFun   = convTestFun;
     195           7 :   primme->monitorFun    = monitorFun;
     196             : #endif
     197           7 :   if (ops->bs > 0) primme->maxBlockSize = (int)ops->bs;
     198             : 
     199           7 :   switch (svd->which) {
     200           6 :     case SVD_LARGEST:
     201           6 :       primme->target = primme_svds_largest;
     202           6 :       break;
     203           1 :     case SVD_SMALLEST:
     204           1 :       primme->target = primme_svds_smallest;
     205           1 :       break;
     206             :   }
     207             : 
     208             :   /* If user sets mpd or ncv, maxBasisSize is modified */
     209           7 :   if (svd->mpd!=PETSC_DETERMINE) {
     210           7 :     primme->maxBasisSize = (int)svd->mpd;
     211           7 :     if (svd->ncv!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: 'ncv' is ignored by PRIMME\n"));
     212           0 :   } else if (svd->ncv!=PETSC_DETERMINE) primme->maxBasisSize = (int)svd->ncv;
     213             : 
     214           7 :   PetscCheck(primme_svds_set_method(ops->method,(primme_preset_method)EPS_PRIMME_DEFAULT_MIN_TIME,PRIMME_DEFAULT_METHOD,primme)>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"PRIMME method not valid");
     215             : 
     216           7 :   svd->mpd = (PetscInt)primme->maxBasisSize;
     217           7 :   svd->ncv = (PetscInt)(primme->locking?svd->nsv:0)+primme->maxBasisSize;
     218           7 :   ops->bs  = (PetscInt)primme->maxBlockSize;
     219             : 
     220             :   /* Set workspace */
     221           7 :   PetscCall(SVDAllocateSolution(svd,0));
     222             : 
     223             :   /* Prepare auxiliary vectors */
     224           7 :   if (!ops->x) PetscCall(MatCreateVecsEmpty(svd->A,&ops->x,&ops->y));
     225           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     226             : }
     227             : 
     228           7 : static PetscErrorCode SVDSolve_PRIMME(SVD svd)
     229             : {
     230           7 :   SVD_PRIMME     *ops = (SVD_PRIMME*)svd->data;
     231           7 :   PetscScalar    *svecs, *a;
     232           7 :   PetscInt       i,ierrprimme,ld;
     233           7 :   PetscReal      *svals,*rnorms;
     234             : 
     235           7 :   PetscFunctionBegin;
     236             :   /* Reset some parameters left from previous runs */
     237           7 :   ops->primme.aNorm    = 0.0;
     238           7 :   ops->primme.initSize = (int)svd->nini;
     239           7 :   ops->primme.iseed[0] = -1;
     240           7 :   ops->primme.iseed[1] = -1;
     241           7 :   ops->primme.iseed[2] = -1;
     242           7 :   ops->primme.iseed[3] = -1;
     243             : 
     244             :   /* Allocating left and right singular vectors contiguously */
     245           7 :   PetscCall(PetscCalloc1(ops->primme.numSvals*(ops->primme.mLocal+ops->primme.nLocal),&svecs));
     246             : 
     247             :   /* Call PRIMME solver */
     248           7 :   PetscCall(PetscMalloc2(svd->ncv,&svals,svd->ncv,&rnorms));
     249           7 :   ierrprimme = PRIMME_DRIVER(svals,svecs,rnorms,&ops->primme);
     250         114 :   for (i=0;i<svd->ncv;i++) svd->sigma[i] = svals[i];
     251         107 :   for (i=0;i<svd->ncv;i++) svd->errest[i] = rnorms[i];
     252           7 :   PetscCall(PetscFree2(svals,rnorms));
     253             : 
     254             :   /* Copy left and right singular vectors into svd */
     255           7 :   PetscCall(BVGetLeadingDimension(svd->U,&ld));
     256           7 :   PetscCall(BVGetArray(svd->U,&a));
     257          35 :   for (i=0;i<ops->primme.initSize;i++) PetscCall(PetscArraycpy(a+i*ld,svecs+i*ops->primme.mLocal,ops->primme.mLocal));
     258           7 :   PetscCall(BVRestoreArray(svd->U,&a));
     259             : 
     260           7 :   PetscCall(BVGetLeadingDimension(svd->V,&ld));
     261           7 :   PetscCall(BVGetArray(svd->V,&a));
     262          35 :   for (i=0;i<ops->primme.initSize;i++) PetscCall(PetscArraycpy(a+i*ld,svecs+ops->primme.mLocal*ops->primme.initSize+i*ops->primme.nLocal,ops->primme.nLocal));
     263           7 :   PetscCall(BVRestoreArray(svd->V,&a));
     264             : 
     265           7 :   PetscCall(PetscFree(svecs));
     266             : 
     267           7 :   svd->nconv  = ops->primme.initSize >= 0 ? (PetscInt)ops->primme.initSize : 0;
     268           7 :   svd->reason = svd->nconv >= svd->nsv ? SVD_CONVERGED_TOL: SVD_DIVERGED_ITS;
     269           7 :   PetscCall(PetscIntCast(ops->primme.stats.numOuterIterations,&svd->its));
     270             : 
     271             :   /* Process PRIMME error code */
     272           7 :   if (ierrprimme != 0) {
     273           0 :     switch (ierrprimme%100) {
     274           0 :       case -1:
     275           0 :         SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": unexpected error",ierrprimme);
     276           0 :       case -2:
     277           0 :         SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": allocation error",ierrprimme);
     278             :       case -3: /* stop due to maximum number of iterations or matvecs */
     279             :         break;
     280           0 :       default:
     281           0 :         PetscCheck(ierrprimme<-39,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": configuration error; check PRIMME's manual",ierrprimme);
     282           0 :         PetscCheck(ierrprimme>=-39,PetscObjectComm((PetscObject)svd),PETSC_ERR_LIB,"PRIMME library failed with error code=%" PetscInt_FMT ": runtime error; check PRIMME's manual",ierrprimme);
     283             :     }
     284             :   }
     285           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     286             : }
     287             : 
     288           6 : static PetscErrorCode SVDReset_PRIMME(SVD svd)
     289             : {
     290           6 :   SVD_PRIMME     *ops = (SVD_PRIMME*)svd->data;
     291             : 
     292           6 :   PetscFunctionBegin;
     293           6 :   primme_svds_free(&ops->primme);
     294           6 :   PetscCall(VecDestroy(&ops->x));
     295           6 :   PetscCall(VecDestroy(&ops->y));
     296           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     297             : }
     298             : 
     299           6 : static PetscErrorCode SVDDestroy_PRIMME(SVD svd)
     300             : {
     301           6 :   PetscFunctionBegin;
     302           6 :   PetscCall(PetscFree(svd->data));
     303           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetBlockSize_C",NULL));
     304           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetBlockSize_C",NULL));
     305           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetMethod_C",NULL));
     306           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetMethod_C",NULL));
     307           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     308             : }
     309             : 
     310           0 : static PetscErrorCode SVDView_PRIMME(SVD svd,PetscViewer viewer)
     311             : {
     312           0 :   PetscBool      isascii;
     313           0 :   SVD_PRIMME     *ctx = (SVD_PRIMME*)svd->data;
     314           0 :   PetscMPIInt    rank;
     315             : 
     316           0 :   PetscFunctionBegin;
     317           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     318           0 :   if (isascii) {
     319           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  block size=%" PetscInt_FMT "\n",ctx->bs));
     320           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  solver method: %s\n",SVDPRIMMEMethods[(SVDPRIMMEMethod)ctx->method]));
     321             : 
     322             :     /* Display PRIMME params */
     323           0 :     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)svd),&rank));
     324           0 :     if (!rank) primme_svds_display_params(ctx->primme);
     325             :   }
     326           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     327             : }
     328             : 
     329           5 : static PetscErrorCode SVDSetFromOptions_PRIMME(SVD svd,PetscOptionItems *PetscOptionsObject)
     330             : {
     331           5 :   SVD_PRIMME      *ctx = (SVD_PRIMME*)svd->data;
     332           5 :   PetscInt        bs;
     333           5 :   SVDPRIMMEMethod meth;
     334           5 :   PetscBool       flg;
     335             : 
     336           5 :   PetscFunctionBegin;
     337           5 :   PetscOptionsHeadBegin(PetscOptionsObject,"SVD PRIMME Options");
     338             : 
     339           5 :     PetscCall(PetscOptionsInt("-svd_primme_blocksize","Maximum block size","SVDPRIMMESetBlockSize",ctx->bs,&bs,&flg));
     340           5 :     if (flg) PetscCall(SVDPRIMMESetBlockSize(svd,bs));
     341             : 
     342           5 :     PetscCall(PetscOptionsEnum("-svd_primme_method","Method for solving the singular value problem","SVDPRIMMESetMethod",SVDPRIMMEMethods,(PetscEnum)ctx->method,(PetscEnum*)&meth,&flg));
     343           5 :     if (flg) PetscCall(SVDPRIMMESetMethod(svd,meth));
     344             : 
     345           5 :   PetscOptionsHeadEnd();
     346           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     347             : }
     348             : 
     349           3 : static PetscErrorCode SVDPRIMMESetBlockSize_PRIMME(SVD svd,PetscInt bs)
     350             : {
     351           3 :   SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
     352             : 
     353           3 :   PetscFunctionBegin;
     354           3 :   if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) ops->bs = 0;
     355             :   else {
     356           3 :     PetscCheck(bs>0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"PRIMME: block size must be positive");
     357           3 :     ops->bs = bs;
     358             :   }
     359           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     360             : }
     361             : 
     362             : /*@
     363             :    SVDPRIMMESetBlockSize - The maximum block size that PRIMME will try to use.
     364             : 
     365             :    Logically Collective
     366             : 
     367             :    Input Parameters:
     368             : +  svd - the singular value solver context
     369             : -  bs - block size
     370             : 
     371             :    Options Database Key:
     372             : .  -svd_primme_blocksize - Sets the max allowed block size value
     373             : 
     374             :    Notes:
     375             :    If the block size is not set, the value established by primme_svds_initialize
     376             :    is used.
     377             : 
     378             :    The user should set the block size based on the architecture specifics
     379             :    of the target computer, as well as any a priori knowledge of multiplicities.
     380             :    The code does NOT require bs > 1 to find multiple eigenvalues. For some
     381             :    methods, keeping bs = 1 yields the best overall performance.
     382             : 
     383             :    Level: advanced
     384             : 
     385             : .seealso: SVDPRIMMEGetBlockSize()
     386             : @*/
     387           3 : PetscErrorCode SVDPRIMMESetBlockSize(SVD svd,PetscInt bs)
     388             : {
     389           3 :   PetscFunctionBegin;
     390           3 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     391           9 :   PetscValidLogicalCollectiveInt(svd,bs,2);
     392           3 :   PetscTryMethod(svd,"SVDPRIMMESetBlockSize_C",(SVD,PetscInt),(svd,bs));
     393           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     394             : }
     395             : 
     396           2 : static PetscErrorCode SVDPRIMMEGetBlockSize_PRIMME(SVD svd,PetscInt *bs)
     397             : {
     398           2 :   SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
     399             : 
     400           2 :   PetscFunctionBegin;
     401           2 :   *bs = ops->bs;
     402           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     403             : }
     404             : 
     405             : /*@
     406             :    SVDPRIMMEGetBlockSize - Get the maximum block size the code will try to use.
     407             : 
     408             :    Not Collective
     409             : 
     410             :    Input Parameter:
     411             : .  svd - the singular value solver context
     412             : 
     413             :    Output Parameter:
     414             : .  bs - returned block size
     415             : 
     416             :    Level: advanced
     417             : 
     418             : .seealso: SVDPRIMMESetBlockSize()
     419             : @*/
     420           2 : PetscErrorCode SVDPRIMMEGetBlockSize(SVD svd,PetscInt *bs)
     421             : {
     422           2 :   PetscFunctionBegin;
     423           2 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     424           2 :   PetscAssertPointer(bs,2);
     425           2 :   PetscUseMethod(svd,"SVDPRIMMEGetBlockSize_C",(SVD,PetscInt*),(svd,bs));
     426           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     427             : }
     428             : 
     429           3 : static PetscErrorCode SVDPRIMMESetMethod_PRIMME(SVD svd,SVDPRIMMEMethod method)
     430             : {
     431           3 :   SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
     432             : 
     433           3 :   PetscFunctionBegin;
     434           3 :   ops->method = (primme_svds_preset_method)method;
     435           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     436             : }
     437             : 
     438             : /*@
     439             :    SVDPRIMMESetMethod - Sets the method for the PRIMME SVD solver.
     440             : 
     441             :    Logically Collective
     442             : 
     443             :    Input Parameters:
     444             : +  svd - the singular value solver context
     445             : -  method - method that will be used by PRIMME
     446             : 
     447             :    Options Database Key:
     448             : .  -svd_primme_method - Sets the method for the PRIMME SVD solver
     449             : 
     450             :    Note:
     451             :    If not set, the method defaults to SVD_PRIMME_HYBRID.
     452             : 
     453             :    Level: advanced
     454             : 
     455             : .seealso: SVDPRIMMEGetMethod(), SVDPRIMMEMethod
     456             : @*/
     457           3 : PetscErrorCode SVDPRIMMESetMethod(SVD svd,SVDPRIMMEMethod method)
     458             : {
     459           3 :   PetscFunctionBegin;
     460           3 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     461           9 :   PetscValidLogicalCollectiveEnum(svd,method,2);
     462           3 :   PetscTryMethod(svd,"SVDPRIMMESetMethod_C",(SVD,SVDPRIMMEMethod),(svd,method));
     463           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     464             : }
     465             : 
     466           2 : static PetscErrorCode SVDPRIMMEGetMethod_PRIMME(SVD svd,SVDPRIMMEMethod *method)
     467             : {
     468           2 :   SVD_PRIMME *ops = (SVD_PRIMME*)svd->data;
     469             : 
     470           2 :   PetscFunctionBegin;
     471           2 :   *method = (SVDPRIMMEMethod)ops->method;
     472           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     473             : }
     474             : 
     475             : /*@
     476             :    SVDPRIMMEGetMethod - Gets the method for the PRIMME SVD solver.
     477             : 
     478             :    Not Collective
     479             : 
     480             :    Input Parameter:
     481             : .  svd - the singular value solver context
     482             : 
     483             :    Output Parameter:
     484             : .  method - method that will be used by PRIMME
     485             : 
     486             :    Level: advanced
     487             : 
     488             : .seealso: SVDPRIMMESetMethod(), SVDPRIMMEMethod
     489             : @*/
     490           2 : PetscErrorCode SVDPRIMMEGetMethod(SVD svd,SVDPRIMMEMethod *method)
     491             : {
     492           2 :   PetscFunctionBegin;
     493           2 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     494           2 :   PetscAssertPointer(method,2);
     495           2 :   PetscUseMethod(svd,"SVDPRIMMEGetMethod_C",(SVD,SVDPRIMMEMethod*),(svd,method));
     496           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     497             : }
     498             : 
     499           6 : SLEPC_EXTERN PetscErrorCode SVDCreate_PRIMME(SVD svd)
     500             : {
     501           6 :   SVD_PRIMME     *primme;
     502             : 
     503           6 :   PetscFunctionBegin;
     504           6 :   PetscCall(PetscNew(&primme));
     505           6 :   svd->data = (void*)primme;
     506             : 
     507           6 :   primme_svds_initialize(&primme->primme);
     508           6 :   primme->bs = 0;
     509           6 :   primme->method = (primme_svds_preset_method)SVD_PRIMME_HYBRID;
     510           6 :   primme->svd = svd;
     511             : 
     512           6 :   svd->ops->solve          = SVDSolve_PRIMME;
     513           6 :   svd->ops->setup          = SVDSetUp_PRIMME;
     514           6 :   svd->ops->setfromoptions = SVDSetFromOptions_PRIMME;
     515           6 :   svd->ops->destroy        = SVDDestroy_PRIMME;
     516           6 :   svd->ops->reset          = SVDReset_PRIMME;
     517           6 :   svd->ops->view           = SVDView_PRIMME;
     518             : 
     519           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetBlockSize_C",SVDPRIMMESetBlockSize_PRIMME));
     520           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetBlockSize_C",SVDPRIMMEGetBlockSize_PRIMME));
     521           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMESetMethod_C",SVDPRIMMESetMethod_PRIMME));
     522           6 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDPRIMMEGetMethod_C",SVDPRIMMEGetMethod_PRIMME));
     523           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     524             : }