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: "lanczos"
      12             : 
      13             :    Method: Explicitly restarted Lanczos
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Golub-Kahan-Lanczos bidiagonalization with explicit restart.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] G.H. Golub and W. Kahan, "Calculating the singular values
      22             :            and pseudo-inverse of a matrix", SIAM J. Numer. Anal. Ser.
      23             :            B 2:205-224, 1965.
      24             : 
      25             :        [2] V. Hernandez, J.E. Roman, and A. Tomas, "A robust and
      26             :            efficient parallel SVD solver based on restarted Lanczos
      27             :            bidiagonalization", Elec. Trans. Numer. Anal. 31:68-85,
      28             :            2008.
      29             : */
      30             : 
      31             : #include <slepc/private/svdimpl.h>                /*I "slepcsvd.h" I*/
      32             : 
      33             : typedef struct {
      34             :   PetscBool oneside;
      35             : } SVD_LANCZOS;
      36             : 
      37          21 : static PetscErrorCode SVDSetUp_Lanczos(SVD svd)
      38             : {
      39          21 :   SVD_LANCZOS    *lanczos = (SVD_LANCZOS*)svd->data;
      40          21 :   PetscInt       N;
      41             : 
      42          21 :   PetscFunctionBegin;
      43          21 :   SVDCheckStandard(svd);
      44          21 :   SVDCheckDefinite(svd);
      45          21 :   PetscCall(MatGetSize(svd->A,NULL,&N));
      46          21 :   PetscCall(SVDSetDimensions_Default(svd));
      47          21 :   PetscCheck(svd->ncv<=svd->nsv+svd->mpd,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
      48          30 :   if (svd->max_it==PETSC_DETERMINE) svd->max_it = PetscMax(N/svd->ncv,100)*((svd->stop==SVD_STOP_THRESHOLD)?10:1);
      49          21 :   svd->leftbasis = PetscNot(lanczos->oneside);
      50          21 :   PetscCall(SVDAllocateSolution(svd,1));
      51          21 :   PetscCall(DSSetType(svd->ds,DSSVD));
      52          21 :   PetscCall(DSSetCompact(svd->ds,PETSC_TRUE));
      53          21 :   PetscCall(DSSetExtraRow(svd->ds,PETSC_TRUE));
      54          21 :   PetscCall(DSAllocate(svd->ds,svd->ncv+1));
      55          21 :   PetscFunctionReturn(PETSC_SUCCESS);
      56             : }
      57             : 
      58        1082 : PetscErrorCode SVDTwoSideLanczos(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt k,PetscInt *n,PetscBool *breakdown)
      59             : {
      60        1082 :   PetscInt       i;
      61        1082 :   Vec            u,v;
      62        1082 :   PetscBool      lindep=PETSC_FALSE;
      63             : 
      64        1082 :   PetscFunctionBegin;
      65        1082 :   PetscCall(BVGetColumn(svd->V,k,&v));
      66        1082 :   PetscCall(BVGetColumn(svd->U,k,&u));
      67        1082 :   PetscCall(MatMult(svd->A,v,u));
      68        1082 :   PetscCall(BVRestoreColumn(svd->V,k,&v));
      69        1082 :   PetscCall(BVRestoreColumn(svd->U,k,&u));
      70        1082 :   PetscCall(BVOrthonormalizeColumn(svd->U,k,PETSC_FALSE,alpha+k,&lindep));
      71        1082 :   if (PetscUnlikely(lindep)) {
      72           0 :     *n = k;
      73           0 :     if (breakdown) *breakdown = lindep;
      74           0 :     PetscFunctionReturn(PETSC_SUCCESS);
      75             :   }
      76             : 
      77        8370 :   for (i=k+1;i<*n;i++) {
      78        7288 :     PetscCall(BVGetColumn(svd->V,i,&v));
      79        7288 :     PetscCall(BVGetColumn(svd->U,i-1,&u));
      80        7288 :     PetscCall(MatMult(svd->AT,u,v));
      81        7288 :     PetscCall(BVRestoreColumn(svd->V,i,&v));
      82        7288 :     PetscCall(BVRestoreColumn(svd->U,i-1,&u));
      83        7288 :     PetscCall(BVOrthonormalizeColumn(svd->V,i,PETSC_FALSE,beta+i-1,&lindep));
      84        7288 :     if (PetscUnlikely(lindep)) {
      85           0 :       *n = i;
      86           0 :       break;
      87             :     }
      88        7288 :     PetscCall(BVGetColumn(svd->V,i,&v));
      89        7288 :     PetscCall(BVGetColumn(svd->U,i,&u));
      90        7288 :     PetscCall(MatMult(svd->A,v,u));
      91        7288 :     PetscCall(BVRestoreColumn(svd->V,i,&v));
      92        7288 :     PetscCall(BVRestoreColumn(svd->U,i,&u));
      93        7288 :     PetscCall(BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,alpha+i,&lindep));
      94        7288 :     if (PetscUnlikely(lindep)) {
      95           0 :       *n = i;
      96           0 :       break;
      97             :     }
      98             :   }
      99             : 
     100        1082 :   if (!lindep) {
     101        1082 :     PetscCall(BVGetColumn(svd->V,*n,&v));
     102        1082 :     PetscCall(BVGetColumn(svd->U,*n-1,&u));
     103        1082 :     PetscCall(MatMult(svd->AT,u,v));
     104        1082 :     PetscCall(BVRestoreColumn(svd->V,*n,&v));
     105        1082 :     PetscCall(BVRestoreColumn(svd->U,*n-1,&u));
     106        1082 :     PetscCall(BVOrthogonalizeColumn(svd->V,*n,NULL,beta+*n-1,&lindep));
     107             :   }
     108        1082 :   if (breakdown) *breakdown = lindep;
     109        1082 :   PetscFunctionReturn(PETSC_SUCCESS);
     110             : }
     111             : 
     112          54 : static PetscErrorCode SVDOneSideLanczos(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,Vec u,Vec u_1,PetscInt k,PetscInt n,PetscScalar* work)
     113             : {
     114          54 :   PetscInt       i,bvl,bvk;
     115          54 :   PetscReal      a,b;
     116          54 :   Vec            z,temp;
     117             : 
     118          54 :   PetscFunctionBegin;
     119          54 :   PetscCall(BVGetActiveColumns(V,&bvl,&bvk));
     120          54 :   PetscCall(BVGetColumn(V,k,&z));
     121          54 :   PetscCall(MatMult(svd->A,z,u));
     122          54 :   PetscCall(BVRestoreColumn(V,k,&z));
     123             : 
     124         468 :   for (i=k+1;i<n;i++) {
     125         414 :     PetscCall(BVGetColumn(V,i,&z));
     126         414 :     PetscCall(MatMult(svd->AT,u,z));
     127         414 :     PetscCall(BVRestoreColumn(V,i,&z));
     128         414 :     PetscCall(VecNormBegin(u,NORM_2,&a));
     129         414 :     PetscCall(BVSetActiveColumns(V,0,i));
     130         414 :     PetscCall(BVDotColumnBegin(V,i,work));
     131         414 :     PetscCall(VecNormEnd(u,NORM_2,&a));
     132         414 :     PetscCall(BVDotColumnEnd(V,i,work));
     133         414 :     PetscCall(VecScale(u,1.0/a));
     134         414 :     PetscCall(BVMultColumn(V,-1.0/a,1.0/a,i,work));
     135             : 
     136             :     /* h = V^* z, z = z - V h  */
     137         414 :     PetscCall(BVDotColumn(V,i,work));
     138         414 :     PetscCall(BVMultColumn(V,-1.0,1.0,i,work));
     139         414 :     PetscCall(BVNormColumn(V,i,NORM_2,&b));
     140         414 :     PetscCheck(PetscAbsReal(b)>10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_PLIB,"Recurrence generated a zero vector; use a two-sided variant");
     141         414 :     PetscCall(BVScaleColumn(V,i,1.0/b));
     142             : 
     143         414 :     PetscCall(BVGetColumn(V,i,&z));
     144         414 :     PetscCall(MatMult(svd->A,z,u_1));
     145         414 :     PetscCall(BVRestoreColumn(V,i,&z));
     146         414 :     PetscCall(VecAXPY(u_1,-b,u));
     147         414 :     alpha[i-1] = a;
     148         414 :     beta[i-1] = b;
     149         414 :     temp = u;
     150         414 :     u = u_1;
     151         414 :     u_1 = temp;
     152             :   }
     153             : 
     154          54 :   PetscCall(BVGetColumn(V,n,&z));
     155          54 :   PetscCall(MatMult(svd->AT,u,z));
     156          54 :   PetscCall(BVRestoreColumn(V,n,&z));
     157          54 :   PetscCall(VecNormBegin(u,NORM_2,&a));
     158          54 :   PetscCall(BVDotColumnBegin(V,n,work));
     159          54 :   PetscCall(VecNormEnd(u,NORM_2,&a));
     160          54 :   PetscCall(BVDotColumnEnd(V,n,work));
     161          54 :   PetscCall(VecScale(u,1.0/a));
     162          54 :   PetscCall(BVMultColumn(V,-1.0/a,1.0/a,n,work));
     163             : 
     164             :   /* h = V^* z, z = z - V h  */
     165          54 :   PetscCall(BVDotColumn(V,n,work));
     166          54 :   PetscCall(BVMultColumn(V,-1.0,1.0,n,work));
     167          54 :   PetscCall(BVNormColumn(V,i,NORM_2,&b));
     168             : 
     169          54 :   alpha[n-1] = a;
     170          54 :   beta[n-1] = b;
     171          54 :   PetscCall(BVSetActiveColumns(V,bvl,bvk));
     172          54 :   PetscFunctionReturn(PETSC_SUCCESS);
     173             : }
     174             : 
     175             : /*
     176             :    SVDKrylovConvergence - Implements the loop that checks for convergence
     177             :    in Krylov methods.
     178             : 
     179             :    Input Parameters:
     180             :      svd     - the solver; some error estimates are updated in svd->errest
     181             :      getall  - whether all residuals must be computed
     182             :      kini    - initial value of k (the loop variable)
     183             :      nits    - number of iterations of the loop
     184             :      normr   - norm of triangular factor of qr([A;B]), used only in GSVD
     185             : 
     186             :    Output Parameter:
     187             :      kout  - the first index where the convergence test failed
     188             : */
     189        2580 : PetscErrorCode SVDKrylovConvergence(SVD svd,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal normr,PetscInt *kout)
     190             : {
     191        2580 :   PetscInt       k,marker,ld;
     192        2580 :   PetscReal      *alpha,*beta,*betah,resnorm;
     193        2580 :   PetscBool      extra;
     194             : 
     195        2580 :   PetscFunctionBegin;
     196        2580 :   if (PetscUnlikely(svd->conv == SVD_CONV_MAXIT && svd->its >= svd->max_it)) *kout = svd->nsv;
     197             :   else {
     198        2580 :     PetscCall(DSGetLeadingDimension(svd->ds,&ld));
     199        2580 :     PetscCall(DSGetExtraRow(svd->ds,&extra));
     200        2580 :     PetscCheck(extra,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Only implemented for DS with extra row");
     201        2580 :     marker = -1;
     202        2580 :     if (svd->trackall) getall = PETSC_TRUE;
     203        2580 :     PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
     204        2580 :     beta = alpha + ld;
     205        2580 :     betah = alpha + 2*ld;
     206        3401 :     for (k=kini;k<kini+nits;k++) {
     207        3396 :       if (svd->isgeneralized) resnorm = SlepcAbs(beta[k],betah[k])*normr;
     208        2156 :       else resnorm = PetscAbsReal(beta[k]);
     209        3396 :       PetscCall((*svd->converged)(svd,svd->sigma[k],resnorm,&svd->errest[k],svd->convergedctx));
     210        3396 :       if (marker==-1 && svd->errest[k] >= svd->tol) marker = k;
     211        3396 :       if (marker!=-1 && !getall) break;
     212             :     }
     213        2580 :     PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
     214        2580 :     if (marker!=-1) k = marker;
     215        2580 :     *kout = k;
     216             :   }
     217        2580 :   PetscFunctionReturn(PETSC_SUCCESS);
     218             : }
     219             : 
     220          21 : static PetscErrorCode SVDSolve_Lanczos(SVD svd)
     221             : {
     222          21 :   SVD_LANCZOS    *lanczos = (SVD_LANCZOS*)svd->data;
     223          21 :   PetscReal      *alpha,*beta;
     224          21 :   PetscScalar    *swork,*w,*P,*aux1,*aux2;
     225          21 :   PetscInt       i,k,j,nv,ld;
     226          21 :   Vec            u=NULL,u_1=NULL;
     227          21 :   Mat            U,V;
     228             : 
     229          21 :   PetscFunctionBegin;
     230             :   /* allocate working space */
     231          21 :   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
     232          21 :   PetscCall(PetscMalloc2(ld,&w,svd->ncv,&swork));
     233             : 
     234          21 :   if (lanczos->oneside) {
     235           2 :     PetscCall(MatCreateVecs(svd->A,NULL,&u));
     236           2 :     PetscCall(MatCreateVecs(svd->A,NULL,&u_1));
     237             :   }
     238             : 
     239             :   /* normalize start vector */
     240          21 :   if (!svd->nini) {
     241          17 :     PetscCall(BVSetRandomColumn(svd->V,0));
     242          17 :     PetscCall(BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL));
     243             :   }
     244             : 
     245         683 :   while (svd->reason == SVD_CONVERGED_ITERATING) {
     246         662 :     svd->its++;
     247             : 
     248             :     /* inner loop */
     249         662 :     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
     250         662 :     PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
     251         662 :     beta = alpha + ld;
     252         662 :     if (lanczos->oneside) PetscCall(SVDOneSideLanczos(svd,alpha,beta,svd->V,u,u_1,svd->nconv,nv,swork));
     253             :     else {
     254         608 :       PetscCall(SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv,&nv,NULL));
     255         608 :       PetscCall(BVSetActiveColumns(svd->U,svd->nconv,nv));
     256             :     }
     257         662 :     PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
     258         662 :     PetscCall(BVSetActiveColumns(svd->V,svd->nconv,nv));
     259             : 
     260             :     /* compute SVD of bidiagonal matrix */
     261         662 :     PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,0));
     262         662 :     PetscCall(DSSVDSetDimensions(svd->ds,nv));
     263         662 :     PetscCall(DSSetState(svd->ds,DS_STATE_INTERMEDIATE));
     264         662 :     PetscCall(DSSolve(svd->ds,w,NULL));
     265         662 :     PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
     266         662 :     PetscCall(DSUpdateExtraRow(svd->ds));
     267         662 :     PetscCall(DSSynchronize(svd->ds,w,NULL));
     268        6965 :     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
     269             : 
     270             :     /* check convergence */
     271         662 :     PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
     272         662 :     SVDSetCtxThreshold(svd,svd->sigma,k);
     273         662 :     PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
     274             : 
     275             :     /* compute restart vector */
     276         662 :     if (svd->reason == SVD_CONVERGED_ITERATING) {
     277         641 :       if (k<nv) {
     278         641 :         PetscCall(DSGetArray(svd->ds,DS_MAT_V,&P));
     279        6749 :         for (j=svd->nconv;j<nv;j++) swork[j-svd->nconv] = PetscConj(P[j+k*ld]);
     280         641 :         PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&P));
     281         641 :         PetscCall(BVMultColumn(svd->V,1.0,0.0,nv,swork));
     282             :       } else {
     283             :         /* all approximations have converged, generate a new initial vector */
     284           0 :         PetscCall(BVSetRandomColumn(svd->V,nv));
     285           0 :         PetscCall(BVOrthonormalizeColumn(svd->V,nv,PETSC_FALSE,NULL,NULL));
     286             :       }
     287             :     }
     288             : 
     289             :     /* compute converged singular vectors */
     290         662 :     PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
     291         662 :     PetscCall(BVMultInPlace(svd->V,V,svd->nconv,k));
     292         662 :     PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
     293         662 :     if (!lanczos->oneside) {
     294         608 :       PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
     295         608 :       PetscCall(BVMultInPlace(svd->U,U,svd->nconv,k));
     296         608 :       PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
     297             :     }
     298             : 
     299         662 :     if (svd->reason == SVD_CONVERGED_ITERATING) {
     300         641 :       PetscCall(BVCopyColumn(svd->V,nv,k));  /* copy restart vector from the last column */
     301         641 :       if (svd->stop==SVD_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
     302           1 :         svd->ncv = svd->mpd+k;
     303           1 :         PetscCall(SVDReallocateSolution(svd,svd->ncv+1));
     304           9 :         for (i=nv;i<svd->ncv;i++) svd->perm[i] = i;
     305           1 :         PetscCall(DSReallocate(svd->ds,svd->ncv+1));
     306           1 :         aux1 = w;
     307           1 :         aux2 = swork;
     308           1 :         PetscCall(PetscMalloc2(svd->ncv+1,&w,svd->ncv,&swork));
     309           1 :         PetscCall(PetscArraycpy(w,aux1,ld));
     310           1 :         PetscCall(PetscFree2(aux1,aux2));
     311           1 :         PetscCall(DSGetLeadingDimension(svd->ds,&ld));
     312             :       }
     313             :     }
     314             : 
     315         662 :     svd->nconv = k;
     316         683 :     PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
     317             :   }
     318             : 
     319             :   /* free working space */
     320          21 :   PetscCall(VecDestroy(&u));
     321          21 :   PetscCall(VecDestroy(&u_1));
     322          21 :   PetscCall(PetscFree2(w,swork));
     323          21 :   PetscFunctionReturn(PETSC_SUCCESS);
     324             : }
     325             : 
     326          15 : static PetscErrorCode SVDSetFromOptions_Lanczos(SVD svd,PetscOptionItems *PetscOptionsObject)
     327             : {
     328          15 :   PetscBool      set,val;
     329          15 :   SVD_LANCZOS    *lanczos = (SVD_LANCZOS*)svd->data;
     330             : 
     331          15 :   PetscFunctionBegin;
     332          15 :   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Lanczos Options");
     333             : 
     334          15 :     PetscCall(PetscOptionsBool("-svd_lanczos_oneside","Use one-side reorthogonalization","SVDLanczosSetOneSide",lanczos->oneside,&val,&set));
     335          15 :     if (set) PetscCall(SVDLanczosSetOneSide(svd,val));
     336             : 
     337          15 :   PetscOptionsHeadEnd();
     338          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     339             : }
     340             : 
     341           2 : static PetscErrorCode SVDLanczosSetOneSide_Lanczos(SVD svd,PetscBool oneside)
     342             : {
     343           2 :   SVD_LANCZOS *lanczos = (SVD_LANCZOS*)svd->data;
     344             : 
     345           2 :   PetscFunctionBegin;
     346           2 :   if (lanczos->oneside != oneside) {
     347           2 :     lanczos->oneside = oneside;
     348           2 :     svd->state = SVD_STATE_INITIAL;
     349             :   }
     350           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     351             : }
     352             : 
     353             : /*@
     354             :    SVDLanczosSetOneSide - Indicate if the variant of the Lanczos method
     355             :    to be used is one-sided or two-sided.
     356             : 
     357             :    Logically Collective
     358             : 
     359             :    Input Parameters:
     360             : +  svd     - singular value solver
     361             : -  oneside - boolean flag indicating if the method is one-sided or not
     362             : 
     363             :    Options Database Key:
     364             : .  -svd_lanczos_oneside <boolean> - Indicates the boolean flag
     365             : 
     366             :    Note:
     367             :    By default, a two-sided variant is selected, which is sometimes slightly
     368             :    more robust. However, the one-sided variant is faster because it avoids
     369             :    the orthogonalization associated to left singular vectors. It also saves
     370             :    the memory required for storing such vectors.
     371             : 
     372             :    Level: advanced
     373             : 
     374             : .seealso: SVDTRLanczosSetOneSide()
     375             : @*/
     376           2 : PetscErrorCode SVDLanczosSetOneSide(SVD svd,PetscBool oneside)
     377             : {
     378           2 :   PetscFunctionBegin;
     379           2 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     380           6 :   PetscValidLogicalCollectiveBool(svd,oneside,2);
     381           2 :   PetscTryMethod(svd,"SVDLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
     382           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     383             : }
     384             : 
     385           2 : static PetscErrorCode SVDLanczosGetOneSide_Lanczos(SVD svd,PetscBool *oneside)
     386             : {
     387           2 :   SVD_LANCZOS *lanczos = (SVD_LANCZOS*)svd->data;
     388             : 
     389           2 :   PetscFunctionBegin;
     390           2 :   *oneside = lanczos->oneside;
     391           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     392             : }
     393             : 
     394             : /*@
     395             :    SVDLanczosGetOneSide - Gets if the variant of the Lanczos method
     396             :    to be used is one-sided or two-sided.
     397             : 
     398             :    Not Collective
     399             : 
     400             :    Input Parameters:
     401             : .  svd     - singular value solver
     402             : 
     403             :    Output Parameters:
     404             : .  oneside - boolean flag indicating if the method is one-sided or not
     405             : 
     406             :    Level: advanced
     407             : 
     408             : .seealso: SVDLanczosSetOneSide()
     409             : @*/
     410           2 : PetscErrorCode SVDLanczosGetOneSide(SVD svd,PetscBool *oneside)
     411             : {
     412           2 :   PetscFunctionBegin;
     413           2 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     414           2 :   PetscAssertPointer(oneside,2);
     415           2 :   PetscUseMethod(svd,"SVDLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
     416           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     417             : }
     418             : 
     419          17 : static PetscErrorCode SVDDestroy_Lanczos(SVD svd)
     420             : {
     421          17 :   PetscFunctionBegin;
     422          17 :   PetscCall(PetscFree(svd->data));
     423          17 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDLanczosSetOneSide_C",NULL));
     424          17 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDLanczosGetOneSide_C",NULL));
     425          17 :   PetscFunctionReturn(PETSC_SUCCESS);
     426             : }
     427             : 
     428           0 : static PetscErrorCode SVDView_Lanczos(SVD svd,PetscViewer viewer)
     429             : {
     430           0 :   SVD_LANCZOS    *lanczos = (SVD_LANCZOS*)svd->data;
     431           0 :   PetscBool      isascii;
     432             : 
     433           0 :   PetscFunctionBegin;
     434           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     435           0 :   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer,"  %s-sided reorthogonalization\n",lanczos->oneside? "one": "two"));
     436           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     437             : }
     438             : 
     439          17 : SLEPC_EXTERN PetscErrorCode SVDCreate_Lanczos(SVD svd)
     440             : {
     441          17 :   SVD_LANCZOS    *ctx;
     442             : 
     443          17 :   PetscFunctionBegin;
     444          17 :   PetscCall(PetscNew(&ctx));
     445          17 :   svd->data = (void*)ctx;
     446             : 
     447          17 :   svd->ops->setup          = SVDSetUp_Lanczos;
     448          17 :   svd->ops->solve          = SVDSolve_Lanczos;
     449          17 :   svd->ops->destroy        = SVDDestroy_Lanczos;
     450          17 :   svd->ops->setfromoptions = SVDSetFromOptions_Lanczos;
     451          17 :   svd->ops->view           = SVDView_Lanczos;
     452          17 :   svd->ops->computevectors = SVDComputeVectors_Left;
     453          17 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDLanczosSetOneSide_C",SVDLanczosSetOneSide_Lanczos));
     454          17 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDLanczosGetOneSide_C",SVDLanczosGetOneSide_Lanczos));
     455          17 :   PetscFunctionReturn(PETSC_SUCCESS);
     456             : }