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             :    SLEPc singular value solver: "cross"
      12             : 
      13             :    Method: Uses a Hermitian eigensolver for A^T*A
      14             : */
      15             : 
      16             : #include <slepc/private/svdimpl.h>                /*I "slepcsvd.h" I*/
      17             : 
      18             : typedef struct {
      19             :   PetscBool explicitmatrix;
      20             :   EPS       eps;
      21             :   PetscBool usereps;
      22             :   Mat       C,D;
      23             : } SVD_CROSS;
      24             : 
      25             : typedef struct {
      26             :   Mat       A,AT;
      27             :   Vec       w,diag,omega;
      28             :   PetscBool swapped;
      29             : } SVD_CROSS_SHELL;
      30             : 
      31       29920 : static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
      32             : {
      33       29920 :   SVD_CROSS_SHELL *ctx;
      34             : 
      35       29920 :   PetscFunctionBegin;
      36       29920 :   PetscCall(MatShellGetContext(B,&ctx));
      37       29920 :   PetscCall(MatMult(ctx->A,x,ctx->w));
      38       29920 :   if (ctx->omega && !ctx->swapped) PetscCall(VecPointwiseMult(ctx->w,ctx->w,ctx->omega));
      39       29920 :   PetscCall(MatMult(ctx->AT,ctx->w,y));
      40       29920 :   PetscFunctionReturn(PETSC_SUCCESS);
      41             : }
      42             : 
      43           9 : static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
      44             : {
      45           9 :   SVD_CROSS_SHELL   *ctx;
      46           9 :   PetscMPIInt       len;
      47           9 :   PetscInt          N,n,i,j,start,end,ncols;
      48           9 :   PetscScalar       *work1,*work2,*diag;
      49           9 :   const PetscInt    *cols;
      50           9 :   const PetscScalar *vals;
      51             : 
      52           9 :   PetscFunctionBegin;
      53           9 :   PetscCall(MatShellGetContext(B,&ctx));
      54           9 :   if (!ctx->diag) {
      55             :     /* compute diagonal from rows and store in ctx->diag */
      56           9 :     PetscCall(VecDuplicate(d,&ctx->diag));
      57           9 :     PetscCall(MatGetSize(ctx->A,NULL,&N));
      58           9 :     PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
      59           9 :     PetscCall(PetscCalloc2(N,&work1,N,&work2));
      60           9 :     if (ctx->swapped) {
      61           2 :       PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
      62          42 :       for (i=start;i<end;i++) {
      63          40 :         PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
      64         120 :         for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
      65          40 :         PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
      66             :       }
      67             :     } else {
      68           7 :       PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
      69         338 :       for (i=start;i<end;i++) {
      70         331 :         PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
      71        2214 :         for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
      72         331 :         PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
      73             :       }
      74             :     }
      75           9 :     PetscCall(PetscMPIIntCast(N,&len));
      76          27 :     PetscCallMPI(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
      77           9 :     PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
      78           9 :     PetscCall(VecGetArrayWrite(ctx->diag,&diag));
      79         364 :     for (i=start;i<end;i++) diag[i-start] = work2[i];
      80           9 :     PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
      81           9 :     PetscCall(PetscFree2(work1,work2));
      82             :   }
      83           9 :   PetscCall(VecCopy(ctx->diag,d));
      84           9 :   PetscFunctionReturn(PETSC_SUCCESS);
      85             : }
      86             : 
      87          48 : static PetscErrorCode MatDestroy_Cross(Mat B)
      88             : {
      89          48 :   SVD_CROSS_SHELL *ctx;
      90             : 
      91          48 :   PetscFunctionBegin;
      92          48 :   PetscCall(MatShellGetContext(B,&ctx));
      93          48 :   PetscCall(VecDestroy(&ctx->w));
      94          48 :   PetscCall(VecDestroy(&ctx->diag));
      95          48 :   PetscCall(PetscFree(ctx));
      96          48 :   PetscFunctionReturn(PETSC_SUCCESS);
      97             : }
      98             : 
      99          77 : static PetscErrorCode SVDCrossGetProductMat(SVD svd,Mat A,Mat AT,Mat *C)
     100             : {
     101          77 :   SVD_CROSS       *cross = (SVD_CROSS*)svd->data;
     102          77 :   SVD_CROSS_SHELL *ctx;
     103          77 :   PetscInt        n;
     104          77 :   VecType         vtype;
     105          77 :   Mat             B;
     106             : 
     107          77 :   PetscFunctionBegin;
     108          77 :   if (cross->explicitmatrix) {
     109          54 :     if (!svd->ishyperbolic || svd->swapped) B = (!svd->expltrans && svd->swapped)? AT: A;
     110             :     else {  /* duplicate A and scale by signature */
     111           2 :       PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
     112           2 :       PetscCall(MatDiagonalScale(B,svd->omega,NULL));
     113             :     }
     114          29 :     if (svd->expltrans) {  /* explicit transpose */
     115          25 :       PetscCall(MatProductCreate(AT,B,NULL,C));
     116          25 :       PetscCall(MatProductSetType(*C,MATPRODUCT_AB));
     117             :     } else {  /* implicit transpose */
     118             : #if defined(PETSC_USE_COMPLEX)
     119             :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Must use explicit transpose with complex scalars");
     120             : #else
     121           4 :       if (!svd->swapped) {
     122           2 :         PetscCall(MatProductCreate(A,B,NULL,C));
     123           2 :         PetscCall(MatProductSetType(*C,MATPRODUCT_AtB));
     124             :       } else {
     125           2 :         PetscCall(MatProductCreate(B,AT,NULL,C));
     126           2 :         PetscCall(MatProductSetType(*C,MATPRODUCT_ABt));
     127             :       }
     128             : #endif
     129             :     }
     130          29 :     PetscCall(MatProductSetFromOptions(*C));
     131          29 :     PetscCall(MatProductSymbolic(*C));
     132          29 :     PetscCall(MatProductNumeric(*C));
     133          29 :     if (svd->ishyperbolic && !svd->swapped) PetscCall(MatDestroy(&B));
     134             :   } else {
     135          48 :     PetscCall(PetscNew(&ctx));
     136          48 :     ctx->A       = A;
     137          48 :     ctx->AT      = AT;
     138          48 :     ctx->omega   = svd->omega;
     139          48 :     ctx->swapped = svd->swapped;
     140          48 :     PetscCall(MatCreateVecs(A,NULL,&ctx->w));
     141          48 :     PetscCall(MatGetLocalSize(A,NULL,&n));
     142          48 :     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)ctx,C));
     143          48 :     PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cross));
     144          48 :     if (!svd->ishyperbolic || svd->swapped) PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cross));
     145          48 :     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cross));
     146          48 :     PetscCall(MatGetVecType(A,&vtype));
     147          48 :     PetscCall(MatSetVecType(*C,vtype));
     148             :   }
     149          77 :   PetscFunctionReturn(PETSC_SUCCESS);
     150             : }
     151             : 
     152             : /* Convergence test relative to the norm of R (used in GSVD only) */
     153         482 : static PetscErrorCode EPSConv_Cross(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     154             : {
     155         482 :   SVD svd = (SVD)ctx;
     156             : 
     157         482 :   PetscFunctionBegin;
     158         482 :   *errest = res/PetscMax(svd->nrma,svd->nrmb);
     159         482 :   PetscFunctionReturn(PETSC_SUCCESS);
     160             : }
     161             : 
     162          62 : static PetscErrorCode SVDSetUp_Cross(SVD svd)
     163             : {
     164          62 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     165          62 :   ST             st;
     166          62 :   PetscBool      trackall,issinv,isks;
     167          62 :   EPSProblemType ptype;
     168          62 :   EPSWhich       which;
     169          62 :   Mat            Omega;
     170          62 :   MatType        Atype;
     171          62 :   PetscInt       n,N;
     172             : 
     173          62 :   PetscFunctionBegin;
     174          62 :   if (svd->nsv==0 && svd->stop!=SVD_STOP_THRESHOLD) svd->nsv = 1;
     175          62 :   if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
     176          62 :   PetscCall(MatDestroy(&cross->C));
     177          62 :   PetscCall(MatDestroy(&cross->D));
     178          62 :   PetscCall(SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C));
     179          62 :   if (svd->isgeneralized) {
     180          15 :     PetscCall(SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D));
     181          15 :     PetscCall(EPSSetOperators(cross->eps,cross->C,cross->D));
     182          15 :     PetscCall(EPSGetProblemType(cross->eps,&ptype));
     183          15 :     if (!ptype) PetscCall(EPSSetProblemType(cross->eps,EPS_GHEP));
     184          47 :   } else if (svd->ishyperbolic && svd->swapped) {
     185           4 :     PetscCall(MatGetType(svd->OP,&Atype));
     186           4 :     PetscCall(MatGetSize(svd->A,NULL,&N));
     187           4 :     PetscCall(MatGetLocalSize(svd->A,NULL,&n));
     188           4 :     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
     189           4 :     PetscCall(MatSetSizes(Omega,n,n,N,N));
     190           4 :     PetscCall(MatSetType(Omega,Atype));
     191           4 :     PetscCall(MatDiagonalSet(Omega,svd->omega,INSERT_VALUES));
     192           4 :     PetscCall(EPSSetOperators(cross->eps,cross->C,Omega));
     193           4 :     PetscCall(EPSSetProblemType(cross->eps,EPS_GHIEP));
     194           4 :     PetscCall(MatDestroy(&Omega));
     195             :   } else {
     196          43 :     PetscCall(EPSSetOperators(cross->eps,cross->C,NULL));
     197          43 :     PetscCall(EPSSetProblemType(cross->eps,EPS_HEP));
     198             :   }
     199          62 :   if (!cross->usereps) {
     200          60 :     PetscCall(EPSGetST(cross->eps,&st));
     201          60 :     PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
     202          60 :     PetscCall(PetscObjectTypeCompare((PetscObject)cross->eps,EPSKRYLOVSCHUR,&isks));
     203          60 :     if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
     204           4 :       if (cross->explicitmatrix && isks && !issinv) {  /* default to shift-and-invert */
     205           2 :         PetscCall(STSetType(st,STSINVERT));
     206           2 :         PetscCall(EPSSetTarget(cross->eps,0.0));
     207             :         which = EPS_TARGET_REAL;
     208           2 :       } else which = issinv?EPS_TARGET_REAL:EPS_SMALLEST_REAL;
     209             :     } else {
     210          56 :       if (issinv) which = EPS_TARGET_MAGNITUDE;
     211          55 :       else if (svd->ishyperbolic) which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
     212          42 :       else which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
     213             :     }
     214          60 :     PetscCall(EPSSetWhichEigenpairs(cross->eps,which));
     215         116 :     PetscCall(EPSSetDimensions(cross->eps,svd->nsv?svd->nsv:PETSC_CURRENT,svd->ncv,svd->mpd));
     216          60 :     if (svd->stop==SVD_STOP_THRESHOLD) PetscCall(EPSSetThreshold(cross->eps,svd->thres*svd->thres,svd->threlative));
     217          78 :     PetscCall(EPSSetTolerances(cross->eps,svd->tol==(PetscReal)PETSC_DETERMINE?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it));
     218          60 :     switch (svd->conv) {
     219           3 :     case SVD_CONV_ABS:
     220           3 :       PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS));break;
     221          42 :     case SVD_CONV_REL:
     222          42 :       PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_REL));break;
     223          15 :     case SVD_CONV_NORM:
     224          15 :       if (svd->isgeneralized) {
     225          15 :         if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
     226          15 :         if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
     227          15 :         PetscCall(EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL));
     228             :       } else {
     229           0 :         PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM));break;
     230             :       }
     231             :       break;
     232           0 :     case SVD_CONV_MAXIT:
     233           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
     234           0 :     case SVD_CONV_USER:
     235           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
     236             :     }
     237           2 :   }
     238             :   /* Transfer the trackall option from svd to eps */
     239          62 :   PetscCall(SVDGetTrackAll(svd,&trackall));
     240          62 :   PetscCall(EPSSetTrackAll(cross->eps,trackall));
     241             :   /* Transfer the initial space from svd to eps */
     242          62 :   if (svd->nini<0) {
     243           5 :     PetscCall(EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS));
     244           5 :     PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
     245             :   }
     246          62 :   PetscCall(EPSSetUp(cross->eps));
     247          62 :   PetscCall(EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd));
     248          62 :   PetscCall(EPSGetTolerances(cross->eps,NULL,&svd->max_it));
     249          62 :   if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
     250             : 
     251          62 :   svd->leftbasis = PETSC_FALSE;
     252          62 :   PetscCall(SVDAllocateSolution(svd,0));
     253          62 :   PetscFunctionReturn(PETSC_SUCCESS);
     254             : }
     255             : 
     256          62 : static PetscErrorCode SVDSolve_Cross(SVD svd)
     257             : {
     258          62 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     259          62 :   PetscInt       i;
     260          62 :   PetscScalar    lambda;
     261          62 :   PetscReal      sigma;
     262             : 
     263          62 :   PetscFunctionBegin;
     264          62 :   PetscCall(EPSSolve(cross->eps));
     265          62 :   PetscCall(EPSGetConverged(cross->eps,&svd->nconv));
     266          62 :   PetscCall(EPSGetIterationNumber(cross->eps,&svd->its));
     267          62 :   PetscCall(EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason));
     268         912 :   for (i=0;i<svd->nconv;i++) {
     269         850 :     PetscCall(EPSGetEigenvalue(cross->eps,i,&lambda,NULL));
     270         850 :     sigma = PetscRealPart(lambda);
     271         850 :     if (svd->ishyperbolic) svd->sigma[i] = PetscSqrtReal(PetscAbsReal(sigma));
     272             :     else {
     273         369 :       PetscCheck(sigma>-10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_FP,"Negative eigenvalue computed by EPS: %g",(double)sigma);
     274         369 :       if (sigma<0.0) {
     275           2 :         PetscCall(PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma));
     276             :         sigma = 0.0;
     277             :       }
     278         369 :       svd->sigma[i] = PetscSqrtReal(sigma);
     279             :     }
     280             :   }
     281          62 :   PetscFunctionReturn(PETSC_SUCCESS);
     282             : }
     283             : 
     284          56 : static PetscErrorCode SVDComputeVectors_Cross(SVD svd)
     285             : {
     286          56 :   SVD_CROSS         *cross = (SVD_CROSS*)svd->data;
     287          56 :   PetscInt          i,mloc,ploc;
     288          56 :   Vec               u,v,x,uv,w,omega2=NULL;
     289          56 :   Mat               Omega;
     290          56 :   PetscScalar       *dst,alpha,lambda,*varray;
     291          56 :   const PetscScalar *src;
     292          56 :   PetscReal         nrm;
     293             : 
     294          56 :   PetscFunctionBegin;
     295          56 :   if (svd->isgeneralized) {
     296          15 :     PetscCall(MatCreateVecs(svd->A,NULL,&u));
     297          15 :     PetscCall(VecGetLocalSize(u,&mloc));
     298          15 :     PetscCall(MatCreateVecs(svd->B,NULL,&v));
     299          15 :     PetscCall(VecGetLocalSize(v,&ploc));
     300         136 :     for (i=0;i<svd->nconv;i++) {
     301         121 :       PetscCall(BVGetColumn(svd->V,i,&x));
     302         121 :       PetscCall(EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL));
     303         121 :       PetscCall(MatMult(svd->A,x,u));     /* u_i*c_i/alpha = A*x_i */
     304         121 :       PetscCall(VecNormalize(u,NULL));
     305         121 :       PetscCall(MatMult(svd->B,x,v));     /* v_i*s_i/alpha = B*x_i */
     306         121 :       PetscCall(VecNormalize(v,&nrm));    /* ||v||_2 = s_i/alpha   */
     307         121 :       alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm);    /* alpha=s_i/||v||_2 */
     308         121 :       PetscCall(VecScale(x,alpha));
     309         121 :       PetscCall(BVRestoreColumn(svd->V,i,&x));
     310             :       /* copy [u;v] to U[i] */
     311         121 :       PetscCall(BVGetColumn(svd->U,i,&uv));
     312         121 :       PetscCall(VecGetArrayWrite(uv,&dst));
     313         121 :       PetscCall(VecGetArrayRead(u,&src));
     314         121 :       PetscCall(PetscArraycpy(dst,src,mloc));
     315         121 :       PetscCall(VecRestoreArrayRead(u,&src));
     316         121 :       PetscCall(VecGetArrayRead(v,&src));
     317         121 :       PetscCall(PetscArraycpy(dst+mloc,src,ploc));
     318         121 :       PetscCall(VecRestoreArrayRead(v,&src));
     319         121 :       PetscCall(VecRestoreArrayWrite(uv,&dst));
     320         121 :       PetscCall(BVRestoreColumn(svd->U,i,&uv));
     321             :     }
     322          15 :     PetscCall(VecDestroy(&v));
     323          15 :     PetscCall(VecDestroy(&u));
     324          41 :   } else if (svd->ishyperbolic && svd->swapped) {  /* was solved as GHIEP, set u=Omega*u and normalize */
     325           4 :     PetscCall(EPSGetOperators(cross->eps,NULL,&Omega));
     326           4 :     PetscCall(MatCreateVecs(Omega,&w,NULL));
     327           4 :     PetscCall(VecCreateSeq(PETSC_COMM_SELF,svd->ncv,&omega2));
     328           4 :     PetscCall(VecGetArrayWrite(omega2,&varray));
     329         120 :     for (i=0;i<svd->nconv;i++) {
     330         116 :       PetscCall(BVGetColumn(svd->V,i,&v));
     331         116 :       PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
     332         116 :       PetscCall(MatMult(Omega,v,w));
     333         116 :       PetscCall(VecDot(v,w,&alpha));
     334         116 :       svd->sign[i] = PetscSign(PetscRealPart(alpha));
     335         116 :       varray[i] = svd->sign[i];
     336         116 :       alpha = 1.0/PetscSqrtScalar(PetscAbsScalar(alpha));
     337         116 :       PetscCall(VecScale(w,alpha));
     338         116 :       PetscCall(VecCopy(w,v));
     339         116 :       PetscCall(BVRestoreColumn(svd->V,i,&v));
     340             :     }
     341           4 :     PetscCall(BVSetSignature(svd->V,omega2));
     342           4 :     PetscCall(VecRestoreArrayWrite(omega2,&varray));
     343           4 :     PetscCall(VecDestroy(&omega2));
     344           4 :     PetscCall(VecDestroy(&w));
     345           4 :     PetscCall(SVDComputeVectors_Left(svd));
     346             :   } else {
     347         635 :     for (i=0;i<svd->nconv;i++) {
     348         598 :       PetscCall(BVGetColumn(svd->V,i,&v));
     349         598 :       PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
     350         598 :       PetscCall(BVRestoreColumn(svd->V,i,&v));
     351             :     }
     352          37 :     PetscCall(SVDComputeVectors_Left(svd));
     353             :   }
     354          56 :   PetscFunctionReturn(PETSC_SUCCESS);
     355             : }
     356             : 
     357        1545 : static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
     358             : {
     359        1545 :   PetscInt       i,ncv;
     360        1545 :   SVD            svd = (SVD)ctx;
     361        1545 :   SVD_CROSS      *cross;
     362        1545 :   PetscScalar    er,ei;
     363        1545 :   ST             st;
     364             : 
     365        1545 :   PetscFunctionBegin;
     366        1545 :   if (svd->stop==SVD_STOP_THRESHOLD) {
     367         115 :     cross = (SVD_CROSS*)svd->data;
     368         115 :     PetscCall(EPSGetDimensions(cross->eps,NULL,&ncv,NULL));
     369         115 :     if (ncv!=svd->ncv) {  /* reallocate */
     370           1 :       PetscCall(SVDReallocateSolution(svd,ncv));
     371          13 :       for (i=svd->ncv;i<ncv;i++) svd->perm[i] = i;
     372           1 :       svd->ncv = ncv;
     373             :     }
     374             :   }
     375        1545 :   PetscCall(EPSGetST(eps,&st));
     376       24919 :   for (i=0;i<PetscMin(nest,svd->ncv);i++) {
     377       23374 :     er = eigr[i]; ei = eigi[i];
     378       23374 :     PetscCall(STBackTransform(st,1,&er,&ei));
     379       23374 :     svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
     380       23374 :     svd->errest[i] = errest[i];
     381             :   }
     382        1545 :   PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
     383        1545 :   PetscFunctionReturn(PETSC_SUCCESS);
     384             : }
     385             : 
     386          49 : static PetscErrorCode SVDSetFromOptions_Cross(SVD svd,PetscOptionItems *PetscOptionsObject)
     387             : {
     388          49 :   PetscBool      set,val;
     389          49 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     390          49 :   ST             st;
     391             : 
     392          49 :   PetscFunctionBegin;
     393          49 :   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cross Options");
     394             : 
     395          49 :     PetscCall(PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set));
     396          49 :     if (set) PetscCall(SVDCrossSetExplicitMatrix(svd,val));
     397             : 
     398          49 :   PetscOptionsHeadEnd();
     399             : 
     400          49 :   if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
     401          49 :   if (!cross->explicitmatrix && !cross->usereps) {
     402             :     /* use as default an ST with shell matrix and Jacobi */
     403          29 :     PetscCall(EPSGetST(cross->eps,&st));
     404          29 :     PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
     405             :   }
     406          49 :   PetscCall(EPSSetFromOptions(cross->eps));
     407          49 :   PetscFunctionReturn(PETSC_SUCCESS);
     408             : }
     409             : 
     410          27 : static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
     411             : {
     412          27 :   SVD_CROSS *cross = (SVD_CROSS*)svd->data;
     413             : 
     414          27 :   PetscFunctionBegin;
     415          27 :   if (cross->explicitmatrix != explicitmatrix) {
     416          18 :     cross->explicitmatrix = explicitmatrix;
     417          18 :     svd->state = SVD_STATE_INITIAL;
     418             :   }
     419          27 :   PetscFunctionReturn(PETSC_SUCCESS);
     420             : }
     421             : 
     422             : /*@
     423             :    SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
     424             :    be computed explicitly.
     425             : 
     426             :    Logically Collective
     427             : 
     428             :    Input Parameters:
     429             : +  svd         - singular value solver
     430             : -  explicitmat - boolean flag indicating if A^T*A is built explicitly
     431             : 
     432             :    Options Database Key:
     433             : .  -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag
     434             : 
     435             :    Level: advanced
     436             : 
     437             : .seealso: SVDCrossGetExplicitMatrix()
     438             : @*/
     439          27 : PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
     440             : {
     441          27 :   PetscFunctionBegin;
     442          27 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     443          81 :   PetscValidLogicalCollectiveBool(svd,explicitmat,2);
     444          27 :   PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
     445          27 :   PetscFunctionReturn(PETSC_SUCCESS);
     446             : }
     447             : 
     448           2 : static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
     449             : {
     450           2 :   SVD_CROSS *cross = (SVD_CROSS*)svd->data;
     451             : 
     452           2 :   PetscFunctionBegin;
     453           2 :   *explicitmat = cross->explicitmatrix;
     454           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     455             : }
     456             : 
     457             : /*@
     458             :    SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.
     459             : 
     460             :    Not Collective
     461             : 
     462             :    Input Parameter:
     463             : .  svd  - singular value solver
     464             : 
     465             :    Output Parameter:
     466             : .  explicitmat - the mode flag
     467             : 
     468             :    Level: advanced
     469             : 
     470             : .seealso: SVDCrossSetExplicitMatrix()
     471             : @*/
     472           2 : PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
     473             : {
     474           2 :   PetscFunctionBegin;
     475           2 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     476           2 :   PetscAssertPointer(explicitmat,2);
     477           2 :   PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
     478           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     479             : }
     480             : 
     481           2 : static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
     482             : {
     483           2 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     484             : 
     485           2 :   PetscFunctionBegin;
     486           2 :   PetscCall(PetscObjectReference((PetscObject)eps));
     487           2 :   PetscCall(EPSDestroy(&cross->eps));
     488           2 :   cross->eps     = eps;
     489           2 :   cross->usereps = PETSC_TRUE;
     490           2 :   svd->state     = SVD_STATE_INITIAL;
     491           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     492             : }
     493             : 
     494             : /*@
     495             :    SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
     496             :    singular value solver.
     497             : 
     498             :    Collective
     499             : 
     500             :    Input Parameters:
     501             : +  svd - singular value solver
     502             : -  eps - the eigensolver object
     503             : 
     504             :    Level: advanced
     505             : 
     506             : .seealso: SVDCrossGetEPS()
     507             : @*/
     508           2 : PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
     509             : {
     510           2 :   PetscFunctionBegin;
     511           2 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     512           2 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
     513           2 :   PetscCheckSameComm(svd,1,eps,2);
     514           2 :   PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
     515           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     516             : }
     517             : 
     518          49 : static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
     519             : {
     520          49 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     521             : 
     522          49 :   PetscFunctionBegin;
     523          49 :   if (!cross->eps) {
     524          49 :     PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps));
     525          49 :     PetscCall(PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1));
     526          49 :     PetscCall(EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix));
     527          49 :     PetscCall(EPSAppendOptionsPrefix(cross->eps,"svd_cross_"));
     528          49 :     PetscCall(PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options));
     529          49 :     PetscCall(EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_MAGNITUDE));
     530          49 :     PetscCall(EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL));
     531             :   }
     532          49 :   *eps = cross->eps;
     533          49 :   PetscFunctionReturn(PETSC_SUCCESS);
     534             : }
     535             : 
     536             : /*@
     537             :    SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
     538             :    to the singular value solver.
     539             : 
     540             :    Collective
     541             : 
     542             :    Input Parameter:
     543             : .  svd - singular value solver
     544             : 
     545             :    Output Parameter:
     546             : .  eps - the eigensolver object
     547             : 
     548             :    Level: advanced
     549             : 
     550             : .seealso: SVDCrossSetEPS()
     551             : @*/
     552          49 : PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
     553             : {
     554          49 :   PetscFunctionBegin;
     555          49 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     556          49 :   PetscAssertPointer(eps,2);
     557          49 :   PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
     558          49 :   PetscFunctionReturn(PETSC_SUCCESS);
     559             : }
     560             : 
     561           1 : static PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
     562             : {
     563           1 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     564           1 :   PetscBool      isascii;
     565             : 
     566           1 :   PetscFunctionBegin;
     567           1 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     568           1 :   if (isascii) {
     569           1 :     if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
     570           2 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cross->explicitmatrix?"explicit":"implicit"));
     571           1 :     PetscCall(PetscViewerASCIIPushTab(viewer));
     572           1 :     PetscCall(EPSView(cross->eps,viewer));
     573           1 :     PetscCall(PetscViewerASCIIPopTab(viewer));
     574             :   }
     575           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     576             : }
     577             : 
     578          52 : static PetscErrorCode SVDReset_Cross(SVD svd)
     579             : {
     580          52 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     581             : 
     582          52 :   PetscFunctionBegin;
     583          52 :   PetscCall(EPSReset(cross->eps));
     584          52 :   PetscCall(MatDestroy(&cross->C));
     585          52 :   PetscCall(MatDestroy(&cross->D));
     586          52 :   PetscFunctionReturn(PETSC_SUCCESS);
     587             : }
     588             : 
     589          51 : static PetscErrorCode SVDDestroy_Cross(SVD svd)
     590             : {
     591          51 :   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
     592             : 
     593          51 :   PetscFunctionBegin;
     594          51 :   PetscCall(EPSDestroy(&cross->eps));
     595          51 :   PetscCall(PetscFree(svd->data));
     596          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL));
     597          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL));
     598          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL));
     599          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL));
     600          51 :   PetscFunctionReturn(PETSC_SUCCESS);
     601             : }
     602             : 
     603          51 : SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
     604             : {
     605          51 :   SVD_CROSS      *cross;
     606             : 
     607          51 :   PetscFunctionBegin;
     608          51 :   PetscCall(PetscNew(&cross));
     609          51 :   svd->data = (void*)cross;
     610             : 
     611          51 :   svd->ops->solve          = SVDSolve_Cross;
     612          51 :   svd->ops->solveg         = SVDSolve_Cross;
     613          51 :   svd->ops->solveh         = SVDSolve_Cross;
     614          51 :   svd->ops->setup          = SVDSetUp_Cross;
     615          51 :   svd->ops->setfromoptions = SVDSetFromOptions_Cross;
     616          51 :   svd->ops->destroy        = SVDDestroy_Cross;
     617          51 :   svd->ops->reset          = SVDReset_Cross;
     618          51 :   svd->ops->view           = SVDView_Cross;
     619          51 :   svd->ops->computevectors = SVDComputeVectors_Cross;
     620          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross));
     621          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross));
     622          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross));
     623          51 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross));
     624          51 :   PetscFunctionReturn(PETSC_SUCCESS);
     625             : }