Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Test BVGetSplit().\n\n";
12 :
13 : #include <slepcbv.h>
14 :
15 : /*
16 : Print the first row of a BV
17 : */
18 48 : PetscErrorCode PrintFirstRow(BV X)
19 : {
20 48 : PetscMPIInt rank;
21 48 : PetscInt i,ld,k,nc;
22 48 : const PetscScalar *pX;
23 48 : const char *name;
24 :
25 48 : PetscFunctionBeginUser;
26 48 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)X),&rank));
27 48 : if (!rank) {
28 24 : PetscCall(BVGetActiveColumns(X,NULL,&k));
29 24 : PetscCall(BVGetLeadingDimension(X,&ld));
30 24 : PetscCall(BVGetNumConstraints(X,&nc));
31 24 : PetscCall(PetscObjectGetName((PetscObject)X,&name));
32 24 : PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"First row of %s =\n",name));
33 24 : PetscCall(BVGetArrayRead(X,&pX));
34 168 : for (i=0;i<nc+k;i++) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"%g ",(double)PetscRealPart(pX[i*ld])));
35 24 : PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"\n"));
36 24 : PetscCall(BVRestoreArrayRead(X,&pX));
37 : }
38 48 : PetscFunctionReturn(PETSC_SUCCESS);
39 : }
40 :
41 16 : int main(int argc,char **argv)
42 : {
43 16 : Vec t,v,*C;
44 16 : BV X,L,R;
45 16 : PetscInt i,j,n=10,k=5,l=3,nc=0;
46 16 : PetscReal norm;
47 16 : PetscScalar alpha;
48 16 : PetscViewer view;
49 16 : PetscBool verbose;
50 :
51 16 : PetscFunctionBeginUser;
52 16 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
53 16 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
54 16 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
55 16 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL));
56 16 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL));
57 16 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
58 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVGetSplit (length %" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT ", nc=%" PetscInt_FMT ").\n",n,l,k,nc));
59 :
60 : /* Create template vector */
61 16 : PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
62 16 : PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
63 16 : PetscCall(VecSetFromOptions(t));
64 :
65 : /* Create BV object X */
66 16 : PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
67 16 : PetscCall(PetscObjectSetName((PetscObject)X,"X"));
68 16 : PetscCall(BVSetSizesFromVec(X,t,k));
69 16 : PetscCall(BVSetFromOptions(X));
70 :
71 : /* Generate constraints and attach them to X */
72 16 : if (nc>0) {
73 8 : PetscCall(VecDuplicateVecs(t,nc,&C));
74 24 : for (j=0;j<nc;j++) {
75 40 : for (i=0;i<=j;i++) PetscCall(VecSetValue(C[j],nc-i+1,1.0,INSERT_VALUES));
76 16 : PetscCall(VecAssemblyBegin(C[j]));
77 16 : PetscCall(VecAssemblyEnd(C[j]));
78 : }
79 8 : PetscCall(BVInsertConstraints(X,&nc,C));
80 8 : PetscCall(VecDestroyVecs(nc,&C));
81 : }
82 :
83 : /* Set up viewer */
84 16 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
85 16 : if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
86 :
87 : /* Fill X entries */
88 96 : for (j=0;j<k;j++) {
89 80 : PetscCall(BVGetColumn(X,j,&v));
90 80 : PetscCall(VecSet(v,0.0));
91 400 : for (i=0;i<4;i++) {
92 320 : if (i+j<n) PetscCall(VecSetValue(v,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES));
93 : }
94 80 : PetscCall(VecAssemblyBegin(v));
95 80 : PetscCall(VecAssemblyEnd(v));
96 80 : PetscCall(BVRestoreColumn(X,j,&v));
97 : }
98 16 : if (verbose) PetscCall(BVView(X,view));
99 :
100 : /* Get split BVs */
101 16 : PetscCall(BVSetActiveColumns(X,l,k));
102 16 : PetscCall(BVGetSplit(X,&L,&R));
103 16 : PetscCall(PetscObjectSetName((PetscObject)L,"L"));
104 16 : PetscCall(PetscObjectSetName((PetscObject)R,"R"));
105 :
106 16 : if (verbose) {
107 0 : PetscCall(BVView(L,view));
108 0 : PetscCall(BVView(R,view));
109 : }
110 :
111 : /* Modify first column of R */
112 16 : PetscCall(BVGetColumn(R,0,&v));
113 16 : PetscCall(VecSet(v,-1.0));
114 16 : PetscCall(BVRestoreColumn(R,0,&v));
115 :
116 : /* Finished using the split BVs */
117 16 : PetscCall(BVRestoreSplit(X,&L,&R));
118 16 : PetscCall(PrintFirstRow(X));
119 16 : if (verbose) PetscCall(BVView(X,view));
120 :
121 : /* Get the left split BV only */
122 16 : PetscCall(BVGetSplit(X,&L,NULL));
123 64 : for (j=0;j<l;j++) {
124 48 : PetscCall(BVOrthogonalizeColumn(L,j,NULL,&norm,NULL));
125 48 : alpha = 1.0/norm;
126 48 : PetscCall(BVScaleColumn(L,j,alpha));
127 : }
128 16 : PetscCall(BVRestoreSplit(X,&L,NULL));
129 16 : PetscCall(PrintFirstRow(X));
130 16 : if (verbose) PetscCall(BVView(X,view));
131 :
132 : /* Now get the right split BV after changing the number of leading columns */
133 16 : PetscCall(BVSetActiveColumns(X,l-1,k));
134 16 : PetscCall(BVGetSplit(X,NULL,&R));
135 16 : PetscCall(BVGetColumn(R,0,&v));
136 16 : PetscCall(BVInsertVec(X,0,v));
137 16 : PetscCall(BVRestoreColumn(R,0,&v));
138 16 : PetscCall(BVRestoreSplit(X,NULL,&R));
139 16 : PetscCall(PrintFirstRow(X));
140 16 : if (verbose) PetscCall(BVView(X,view));
141 :
142 16 : PetscCall(BVDestroy(&X));
143 16 : PetscCall(VecDestroy(&t));
144 16 : PetscCall(SlepcFinalize());
145 : return 0;
146 : }
147 :
148 : /*TEST
149 :
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 :
157 : test:
158 : suffix: 1_cuda
159 : args: -bv_type {{svec mat}} -vec_type cuda
160 : requires: cuda
161 :
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
172 :
173 : TEST*/
|
