Actual source code: test15.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 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,nloc,k,nc;
 22:   const PetscScalar *pX;
 23:   const char        *name;

 26:   MPI_Comm_rank(PetscObjectComm((PetscObject)X),&rank);
 27:   if (!rank) {
 28:     BVGetActiveColumns(X,NULL,&k);
 29:     BVGetSizes(X,&nloc,NULL,NULL);
 30:     BVGetNumConstraints(X,&nc);
 31:     PetscObjectGetName((PetscObject)X,&name);
 32:     PetscPrintf(PetscObjectComm((PetscObject)X),"First row of %s =\n",name);
 33:     BVGetArrayRead(X,&pX);
 34:     for (i=0;i<nc+k;i++) PetscPrintf(PetscObjectComm((PetscObject)X),"%g ",(double)PetscRealPart(pX[i*nloc]));
 35:     PetscPrintf(PetscObjectComm((PetscObject)X),"\n");
 36:     BVRestoreArrayRead(X,&pX);
 37:   }
 38:   return 0;
 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,nloc;
 46:   PetscReal      norm;
 47:   PetscScalar    alpha;
 48:   PetscViewer    view;
 49:   PetscBool      verbose;

 52:   SlepcInitialize(&argc,&argv,(char*)0,help);
 53:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 54:   PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
 55:   PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL);
 56:   PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
 57:   PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
 58:   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:   VecCreate(PETSC_COMM_WORLD,&t);
 62:   VecSetSizes(t,PETSC_DECIDE,n);
 63:   VecSetFromOptions(t);
 64:   VecGetLocalSize(t,&nloc);

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

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

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

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

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

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

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

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

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

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

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

149: /*TEST

151:    testset:
152:       nsize: 2
153:       output_file: output/test15_1.out
154:       test:
155:          suffix: 1
156:          args: -bv_type {{vecs contiguous svec mat}shared output}

158:       test:
159:          suffix: 1_cuda
160:          args: -bv_type svec -vec_type cuda
161:          requires: cuda

163:    testset:
164:       nsize: 2
165:       output_file: output/test15_2.out
166:       test:
167:          suffix: 2
168:          args: -nc 2 -bv_type {{vecs contiguous svec mat}shared output}
169:       test:
170:          suffix: 2_cuda
171:          args: -nc 2 -bv_type svec -vec_type cuda
172:          requires: cuda

174: TEST*/