Actual source code: test15.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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:   PetscErrorCode    ierr;
 21:   PetscMPIInt       rank;
 22:   PetscInt          i,nloc,k,nc;
 23:   const PetscScalar *pX;
 24:   const char        *name;

 27:   MPI_Comm_rank(PetscObjectComm((PetscObject)X),&rank);
 28:   if (!rank) {
 29:     BVGetActiveColumns(X,NULL,&k);
 30:     BVGetSizes(X,&nloc,NULL,NULL);
 31:     BVGetNumConstraints(X,&nc);
 32:     PetscObjectGetName((PetscObject)X,&name);
 33:     PetscPrintf(PetscObjectComm((PetscObject)X),"First row of %s =\n",name);
 34:     BVGetArrayRead(X,&pX);
 35:     for (i=0;i<nc+k;i++) {
 36:       PetscPrintf(PetscObjectComm((PetscObject)X),"%g ",(double)PetscRealPart(pX[i*nloc]));
 37:     }
 38:     PetscPrintf(PetscObjectComm((PetscObject)X),"\n");
 39:     BVRestoreArrayRead(X,&pX);
 40:   }
 41:   return(0);
 42: }

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

 55:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 56:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 57:   PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
 58:   PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL);
 59:   PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
 60:   PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
 61:   PetscPrintf(PETSC_COMM_WORLD,"Test BVGetSplit (length %D, l=%D, k=%D, nc=%D).\n",n,l,k,nc);

 63:   /* Create template vector */
 64:   VecCreate(PETSC_COMM_WORLD,&t);
 65:   VecSetSizes(t,PETSC_DECIDE,n);
 66:   VecSetFromOptions(t);
 67:   VecGetLocalSize(t,&nloc);

 69:   /* Create BV object X */
 70:   BVCreate(PETSC_COMM_WORLD,&X);
 71:   PetscObjectSetName((PetscObject)X,"X");
 72:   BVSetSizesFromVec(X,t,k);
 73:   BVSetFromOptions(X);

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

 89:   /* Set up viewer */
 90:   PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
 91:   if (verbose) {
 92:     PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
 93:   }

 95:   /* Fill X entries */
 96:   for (j=0;j<k;j++) {
 97:     BVGetColumn(X,j,&v);
 98:     VecSet(v,0.0);
 99:     for (i=0;i<4;i++) {
100:       if (i+j<n) {
101:         VecSetValue(v,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES);
102:       }
103:     }
104:     VecAssemblyBegin(v);
105:     VecAssemblyEnd(v);
106:     BVRestoreColumn(X,j,&v);
107:   }
108:   if (verbose) {
109:     BVView(X,view);
110:   }

112:   /* Get split BVs */
113:   BVSetActiveColumns(X,l,k);
114:   BVGetSplit(X,&L,&R);
115:   PetscObjectSetName((PetscObject)L,"L");
116:   PetscObjectSetName((PetscObject)R,"R");

118:   if (verbose) {
119:     BVView(L,view);
120:     BVView(R,view);
121:   }

123:   /* Modify first column of R */
124:   BVGetColumn(R,0,&v);
125:   VecSet(v,-1.0);
126:   BVRestoreColumn(R,0,&v);

128:   /* Finished using the split BVs */
129:   BVRestoreSplit(X,&L,&R);
130:   PrintFirstRow(X);
131:   if (verbose) {
132:     BVView(X,view);
133:   }

135:   /* Get the left split BV only */
136:   BVGetSplit(X,&L,NULL);
137:   for (j=0;j<l;j++) {
138:     BVOrthogonalizeColumn(L,j,NULL,&norm,NULL);
139:     alpha = 1.0/norm;
140:     BVScaleColumn(L,j,alpha);
141:   }
142:   BVRestoreSplit(X,&L,NULL);
143:   PrintFirstRow(X);
144:   if (verbose) {
145:     BVView(X,view);
146:   }

148:   /* Now get the right split BV after changing the number of leading columns */
149:   BVSetActiveColumns(X,l-1,k);
150:   BVGetSplit(X,NULL,&R);
151:   BVGetColumn(R,0,&v);
152:   BVInsertVec(X,0,v);
153:   BVRestoreColumn(R,0,&v);
154:   BVRestoreSplit(X,NULL,&R);
155:   PrintFirstRow(X);
156:   if (verbose) {
157:     BVView(X,view);
158:   }

160:   BVDestroy(&X);
161:   VecDestroy(&t);
162:   SlepcFinalize();
163:   return ierr;
164: }

166: /*TEST

168:    test:
169:       suffix: 1
170:       nsize: 2
171:       args: -bv_type {{vecs contiguous svec mat}shared output}
172:       output_file: output/test15_1.out

174:    test:
175:       suffix: 1_cuda
176:       nsize: 2
177:       args: -bv_type svec -vec_type cuda
178:       requires: cuda
179:       output_file: output/test15_1.out

181:    test:
182:       suffix: 2
183:       nsize: 2
184:       args: -nc 2 -bv_type {{vecs contiguous svec mat}shared output}
185:       output_file: output/test15_2.out

187:    test:
188:       suffix: 2_cuda
189:       nsize: 2
190:       args: -nc 2 -bv_type svec -vec_type cuda
191:       requires: cuda
192:       output_file: output/test15_2.out

194: TEST*/