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: "cyclic"
      12             : 
      13             :    Method: Uses a Hermitian eigensolver for H(A) = [ 0  A ; A^T 0 ]
      14             : */
      15             : 
      16             : #include <slepc/private/svdimpl.h>                /*I "slepcsvd.h" I*/
      17             : #include <slepc/private/bvimpl.h>
      18             : #include "cyclic.h"
      19             : 
      20       43800 : static PetscErrorCode MatMult_Cyclic(Mat B,Vec x,Vec y)
      21             : {
      22       43800 :   SVD_CYCLIC_SHELL  *ctx;
      23       43800 :   const PetscScalar *px;
      24       43800 :   PetscScalar       *py;
      25       43800 :   PetscInt          m;
      26             : 
      27       43800 :   PetscFunctionBegin;
      28       43800 :   PetscCall(MatShellGetContext(B,&ctx));
      29       43800 :   PetscCall(MatGetLocalSize(ctx->A,&m,NULL));
      30       43800 :   PetscCall(VecGetArrayRead(x,&px));
      31       43800 :   PetscCall(VecGetArrayWrite(y,&py));
      32       43800 :   PetscCall(VecPlaceArray(ctx->x1,px));
      33       43800 :   PetscCall(VecPlaceArray(ctx->x2,px+m));
      34       43800 :   PetscCall(VecPlaceArray(ctx->y1,py));
      35       43800 :   PetscCall(VecPlaceArray(ctx->y2,py+m));
      36       43800 :   PetscCall(MatMult(ctx->A,ctx->x2,ctx->y1));
      37       43800 :   PetscCall(MatMult(ctx->AT,ctx->x1,ctx->y2));
      38       43800 :   PetscCall(VecResetArray(ctx->x1));
      39       43800 :   PetscCall(VecResetArray(ctx->x2));
      40       43800 :   PetscCall(VecResetArray(ctx->y1));
      41       43800 :   PetscCall(VecResetArray(ctx->y2));
      42       43800 :   PetscCall(VecRestoreArrayRead(x,&px));
      43       43800 :   PetscCall(VecRestoreArrayWrite(y,&py));
      44       43800 :   PetscFunctionReturn(PETSC_SUCCESS);
      45             : }
      46             : 
      47           0 : static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
      48             : {
      49           0 :   PetscFunctionBegin;
      50           0 :   PetscCall(VecSet(diag,0.0));
      51           0 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54          29 : static PetscErrorCode MatDestroy_Cyclic(Mat B)
      55             : {
      56          29 :   SVD_CYCLIC_SHELL *ctx;
      57             : 
      58          29 :   PetscFunctionBegin;
      59          29 :   PetscCall(MatShellGetContext(B,&ctx));
      60          29 :   PetscCall(VecDestroy(&ctx->x1));
      61          29 :   PetscCall(VecDestroy(&ctx->x2));
      62          29 :   PetscCall(VecDestroy(&ctx->y1));
      63          29 :   PetscCall(VecDestroy(&ctx->y2));
      64          29 :   if (ctx->misaligned) {
      65           0 :     PetscCall(VecDestroy(&ctx->wx2));
      66           0 :     PetscCall(VecDestroy(&ctx->wy2));
      67             :   }
      68          29 :   PetscCall(PetscFree(ctx));
      69          29 :   PetscFunctionReturn(PETSC_SUCCESS);
      70             : }
      71             : 
      72             : /*
      73             :    Builds cyclic matrix   C = | 0   A |
      74             :                               | AT  0 |
      75             : */
      76          46 : static PetscErrorCode SVDCyclicGetCyclicMat(SVD svd,Mat A,Mat AT,Mat *C)
      77             : {
      78          46 :   SVD_CYCLIC       *cyclic = (SVD_CYCLIC*)svd->data;
      79          46 :   SVD_CYCLIC_SHELL *ctx;
      80          46 :   PetscInt         i,M,N,m,n,Istart,Iend;
      81          46 :   VecType          vtype;
      82          46 :   Mat              Zm,Zn;
      83             : #if defined(PETSC_HAVE_CUDA)
      84             :   PetscBool        cuda;
      85             :   const PetscInt   *ranges;
      86             :   PetscMPIInt      size;
      87             : #endif
      88             : 
      89          46 :   PetscFunctionBegin;
      90          46 :   PetscCall(MatGetSize(A,&M,&N));
      91          46 :   PetscCall(MatGetLocalSize(A,&m,&n));
      92             : 
      93          46 :   if (cyclic->explicitmatrix) {
      94          17 :     PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
      95          17 :     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
      96          17 :     PetscCall(MatSetSizes(Zm,m,m,M,M));
      97          17 :     PetscCall(MatSetFromOptions(Zm));
      98          17 :     PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
      99        2794 :     for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
     100          17 :     PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
     101          17 :     PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
     102          17 :     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
     103          17 :     PetscCall(MatSetSizes(Zn,n,n,N,N));
     104          17 :     PetscCall(MatSetFromOptions(Zn));
     105          17 :     PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
     106        1360 :     for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
     107          17 :     PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
     108          17 :     PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
     109          17 :     PetscCall(MatCreateTile(1.0,Zm,1.0,A,1.0,AT,1.0,Zn,C));
     110          17 :     PetscCall(MatDestroy(&Zm));
     111          17 :     PetscCall(MatDestroy(&Zn));
     112             :   } else {
     113          29 :     PetscCall(PetscNew(&ctx));
     114          29 :     ctx->A       = A;
     115          29 :     ctx->AT      = AT;
     116          29 :     ctx->swapped = svd->swapped;
     117          29 :     PetscCall(MatCreateVecsEmpty(A,&ctx->x2,&ctx->x1));
     118          29 :     PetscCall(MatCreateVecsEmpty(A,&ctx->y2,&ctx->y1));
     119          29 :     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
     120          29 :     PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic));
     121          29 :     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic));
     122             : #if defined(PETSC_HAVE_CUDA)
     123             :     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
     124             :     if (cuda) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA));
     125             :     else
     126             : #endif
     127          29 :       PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic));
     128          29 :     PetscCall(MatGetVecType(A,&vtype));
     129          29 :     PetscCall(MatSetVecType(*C,vtype));
     130             : #if defined(PETSC_HAVE_CUDA)
     131             :     if (cuda) {
     132             :       /* check alignment of bottom block */
     133             :       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
     134             :       PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
     135             :       for (i=0;i<size;i++) {
     136             :         ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
     137             :         if (ctx->misaligned) break;
     138             :       }
     139             :       if (ctx->misaligned) {  /* create work vectors for MatMult */
     140             :         PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
     141             :         PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
     142             :       }
     143             :     }
     144             : #endif
     145             :   }
     146          46 :   PetscFunctionReturn(PETSC_SUCCESS);
     147             : }
     148             : 
     149        7649 : static PetscErrorCode MatMult_ECross(Mat B,Vec x,Vec y)
     150             : {
     151        7649 :   SVD_CYCLIC_SHELL  *ctx;
     152        7649 :   const PetscScalar *px;
     153        7649 :   PetscScalar       *py;
     154        7649 :   PetscInt          mn,m,n;
     155             : 
     156        7649 :   PetscFunctionBegin;
     157        7649 :   PetscCall(MatShellGetContext(B,&ctx));
     158        7649 :   PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
     159        7649 :   PetscCall(VecGetLocalSize(y,&mn));
     160        7649 :   m = mn-n;
     161        7649 :   PetscCall(VecGetArrayRead(x,&px));
     162        7649 :   PetscCall(VecGetArrayWrite(y,&py));
     163        7649 :   PetscCall(VecPlaceArray(ctx->x1,px));
     164        7649 :   PetscCall(VecPlaceArray(ctx->x2,px+m));
     165        7649 :   PetscCall(VecPlaceArray(ctx->y1,py));
     166        7649 :   PetscCall(VecPlaceArray(ctx->y2,py+m));
     167        7649 :   PetscCall(VecCopy(ctx->x1,ctx->y1));
     168        7649 :   PetscCall(MatMult(ctx->A,ctx->x2,ctx->w));
     169        7649 :   PetscCall(MatMult(ctx->AT,ctx->w,ctx->y2));
     170        7649 :   PetscCall(VecResetArray(ctx->x1));
     171        7649 :   PetscCall(VecResetArray(ctx->x2));
     172        7649 :   PetscCall(VecResetArray(ctx->y1));
     173        7649 :   PetscCall(VecResetArray(ctx->y2));
     174        7649 :   PetscCall(VecRestoreArrayRead(x,&px));
     175        7649 :   PetscCall(VecRestoreArrayWrite(y,&py));
     176        7649 :   PetscFunctionReturn(PETSC_SUCCESS);
     177             : }
     178             : 
     179           6 : static PetscErrorCode MatGetDiagonal_ECross(Mat B,Vec d)
     180             : {
     181           6 :   SVD_CYCLIC_SHELL  *ctx;
     182           6 :   PetscScalar       *pd;
     183           6 :   PetscMPIInt       len;
     184           6 :   PetscInt          mn,m,n,N,i,j,start,end,ncols;
     185           6 :   PetscScalar       *work1,*work2,*diag;
     186           6 :   const PetscInt    *cols;
     187           6 :   const PetscScalar *vals;
     188             : 
     189           6 :   PetscFunctionBegin;
     190           6 :   PetscCall(MatShellGetContext(B,&ctx));
     191           6 :   PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
     192           6 :   PetscCall(VecGetLocalSize(d,&mn));
     193           6 :   m = mn-n;
     194           6 :   PetscCall(VecGetArrayWrite(d,&pd));
     195           6 :   PetscCall(VecPlaceArray(ctx->y1,pd));
     196           6 :   PetscCall(VecSet(ctx->y1,1.0));
     197           6 :   PetscCall(VecResetArray(ctx->y1));
     198           6 :   PetscCall(VecPlaceArray(ctx->y2,pd+m));
     199           6 :   if (!ctx->diag) {
     200             :     /* compute diagonal from rows and store in ctx->diag */
     201           6 :     PetscCall(VecDuplicate(ctx->y2,&ctx->diag));
     202           6 :     PetscCall(MatGetSize(ctx->A,NULL,&N));
     203           6 :     PetscCall(PetscCalloc2(N,&work1,N,&work2));
     204           6 :     if (ctx->swapped) {
     205           0 :       PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
     206           0 :       for (i=start;i<end;i++) {
     207           0 :         PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
     208           0 :         for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
     209           0 :         PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
     210             :       }
     211             :     } else {
     212           6 :       PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
     213         137 :       for (i=start;i<end;i++) {
     214         131 :         PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
     215         894 :         for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
     216         131 :         PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
     217             :       }
     218             :     }
     219           6 :     PetscCall(PetscMPIIntCast(N,&len));
     220          24 :     PetscCall(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
     221           6 :     PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
     222           6 :     PetscCall(VecGetArrayWrite(ctx->diag,&diag));
     223         121 :     for (i=start;i<end;i++) diag[i-start] = work2[i];
     224           6 :     PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
     225           6 :     PetscCall(PetscFree2(work1,work2));
     226             :   }
     227           6 :   PetscCall(VecCopy(ctx->diag,ctx->y2));
     228           6 :   PetscCall(VecResetArray(ctx->y2));
     229           6 :   PetscCall(VecRestoreArrayWrite(d,&pd));
     230           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     231             : }
     232             : 
     233           6 : static PetscErrorCode MatDestroy_ECross(Mat B)
     234             : {
     235           6 :   SVD_CYCLIC_SHELL *ctx;
     236             : 
     237           6 :   PetscFunctionBegin;
     238           6 :   PetscCall(MatShellGetContext(B,&ctx));
     239           6 :   PetscCall(VecDestroy(&ctx->x1));
     240           6 :   PetscCall(VecDestroy(&ctx->x2));
     241           6 :   PetscCall(VecDestroy(&ctx->y1));
     242           6 :   PetscCall(VecDestroy(&ctx->y2));
     243           6 :   PetscCall(VecDestroy(&ctx->diag));
     244           6 :   PetscCall(VecDestroy(&ctx->w));
     245           6 :   if (ctx->misaligned) {
     246           0 :     PetscCall(VecDestroy(&ctx->wx2));
     247           0 :     PetscCall(VecDestroy(&ctx->wy2));
     248             :   }
     249           6 :   PetscCall(PetscFree(ctx));
     250           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     251             : }
     252             : 
     253             : /*
     254             :    Builds extended cross product matrix   C = | I_m   0  |
     255             :                                               |  0  AT*A |
     256             :    t is an auxiliary Vec used to take the dimensions of the upper block
     257             : */
     258          13 : static PetscErrorCode SVDCyclicGetECrossMat(SVD svd,Mat A,Mat AT,Mat *C,Vec t)
     259             : {
     260          13 :   SVD_CYCLIC       *cyclic = (SVD_CYCLIC*)svd->data;
     261          13 :   SVD_CYCLIC_SHELL *ctx;
     262          13 :   PetscInt         i,M,N,m,n,Istart,Iend;
     263          13 :   VecType          vtype;
     264          13 :   Mat              Id,Zm,Zn,ATA;
     265             : #if defined(PETSC_HAVE_CUDA)
     266             :   PetscBool        cuda;
     267             :   const PetscInt   *ranges;
     268             :   PetscMPIInt      size;
     269             : #endif
     270             : 
     271          13 :   PetscFunctionBegin;
     272          13 :   PetscCall(MatGetSize(A,NULL,&N));
     273          13 :   PetscCall(MatGetLocalSize(A,NULL,&n));
     274          13 :   PetscCall(VecGetSize(t,&M));
     275          13 :   PetscCall(VecGetLocalSize(t,&m));
     276             : 
     277          13 :   if (cyclic->explicitmatrix) {
     278           7 :     PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
     279           7 :     PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id));
     280           7 :     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
     281           7 :     PetscCall(MatSetSizes(Zm,m,n,M,N));
     282           7 :     PetscCall(MatSetFromOptions(Zm));
     283           7 :     PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
     284         355 :     for (i=Istart;i<Iend;i++) {
     285         348 :       if (i<N) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
     286             :     }
     287           7 :     PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
     288           7 :     PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
     289           7 :     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
     290           7 :     PetscCall(MatSetSizes(Zn,n,m,N,M));
     291           7 :     PetscCall(MatSetFromOptions(Zn));
     292           7 :     PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
     293         364 :     for (i=Istart;i<Iend;i++) {
     294         357 :       if (i<m) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
     295             :     }
     296           7 :     PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
     297           7 :     PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
     298           7 :     PetscCall(MatProductCreate(AT,A,NULL,&ATA));
     299           7 :     PetscCall(MatProductSetType(ATA,MATPRODUCT_AB));
     300           7 :     PetscCall(MatProductSetFromOptions(ATA));
     301           7 :     PetscCall(MatProductSymbolic(ATA));
     302           7 :     PetscCall(MatProductNumeric(ATA));
     303           7 :     PetscCall(MatCreateTile(1.0,Id,1.0,Zm,1.0,Zn,1.0,ATA,C));
     304           7 :     PetscCall(MatDestroy(&Id));
     305           7 :     PetscCall(MatDestroy(&Zm));
     306           7 :     PetscCall(MatDestroy(&Zn));
     307           7 :     PetscCall(MatDestroy(&ATA));
     308             :   } else {
     309           6 :     PetscCall(PetscNew(&ctx));
     310           6 :     ctx->A       = A;
     311           6 :     ctx->AT      = AT;
     312           6 :     ctx->swapped = svd->swapped;
     313           6 :     PetscCall(VecDuplicateEmpty(t,&ctx->x1));
     314           6 :     PetscCall(VecDuplicateEmpty(t,&ctx->y1));
     315           6 :     PetscCall(MatCreateVecsEmpty(A,&ctx->x2,NULL));
     316           6 :     PetscCall(MatCreateVecsEmpty(A,&ctx->y2,NULL));
     317           6 :     PetscCall(MatCreateVecs(A,NULL,&ctx->w));
     318           6 :     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
     319           6 :     PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ECross));
     320           6 :     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_ECross));
     321             : #if defined(PETSC_HAVE_CUDA)
     322             :     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
     323             :     if (cuda) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA));
     324             :     else
     325             : #endif
     326           6 :       PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross));
     327           6 :     PetscCall(MatGetVecType(A,&vtype));
     328           6 :     PetscCall(MatSetVecType(*C,vtype));
     329             : #if defined(PETSC_HAVE_CUDA)
     330             :     if (cuda) {
     331             :       /* check alignment of bottom block */
     332             :       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
     333             :       PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
     334             :       for (i=0;i<size;i++) {
     335             :         ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
     336             :         if (ctx->misaligned) break;
     337             :       }
     338             :       if (ctx->misaligned) {  /* create work vectors for MatMult */
     339             :         PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
     340             :         PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
     341             :       }
     342             :     }
     343             : #endif
     344             :   }
     345          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     346             : }
     347             : 
     348             : /* Convergence test relative to the norm of R (used in GSVD only) */
     349         175 : static PetscErrorCode EPSConv_Cyclic(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     350             : {
     351         175 :   SVD svd = (SVD)ctx;
     352             : 
     353         175 :   PetscFunctionBegin;
     354         175 :   *errest = res/PetscMax(svd->nrma,svd->nrmb);
     355         175 :   PetscFunctionReturn(PETSC_SUCCESS);
     356             : }
     357             : 
     358          46 : static PetscErrorCode SVDSetUp_Cyclic(SVD svd)
     359             : {
     360          46 :   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
     361          46 :   PetscInt          M,N,m,n,p,k,i,isl,offset,nev,ncv,mpd,maxit;
     362          46 :   PetscReal         tol;
     363          46 :   const PetscScalar *isa,*oa;
     364          46 :   PetscScalar       *va;
     365          46 :   EPSProblemType    ptype;
     366          46 :   PetscBool         trackall,issinv;
     367          46 :   Vec               v,t;
     368          46 :   ST                st;
     369          46 :   Mat               Omega;
     370          46 :   MatType           Atype;
     371             : 
     372          46 :   PetscFunctionBegin;
     373          46 :   PetscCall(MatGetSize(svd->A,&M,&N));
     374          46 :   PetscCall(MatGetLocalSize(svd->A,&m,&n));
     375          46 :   if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
     376          46 :   PetscCall(MatDestroy(&cyclic->C));
     377          46 :   PetscCall(MatDestroy(&cyclic->D));
     378          46 :   if (svd->isgeneralized) {
     379          13 :     if (svd->which==SVD_SMALLEST) {  /* alternative pencil */
     380           3 :       PetscCall(MatCreateVecs(svd->B,NULL,&t));
     381           3 :       PetscCall(SVDCyclicGetCyclicMat(svd,svd->B,svd->BT,&cyclic->C));
     382           3 :       PetscCall(SVDCyclicGetECrossMat(svd,svd->A,svd->AT,&cyclic->D,t));
     383             :     } else {
     384          10 :       PetscCall(MatCreateVecs(svd->A,NULL,&t));
     385          10 :       PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
     386          10 :       PetscCall(SVDCyclicGetECrossMat(svd,svd->B,svd->BT,&cyclic->D,t));
     387             :     }
     388          13 :     PetscCall(VecDestroy(&t));
     389          13 :     PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,cyclic->D));
     390          13 :     PetscCall(EPSGetProblemType(cyclic->eps,&ptype));
     391          13 :     if (!ptype) PetscCall(EPSSetProblemType(cyclic->eps,EPS_GHEP));
     392          33 :   } else if (svd->ishyperbolic) {
     393           8 :     PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
     394           8 :     PetscCall(MatCreateVecs(cyclic->C,&v,NULL));
     395           8 :     PetscCall(VecSet(v,1.0));
     396           8 :     PetscCall(VecGetArrayRead(svd->omega,&oa));
     397           8 :     PetscCall(VecGetArray(v,&va));
     398           8 :     if (svd->swapped) PetscCall(PetscArraycpy(va+m,oa,n));
     399           6 :     else PetscCall(PetscArraycpy(va,oa,m));
     400           8 :     PetscCall(VecRestoreArrayRead(svd->omega,&oa));
     401           8 :     PetscCall(VecRestoreArray(v,&va));
     402           8 :     PetscCall(MatGetType(svd->OP,&Atype));
     403           8 :     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
     404           8 :     PetscCall(MatSetSizes(Omega,m+n,m+n,M+N,M+N));
     405           8 :     PetscCall(MatSetType(Omega,Atype));
     406           8 :     PetscCall(MatDiagonalSet(Omega,v,INSERT_VALUES));
     407           8 :     PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,Omega));
     408           8 :     PetscCall(EPSSetProblemType(cyclic->eps,EPS_GHIEP));
     409           8 :     PetscCall(MatDestroy(&Omega));
     410           8 :     PetscCall(VecDestroy(&v));
     411             :   } else {
     412          25 :     PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
     413          25 :     PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,NULL));
     414          25 :     PetscCall(EPSSetProblemType(cyclic->eps,EPS_HEP));
     415             :   }
     416          46 :   if (!cyclic->usereps) {
     417          45 :     if (svd->which == SVD_LARGEST) {
     418          39 :       PetscCall(EPSGetST(cyclic->eps,&st));
     419          39 :       PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
     420          39 :       if (issinv) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE));
     421          38 :       else if (svd->ishyperbolic) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_MAGNITUDE));
     422          31 :       else PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
     423             :     } else {
     424           6 :       if (svd->isgeneralized) {  /* computes sigma^{-1} via alternative pencil */
     425           3 :         PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
     426             :       } else {
     427           3 :         if (svd->ishyperbolic) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE));
     428           2 :         else PetscCall(EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL));
     429           3 :         PetscCall(EPSSetTarget(cyclic->eps,0.0));
     430             :       }
     431             :     }
     432          45 :     PetscCall(EPSGetDimensions(cyclic->eps,&nev,&ncv,&mpd));
     433          45 :     PetscCheck(nev==1 || nev>=2*svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"The number of requested eigenvalues %" PetscInt_FMT " must be at least 2*%" PetscInt_FMT,nev,svd->nsv);
     434          45 :     nev = PetscMax(nev,2*svd->nsv);
     435          45 :     if (ncv==PETSC_DEFAULT && svd->ncv!=PETSC_DEFAULT) ncv = PetscMax(3*svd->nsv,svd->ncv);
     436          45 :     if (mpd==PETSC_DEFAULT && svd->mpd!=PETSC_DEFAULT) mpd = svd->mpd;
     437          45 :     PetscCall(EPSSetDimensions(cyclic->eps,nev,ncv,mpd));
     438          45 :     PetscCall(EPSGetTolerances(cyclic->eps,&tol,&maxit));
     439          51 :     if (tol==(PetscReal)PETSC_DEFAULT) tol = svd->tol==(PetscReal)PETSC_DEFAULT? SLEPC_DEFAULT_TOL/10.0: svd->tol;
     440          45 :     if (maxit==PETSC_DEFAULT && svd->max_it!=PETSC_DEFAULT) maxit = svd->max_it;
     441          45 :     PetscCall(EPSSetTolerances(cyclic->eps,tol,maxit));
     442          45 :     switch (svd->conv) {
     443           3 :     case SVD_CONV_ABS:
     444           3 :       PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS));break;
     445          29 :     case SVD_CONV_REL:
     446          29 :       PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL));break;
     447          13 :     case SVD_CONV_NORM:
     448          13 :       if (svd->isgeneralized) {
     449          13 :         if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
     450          13 :         if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
     451          13 :         PetscCall(EPSSetConvergenceTestFunction(cyclic->eps,EPSConv_Cyclic,svd,NULL));
     452             :       } else {
     453           0 :         PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_NORM));break;
     454             :       }
     455             :       break;
     456           0 :     case SVD_CONV_MAXIT:
     457           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
     458           0 :     case SVD_CONV_USER:
     459           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
     460             :     }
     461           1 :   }
     462          46 :   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
     463             :   /* Transfer the trackall option from svd to eps */
     464          46 :   PetscCall(SVDGetTrackAll(svd,&trackall));
     465          46 :   PetscCall(EPSSetTrackAll(cyclic->eps,trackall));
     466             :   /* Transfer the initial subspace from svd to eps */
     467          46 :   if (svd->nini<0 || svd->ninil<0) {
     468           8 :     for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
     469           4 :       PetscCall(MatCreateVecs(cyclic->C,&v,NULL));
     470           4 :       PetscCall(VecGetArrayWrite(v,&va));
     471           4 :       if (svd->isgeneralized) PetscCall(MatGetLocalSize(svd->B,&p,NULL));
     472           4 :       k = (svd->isgeneralized && svd->which==SVD_SMALLEST)? p: m;  /* size of upper block row */
     473           4 :       if (i<-svd->ninil) {
     474           4 :         PetscCall(VecGetArrayRead(svd->ISL[i],&isa));
     475           4 :         if (svd->isgeneralized) {
     476           1 :           PetscCall(VecGetLocalSize(svd->ISL[i],&isl));
     477           1 :           PetscCheck(isl==m+p,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
     478           1 :           offset = (svd->which==SVD_SMALLEST)? m: 0;
     479           1 :           PetscCall(PetscArraycpy(va,isa+offset,k));
     480             :         } else {
     481           3 :           PetscCall(VecGetLocalSize(svd->ISL[i],&isl));
     482           3 :           PetscCheck(isl==k,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
     483           3 :           PetscCall(PetscArraycpy(va,isa,k));
     484             :         }
     485           4 :         PetscCall(VecRestoreArrayRead(svd->IS[i],&isa));
     486           0 :       } else PetscCall(PetscArrayzero(&va,k));
     487           4 :       if (i<-svd->nini) {
     488           4 :         PetscCall(VecGetLocalSize(svd->IS[i],&isl));
     489           4 :         PetscCheck(isl==n,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for right initial vector");
     490           4 :         PetscCall(VecGetArrayRead(svd->IS[i],&isa));
     491           4 :         PetscCall(PetscArraycpy(va+k,isa,n));
     492           4 :         PetscCall(VecRestoreArrayRead(svd->IS[i],&isa));
     493           0 :       } else PetscCall(PetscArrayzero(va+k,n));
     494           4 :       PetscCall(VecRestoreArrayWrite(v,&va));
     495           4 :       PetscCall(VecDestroy(&svd->IS[i]));
     496           4 :       svd->IS[i] = v;
     497             :     }
     498           4 :     svd->nini = PetscMin(svd->nini,svd->ninil);
     499           4 :     PetscCall(EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS));
     500           4 :     PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
     501           4 :     PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
     502             :   }
     503          46 :   PetscCall(EPSSetUp(cyclic->eps));
     504          46 :   PetscCall(EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd));
     505          46 :   svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
     506          46 :   PetscCall(EPSGetTolerances(cyclic->eps,NULL,&svd->max_it));
     507          46 :   if (svd->tol==(PetscReal)PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
     508             : 
     509          46 :   svd->leftbasis = PETSC_TRUE;
     510          46 :   PetscCall(SVDAllocateSolution(svd,0));
     511          46 :   PetscFunctionReturn(PETSC_SUCCESS);
     512             : }
     513             : 
     514       96350 : static PetscErrorCode SVDCyclicCheckEigenvalue(SVD svd,PetscScalar er,PetscScalar ei,PetscReal *sigma,PetscBool *isreal)
     515             : {
     516       96350 :   PetscFunctionBegin;
     517       96350 :   if (svd->ishyperbolic && PetscDefined(USE_COMPLEX) && PetscAbsReal(PetscImaginaryPart(er))>10*PetscAbsReal(PetscRealPart(er))) {
     518             :     *sigma = PetscImaginaryPart(er);
     519             :     if (isreal) *isreal = PETSC_FALSE;
     520       96350 :   } else if (svd->ishyperbolic && !PetscDefined(USE_COMPLEX) && PetscAbsScalar(ei)>10*PetscAbsScalar(er)) {
     521        6048 :     *sigma = PetscRealPart(ei);
     522        6048 :     if (isreal) *isreal = PETSC_FALSE;
     523             :   } else {
     524       90302 :     *sigma = PetscRealPart(er);
     525       90302 :     if (isreal) *isreal = PETSC_TRUE;
     526             :   }
     527       96350 :   PetscFunctionReturn(PETSC_SUCCESS);
     528             : }
     529             : 
     530          46 : static PetscErrorCode SVDSolve_Cyclic(SVD svd)
     531             : {
     532          46 :   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
     533          46 :   PetscInt       i,j,nconv;
     534          46 :   PetscScalar    er,ei;
     535          46 :   PetscReal      sigma;
     536             : 
     537          46 :   PetscFunctionBegin;
     538          46 :   PetscCall(EPSSolve(cyclic->eps));
     539          46 :   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
     540          46 :   PetscCall(EPSGetIterationNumber(cyclic->eps,&svd->its));
     541          46 :   PetscCall(EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason));
     542         744 :   for (i=0,j=0;i<nconv;i++) {
     543         698 :     PetscCall(EPSGetEigenvalue(cyclic->eps,i,&er,&ei));
     544         698 :     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
     545         698 :     if (sigma>0.0) {
     546         469 :       if (svd->isgeneralized && svd->which==SVD_SMALLEST) svd->sigma[j] = 1.0/sigma;
     547         453 :       else svd->sigma[j] = sigma;
     548         469 :       j++;
     549             :     }
     550             :   }
     551          46 :   svd->nconv = j;
     552          46 :   PetscFunctionReturn(PETSC_SUCCESS);
     553             : }
     554             : 
     555          21 : static PetscErrorCode SVDComputeVectors_Cyclic_Standard(SVD svd)
     556             : {
     557          21 :   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
     558          21 :   PetscInt          i,j,m,nconv;
     559          21 :   PetscScalar       er,ei;
     560          21 :   PetscReal         sigma;
     561          21 :   const PetscScalar *px;
     562          21 :   Vec               x,x1,x2;
     563             : 
     564          21 :   PetscFunctionBegin;
     565          21 :   PetscCall(MatCreateVecs(cyclic->C,&x,NULL));
     566          21 :   PetscCall(MatGetLocalSize(svd->A,&m,NULL));
     567          21 :   PetscCall(MatCreateVecsEmpty(svd->A,&x2,&x1));
     568          21 :   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
     569         186 :   for (i=0,j=0;i<nconv;i++) {
     570         165 :     PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL));
     571         165 :     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
     572         165 :     if (sigma<0.0) continue;
     573         153 :     PetscCall(VecGetArrayRead(x,&px));
     574         153 :     PetscCall(VecPlaceArray(x1,px));
     575         153 :     PetscCall(VecPlaceArray(x2,px+m));
     576         153 :     PetscCall(BVInsertVec(svd->U,j,x1));
     577         153 :     PetscCall(BVScaleColumn(svd->U,j,PETSC_SQRT2));
     578         153 :     PetscCall(BVInsertVec(svd->V,j,x2));
     579         153 :     PetscCall(BVScaleColumn(svd->V,j,PETSC_SQRT2));
     580         153 :     PetscCall(VecResetArray(x1));
     581         153 :     PetscCall(VecResetArray(x2));
     582         153 :     PetscCall(VecRestoreArrayRead(x,&px));
     583         153 :     j++;
     584             :   }
     585          21 :   PetscCall(VecDestroy(&x));
     586          21 :   PetscCall(VecDestroy(&x1));
     587          21 :   PetscCall(VecDestroy(&x2));
     588          21 :   PetscFunctionReturn(PETSC_SUCCESS);
     589             : }
     590             : 
     591          13 : static PetscErrorCode SVDComputeVectors_Cyclic_Generalized(SVD svd)
     592             : {
     593          13 :   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
     594          13 :   PetscInt          i,j,m,p,nconv;
     595          13 :   PetscScalar       *dst,er,ei;
     596          13 :   PetscReal         sigma;
     597          13 :   const PetscScalar *src,*px;
     598          13 :   Vec               u,v,x,x1,x2,uv;
     599             : 
     600          13 :   PetscFunctionBegin;
     601          13 :   PetscCall(MatGetLocalSize(svd->A,&m,NULL));
     602          13 :   PetscCall(MatGetLocalSize(svd->B,&p,NULL));
     603          13 :   PetscCall(MatCreateVecs(cyclic->C,&x,NULL));
     604          13 :   if (svd->which==SVD_SMALLEST) PetscCall(MatCreateVecsEmpty(svd->B,&x1,&x2));
     605          10 :   else PetscCall(MatCreateVecsEmpty(svd->A,&x2,&x1));
     606          13 :   PetscCall(MatCreateVecs(svd->A,NULL,&u));
     607          13 :   PetscCall(MatCreateVecs(svd->B,NULL,&v));
     608          13 :   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
     609         100 :   for (i=0,j=0;i<nconv;i++) {
     610          87 :     PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL));
     611          87 :     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
     612          87 :     if (sigma<0.0) continue;
     613          87 :     if (svd->which==SVD_SMALLEST) {
     614             :       /* evec_i = 1/sqrt(2)*[ v_i; w_i ],  w_i = x_i/c_i */
     615          16 :       PetscCall(VecGetArrayRead(x,&px));
     616          16 :       PetscCall(VecPlaceArray(x2,px));
     617          16 :       PetscCall(VecPlaceArray(x1,px+p));
     618          16 :       PetscCall(VecCopy(x2,v));
     619          16 :       PetscCall(VecScale(v,PETSC_SQRT2));  /* v_i = sqrt(2)*evec_i_1 */
     620          16 :       PetscCall(VecScale(x1,PETSC_SQRT2)); /* w_i = sqrt(2)*evec_i_2 */
     621          16 :       PetscCall(MatMult(svd->A,x1,u));     /* A*w_i = u_i */
     622          16 :       PetscCall(VecScale(x1,1.0/PetscSqrtScalar(1.0+sigma*sigma)));  /* x_i = w_i*c_i */
     623          16 :       PetscCall(BVInsertVec(svd->V,j,x1));
     624          16 :       PetscCall(VecResetArray(x2));
     625          16 :       PetscCall(VecResetArray(x1));
     626          16 :       PetscCall(VecRestoreArrayRead(x,&px));
     627             :     } else {
     628             :       /* evec_i = 1/sqrt(2)*[ u_i; w_i ],  w_i = x_i/s_i */
     629          71 :       PetscCall(VecGetArrayRead(x,&px));
     630          71 :       PetscCall(VecPlaceArray(x1,px));
     631          71 :       PetscCall(VecPlaceArray(x2,px+m));
     632          71 :       PetscCall(VecCopy(x1,u));
     633          71 :       PetscCall(VecScale(u,PETSC_SQRT2));  /* u_i = sqrt(2)*evec_i_1 */
     634          71 :       PetscCall(VecScale(x2,PETSC_SQRT2)); /* w_i = sqrt(2)*evec_i_2 */
     635          71 :       PetscCall(MatMult(svd->B,x2,v));     /* B*w_i = v_i */
     636          71 :       PetscCall(VecScale(x2,1.0/PetscSqrtScalar(1.0+sigma*sigma)));  /* x_i = w_i*s_i */
     637          71 :       PetscCall(BVInsertVec(svd->V,j,x2));
     638          71 :       PetscCall(VecResetArray(x1));
     639          71 :       PetscCall(VecResetArray(x2));
     640          71 :       PetscCall(VecRestoreArrayRead(x,&px));
     641             :     }
     642             :     /* copy [u;v] to U[j] */
     643          87 :     PetscCall(BVGetColumn(svd->U,j,&uv));
     644          87 :     PetscCall(VecGetArrayWrite(uv,&dst));
     645          87 :     PetscCall(VecGetArrayRead(u,&src));
     646          87 :     PetscCall(PetscArraycpy(dst,src,m));
     647          87 :     PetscCall(VecRestoreArrayRead(u,&src));
     648          87 :     PetscCall(VecGetArrayRead(v,&src));
     649          87 :     PetscCall(PetscArraycpy(dst+m,src,p));
     650          87 :     PetscCall(VecRestoreArrayRead(v,&src));
     651          87 :     PetscCall(VecRestoreArrayWrite(uv,&dst));
     652          87 :     PetscCall(BVRestoreColumn(svd->U,j,&uv));
     653          87 :     j++;
     654             :   }
     655          13 :   PetscCall(VecDestroy(&x));
     656          13 :   PetscCall(VecDestroy(&x1));
     657          13 :   PetscCall(VecDestroy(&x2));
     658          13 :   PetscCall(VecDestroy(&u));
     659          13 :   PetscCall(VecDestroy(&v));
     660          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     661             : }
     662             : 
     663             : #if defined(PETSC_USE_COMPLEX)
     664             : /* VecMaxAbs: returns the entry of x that has max(abs(x(i))), using w as a workspace vector */
     665             : static PetscErrorCode VecMaxAbs(Vec x,Vec w,PetscScalar *v)
     666             : {
     667             :   PetscMPIInt       size,rank,root;
     668             :   const PetscScalar *xx;
     669             :   const PetscInt    *ranges;
     670             :   PetscReal         val;
     671             :   PetscInt          p;
     672             : 
     673             :   PetscFunctionBegin;
     674             :   PetscCall(VecCopy(x,w));
     675             :   PetscCall(VecAbs(w));
     676             :   PetscCall(VecMax(w,&p,&val));
     677             :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)x),&size));
     678             :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)x),&rank));
     679             :   PetscCall(VecGetOwnershipRanges(x,&ranges));
     680             :   for (root=0;root<size;root++) if (p>=ranges[root] && p<ranges[root+1]) break;
     681             :   if (rank==root) {
     682             :     PetscCall(VecGetArrayRead(x,&xx));
     683             :     *v = xx[p-ranges[root]];
     684             :     PetscCall(VecRestoreArrayRead(x,&xx));
     685             :   }
     686             :   PetscCallMPI(MPI_Bcast(v,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)x)));
     687             :   PetscFunctionReturn(PETSC_SUCCESS);
     688             : }
     689             : #endif
     690             : 
     691           8 : static PetscErrorCode SVDComputeVectors_Cyclic_Hyperbolic(SVD svd)
     692             : {
     693           8 :   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
     694           8 :   PetscInt          i,j,m,n,nconv;
     695           8 :   PetscScalar       er,ei;
     696           8 :   PetscReal         sigma,nrm;
     697           8 :   PetscBool         isreal;
     698           8 :   const PetscScalar *px;
     699           8 :   Vec               u,x,xi=NULL,x1,x2,x1i=NULL,x2i;
     700           8 :   BV                U=NULL,V=NULL;
     701             : #if !defined(PETSC_USE_COMPLEX)
     702           8 :   const PetscScalar *pxi;
     703           8 :   PetscReal         nrmr,nrmi;
     704             : #else
     705             :   PetscScalar       alpha;
     706             : #endif
     707             : 
     708           8 :   PetscFunctionBegin;
     709           8 :   PetscCall(MatCreateVecs(cyclic->C,&x,svd->ishyperbolic?&xi:NULL));
     710           8 :   PetscCall(MatGetLocalSize(svd->A,&m,NULL));
     711           8 :   PetscCall(MatCreateVecsEmpty(svd->OP,&x2,&x1));
     712             : #if defined(PETSC_USE_COMPLEX)
     713             :   PetscCall(MatCreateVecs(svd->OP,&x2i,&x1i));
     714             : #else
     715           8 :   PetscCall(MatCreateVecsEmpty(svd->OP,&x2i,&x1i));
     716             : #endif
     717             :   /* set-up Omega-normalization of U */
     718           8 :   U = svd->swapped? svd->V: svd->U;
     719           8 :   V = svd->swapped? svd->U: svd->V;
     720           8 :   PetscCall(BVGetSizes(U,&n,NULL,NULL));
     721           8 :   PetscCall(BV_SetMatrixDiagonal(U,svd->omega,svd->A));
     722           8 :   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
     723         442 :   for (i=0,j=0;i<nconv;i++) {
     724         434 :     PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,xi));
     725         434 :     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,&isreal));
     726         434 :     if (sigma<0.0) continue;
     727         217 :     PetscCall(VecGetArrayRead(x,&px));
     728         217 :     if (svd->swapped) {
     729          12 :       PetscCall(VecPlaceArray(x2,px));
     730          12 :       PetscCall(VecPlaceArray(x1,px+m));
     731             :     } else {
     732         205 :       PetscCall(VecPlaceArray(x1,px));
     733         205 :       PetscCall(VecPlaceArray(x2,px+n));
     734             :     }
     735             : #if defined(PETSC_USE_COMPLEX)
     736             :     PetscCall(BVInsertVec(U,j,x1));
     737             :     PetscCall(BVInsertVec(V,j,x2));
     738             :     if (!isreal) {
     739             :       PetscCall(VecMaxAbs(x1,x1i,&alpha));
     740             :       PetscCall(BVScaleColumn(U,j,PetscAbsScalar(alpha)/alpha));
     741             :       PetscCall(BVScaleColumn(V,j,PetscAbsScalar(alpha)/(alpha*PETSC_i)));
     742             :     }
     743             : #else
     744         217 :     PetscCall(VecGetArrayRead(xi,&pxi));
     745         217 :     if (svd->swapped) {
     746          12 :       PetscCall(VecPlaceArray(x2i,pxi));
     747          12 :       PetscCall(VecPlaceArray(x1i,pxi+m));
     748             :     } else {
     749         205 :       PetscCall(VecPlaceArray(x1i,pxi));
     750         205 :       PetscCall(VecPlaceArray(x2i,pxi+n));
     751             :     }
     752         217 :     PetscCall(VecNorm(x2,NORM_2,&nrmr));
     753         217 :     PetscCall(VecNorm(x2i,NORM_2,&nrmi));
     754         217 :     if (nrmi>nrmr) {
     755           3 :       if (isreal) {
     756           0 :         PetscCall(BVInsertVec(U,j,x1i));
     757           0 :         PetscCall(BVInsertVec(V,j,x2i));
     758             :       } else {
     759           3 :         PetscCall(BVInsertVec(U,j,x1));
     760           3 :         PetscCall(BVInsertVec(V,j,x2i));
     761             :       }
     762             :     } else {
     763         214 :       if (isreal) {
     764         213 :         PetscCall(BVInsertVec(U,j,x1));
     765         213 :         PetscCall(BVInsertVec(V,j,x2));
     766             :       } else {
     767           1 :         PetscCall(BVInsertVec(U,j,x1i));
     768           1 :         PetscCall(BVScaleColumn(U,j,-1.0));
     769           1 :         PetscCall(BVInsertVec(V,j,x2));
     770             :       }
     771             :     }
     772         217 :     PetscCall(VecResetArray(x1i));
     773         217 :     PetscCall(VecResetArray(x2i));
     774         217 :     PetscCall(VecRestoreArrayRead(xi,&pxi));
     775             : #endif
     776         217 :     PetscCall(VecResetArray(x1));
     777         217 :     PetscCall(VecResetArray(x2));
     778         217 :     PetscCall(VecRestoreArrayRead(x,&px));
     779         217 :     PetscCall(BVGetColumn(U,j,&u));
     780         217 :     PetscCall(VecPointwiseMult(u,u,svd->omega));
     781         217 :     PetscCall(BVRestoreColumn(U,j,&u));
     782         217 :     PetscCall(BVNormColumn(U,j,NORM_2,&nrm));
     783         217 :     PetscCall(BVScaleColumn(U,j,1.0/PetscAbs(nrm)));
     784         217 :     svd->sign[j] = PetscSign(nrm);
     785         217 :     PetscCall(BVNormColumn(V,j,NORM_2,&nrm));
     786         217 :     PetscCall(BVScaleColumn(V,j,1.0/nrm));
     787         217 :     j++;
     788             :   }
     789           8 :   PetscCall(VecDestroy(&x));
     790           8 :   PetscCall(VecDestroy(&x1));
     791           8 :   PetscCall(VecDestroy(&x2));
     792           8 :   PetscCall(VecDestroy(&xi));
     793           8 :   PetscCall(VecDestroy(&x1i));
     794           8 :   PetscCall(VecDestroy(&x2i));
     795           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     796             : }
     797             : 
     798          42 : static PetscErrorCode SVDComputeVectors_Cyclic(SVD svd)
     799             : {
     800          42 :   PetscFunctionBegin;
     801          42 :   switch (svd->problem_type) {
     802          21 :     case SVD_STANDARD:
     803          21 :       PetscCall(SVDComputeVectors_Cyclic_Standard(svd));
     804             :       break;
     805          13 :     case SVD_GENERALIZED:
     806          13 :       PetscCall(SVDComputeVectors_Cyclic_Generalized(svd));
     807             :       break;
     808           8 :     case SVD_HYPERBOLIC:
     809           8 :       PetscCall(SVDComputeVectors_Cyclic_Hyperbolic(svd));
     810             :       break;
     811           0 :     default:
     812           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Unknown singular value problem type");
     813             :   }
     814          42 :   PetscFunctionReturn(PETSC_SUCCESS);
     815             : }
     816             : 
     817        3429 : static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
     818             : {
     819        3429 :   PetscInt       i,j;
     820        3429 :   SVD            svd = (SVD)ctx;
     821        3429 :   PetscScalar    er,ei;
     822        3429 :   PetscReal      sigma;
     823        3429 :   ST             st;
     824             : 
     825        3429 :   PetscFunctionBegin;
     826        3429 :   nconv = 0;
     827        3429 :   PetscCall(EPSGetST(eps,&st));
     828       98395 :   for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
     829       94966 :     er = eigr[i]; ei = eigi[i];
     830       94966 :     PetscCall(STBackTransform(st,1,&er,&ei));
     831       94966 :     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
     832       94966 :     if (sigma>0.0) {
     833       42441 :       svd->sigma[j]  = sigma;
     834       42441 :       svd->errest[j] = errest[i];
     835       42441 :       if (errest[i] && errest[i] < svd->tol) nconv++;
     836       42441 :       j++;
     837             :     }
     838             :   }
     839        3429 :   nest = j;
     840        3429 :   PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
     841        3429 :   PetscFunctionReturn(PETSC_SUCCESS);
     842             : }
     843             : 
     844          34 : static PetscErrorCode SVDSetFromOptions_Cyclic(SVD svd,PetscOptionItems *PetscOptionsObject)
     845             : {
     846          34 :   PetscBool      set,val;
     847          34 :   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
     848          34 :   ST             st;
     849             : 
     850          34 :   PetscFunctionBegin;
     851          34 :   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cyclic Options");
     852             : 
     853          34 :     PetscCall(PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set));
     854          34 :     if (set) PetscCall(SVDCyclicSetExplicitMatrix(svd,val));
     855             : 
     856          34 :   PetscOptionsHeadEnd();
     857             : 
     858          34 :   if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
     859          34 :   if (!cyclic->explicitmatrix && !cyclic->usereps) {
     860             :     /* use as default an ST with shell matrix and Jacobi */
     861          20 :     PetscCall(EPSGetST(cyclic->eps,&st));
     862          20 :     PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
     863             :   }
     864          34 :   PetscCall(EPSSetFromOptions(cyclic->eps));
     865          34 :   PetscFunctionReturn(PETSC_SUCCESS);
     866             : }
     867             : 
     868          22 : static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
     869             : {
     870          22 :   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
     871             : 
     872          22 :   PetscFunctionBegin;
     873          22 :   if (cyclic->explicitmatrix != explicitmat) {
     874          14 :     cyclic->explicitmatrix = explicitmat;
     875          14 :     svd->state = SVD_STATE_INITIAL;
     876             :   }
     877          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     878             : }
     879             : 
     880             : /*@
     881             :    SVDCyclicSetExplicitMatrix - Indicate if the eigensolver operator
     882             :    H(A) = [ 0  A ; A^T 0 ] must be computed explicitly.
     883             : 
     884             :    Logically Collective
     885             : 
     886             :    Input Parameters:
     887             : +  svd         - singular value solver
     888             : -  explicitmat - boolean flag indicating if H(A) is built explicitly
     889             : 
     890             :    Options Database Key:
     891             : .  -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag
     892             : 
     893             :    Level: advanced
     894             : 
     895             : .seealso: SVDCyclicGetExplicitMatrix()
     896             : @*/
     897          22 : PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
     898             : {
     899          22 :   PetscFunctionBegin;
     900          22 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     901          88 :   PetscValidLogicalCollectiveBool(svd,explicitmat,2);
     902          22 :   PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
     903          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     904             : }
     905             : 
     906           1 : static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
     907             : {
     908           1 :   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
     909             : 
     910           1 :   PetscFunctionBegin;
     911           1 :   *explicitmat = cyclic->explicitmatrix;
     912           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     913             : }
     914             : 
     915             : /*@
     916             :    SVDCyclicGetExplicitMatrix - Returns the flag indicating if H(A) is built explicitly.
     917             : 
     918             :    Not Collective
     919             : 
     920             :    Input Parameter:
     921             : .  svd  - singular value solver
     922             : 
     923             :    Output Parameter:
     924             : .  explicitmat - the mode flag
     925             : 
     926             :    Level: advanced
     927             : 
     928             : .seealso: SVDCyclicSetExplicitMatrix()
     929             : @*/
     930           1 : PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
     931             : {
     932           1 :   PetscFunctionBegin;
     933           1 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     934           1 :   PetscAssertPointer(explicitmat,2);
     935           1 :   PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
     936           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     937             : }
     938             : 
     939           1 : static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
     940             : {
     941           1 :   SVD_CYCLIC      *cyclic = (SVD_CYCLIC*)svd->data;
     942             : 
     943           1 :   PetscFunctionBegin;
     944           1 :   PetscCall(PetscObjectReference((PetscObject)eps));
     945           1 :   PetscCall(EPSDestroy(&cyclic->eps));
     946           1 :   cyclic->eps     = eps;
     947           1 :   cyclic->usereps = PETSC_TRUE;
     948           1 :   svd->state      = SVD_STATE_INITIAL;
     949           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     950             : }
     951             : 
     952             : /*@
     953             :    SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
     954             :    singular value solver.
     955             : 
     956             :    Collective
     957             : 
     958             :    Input Parameters:
     959             : +  svd - singular value solver
     960             : -  eps - the eigensolver object
     961             : 
     962             :    Level: advanced
     963             : 
     964             : .seealso: SVDCyclicGetEPS()
     965             : @*/
     966           1 : PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
     967             : {
     968           1 :   PetscFunctionBegin;
     969           1 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     970           1 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
     971           1 :   PetscCheckSameComm(svd,1,eps,2);
     972           1 :   PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
     973           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     974             : }
     975             : 
     976          35 : static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
     977             : {
     978          35 :   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
     979             : 
     980          35 :   PetscFunctionBegin;
     981          35 :   if (!cyclic->eps) {
     982          35 :     PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps));
     983          35 :     PetscCall(PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1));
     984          35 :     PetscCall(EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix));
     985          35 :     PetscCall(EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_"));
     986          35 :     PetscCall(PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options));
     987          35 :     PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
     988          35 :     PetscCall(EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL));
     989             :   }
     990          35 :   *eps = cyclic->eps;
     991          35 :   PetscFunctionReturn(PETSC_SUCCESS);
     992             : }
     993             : 
     994             : /*@
     995             :    SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
     996             :    to the singular value solver.
     997             : 
     998             :    Collective
     999             : 
    1000             :    Input Parameter:
    1001             : .  svd - singular value solver
    1002             : 
    1003             :    Output Parameter:
    1004             : .  eps - the eigensolver object
    1005             : 
    1006             :    Level: advanced
    1007             : 
    1008             : .seealso: SVDCyclicSetEPS()
    1009             : @*/
    1010          35 : PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
    1011             : {
    1012          35 :   PetscFunctionBegin;
    1013          35 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
    1014          35 :   PetscAssertPointer(eps,2);
    1015          35 :   PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
    1016          35 :   PetscFunctionReturn(PETSC_SUCCESS);
    1017             : }
    1018             : 
    1019           1 : static PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
    1020             : {
    1021           1 :   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
    1022           1 :   PetscBool      isascii;
    1023             : 
    1024           1 :   PetscFunctionBegin;
    1025           1 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    1026           1 :   if (isascii) {
    1027           1 :     if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
    1028           1 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit"));
    1029           1 :     PetscCall(PetscViewerASCIIPushTab(viewer));
    1030           1 :     PetscCall(EPSView(cyclic->eps,viewer));
    1031           1 :     PetscCall(PetscViewerASCIIPopTab(viewer));
    1032             :   }
    1033           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1034             : }
    1035             : 
    1036          37 : static PetscErrorCode SVDReset_Cyclic(SVD svd)
    1037             : {
    1038          37 :   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
    1039             : 
    1040          37 :   PetscFunctionBegin;
    1041          37 :   PetscCall(EPSReset(cyclic->eps));
    1042          37 :   PetscCall(MatDestroy(&cyclic->C));
    1043          37 :   PetscCall(MatDestroy(&cyclic->D));
    1044          37 :   PetscFunctionReturn(PETSC_SUCCESS);
    1045             : }
    1046             : 
    1047          36 : static PetscErrorCode SVDDestroy_Cyclic(SVD svd)
    1048             : {
    1049          36 :   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
    1050             : 
    1051          36 :   PetscFunctionBegin;
    1052          36 :   PetscCall(EPSDestroy(&cyclic->eps));
    1053          36 :   PetscCall(PetscFree(svd->data));
    1054          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL));
    1055          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL));
    1056          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL));
    1057          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL));
    1058          36 :   PetscFunctionReturn(PETSC_SUCCESS);
    1059             : }
    1060             : 
    1061          36 : SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
    1062             : {
    1063          36 :   SVD_CYCLIC     *cyclic;
    1064             : 
    1065          36 :   PetscFunctionBegin;
    1066          36 :   PetscCall(PetscNew(&cyclic));
    1067          36 :   svd->data                = (void*)cyclic;
    1068          36 :   svd->ops->solve          = SVDSolve_Cyclic;
    1069          36 :   svd->ops->solveg         = SVDSolve_Cyclic;
    1070          36 :   svd->ops->solveh         = SVDSolve_Cyclic;
    1071          36 :   svd->ops->setup          = SVDSetUp_Cyclic;
    1072          36 :   svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
    1073          36 :   svd->ops->destroy        = SVDDestroy_Cyclic;
    1074          36 :   svd->ops->reset          = SVDReset_Cyclic;
    1075          36 :   svd->ops->view           = SVDView_Cyclic;
    1076          36 :   svd->ops->computevectors = SVDComputeVectors_Cyclic;
    1077          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic));
    1078          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic));
    1079          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic));
    1080          36 :   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic));
    1081          36 :   PetscFunctionReturn(PETSC_SUCCESS);
    1082             : }