Actual source code: test7.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[] = "Test multiplication of a Mat times a BV.\n\n";

     13: #include <slepcbv.h>

     15: int main(int argc,char **argv)
     16: {
     17:   Vec            t,r,v;
     18:   Mat            B,Ymat;
     19:   BV             X,Y,Z=NULL,Zcopy=NULL;
     20:   PetscInt       i,j,m=10,n,k=5,rep=1,Istart,Iend,ld;
     21:   PetscScalar    *pZ;
     22:   PetscReal      norm;
     23:   PetscViewer    view;
     24:   PetscBool      flg,verbose,fromfile;
     25:   char           filename[PETSC_MAX_PATH_LEN];
     26:   PetscViewer    viewer;
     27:   BVMatMultType  vmm;

     29:   PetscFunctionBeginUser;
     30:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     31:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
     32:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-rep",&rep,NULL));
     33:   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
     34:   PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&fromfile));
     35:   PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
     36:   PetscCall(PetscObjectSetName((PetscObject)B,"B"));
     37:   if (fromfile) {
     38: #if defined(PETSC_USE_COMPLEX)
     39:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
     40: #else
     41:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
     42: #endif
     43:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
     44:     PetscCall(MatSetFromOptions(B));
     45:     PetscCall(MatLoad(B,viewer));
     46:     PetscCall(PetscViewerDestroy(&viewer));
     47:     PetscCall(MatGetSize(B,&m,&n));
     48:     PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
     49:   } else {
     50:     /* Create 1-D Laplacian matrix */
     51:     PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
     52:     PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg));
     53:     if (!flg) n = m;
     54:     PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n));
     55:     PetscCall(MatSetFromOptions(B));
     56:     PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
     57:     for (i=Istart;i<Iend;i++) {
     58:       if (i>0 && i-1<n) PetscCall(MatSetValue(B,i,i-1,-1.0,INSERT_VALUES));
     59:       if (i+1<n) PetscCall(MatSetValue(B,i,i+1,-1.0,INSERT_VALUES));
     60:       if (i<n) PetscCall(MatSetValue(B,i,i,2.0,INSERT_VALUES));
     61:     }
     62:     PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
     63:     PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
     64:   }

     66:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVMatMult (m=%" PetscInt_FMT ", n=%" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",m,n,k));
     67:   PetscCall(MatCreateVecs(B,&t,&r));

     69:   /* Create BV object X */
     70:   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
     71:   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
     72:   PetscCall(BVSetSizesFromVec(X,t,k));
     73:   PetscCall(BVSetMatMultMethod(X,BV_MATMULT_VECS));
     74:   PetscCall(BVSetFromOptions(X));
     75:   PetscCall(BVGetMatMultMethod(X,&vmm));
     76:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Using method: %s\n",BVMatMultTypes[vmm]));

     78:   /* Set up viewer */
     79:   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
     80:   if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));

     82:   /* Fill X entries */
     83:   for (j=0;j<k;j++) {
     84:     PetscCall(BVGetColumn(X,j,&v));
     85:     PetscCall(VecSet(v,0.0));
     86:     for (i=Istart;i<PetscMin(j+1,Iend);i++) PetscCall(VecSetValue(v,i,1.0,INSERT_VALUES));
     87:     PetscCall(VecAssemblyBegin(v));
     88:     PetscCall(VecAssemblyEnd(v));
     89:     PetscCall(BVRestoreColumn(X,j,&v));
     90:   }
     91:   if (verbose) {
     92:     PetscCall(MatView(B,view));
     93:     PetscCall(BVView(X,view));
     94:   }

     96:   /* Create BV object Y */
     97:   PetscCall(BVCreate(PETSC_COMM_WORLD,&Y));
     98:   PetscCall(PetscObjectSetName((PetscObject)Y,"Y"));
     99:   PetscCall(BVSetSizesFromVec(Y,r,k+4));
    100:   PetscCall(BVSetMatMultMethod(Y,BV_MATMULT_VECS));
    101:   PetscCall(BVSetFromOptions(Y));
    102:   PetscCall(BVSetActiveColumns(Y,2,k+2));

    104:   /* Test BVMatMult */
    105:   for (i=0;i<rep;i++) PetscCall(BVMatMult(X,B,Y));
    106:   if (verbose) PetscCall(BVView(Y,view));

    108:   if (fromfile) {
    109:     /* Test BVMatMultTranspose */
    110:     PetscCall(BVDuplicate(X,&Z));
    111:     PetscCall(BVSetRandom(Z));
    112:     for (i=0;i<rep;i++) PetscCall(BVMatMultTranspose(Z,B,Y));
    113:     if (verbose) {
    114:       PetscCall(BVView(Z,view));
    115:       PetscCall(BVView(Y,view));
    116:     }
    117:     PetscCall(BVDestroy(&Z));
    118:     PetscCall(BVMatMultTransposeColumn(Y,B,2));
    119:     if (verbose) PetscCall(BVView(Y,view));
    120:   }

    122:   /* Test BVGetMat/RestoreMat */
    123:   PetscCall(BVGetMat(Y,&Ymat));
    124:   PetscCall(PetscObjectSetName((PetscObject)Ymat,"Ymat"));
    125:   if (verbose) PetscCall(MatView(Ymat,view));
    126:   PetscCall(BVRestoreMat(Y,&Ymat));

    128:   if (!fromfile) {
    129:     /* Create BV object Z */
    130:     PetscCall(BVDuplicateResize(Y,k,&Z));
    131:     PetscCall(PetscObjectSetName((PetscObject)Z,"Z"));

    133:     /* Fill Z entries */
    134:     for (j=0;j<k;j++) {
    135:       PetscCall(BVGetColumn(Z,j,&v));
    136:       PetscCall(VecSet(v,0.0));
    137:       if (!Istart) PetscCall(VecSetValue(v,0,1.0,ADD_VALUES));
    138:       if (j<n && j>=Istart && j<Iend) PetscCall(VecSetValue(v,j,1.0,ADD_VALUES));
    139:       if (j+1<n && j>=Istart && j<Iend) PetscCall(VecSetValue(v,j+1,-1.0,ADD_VALUES));
    140:       PetscCall(VecAssemblyBegin(v));
    141:       PetscCall(VecAssemblyEnd(v));
    142:       PetscCall(BVRestoreColumn(Z,j,&v));
    143:     }
    144:     if (verbose) PetscCall(BVView(Z,view));

    146:     /* Save a copy of Z */
    147:     PetscCall(BVDuplicate(Z,&Zcopy));
    148:     PetscCall(BVCopy(Z,Zcopy));

    150:     /* Test BVMult, check result of previous operations */
    151:     PetscCall(BVMult(Z,-1.0,1.0,Y,NULL));
    152:     PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
    153:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm));
    154:   }

    156:   /* Test BVMatMultColumn, multiply Y(:,2), result in Y(:,3) */
    157:   if (m==n) {
    158:     PetscCall(BVMatMultColumn(Y,B,2));
    159:     if (verbose) PetscCall(BVView(Y,view));

    161:     if (!fromfile) {
    162:       /* Test BVGetArray, modify Z to match Y */
    163:       PetscCall(BVCopy(Zcopy,Z));
    164:       PetscCall(BVGetArray(Z,&pZ));
    165:       if (Istart==0) {
    166:         PetscCheck(Iend>2,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"First process must have at least 3 rows");
    167:         PetscCall(BVGetLeadingDimension(Z,&ld));
    168:         pZ[ld]   = 5.0;   /* modify 3 first entries of second column */
    169:         pZ[ld+1] = -4.0;
    170:         pZ[ld+2] = 1.0;
    171:       }
    172:       PetscCall(BVRestoreArray(Z,&pZ));
    173:       if (verbose) PetscCall(BVView(Z,view));

    175:       /* Check result again with BVMult */
    176:       PetscCall(BVMult(Z,-1.0,1.0,Y,NULL));
    177:       PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
    178:       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm));
    179:     }
    180:   }

    182:   PetscCall(BVDestroy(&Z));
    183:   PetscCall(BVDestroy(&Zcopy));
    184:   PetscCall(BVDestroy(&X));
    185:   PetscCall(BVDestroy(&Y));
    186:   PetscCall(MatDestroy(&B));
    187:   PetscCall(VecDestroy(&t));
    188:   PetscCall(VecDestroy(&r));
    189:   PetscCall(SlepcFinalize());
    190:   return 0;
    191: }

    193: /*TEST

    195:    testset:
    196:       output_file: output/test7_1.out
    197:       filter: grep -v "Using method"
    198:       test:
    199:          suffix: 1
    200:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
    201:       test:
    202:          suffix: 1_cuda
    203:          args: -bv_type {{svec mat}} -mat_type aijcusparse -bv_matmult vecs
    204:          requires: cuda
    205:       test:
    206:          suffix: 1_hip
    207:          args: -bv_type {{svec mat}} -mat_type aijhipsparse -bv_matmult vecs
    208:          requires: hip
    209:       test:
    210:          suffix: 1_mat
    211:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat

    213:    testset:
    214:       output_file: output/test7_2.out
    215:       filter: grep -v "Using method" | sed -e "s/error: -0\./error: 0./"
    216:       args: -m 34 -n 38 -k 9
    217:       nsize: 2
    218:       test:
    219:          suffix: 2
    220:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
    221:       test:
    222:          suffix: 2_cuda
    223:          args: -bv_type {{svec mat}} -mat_type aijcusparse -bv_matmult vecs
    224:          requires: cuda
    225:       test:
    226:          suffix: 2_hip
    227:          args: -bv_type {{svec mat}} -mat_type aijhipsparse -bv_matmult vecs
    228:          requires: hip
    229:       test:
    230:          suffix: 2_mat
    231:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat

    233:    testset:
    234:       output_file: output/test7_3.out
    235:       filter: grep -v "Using method"
    236:       args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -bv_reproducible_random
    237:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
    238:       nsize: 2
    239:       test:
    240:          suffix: 3
    241:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult {{vecs mat}}

    243: TEST*/