Actual source code: test15.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 BVGetSplit().\n\n";

     13: #include <slepcbv.h>

     15: /*
     16:    Print the first row of a BV
     17:  */
     18: PetscErrorCode PrintFirstRow(BV X)
     19: {
     20:   PetscMPIInt       rank;
     21:   PetscInt          i,ld,k,nc;
     22:   const PetscScalar *pX;
     23:   const char        *name;

     25:   PetscFunctionBeginUser;
     26:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)X),&rank));
     27:   if (!rank) {
     28:     PetscCall(BVGetActiveColumns(X,NULL,&k));
     29:     PetscCall(BVGetLeadingDimension(X,&ld));
     30:     PetscCall(BVGetNumConstraints(X,&nc));
     31:     PetscCall(PetscObjectGetName((PetscObject)X,&name));
     32:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"First row of %s =\n",name));
     33:     PetscCall(BVGetArrayRead(X,&pX));
     34:     for (i=0;i<nc+k;i++) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"%g ",(double)PetscRealPart(pX[i*ld])));
     35:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"\n"));
     36:     PetscCall(BVRestoreArrayRead(X,&pX));
     37:   }
     38:   PetscFunctionReturn(PETSC_SUCCESS);
     39: }

     41: int main(int argc,char **argv)
     42: {
     43:   Vec            t,v,*C;
     44:   BV             X,L,R;
     45:   PetscInt       i,j,n=10,k=5,l=3,nc=0;
     46:   PetscReal      norm;
     47:   PetscScalar    alpha;
     48:   PetscViewer    view;
     49:   PetscBool      verbose;

     51:   PetscFunctionBeginUser;
     52:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     53:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     54:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
     55:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL));
     56:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL));
     57:   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
     58:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVGetSplit (length %" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT ", nc=%" PetscInt_FMT ").\n",n,l,k,nc));

     60:   /* Create template vector */
     61:   PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
     62:   PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
     63:   PetscCall(VecSetFromOptions(t));

     65:   /* Create BV object X */
     66:   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
     67:   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
     68:   PetscCall(BVSetSizesFromVec(X,t,k));
     69:   PetscCall(BVSetFromOptions(X));

     71:   /* Generate constraints and attach them to X */
     72:   if (nc>0) {
     73:     PetscCall(VecDuplicateVecs(t,nc,&C));
     74:     for (j=0;j<nc;j++) {
     75:       for (i=0;i<=j;i++) PetscCall(VecSetValue(C[j],nc-i+1,1.0,INSERT_VALUES));
     76:       PetscCall(VecAssemblyBegin(C[j]));
     77:       PetscCall(VecAssemblyEnd(C[j]));
     78:     }
     79:     PetscCall(BVInsertConstraints(X,&nc,C));
     80:     PetscCall(VecDestroyVecs(nc,&C));
     81:   }

     83:   /* Set up viewer */
     84:   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
     85:   if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));

     87:   /* Fill X entries */
     88:   for (j=0;j<k;j++) {
     89:     PetscCall(BVGetColumn(X,j,&v));
     90:     PetscCall(VecSet(v,0.0));
     91:     for (i=0;i<4;i++) {
     92:       if (i+j<n) PetscCall(VecSetValue(v,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES));
     93:     }
     94:     PetscCall(VecAssemblyBegin(v));
     95:     PetscCall(VecAssemblyEnd(v));
     96:     PetscCall(BVRestoreColumn(X,j,&v));
     97:   }
     98:   if (verbose) PetscCall(BVView(X,view));

    100:   /* Get split BVs */
    101:   PetscCall(BVSetActiveColumns(X,l,k));
    102:   PetscCall(BVGetSplit(X,&L,&R));
    103:   PetscCall(PetscObjectSetName((PetscObject)L,"L"));
    104:   PetscCall(PetscObjectSetName((PetscObject)R,"R"));

    106:   if (verbose) {
    107:     PetscCall(BVView(L,view));
    108:     PetscCall(BVView(R,view));
    109:   }

    111:   /* Modify first column of R */
    112:   PetscCall(BVGetColumn(R,0,&v));
    113:   PetscCall(VecSet(v,-1.0));
    114:   PetscCall(BVRestoreColumn(R,0,&v));

    116:   /* Finished using the split BVs */
    117:   PetscCall(BVRestoreSplit(X,&L,&R));
    118:   PetscCall(PrintFirstRow(X));
    119:   if (verbose) PetscCall(BVView(X,view));

    121:   /* Get the left split BV only */
    122:   PetscCall(BVGetSplit(X,&L,NULL));
    123:   for (j=0;j<l;j++) {
    124:     PetscCall(BVOrthogonalizeColumn(L,j,NULL,&norm,NULL));
    125:     alpha = 1.0/norm;
    126:     PetscCall(BVScaleColumn(L,j,alpha));
    127:   }
    128:   PetscCall(BVRestoreSplit(X,&L,NULL));
    129:   PetscCall(PrintFirstRow(X));
    130:   if (verbose) PetscCall(BVView(X,view));

    132:   /* Now get the right split BV after changing the number of leading columns */
    133:   PetscCall(BVSetActiveColumns(X,l-1,k));
    134:   PetscCall(BVGetSplit(X,NULL,&R));
    135:   PetscCall(BVGetColumn(R,0,&v));
    136:   PetscCall(BVInsertVec(X,0,v));
    137:   PetscCall(BVRestoreColumn(R,0,&v));
    138:   PetscCall(BVRestoreSplit(X,NULL,&R));
    139:   PetscCall(PrintFirstRow(X));
    140:   if (verbose) PetscCall(BVView(X,view));

    142:   PetscCall(BVDestroy(&X));
    143:   PetscCall(VecDestroy(&t));
    144:   PetscCall(SlepcFinalize());
    145:   return 0;
    146: }

    148: /*TEST

    150:    testset:
    151:       nsize: 2
    152:       output_file: output/test15_1.out
    153:       test:
    154:          suffix: 1
    155:          args: -bv_type {{vecs contiguous svec mat}shared output}
    156:       test:
    157:          suffix: 1_cuda
    158:          args: -bv_type {{svec mat}} -vec_type cuda
    159:          requires: cuda
    160:       test:
    161:          suffix: 1_hip
    162:          args: -bv_type {{svec mat}} -vec_type hip
    163:          requires: hip

    165:    testset:
    166:       nsize: 2
    167:       output_file: output/test15_2.out
    168:       test:
    169:          suffix: 2
    170:          args: -nc 2 -bv_type {{vecs contiguous svec mat}shared output}
    171:       test:
    172:          suffix: 2_cuda
    173:          args: -nc 2 -bv_type {{svec mat}} -vec_type cuda
    174:          requires: cuda
    175:       test:
    176:          suffix: 2_hip
    177:          args: -nc 2 -bv_type {{svec mat}} -vec_type hip
    178:          requires: hip

    180: TEST*/