Actual source code: test19.c

slepc-3.18.1 2022-11-02
Report Typos and Errors
  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: */

 11: static char help[] = "Test RSVD on a low-rank matrix.\n\n";

 13: #include <slepcsvd.h>

 15: int main(int argc,char **argv)
 16: {
 17:   Mat            A,Ur,Vr;
 18:   SVD            svd;
 19:   PetscInt       m=80,n=40,rank=3,Istart,Iend,i,j;
 20:   PetscScalar    *u;
 21:   PetscReal      tol=PETSC_SMALL;

 24:   SlepcInitialize(&argc,&argv,(char*)0,help);

 26:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 28:   PetscOptionsGetInt(NULL,NULL,"-rank",&rank,NULL);
 31:   PetscPrintf(PETSC_COMM_WORLD,"\nSVD of low-rank matrix, size %" PetscInt_FMT "x%" PetscInt_FMT ", rank %" PetscInt_FMT "\n\n",m,n,rank);

 33:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 34:                   Create a low-rank matrix A = Ur*Vr'
 35:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 37:   MatCreate(PETSC_COMM_WORLD,&Ur);
 38:   MatSetSizes(Ur,PETSC_DECIDE,PETSC_DECIDE,m,rank);
 39:   MatSetType(Ur,MATDENSE);
 40:   MatSetUp(Ur);
 41:   MatGetOwnershipRange(Ur,&Istart,&Iend);
 42:   MatDenseGetArray(Ur,&u);
 43:   for (i=Istart;i<Iend;i++) {
 44:     if (i<m/2) u[i-Istart] = 1.0;
 45:     if (i==0) u[i+Iend-2*Istart] = -2.0;
 46:     if (i==1) u[i+Iend-2*Istart] = -1.0;
 47:     if (i==2) u[i+Iend-2*Istart] = -1.0;
 48:     for (j=2;j<rank;j++) if (i==j) u[j+j*Iend-(j+1)*Istart] = j;
 49:   }
 50:   MatDenseRestoreArray(Ur,&u);
 51:   MatAssemblyBegin(Ur,MAT_FINAL_ASSEMBLY);
 52:   MatAssemblyEnd(Ur,MAT_FINAL_ASSEMBLY);
 53:   MatCreate(PETSC_COMM_WORLD,&Vr);
 54:   MatSetSizes(Vr,PETSC_DECIDE,PETSC_DECIDE,n,rank);
 55:   MatSetType(Vr,MATDENSE);
 56:   MatSetUp(Vr);
 57:   MatGetOwnershipRange(Vr,&Istart,&Iend);
 58:   MatDenseGetArray(Vr,&u);
 59:   for (i=Istart;i<Iend;i++) {
 60:     if (i>n/2) u[i-Istart] = 1.0;
 61:     if (i==0) u[i+Iend-2*Istart] = 1.0;
 62:     if (i==1) u[i+Iend-2*Istart] = 2.0;
 63:     if (i==2) u[i+Iend-2*Istart] = 2.0;
 64:     for (j=2;j<rank;j++) if (i==j) u[j+j*Iend-(j+1)*Istart] = j;
 65:   }
 66:   MatDenseRestoreArray(Vr,&u);
 67:   MatAssemblyBegin(Vr,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyEnd(Vr,MAT_FINAL_ASSEMBLY);
 69:   MatCreateLRC(NULL,Ur,NULL,Vr,&A);
 70:   MatDestroy(&Ur);
 71:   MatDestroy(&Vr);

 73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:                      Create the SVD solver
 75:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 77:   SVDCreate(PETSC_COMM_WORLD,&svd);
 78:   SVDSetOperators(svd,A,NULL);
 79:   SVDSetType(svd,SVDRANDOMIZED);
 80:   SVDSetDimensions(svd,2,PETSC_DEFAULT,PETSC_DEFAULT);
 81:   SVDSetConvergenceTest(svd,SVD_CONV_MAXIT);
 82:   SVDSetTolerances(svd,tol,1);   /* maxit=1 to disable outer iteration */
 83:   SVDSetFromOptions(svd);

 85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 86:                 Compute the singular triplets and display solution
 87:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 89:   SVDSolve(svd);
 90:   SVDErrorView(svd,SVD_ERROR_RELATIVE,NULL);
 91:   SVDDestroy(&svd);
 92:   MatDestroy(&A);
 93:   SlepcFinalize();
 94:   return 0;
 95: }

 97: /*TEST

 99:    test:
100:       args: -bv_orthog_block tsqr    # currently fail with other block orthogonalization methods
101:       requires: !single

103: TEST*/