Actual source code: ex37.c
1: static char help[] = "Nest vector functionality.\n\n";
3: /*T
4: Concepts: vectors^block operators
5: Concepts: vectors^setting values
6: Concepts: vectors^local access to
7: Processors: n
8: T*/
10: #include <petscvec.h>
12: static PetscErrorCode GetISs(Vec vecs[],IS is[])
13: {
14: PetscInt rstart[2],rend[2];
16: VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
17: VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
18: ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
19: ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
20: return 0;
21: }
23: PetscErrorCode test_view(void)
24: {
25: Vec X, a,b;
26: Vec c,d,e,f;
27: Vec tmp_buf[2];
28: IS tmp_is[2];
29: PetscInt index;
30: PetscReal val;
31: PetscInt list[]={0,1,2};
32: PetscScalar vals[]={0.720032,0.061794,0.0100223};
33: PetscBool explcit = PETSC_FALSE;
35: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
37: VecCreate(PETSC_COMM_WORLD, &c);
38: VecSetSizes(c, PETSC_DECIDE, 3);
39: VecSetFromOptions(c);
40: VecDuplicate(c, &d);
41: VecDuplicate(c, &e);
42: VecDuplicate(c, &f);
44: VecSet(c, 1.0);
45: VecSet(d, 2.0);
46: VecSet(e, 3.0);
47: VecSetValues(f,3,list,vals,INSERT_VALUES);
48: VecAssemblyBegin(f);
49: VecAssemblyEnd(f);
50: VecScale(f, 10.0);
52: tmp_buf[0] = e;
53: tmp_buf[1] = f;
54: PetscOptionsGetBool(NULL,0,"-explicit_is",&explcit,0);
55: GetISs(tmp_buf,tmp_is);
56: VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
57: VecDestroy(&e);
58: VecDestroy(&f);
59: ISDestroy(&tmp_is[0]);
60: ISDestroy(&tmp_is[1]);
62: tmp_buf[0] = c;
63: tmp_buf[1] = d;
64: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
65: VecDestroy(&c)); PetscCall(VecDestroy(&d);
67: tmp_buf[0] = a;
68: tmp_buf[1] = b;
69: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
70: VecDestroy(&a);
72: VecAssemblyBegin(X);
73: VecAssemblyEnd(X);
75: VecMax(b, &index, &val);
76: PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
78: VecMin(b, &index, &val);
79: PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
81: VecDestroy(&b);
83: VecMax(X, &index, &val);
84: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
85: VecMin(X, &index, &val);
86: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
88: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
90: VecDestroy(&X);
91: return 0;
92: }
94: #if 0
95: PetscErrorCode test_vec_ops(void)
96: {
97: Vec X, a,b;
98: Vec c,d,e,f;
99: PetscScalar val;
101: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
103: VecCreate(PETSC_COMM_WORLD, &X);
104: VecSetSizes(X, 2, 2);
105: VecSetType(X, VECNEST);
107: VecCreate(PETSC_COMM_WORLD, &a);
108: VecSetSizes(a, 2, 2);
109: VecSetType(a, VECNEST);
111: VecCreate(PETSC_COMM_WORLD, &b);
112: VecSetSizes(b, 2, 2);
113: VecSetType(b, VECNEST);
115: /* assemble X */
116: VecNestSetSubVec(X, 0, a);
117: VecNestSetSubVec(X, 1, b);
118: VecAssemblyBegin(X);
119: VecAssemblyEnd(X);
121: VecCreate(PETSC_COMM_WORLD, &c);
122: VecSetSizes(c, 3, 3);
123: VecSetType(c, VECSEQ);
124: VecDuplicate(c, &d);
125: VecDuplicate(c, &e);
126: VecDuplicate(c, &f);
128: VecSet(c, 1.0);
129: VecSet(d, 2.0);
130: VecSet(e, 3.0);
131: VecSet(f, 4.0);
133: /* assemble a */
134: VecNestSetSubVec(a, 0, c);
135: VecNestSetSubVec(a, 1, d);
136: VecAssemblyBegin(a);
137: VecAssemblyEnd(a);
139: /* assemble b */
140: VecNestSetSubVec(b, 0, e);
141: VecNestSetSubVec(b, 1, f);
142: VecAssemblyBegin(b);
143: VecAssemblyEnd(b);
145: VecDot(X,X, &val);
146: PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
147: return 0;
148: }
149: #endif
151: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
152: {
153: int size;
154: Vec v;
155: PetscInt i;
156: PetscScalar vx;
158: MPI_Comm_size(comm, &size);
159: VecCreate(comm, &v);
160: VecSetSizes(v, PETSC_DECIDE, length);
161: if (size == 1) VecSetType(v, VECSEQ);
162: else VecSetType(v, VECMPI);
164: for (i=0; i<length; i++) {
165: vx = (PetscScalar)(start_value + i * stride);
166: VecSetValue(v, i, vx, INSERT_VALUES);
167: }
168: VecAssemblyBegin(v);
169: VecAssemblyEnd(v);
171: *_v = v;
172: return 0;
173: }
175: /*
176: X = ([0,1,2,3], [10,12,14,16,18])
177: Y = ([4,7,10,13], [5,6,7,8,9])
179: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
180: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
182: */
183: PetscErrorCode test_axpy_dot_max(void)
184: {
185: Vec x1,y1, x2,y2;
186: Vec tmp_buf[2];
187: Vec X, Y;
188: PetscReal real,real2;
189: PetscScalar scalar;
190: PetscInt index;
192: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
194: gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
195: gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);
197: gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
198: gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);
200: tmp_buf[0] = x1;
201: tmp_buf[1] = x2;
202: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
203: VecAssemblyBegin(X);
204: VecAssemblyEnd(X);
205: VecDestroy(&x1);
206: VecDestroy(&x2);
208: tmp_buf[0] = y1;
209: tmp_buf[1] = y2;
210: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
211: VecAssemblyBegin(Y);
212: VecAssemblyEnd(Y);
213: VecDestroy(&y1);
214: VecDestroy(&y2);
216: PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
217: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
218: VecNestGetSubVec(Y, 0, &y1);
219: VecNestGetSubVec(Y, 1, &y2);
220: PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
221: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
222: PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
223: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
224: VecDot(X,Y, &scalar);
226: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
228: VecDotNorm2(X,Y, &scalar, &real2);
229: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
231: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
232: VecNestGetSubVec(Y, 0, &y1);
233: VecNestGetSubVec(Y, 1, &y2);
234: PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
235: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
236: PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
237: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
238: VecDot(X,Y, &scalar);
240: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
241: VecDotNorm2(X,Y, &scalar, &real2);
242: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
244: VecMax(X, &index, &real);
245: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);
246: VecMin(X, &index, &real);
247: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);
249: VecDestroy(&X);
250: VecDestroy(&Y);
251: return 0;
252: }
254: int main(int argc, char **args)
255: {
257: PetscInitialize(&argc, &args,(char*)0, help);
258: test_view();
259: test_axpy_dot_max();
260: PetscFinalize();
261: return 0;
262: }
264: /*TEST
266: test:
267: args: -explicit_is 0
269: test:
270: suffix: 2
271: args: -explicit_is 1
272: output_file: output/ex37_1.out
274: test:
275: suffix: 3
276: nsize: 2
277: args: -explicit_is 0
279: testset:
280: nsize: 2
281: args: -explicit_is 1
282: output_file: output/ex37_4.out
283: filter: grep -v -e "type: mpi" -e "type=mpi"
285: test:
286: suffix: 4
288: test:
289: requires: kokkos_kernels
290: suffix: kokkos
291: args: -vec_type kokkos
293: test:
294: requires: hip
295: suffix: hip
296: args: -vec_type hip
298: TEST*/