Actual source code: cross.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, 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:    SLEPc singular value solver: "cross"

     13:    Method: Uses a Hermitian eigensolver for A^T*A
     14: */

     16: #include <slepc/private/svdimpl.h>

     18: typedef struct {
     19:   PetscBool explicitmatrix;
     20:   EPS       eps;
     21:   PetscBool usereps;
     22:   Mat       C,D;
     23: } SVD_CROSS;

     25: typedef struct {
     26:   Mat       A,AT;
     27:   Vec       w,diag,omega;
     28:   PetscBool swapped;
     29: } SVD_CROSS_SHELL;

     31: static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
     32: {
     33:   SVD_CROSS_SHELL *ctx;

     35:   PetscFunctionBegin;
     36:   PetscCall(MatShellGetContext(B,&ctx));
     37:   PetscCall(MatMult(ctx->A,x,ctx->w));
     38:   if (ctx->omega && !ctx->swapped) PetscCall(VecPointwiseMult(ctx->w,ctx->w,ctx->omega));
     39:   PetscCall(MatMult(ctx->AT,ctx->w,y));
     40:   PetscFunctionReturn(PETSC_SUCCESS);
     41: }

     43: static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
     44: {
     45:   SVD_CROSS_SHELL   *ctx;
     46:   PetscMPIInt       len;
     47:   PetscInt          N,n,i,j,start,end,ncols;
     48:   PetscScalar       *work1,*work2,*diag;
     49:   const PetscInt    *cols;
     50:   const PetscScalar *vals;

     52:   PetscFunctionBegin;
     53:   PetscCall(MatShellGetContext(B,&ctx));
     54:   if (!ctx->diag) {
     55:     /* compute diagonal from rows and store in ctx->diag */
     56:     PetscCall(VecDuplicate(d,&ctx->diag));
     57:     PetscCall(MatGetSize(ctx->A,NULL,&N));
     58:     PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
     59:     PetscCall(PetscCalloc2(N,&work1,N,&work2));
     60:     if (ctx->swapped) {
     61:       PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
     62:       for (i=start;i<end;i++) {
     63:         PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
     64:         for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
     65:         PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
     66:       }
     67:     } else {
     68:       PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
     69:       for (i=start;i<end;i++) {
     70:         PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
     71:         for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
     72:         PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
     73:       }
     74:     }
     75:     PetscCall(PetscMPIIntCast(N,&len));
     76:     PetscCallMPI(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
     77:     PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
     78:     PetscCall(VecGetArrayWrite(ctx->diag,&diag));
     79:     for (i=start;i<end;i++) diag[i-start] = work2[i];
     80:     PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
     81:     PetscCall(PetscFree2(work1,work2));
     82:   }
     83:   PetscCall(VecCopy(ctx->diag,d));
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: static PetscErrorCode MatDestroy_Cross(Mat B)
     88: {
     89:   SVD_CROSS_SHELL *ctx;

     91:   PetscFunctionBegin;
     92:   PetscCall(MatShellGetContext(B,&ctx));
     93:   PetscCall(VecDestroy(&ctx->w));
     94:   PetscCall(VecDestroy(&ctx->diag));
     95:   PetscCall(PetscFree(ctx));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode SVDCrossGetProductMat(SVD svd,Mat A,Mat AT,Mat *C)
    100: {
    101:   SVD_CROSS       *cross = (SVD_CROSS*)svd->data;
    102:   SVD_CROSS_SHELL *ctx;
    103:   PetscInt        n;
    104:   VecType         vtype;
    105:   Mat             B;

    107:   PetscFunctionBegin;
    108:   if (cross->explicitmatrix) {
    109:     if (!svd->ishyperbolic || svd->swapped) B = (!svd->expltrans && svd->swapped)? AT: A;
    110:     else {  /* duplicate A and scale by signature */
    111:       PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
    112:       PetscCall(MatDiagonalScale(B,svd->omega,NULL));
    113:     }
    114:     if (svd->expltrans) {  /* explicit transpose */
    115:       PetscCall(MatProductCreate(AT,B,NULL,C));
    116:       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:       if (!svd->swapped) {
    122:         PetscCall(MatProductCreate(A,B,NULL,C));
    123:         PetscCall(MatProductSetType(*C,MATPRODUCT_AtB));
    124:       } else {
    125:         PetscCall(MatProductCreate(B,AT,NULL,C));
    126:         PetscCall(MatProductSetType(*C,MATPRODUCT_ABt));
    127:       }
    128: #endif
    129:     }
    130:     PetscCall(MatProductSetFromOptions(*C));
    131:     PetscCall(MatProductSymbolic(*C));
    132:     PetscCall(MatProductNumeric(*C));
    133:     if (svd->ishyperbolic && !svd->swapped) PetscCall(MatDestroy(&B));
    134:   } else {
    135:     PetscCall(PetscNew(&ctx));
    136:     ctx->A       = A;
    137:     ctx->AT      = AT;
    138:     ctx->omega   = svd->omega;
    139:     ctx->swapped = svd->swapped;
    140:     PetscCall(MatCreateVecs(A,NULL,&ctx->w));
    141:     PetscCall(MatGetLocalSize(A,NULL,&n));
    142:     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)ctx,C));
    143:     PetscCall(MatShellSetOperation(*C,MATOP_MULT,(PetscErrorCodeFn*)MatMult_Cross));
    144:     if (!svd->ishyperbolic || svd->swapped) PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_Cross));
    145:     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_Cross));
    146:     PetscCall(MatGetVecType(A,&vtype));
    147:     PetscCall(MatSetVecType(*C,vtype));
    148:   }
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: /* Convergence test relative to the norm of R (used in GSVD only) */
    153: static PetscErrorCode EPSConv_Cross(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
    154: {
    155:   SVD svd = (SVD)ctx;

    157:   PetscFunctionBegin;
    158:   *errest = res/PetscMax(svd->nrma,svd->nrmb);
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: static PetscErrorCode SVDSetUp_Cross(SVD svd)
    163: {
    164:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
    165:   ST             st;
    166:   PetscBool      trackall,issinv,isks;
    167:   EPSProblemType ptype;
    168:   EPSWhich       which;
    169:   Mat            Omega;
    170:   MatType        Atype;
    171:   PetscInt       n,N;

    173:   PetscFunctionBegin;
    174:   if (svd->nsv==0 && svd->stop!=SVD_STOP_THRESHOLD) svd->nsv = 1;
    175:   if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
    176:   PetscCall(MatDestroy(&cross->C));
    177:   PetscCall(MatDestroy(&cross->D));
    178:   PetscCall(SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C));
    179:   if (svd->isgeneralized) {
    180:     PetscCall(SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D));
    181:     PetscCall(EPSSetOperators(cross->eps,cross->C,cross->D));
    182:     PetscCall(EPSGetProblemType(cross->eps,&ptype));
    183:     if (!ptype) PetscCall(EPSSetProblemType(cross->eps,EPS_GHEP));
    184:   } else if (svd->ishyperbolic && svd->swapped) {
    185:     PetscCall(MatGetType(svd->OP,&Atype));
    186:     PetscCall(MatGetSize(svd->A,NULL,&N));
    187:     PetscCall(MatGetLocalSize(svd->A,NULL,&n));
    188:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
    189:     PetscCall(MatSetSizes(Omega,n,n,N,N));
    190:     PetscCall(MatSetType(Omega,Atype));
    191:     PetscCall(MatDiagonalSet(Omega,svd->omega,INSERT_VALUES));
    192:     PetscCall(EPSSetOperators(cross->eps,cross->C,Omega));
    193:     PetscCall(EPSSetProblemType(cross->eps,EPS_GHIEP));
    194:     PetscCall(MatDestroy(&Omega));
    195:   } else {
    196:     PetscCall(EPSSetOperators(cross->eps,cross->C,NULL));
    197:     PetscCall(EPSSetProblemType(cross->eps,EPS_HEP));
    198:   }
    199:   if (!cross->usereps) {
    200:     PetscCall(EPSGetST(cross->eps,&st));
    201:     PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
    202:     PetscCall(PetscObjectTypeCompare((PetscObject)cross->eps,EPSKRYLOVSCHUR,&isks));
    203:     if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
    204:       if (cross->explicitmatrix && isks && !issinv) {  /* default to shift-and-invert */
    205:         PetscCall(STSetType(st,STSINVERT));
    206:         PetscCall(EPSSetTarget(cross->eps,0.0));
    207:         which = EPS_TARGET_REAL;
    208:       } else which = issinv?EPS_TARGET_REAL:EPS_SMALLEST_REAL;
    209:     } else {
    210:       if (issinv) which = EPS_TARGET_MAGNITUDE;
    211:       else if (svd->ishyperbolic) which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
    212:       else which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
    213:     }
    214:     PetscCall(EPSSetWhichEigenpairs(cross->eps,which));
    215:     PetscCall(EPSSetDimensions(cross->eps,svd->nsv?svd->nsv:PETSC_CURRENT,svd->ncv,svd->mpd));
    216:     if (svd->stop==SVD_STOP_THRESHOLD) PetscCall(EPSSetThreshold(cross->eps,svd->thres*svd->thres,svd->threlative));
    217:     PetscCall(EPSSetTolerances(cross->eps,svd->tol==(PetscReal)PETSC_DETERMINE?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it));
    218:     switch (svd->conv) {
    219:     case SVD_CONV_ABS:
    220:       PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS));break;
    221:     case SVD_CONV_REL:
    222:       PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_REL));break;
    223:     case SVD_CONV_NORM:
    224:       if (svd->isgeneralized) {
    225:         if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
    226:         if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
    227:         PetscCall(EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL));
    228:       } else {
    229:         PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM));break;
    230:       }
    231:       break;
    232:     case SVD_CONV_MAXIT:
    233:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
    234:     case SVD_CONV_USER:
    235:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
    236:     }
    237:   }
    238:   /* Transfer the trackall option from svd to eps */
    239:   PetscCall(SVDGetTrackAll(svd,&trackall));
    240:   PetscCall(EPSSetTrackAll(cross->eps,trackall));
    241:   /* Transfer the initial space from svd to eps */
    242:   if (svd->nini<0) {
    243:     PetscCall(EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS));
    244:     PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
    245:   }
    246:   PetscCall(EPSSetUp(cross->eps));
    247:   PetscCall(EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd));
    248:   PetscCall(EPSGetTolerances(cross->eps,NULL,&svd->max_it));
    249:   if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;

    251:   svd->leftbasis = PETSC_FALSE;
    252:   PetscCall(SVDAllocateSolution(svd,0));
    253:   PetscFunctionReturn(PETSC_SUCCESS);
    254: }

    256: static PetscErrorCode SVDSolve_Cross(SVD svd)
    257: {
    258:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
    259:   PetscInt       i;
    260:   PetscScalar    lambda;
    261:   PetscReal      sigma;

    263:   PetscFunctionBegin;
    264:   PetscCall(EPSSolve(cross->eps));
    265:   PetscCall(EPSGetConverged(cross->eps,&svd->nconv));
    266:   PetscCall(EPSGetIterationNumber(cross->eps,&svd->its));
    267:   PetscCall(EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason));
    268:   for (i=0;i<svd->nconv;i++) {
    269:     PetscCall(EPSGetEigenvalue(cross->eps,i,&lambda,NULL));
    270:     sigma = PetscRealPart(lambda);
    271:     if (svd->ishyperbolic) svd->sigma[i] = PetscSqrtReal(PetscAbsReal(sigma));
    272:     else {
    273:       PetscCheck(sigma>-10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_FP,"Negative eigenvalue computed by EPS: %g",(double)sigma);
    274:       if (sigma<0.0) {
    275:         PetscCall(PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma));
    276:         sigma = 0.0;
    277:       }
    278:       svd->sigma[i] = PetscSqrtReal(sigma);
    279:     }
    280:   }
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }

    284: static PetscErrorCode SVDComputeVectors_Cross(SVD svd)
    285: {
    286:   SVD_CROSS         *cross = (SVD_CROSS*)svd->data;
    287:   PetscInt          i,mloc,ploc;
    288:   Vec               u,v,x,uv,w,omega2=NULL;
    289:   Mat               Omega;
    290:   PetscScalar       *dst,alpha,lambda,*varray;
    291:   const PetscScalar *src;
    292:   PetscReal         nrm;

    294:   PetscFunctionBegin;
    295:   if (svd->isgeneralized) {
    296:     PetscCall(MatCreateVecs(svd->A,NULL,&u));
    297:     PetscCall(VecGetLocalSize(u,&mloc));
    298:     PetscCall(MatCreateVecs(svd->B,NULL,&v));
    299:     PetscCall(VecGetLocalSize(v,&ploc));
    300:     for (i=0;i<svd->nconv;i++) {
    301:       PetscCall(BVGetColumn(svd->V,i,&x));
    302:       PetscCall(EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL));
    303:       PetscCall(MatMult(svd->A,x,u));     /* u_i*c_i/alpha = A*x_i */
    304:       PetscCall(VecNormalize(u,NULL));
    305:       PetscCall(MatMult(svd->B,x,v));     /* v_i*s_i/alpha = B*x_i */
    306:       PetscCall(VecNormalize(v,&nrm));    /* ||v||_2 = s_i/alpha   */
    307:       alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm);    /* alpha=s_i/||v||_2 */
    308:       PetscCall(VecScale(x,alpha));
    309:       PetscCall(BVRestoreColumn(svd->V,i,&x));
    310:       /* copy [u;v] to U[i] */
    311:       PetscCall(BVGetColumn(svd->U,i,&uv));
    312:       PetscCall(VecGetArrayWrite(uv,&dst));
    313:       PetscCall(VecGetArrayRead(u,&src));
    314:       PetscCall(PetscArraycpy(dst,src,mloc));
    315:       PetscCall(VecRestoreArrayRead(u,&src));
    316:       PetscCall(VecGetArrayRead(v,&src));
    317:       PetscCall(PetscArraycpy(dst+mloc,src,ploc));
    318:       PetscCall(VecRestoreArrayRead(v,&src));
    319:       PetscCall(VecRestoreArrayWrite(uv,&dst));
    320:       PetscCall(BVRestoreColumn(svd->U,i,&uv));
    321:     }
    322:     PetscCall(VecDestroy(&v));
    323:     PetscCall(VecDestroy(&u));
    324:   } else if (svd->ishyperbolic && svd->swapped) {  /* was solved as GHIEP, set u=Omega*u and normalize */
    325:     PetscCall(EPSGetOperators(cross->eps,NULL,&Omega));
    326:     PetscCall(MatCreateVecs(Omega,&w,NULL));
    327:     PetscCall(VecCreateSeq(PETSC_COMM_SELF,svd->ncv,&omega2));
    328:     PetscCall(VecGetArrayWrite(omega2,&varray));
    329:     for (i=0;i<svd->nconv;i++) {
    330:       PetscCall(BVGetColumn(svd->V,i,&v));
    331:       PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
    332:       PetscCall(MatMult(Omega,v,w));
    333:       PetscCall(VecDot(v,w,&alpha));
    334:       svd->sign[i] = PetscSign(PetscRealPart(alpha));
    335:       varray[i] = svd->sign[i];
    336:       alpha = 1.0/PetscSqrtScalar(PetscAbsScalar(alpha));
    337:       PetscCall(VecScale(w,alpha));
    338:       PetscCall(VecCopy(w,v));
    339:       PetscCall(BVRestoreColumn(svd->V,i,&v));
    340:     }
    341:     PetscCall(BVSetSignature(svd->V,omega2));
    342:     PetscCall(VecRestoreArrayWrite(omega2,&varray));
    343:     PetscCall(VecDestroy(&omega2));
    344:     PetscCall(VecDestroy(&w));
    345:     PetscCall(SVDComputeVectors_Left(svd));
    346:   } else {
    347:     for (i=0;i<svd->nconv;i++) {
    348:       PetscCall(BVGetColumn(svd->V,i,&v));
    349:       PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
    350:       PetscCall(BVRestoreColumn(svd->V,i,&v));
    351:     }
    352:     PetscCall(SVDComputeVectors_Left(svd));
    353:   }
    354:   PetscFunctionReturn(PETSC_SUCCESS);
    355: }

    357: static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
    358: {
    359:   PetscInt       i,ncv;
    360:   SVD            svd = (SVD)ctx;
    361:   SVD_CROSS      *cross;
    362:   PetscScalar    er,ei;
    363:   ST             st;

    365:   PetscFunctionBegin;
    366:   if (svd->stop==SVD_STOP_THRESHOLD) {
    367:     cross = (SVD_CROSS*)svd->data;
    368:     PetscCall(EPSGetDimensions(cross->eps,NULL,&ncv,NULL));
    369:     if (ncv!=svd->ncv) {  /* reallocate */
    370:       PetscCall(SVDReallocateSolution(svd,ncv));
    371:       for (i=svd->ncv;i<ncv;i++) svd->perm[i] = i;
    372:       svd->ncv = ncv;
    373:     }
    374:   }
    375:   PetscCall(EPSGetST(eps,&st));
    376:   for (i=0;i<PetscMin(nest,svd->ncv);i++) {
    377:     er = eigr[i]; ei = eigi[i];
    378:     PetscCall(STBackTransform(st,1,&er,&ei));
    379:     svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
    380:     svd->errest[i] = errest[i];
    381:   }
    382:   PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
    383:   PetscFunctionReturn(PETSC_SUCCESS);
    384: }

    386: static PetscErrorCode SVDSetFromOptions_Cross(SVD svd,PetscOptionItems PetscOptionsObject)
    387: {
    388:   PetscBool      set,val;
    389:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
    390:   ST             st;

    392:   PetscFunctionBegin;
    393:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cross Options");

    395:     PetscCall(PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set));
    396:     if (set) PetscCall(SVDCrossSetExplicitMatrix(svd,val));

    398:   PetscOptionsHeadEnd();

    400:   if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
    401:   if (!cross->explicitmatrix && !cross->usereps) {
    402:     /* use as default an ST with shell matrix and Jacobi */
    403:     PetscCall(EPSGetST(cross->eps,&st));
    404:     PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
    405:   }
    406:   PetscCall(EPSSetFromOptions(cross->eps));
    407:   PetscFunctionReturn(PETSC_SUCCESS);
    408: }

    410: static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
    411: {
    412:   SVD_CROSS *cross = (SVD_CROSS*)svd->data;

    414:   PetscFunctionBegin;
    415:   if (cross->explicitmatrix != explicitmatrix) {
    416:     cross->explicitmatrix = explicitmatrix;
    417:     svd->state = SVD_STATE_INITIAL;
    418:   }
    419:   PetscFunctionReturn(PETSC_SUCCESS);
    420: }

    422: /*@
    423:    SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator $A^*A$ must
    424:    be computed explicitly.

    426:    Logically Collective

    428:    Input Parameters:
    429: +  svd         - the singular value solver context
    430: -  explicitmat - `PETSC_TRUE` if $A^*A$ must be built explicitly

    432:    Options Database Key:
    433: .  -svd_cross_explicitmatrix (true|false) - toggle the explicit construction of the matrix

    435:    Notes:
    436:    In GSVD there are two cross product matrices, $A^*A$ and $B^*B$. In HSVD the
    437:    expression for the cross product matrix is different, $A^*\Omega A$. See
    438:    [](#sec:svdback) for details.

    440:    By default the matrices are not built explicitly, but handled as shell matrices,
    441:    see `MATSHELL`.

    443:    Level: advanced

    445: .seealso: [](ch:svd), [](#sec:svdback), `SVDCROSS`, `SVDCrossGetExplicitMatrix()`, `MATSHELL`
    446: @*/
    447: PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
    448: {
    449:   PetscFunctionBegin;
    452:   PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
    453:   PetscFunctionReturn(PETSC_SUCCESS);
    454: }

    456: static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
    457: {
    458:   SVD_CROSS *cross = (SVD_CROSS*)svd->data;

    460:   PetscFunctionBegin;
    461:   *explicitmat = cross->explicitmatrix;
    462:   PetscFunctionReturn(PETSC_SUCCESS);
    463: }

    465: /*@
    466:    SVDCrossGetExplicitMatrix - Returns the flag indicating if $A^*A$ is built explicitly.

    468:    Not Collective

    470:    Input Parameter:
    471: .  svd  - the singular value solver context

    473:    Output Parameter:
    474: .  explicitmat - the mode flag

    476:    Level: advanced

    478: .seealso: [](ch:svd), `SVDCROSS`, `SVDCrossSetExplicitMatrix()`
    479: @*/
    480: PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
    481: {
    482:   PetscFunctionBegin;
    484:   PetscAssertPointer(explicitmat,2);
    485:   PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
    486:   PetscFunctionReturn(PETSC_SUCCESS);
    487: }

    489: static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
    490: {
    491:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

    493:   PetscFunctionBegin;
    494:   PetscCall(PetscObjectReference((PetscObject)eps));
    495:   PetscCall(EPSDestroy(&cross->eps));
    496:   cross->eps     = eps;
    497:   cross->usereps = PETSC_TRUE;
    498:   svd->state     = SVD_STATE_INITIAL;
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: /*@
    503:    SVDCrossSetEPS - Associate an eigensolver object (`EPS`) to the
    504:    singular value solver.

    506:    Collective

    508:    Input Parameters:
    509: +  svd - the singular value solver context
    510: -  eps - the linear eigensolver context

    512:    Level: advanced

    514: .seealso: [](ch:svd), `SVDCROSS`, `SVDCrossGetEPS()`
    515: @*/
    516: PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
    517: {
    518:   PetscFunctionBegin;
    521:   PetscCheckSameComm(svd,1,eps,2);
    522:   PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
    523:   PetscFunctionReturn(PETSC_SUCCESS);
    524: }

    526: static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
    527: {
    528:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

    530:   PetscFunctionBegin;
    531:   if (!cross->eps) {
    532:     PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps));
    533:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1));
    534:     PetscCall(EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix));
    535:     PetscCall(EPSAppendOptionsPrefix(cross->eps,"svd_cross_"));
    536:     PetscCall(PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options));
    537:     PetscCall(EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_MAGNITUDE));
    538:     PetscCall(EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL));
    539:   }
    540:   *eps = cross->eps;
    541:   PetscFunctionReturn(PETSC_SUCCESS);
    542: }

    544: /*@
    545:    SVDCrossGetEPS - Retrieve the eigensolver object (`EPS`) associated
    546:    to the singular value solver.

    548:    Collective

    550:    Input Parameter:
    551: .  svd - the singular value solver context

    553:    Output Parameter:
    554: .  eps - the linear eigensolver context

    556:    Level: advanced

    558: .seealso: [](ch:svd), `SVDCROSS`, `SVDCrossSetEPS()`
    559: @*/
    560: PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
    561: {
    562:   PetscFunctionBegin;
    564:   PetscAssertPointer(eps,2);
    565:   PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
    566:   PetscFunctionReturn(PETSC_SUCCESS);
    567: }

    569: static PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
    570: {
    571:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
    572:   PetscBool      isascii;

    574:   PetscFunctionBegin;
    575:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    576:   if (isascii) {
    577:     if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
    578:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cross->explicitmatrix?"explicit":"implicit"));
    579:     PetscCall(PetscViewerASCIIPushTab(viewer));
    580:     PetscCall(EPSView(cross->eps,viewer));
    581:     PetscCall(PetscViewerASCIIPopTab(viewer));
    582:   }
    583:   PetscFunctionReturn(PETSC_SUCCESS);
    584: }

    586: static PetscErrorCode SVDReset_Cross(SVD svd)
    587: {
    588:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

    590:   PetscFunctionBegin;
    591:   PetscCall(EPSReset(cross->eps));
    592:   PetscCall(MatDestroy(&cross->C));
    593:   PetscCall(MatDestroy(&cross->D));
    594:   PetscFunctionReturn(PETSC_SUCCESS);
    595: }

    597: static PetscErrorCode SVDDestroy_Cross(SVD svd)
    598: {
    599:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

    601:   PetscFunctionBegin;
    602:   PetscCall(EPSDestroy(&cross->eps));
    603:   PetscCall(PetscFree(svd->data));
    604:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL));
    605:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL));
    606:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL));
    607:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL));
    608:   PetscFunctionReturn(PETSC_SUCCESS);
    609: }

    611: /*MC
    612:    SVDCROSS - SVDCROSS = "cross" - Solve the singular value problem
    613:    via an equivalent eigenvalue problem with the cross product matrix.

    615:    Notes:
    616:    This will do the computation with a subsidiary eigensolver on an
    617:    equivalent eigenvalue problem. For the standard SVD, the eigensolver
    618:    operates with the cross product matrix $A^*A$. See the section
    619:    [](#sec:svdback) for details on the formulation for each SVD type.

    621:    This is the default solver. In the case of standard SVD, the computation
    622:    done by this solver will be almost identical to the one with `SVDTRLANCZOS`.

    624:    To manipulate the internal eigensolvers, use `SVDCrossGetEPS()` or
    625:    use the corresponding command-line options.

    627:    Level: beginner

    629: .seealso: [](ch:svd), [](#sec:svdback), `SVD`, `SVDType`, `SVDSetType()`, `SVDSetProblemType()`, `SVDTRLANCZOS`, `SVDCrossGetEPS()`
    630: M*/
    631: SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
    632: {
    633:   SVD_CROSS      *cross;

    635:   PetscFunctionBegin;
    636:   PetscCall(PetscNew(&cross));
    637:   svd->data = (void*)cross;

    639:   svd->ops->solve          = SVDSolve_Cross;
    640:   svd->ops->solveg         = SVDSolve_Cross;
    641:   svd->ops->solveh         = SVDSolve_Cross;
    642:   svd->ops->setup          = SVDSetUp_Cross;
    643:   svd->ops->setfromoptions = SVDSetFromOptions_Cross;
    644:   svd->ops->destroy        = SVDDestroy_Cross;
    645:   svd->ops->reset          = SVDReset_Cross;
    646:   svd->ops->view           = SVDView_Cross;
    647:   svd->ops->computevectors = SVDComputeVectors_Cross;
    648:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross));
    649:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross));
    650:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross));
    651:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross));
    652:   PetscFunctionReturn(PETSC_SUCCESS);
    653: }