Actual source code: test16.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[] = "Tests multiple calls to SVDSolve with equal matrix size (GSVD).\n\n"
     12:   "The command line options are:\n"
     13:   "  -m <m>, where <m> = number of rows of A.\n"
     14:   "  -n <n>, where <n> = number of columns of A.\n"
     15:   "  -p <p>, where <p> = number of rows of B.\n\n";

     17: #include <slepcsvd.h>

     19: /*
     20:    This example solves two GSVD problems for the bidiagonal matrices

     22:               |  1  2                     |       |  1                        |
     23:               |     1  2                  |       |  2  1                     |
     24:               |        1  2               |       |     2  1                  |
     25:          A1 = |          .  .             |  A2 = |       .  .                |
     26:               |             .  .          |       |          .  .             |
     27:               |                1  2       |       |             2  1          |
     28:               |                   1  2    |       |                2  1       |

     30:    with B = tril(ones(p,n))
     31:  */

     33: int main(int argc,char **argv)
     34: {
     35:   Mat            A1,A2,B;
     36:   SVD            svd;
     37:   PetscInt       m=15,n=20,p=21,Istart,Iend,i,j,d,col[2];
     38:   PetscScalar    valsa[] = { 1, 2 }, valsb[] = { 2, 1 };

     40:   PetscFunctionBeginUser;
     41:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     42:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
     43:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     44:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL));
     45:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized singular value decomposition, (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT "\n\n",m,p,n));

     47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     48:                      Generate the matrices
     49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     51:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A1));
     52:   PetscCall(MatSetSizes(A1,PETSC_DECIDE,PETSC_DECIDE,m,n));
     53:   PetscCall(MatSetFromOptions(A1));
     54:   PetscCall(MatGetOwnershipRange(A1,&Istart,&Iend));
     55:   for (i=Istart;i<Iend;i++) {
     56:     col[0]=i; col[1]=i+1;
     57:     if (i<n-1) PetscCall(MatSetValues(A1,1,&i,2,col,valsa,INSERT_VALUES));
     58:     else if (i==n-1) PetscCall(MatSetValue(A1,i,col[0],valsa[0],INSERT_VALUES));
     59:   }
     60:   PetscCall(MatAssemblyBegin(A1,MAT_FINAL_ASSEMBLY));
     61:   PetscCall(MatAssemblyEnd(A1,MAT_FINAL_ASSEMBLY));

     63:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A2));
     64:   PetscCall(MatSetSizes(A2,PETSC_DECIDE,PETSC_DECIDE,m,n));
     65:   PetscCall(MatSetFromOptions(A2));
     66:   PetscCall(MatGetOwnershipRange(A2,&Istart,&Iend));
     67:   for (i=Istart;i<Iend;i++) {
     68:     col[0]=i-1; col[1]=i;
     69:     if (i==0) PetscCall(MatSetValue(A2,i,col[1],valsb[1],INSERT_VALUES));
     70:     else if (i<n) PetscCall(MatSetValues(A2,1,&i,2,col,valsb,INSERT_VALUES));
     71:     else if (i==n) PetscCall(MatSetValue(A2,i,col[0],valsb[0],INSERT_VALUES));
     72:   }
     73:   PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY));
     74:   PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY));

     76:   PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
     77:   PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
     78:   PetscCall(MatSetFromOptions(B));
     79:   PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
     80:   d = PetscMax(0,n-p);
     81:   for (i=Istart;i<Iend;i++) {
     82:     for (j=PetscMax(0,i-5);j<=PetscMin(i,n-1);j++) PetscCall(MatSetValue(B,i,j+d,1.0,INSERT_VALUES));
     83:   }
     84:   PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
     85:   PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));

     87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     88:          Create the singular value solver, set options and solve
     89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     91:   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
     92:   PetscCall(SVDSetOperators(svd,A1,B));
     93:   PetscCall(SVDSetFromOptions(svd));
     94:   PetscCall(SVDSolve(svd));
     95:   PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));

     97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     98:                        Solve second problem
     99:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    101:   PetscCall(SVDSetOperators(svd,A2,B));
    102:   PetscCall(SVDSolve(svd));
    103:   PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));

    105:   /* Free work space */
    106:   PetscCall(SVDDestroy(&svd));
    107:   PetscCall(MatDestroy(&A1));
    108:   PetscCall(MatDestroy(&A2));
    109:   PetscCall(MatDestroy(&B));
    110:   PetscCall(SlepcFinalize());
    111:   return 0;
    112: }

    114: /*TEST

    116:    testset:
    117:       args: -svd_nsv 3
    118:       requires: !single
    119:       output_file: output/test16_1.out
    120:       test:
    121:          suffix: 1_lapack
    122:          args: -svd_type lapack
    123:       test:
    124:          suffix: 1_cross
    125:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
    126:       test:
    127:          suffix: 1_cyclic
    128:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
    129:       test:
    130:          suffix: 1_trlanczos
    131:          args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single lower}} -svd_trlanczos_ksp_rtol 1e-10
    132:          requires: double
    133:       test:
    134:          suffix: 1_trlanczos_par
    135:          nsize: 2
    136:          args: -svd_type trlanczos -ds_parallel {{redundant synchronized}}

    138:    testset:
    139:       args: -svd_nsv 3 -mat_type aijcusparse
    140:       requires: cuda !single
    141:       output_file: output/test16_1.out
    142:       test:
    143:          suffix: 2_cross
    144:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
    145:       test:
    146:          suffix: 2_cyclic
    147:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
    148:       test:
    149:          suffix: 2_trlanczos
    150:          args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single lower}} -svd_trlanczos_ksp_rtol 1e-10
    151:          requires: double

    153:    testset:
    154:       args: -svd_nsv 3 -mat_type aijhipsparse
    155:       requires: hip !single
    156:       output_file: output/test16_1.out
    157:       test:
    158:          suffix: 2_hip_cross
    159:          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
    160:       test:
    161:          suffix: 2_hip_cyclic
    162:          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
    163:       test:
    164:          suffix: 2_hip_trlanczos
    165:          args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single lower}} -svd_trlanczos_ksp_rtol 1e-10
    166:          requires: double

    168: TEST*/