Actual source code: test15.c
slepc-3.11.2 2019-07-30
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*/