Actual source code: ex52.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: static char help[] = "Partial hyperbolic singular value decomposition (HSVD) from a file.\n\n"
     12:   "The command line options are:\n"
     13:   "  -file <filename>, PETSc binary file containing matrix A.\n"
     14:   "  -p <p>, where <p> = number of -1's in signature.\n"
     15:   "  -transpose, to transpose the matrix before doing the computation.\n\n";

     17: #include <slepcsvd.h>

     19: int main(int argc,char **argv)
     20: {
     21:   Mat            A,Omega;         /* operator matrix, signature matrix */
     22:   SVD            svd;             /* singular value problem solver context */
     23:   Mat            At;
     24:   Vec            u,v,vomega,*U,*V;
     25:   MatType        Atype;
     26:   PetscReal      tol,lev1=0.0,lev2=0.0;
     27:   PetscInt       M,N,p=0,i,Istart,Iend,nconv,nsv;
     28:   char           filename[PETSC_MAX_PATH_LEN];
     29:   PetscViewer    viewer;
     30:   PetscBool      flg,terse,skiporth=PETSC_FALSE,transpose=PETSC_FALSE;

     32:   PetscFunctionBeginUser;
     33:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

     35:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     36:         Load matrix that defines the hyperbolic singular value problem
     37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     39:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHyperbolic singular value problem stored in file.\n\n"));
     40:   PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
     41:   PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -file option");

     43: #if defined(PETSC_USE_COMPLEX)
     44:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
     45: #else
     46:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
     47: #endif
     48:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
     49:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     50:   PetscCall(MatSetFromOptions(A));
     51:   PetscCall(MatLoad(A,viewer));
     52:   PetscCall(PetscViewerDestroy(&viewer));

     54:   /* transpose the matrix if requested */
     55:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-transpose",&transpose,NULL));
     56:   if (transpose) {
     57:     PetscCall(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&At));
     58:     PetscCall(MatDestroy(&A));
     59:     A = At;
     60:   }

     62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     63:                           Create the signature
     64:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     66:   PetscCall(MatGetSize(A,&M,&N));
     67:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
     68:   PetscCheck(p>=0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p cannot be negative");
     69:   PetscCheck(p<M,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p cannot be larger than the number of rows of A");
     70:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Matrix dimensions: %" PetscInt_FMT "x%" PetscInt_FMT,M,N));
     71:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,", using signature Omega=blkdiag(-I_%" PetscInt_FMT ",I_%" PetscInt_FMT ")\n\n",p,M-p));

     73:   PetscCall(MatCreateVecs(A,NULL,&vomega));
     74:   PetscCall(VecSet(vomega,1.0));
     75:   PetscCall(VecGetOwnershipRange(vomega,&Istart,&Iend));
     76:   for (i=Istart;i<Iend;i++) {
     77:     if (i<p) PetscCall(VecSetValue(vomega,i,-1.0,INSERT_VALUES));
     78:   }
     79:   PetscCall(VecAssemblyBegin(vomega));
     80:   PetscCall(VecAssemblyEnd(vomega));

     82:   PetscCall(MatGetType(A,&Atype));
     83:   PetscCall(MatCreate(PETSC_COMM_WORLD,&Omega));
     84:   PetscCall(MatSetSizes(Omega,PETSC_DECIDE,PETSC_DECIDE,M,M));
     85:   PetscCall(MatSetType(Omega,Atype));
     86:   PetscCall(MatDiagonalSet(Omega,vomega,INSERT_VALUES));

     88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     89:           Create the singular value solver and set various options
     90:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     92:   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));

     94:   PetscCall(SVDSetOperators(svd,A,NULL));
     95:   PetscCall(SVDSetSignature(svd,vomega));
     96:   PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));

     98:   PetscCall(SVDSetFromOptions(svd));

    100:   PetscCall(SVDIsHyperbolic(svd,&flg));
    101:   PetscCheck(flg,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_COR,"Problem should be hyperbolic");

    103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    104:                  Solve the problem, display solution
    105:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    107:   PetscCall(MatCreateVecs(A,&v,&u));
    108:   PetscCall(VecSet(u,1.0));
    109:   PetscCall(VecSet(v,1.0));
    110:   PetscCall(SVDSolve(svd));

    112:   /* show detailed info unless -terse option is given by user */
    113:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
    114:   if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
    115:   else {
    116:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
    117:     PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
    118:     PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
    119:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
    120:   }

    122:   /* check orthogonality */
    123:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-skiporth",&skiporth,NULL));
    124:   PetscCall(SVDGetConverged(svd,&nconv));
    125:   PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
    126:   if (nsv) nconv = PetscMin(nconv,nsv);
    127:   if (nconv>0 && !skiporth) {
    128:     PetscCall(SVDGetTolerances(svd,&tol,NULL));
    129:     PetscCall(VecDuplicateVecs(u,nconv,&U));
    130:     PetscCall(VecDuplicateVecs(v,nconv,&V));
    131:     for (i=0;i<nconv;i++) PetscCall(SVDGetSingularTriplet(svd,i,NULL,U[i],V[i]));
    132:     PetscCall(VecCheckOrthonormality(U,nconv,NULL,nconv,Omega,NULL,&lev1));
    133:     PetscCall(VecCheckOrthonormality(V,nconv,NULL,nconv,NULL,NULL,&lev2));
    134:     if (lev1+lev2<20*tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n"));
    135:     else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g (U) %g (V)\n",(double)lev1,(double)lev2));
    136:     PetscCall(VecDestroyVecs(nconv,&U));
    137:     PetscCall(VecDestroyVecs(nconv,&V));
    138:   }
    139:   PetscCall(VecDestroy(&u));
    140:   PetscCall(VecDestroy(&v));

    142:   /* free work space */
    143:   PetscCall(SVDDestroy(&svd));
    144:   PetscCall(MatDestroy(&A));
    145:   PetscCall(MatDestroy(&Omega));
    146:   PetscCall(VecDestroy(&vomega));
    147:   PetscCall(SlepcFinalize());
    148:   return 0;
    149: }

    151: /*TEST

    153:    testset:
    154:       args: -file ${DATAFILESPATH}/matrices/real/illc1033.petsc -svd_nsv 62 -p 40 -terse
    155:       requires: double !complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
    156:       filter: grep -v Reading
    157:       output_file: output/ex52_1.out
    158:       test:
    159:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}} -svd_implicittranspose {{0 1}}
    160:          suffix: 1_cross
    161:       test:
    162:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_ncv 300
    163:          suffix: 1_cyclic
    164:       test:
    165:          args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}}
    166:          suffix: 1_trlanczos
    167:       test:
    168:          args: -svd_type lapack
    169:          suffix: 1_lapack

    171:    testset:
    172:       args: -file ${DATAFILESPATH}/matrices/real/illc1033.petsc -transpose -svd_nsv 6 -p 130 -terse
    173:       requires: double !complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
    174:       filter: grep -v Reading
    175:       output_file: output/ex52_2.out
    176:       test:
    177:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}} -svd_implicittranspose {{0 1}} -svd_ncv 100
    178:          suffix: 2_cross
    179:       test:
    180:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_ncv 25
    181:          suffix: 2_cyclic
    182:       test:
    183:          args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}} -svd_implicittranspose {{0 1}}
    184:          suffix: 2_trlanczos

    186:    testset:
    187:       args: -file ${DATAFILESPATH}/matrices/complex/illc1033.petsc -svd_nsv 62 -p 40 -terse
    188:       requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
    189:       filter: grep -v Reading
    190:       output_file: output/ex52_1.out
    191:       test:
    192:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
    193:          suffix: 3_cross
    194:       test:
    195:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_cyclic_bv_definite_tol 1e-13 -svd_cyclic_st_ksp_type gcr -svd_cyclic_st_pc_type jacobi -svd_ncv 230
    196:          suffix: 3_cyclic
    197:       test:
    198:          args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}} -bv_definite_tol 1e-13
    199:          suffix: 3_trlanczos
    200:       test:
    201:          args: -svd_type lapack
    202:          suffix: 3_lapack

    204:    testset:
    205:       args: -file ${DATAFILESPATH}/matrices/complex/illc1033.petsc -transpose -svd_nsv 6 -p 130 -terse
    206:       requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
    207:       filter: grep -v Reading
    208:       output_file: output/ex52_2.out
    209:       test:
    210:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}} -svd_ncv 100 -svd_cross_bv_definite_tol 1e-14
    211:          suffix: 4_cross
    212:       test:
    213:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_ncv 26
    214:          suffix: 4_cyclic
    215:       test:
    216:          args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}}
    217:          suffix: 4_trlanczos

    219:    testset:
    220:       args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/rdb200.petsc -svd_smallest -svd_nsv 3 -p 1 -terse
    221:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
    222:       filter: grep -v Reading
    223:       output_file: output/ex52_5.out
    224:       test:
    225:          args: -svd_type cross -svd_max_it 1000 -svd_cross_bv_definite_tol 1e-14
    226:          suffix: 5_cross
    227:       test:
    228:          args: -svd_type cyclic -svd_max_it 4000 -svd_ncv 36 -svd_cyclic_st_ksp_type preonly -svd_cyclic_st_pc_type jacobi -svd_cyclic_bv_definite_tol 1e-14 -svd_cyclic_eps_krylovschur_restart .7
    229:          suffix: 5_cyclic
    230:       test:
    231:          args: -svd_type trlanczos -svd_max_it 4000 -svd_ncv 28 -bv_definite_tol 1e-14
    232:          suffix: 5_trlanczos
    233:       test:
    234:          args: -svd_type lapack
    235:          suffix: 5_lapack

    237:    testset:
    238:       args: -svd_type {{trlanczos cross}} -terse -p 100
    239:       test:
    240:          suffix: 6
    241:          filter: sed -e "s/27.29445, 27.29445/27.29445/"
    242:          args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/rdb200.petsc -svd_threshold_relative 0.8
    243:          requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
    244:       test:
    245:          suffix: 6_complex
    246:          args: -file ${DATAFILESPATH}/matrices/complex/qc324.petsc -svd_threshold_relative 0.6
    247:          requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)

    249: TEST*/