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,(char*)0,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}

157:       test:
158:          suffix: 1_cuda
159:          args: -bv_type {{svec mat}} -vec_type cuda
160:          requires: cuda

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

173: TEST*/