Actual source code: ex37.c
1: static char help[] = "Nest vector functionality.\n\n";
3: #include <petscvec.h>
5: static PetscErrorCode GetISs(Vec vecs[],IS is[])
6: {
7: PetscInt rstart[2],rend[2];
9: VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
10: VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
11: ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
12: ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
13: return 0;
14: }
16: PetscErrorCode test_view(void)
17: {
18: Vec X, a,b;
19: Vec c,d,e,f;
20: Vec tmp_buf[2];
21: IS tmp_is[2];
22: PetscInt index;
23: PetscReal val;
24: PetscInt list[]={0,1,2};
25: PetscScalar vals[]={0.720032,0.061794,0.0100223};
26: PetscBool explcit = PETSC_FALSE;
28: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
30: VecCreate(PETSC_COMM_WORLD, &c);
31: VecSetSizes(c, PETSC_DECIDE, 3);
32: VecSetFromOptions(c);
33: VecDuplicate(c, &d);
34: VecDuplicate(c, &e);
35: VecDuplicate(c, &f);
37: VecSet(c, 1.0);
38: VecSet(d, 2.0);
39: VecSet(e, 3.0);
40: VecSetValues(f,3,list,vals,INSERT_VALUES);
41: VecAssemblyBegin(f);
42: VecAssemblyEnd(f);
43: VecScale(f, 10.0);
45: tmp_buf[0] = e;
46: tmp_buf[1] = f;
47: PetscOptionsGetBool(NULL,0,"-explicit_is",&explcit,0);
48: GetISs(tmp_buf,tmp_is);
49: VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
50: VecDestroy(&e);
51: VecDestroy(&f);
52: ISDestroy(&tmp_is[0]);
53: ISDestroy(&tmp_is[1]);
55: tmp_buf[0] = c;
56: tmp_buf[1] = d;
57: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
58: VecDestroy(&c)); PetscCall(VecDestroy(&d);
60: tmp_buf[0] = a;
61: tmp_buf[1] = b;
62: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
63: VecDestroy(&a);
65: VecAssemblyBegin(X);
66: VecAssemblyEnd(X);
68: VecMax(b, &index, &val);
69: PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
71: VecMin(b, &index, &val);
72: PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
74: VecDestroy(&b);
76: VecMax(X, &index, &val);
77: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
78: VecMin(X, &index, &val);
79: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
81: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
83: VecDestroy(&X);
84: return 0;
85: }
87: #if 0
88: PetscErrorCode test_vec_ops(void)
89: {
90: Vec X, a,b;
91: Vec c,d,e,f;
92: PetscScalar val;
94: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
96: VecCreate(PETSC_COMM_WORLD, &X);
97: VecSetSizes(X, 2, 2);
98: VecSetType(X, VECNEST);
100: VecCreate(PETSC_COMM_WORLD, &a);
101: VecSetSizes(a, 2, 2);
102: VecSetType(a, VECNEST);
104: VecCreate(PETSC_COMM_WORLD, &b);
105: VecSetSizes(b, 2, 2);
106: VecSetType(b, VECNEST);
108: /* assemble X */
109: VecNestSetSubVec(X, 0, a);
110: VecNestSetSubVec(X, 1, b);
111: VecAssemblyBegin(X);
112: VecAssemblyEnd(X);
114: VecCreate(PETSC_COMM_WORLD, &c);
115: VecSetSizes(c, 3, 3);
116: VecSetType(c, VECSEQ);
117: VecDuplicate(c, &d);
118: VecDuplicate(c, &e);
119: VecDuplicate(c, &f);
121: VecSet(c, 1.0);
122: VecSet(d, 2.0);
123: VecSet(e, 3.0);
124: VecSet(f, 4.0);
126: /* assemble a */
127: VecNestSetSubVec(a, 0, c);
128: VecNestSetSubVec(a, 1, d);
129: VecAssemblyBegin(a);
130: VecAssemblyEnd(a);
132: /* assemble b */
133: VecNestSetSubVec(b, 0, e);
134: VecNestSetSubVec(b, 1, f);
135: VecAssemblyBegin(b);
136: VecAssemblyEnd(b);
138: VecDot(X,X, &val);
139: PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
140: return 0;
141: }
142: #endif
144: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
145: {
146: int size;
147: Vec v;
148: PetscInt i;
149: PetscScalar vx;
151: MPI_Comm_size(comm, &size);
152: VecCreate(comm, &v);
153: VecSetSizes(v, PETSC_DECIDE, length);
154: if (size == 1) VecSetType(v, VECSEQ);
155: else VecSetType(v, VECMPI);
157: for (i=0; i<length; i++) {
158: vx = (PetscScalar)(start_value + i * stride);
159: VecSetValue(v, i, vx, INSERT_VALUES);
160: }
161: VecAssemblyBegin(v);
162: VecAssemblyEnd(v);
164: *_v = v;
165: return 0;
166: }
168: /*
169: X = ([0,1,2,3], [10,12,14,16,18])
170: Y = ([4,7,10,13], [5,6,7,8,9])
172: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
173: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
175: */
176: PetscErrorCode test_axpy_dot_max(void)
177: {
178: Vec x1,y1, x2,y2;
179: Vec tmp_buf[2];
180: Vec X, Y;
181: PetscReal real,real2;
182: PetscScalar scalar;
183: PetscInt index;
185: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
187: gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
188: gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);
190: gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
191: gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);
193: tmp_buf[0] = x1;
194: tmp_buf[1] = x2;
195: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
196: VecAssemblyBegin(X);
197: VecAssemblyEnd(X);
198: VecDestroy(&x1);
199: VecDestroy(&x2);
201: tmp_buf[0] = y1;
202: tmp_buf[1] = y2;
203: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
204: VecAssemblyBegin(Y);
205: VecAssemblyEnd(Y);
206: VecDestroy(&y1);
207: VecDestroy(&y2);
209: PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
210: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
211: VecNestGetSubVec(Y, 0, &y1);
212: VecNestGetSubVec(Y, 1, &y2);
213: PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
214: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
215: PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
216: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
217: VecDot(X,Y, &scalar);
219: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
221: VecDotNorm2(X,Y, &scalar, &real2);
222: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
224: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
225: VecNestGetSubVec(Y, 0, &y1);
226: VecNestGetSubVec(Y, 1, &y2);
227: PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
228: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
229: PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
230: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
231: VecDot(X,Y, &scalar);
233: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
234: VecDotNorm2(X,Y, &scalar, &real2);
235: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
237: VecMax(X, &index, &real);
238: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);
239: VecMin(X, &index, &real);
240: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);
242: VecDestroy(&X);
243: VecDestroy(&Y);
244: return 0;
245: }
247: int main(int argc, char **args)
248: {
250: PetscInitialize(&argc, &args,(char*)0, help);
251: test_view();
252: test_axpy_dot_max();
253: PetscFinalize();
254: return 0;
255: }
257: /*TEST
259: test:
260: args: -explicit_is 0
262: test:
263: suffix: 2
264: args: -explicit_is 1
265: output_file: output/ex37_1.out
267: test:
268: suffix: 3
269: nsize: 2
270: args: -explicit_is 0
272: testset:
273: nsize: 2
274: args: -explicit_is 1
275: output_file: output/ex37_4.out
276: filter: grep -v -e "type: mpi" -e "type=mpi"
278: test:
279: suffix: 4
281: test:
282: requires: kokkos_kernels
283: suffix: kokkos
284: args: -vec_type kokkos
286: test:
287: requires: hip
288: suffix: hip
289: args: -vec_type hip
291: TEST*/