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: {
15: PetscInt rstart[2],rend[2];
18: VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
19: VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
20: ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
21: ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
22: return(0);
23: }
25: PetscErrorCode test_view(void)
26: {
27: Vec X, a,b;
28: Vec c,d,e,f;
29: Vec tmp_buf[2];
30: IS tmp_is[2];
31: PetscInt index;
32: PetscReal val;
33: PetscInt list[]={0,1,2};
34: PetscScalar vals[]={0.720032,0.061794,0.0100223};
36: PetscBool explcit = PETSC_FALSE;
39: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
41: VecCreate(PETSC_COMM_WORLD, &c);
42: VecSetSizes(c, PETSC_DECIDE, 3);
43: VecSetFromOptions(c);
44: VecDuplicate(c, &d);
45: VecDuplicate(c, &e);
46: VecDuplicate(c, &f);
48: VecSet(c, 1.0);
49: VecSet(d, 2.0);
50: VecSet(e, 3.0);
51: VecSetValues(f,3,list,vals,INSERT_VALUES);
52: VecAssemblyBegin(f);
53: VecAssemblyEnd(f);
54: VecScale(f, 10.0);
56: tmp_buf[0] = e;
57: tmp_buf[1] = f;
58: PetscOptionsGetBool(NULL,0,"-explicit_is",&explcit,0);
59: GetISs(tmp_buf,tmp_is);
60: VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
61: VecDestroy(&e);
62: VecDestroy(&f);
63: ISDestroy(&tmp_is[0]);
64: ISDestroy(&tmp_is[1]);
66: tmp_buf[0] = c;
67: tmp_buf[1] = d;
68: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
69: VecDestroy(&c); VecDestroy(&d);
71: tmp_buf[0] = a;
72: tmp_buf[1] = b;
73: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
74: VecDestroy(&a);
76: VecAssemblyBegin(X);
77: VecAssemblyEnd(X);
79: VecMax(b, &index, &val);
80: PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %D \n",(double) val, index);
82: VecMin(b, &index, &val);
83: PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %D \n",(double) val, index);
85: VecDestroy(&b);
87: VecMax(X, &index, &val);
88: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %D \n",(double) val, index);
89: VecMin(X, &index, &val);
90: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %D \n",(double) val, index);
92: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
94: VecDestroy(&X);
95: return(0);
96: }
98: #if 0
99: PetscErrorCode test_vec_ops(void)
100: {
101: Vec X, a,b;
102: Vec c,d,e,f;
103: PetscScalar val;
107: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
109: VecCreate(PETSC_COMM_WORLD, &X);
110: VecSetSizes(X, 2, 2);
111: VecSetType(X, VECNEST);
113: VecCreate(PETSC_COMM_WORLD, &a);
114: VecSetSizes(a, 2, 2);
115: VecSetType(a, VECNEST);
117: VecCreate(PETSC_COMM_WORLD, &b);
118: VecSetSizes(b, 2, 2);
119: VecSetType(b, VECNEST);
121: /* assemble X */
122: VecNestSetSubVec(X, 0, a);
123: VecNestSetSubVec(X, 1, b);
124: VecAssemblyBegin(X);
125: VecAssemblyEnd(X);
127: VecCreate(PETSC_COMM_WORLD, &c);
128: VecSetSizes(c, 3, 3);
129: VecSetType(c, VECSEQ);
130: VecDuplicate(c, &d);
131: VecDuplicate(c, &e);
132: VecDuplicate(c, &f);
134: VecSet(c, 1.0);
135: VecSet(d, 2.0);
136: VecSet(e, 3.0);
137: VecSet(f, 4.0);
139: /* assemble a */
140: VecNestSetSubVec(a, 0, c);
141: VecNestSetSubVec(a, 1, d);
142: VecAssemblyBegin(a);
143: VecAssemblyEnd(a);
145: /* assemble b */
146: VecNestSetSubVec(b, 0, e);
147: VecNestSetSubVec(b, 1, f);
148: VecAssemblyBegin(b);
149: VecAssemblyEnd(b);
151: VecDot(X,X, &val);
152: PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
153: return(0);
154: }
155: #endif
157: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
158: {
159: int size;
160: Vec v;
161: PetscInt i;
162: PetscScalar vx;
166: MPI_Comm_size(comm, &size);
167: VecCreate(comm, &v);
168: VecSetSizes(v, PETSC_DECIDE, length);
169: if (size == 1) { VecSetType(v, VECSEQ); }
170: else { VecSetType(v, VECMPI); }
172: for (i=0; i<length; i++) {
173: vx = (PetscScalar)(start_value + i * stride);
174: VecSetValue(v, i, vx, INSERT_VALUES);
175: }
176: VecAssemblyBegin(v);
177: VecAssemblyEnd(v);
179: *_v = v;
180: return(0);
181: }
183: /*
184: X = ([0,1,2,3], [10,12,14,16,18])
185: Y = ([4,7,10,13], [5,6,7,8,9])
187: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
188: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
190: */
191: PetscErrorCode test_axpy_dot_max(void)
192: {
193: Vec x1,y1, x2,y2;
194: Vec tmp_buf[2];
195: Vec X, Y;
196: PetscReal real,real2;
197: PetscScalar scalar;
198: PetscInt index;
202: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
204: gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
205: gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);
207: gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
208: gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);
210: tmp_buf[0] = x1;
211: tmp_buf[1] = x2;
212: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
213: VecAssemblyBegin(X);
214: VecAssemblyEnd(X);
215: VecDestroy(&x1);
216: VecDestroy(&x2);
218: tmp_buf[0] = y1;
219: tmp_buf[1] = y2;
220: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
221: VecAssemblyBegin(Y);
222: VecAssemblyEnd(Y);
223: VecDestroy(&y1);
224: VecDestroy(&y2);
226: PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
227: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
228: VecNestGetSubVec(Y, 0, &y1);
229: VecNestGetSubVec(Y, 1, &y2);
230: PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
231: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
232: PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
233: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
234: VecDot(X,Y, &scalar);
236: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
238: VecDotNorm2(X,Y, &scalar, &real2);
239: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
241: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
242: VecNestGetSubVec(Y, 0, &y1);
243: VecNestGetSubVec(Y, 1, &y2);
244: PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
245: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
246: PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
247: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
248: VecDot(X,Y, &scalar);
250: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
251: VecDotNorm2(X,Y, &scalar, &real2);
252: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
254: VecMax(X, &index, &real);
255: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %D \n",(double) real, index);
256: VecMin(X, &index, &real);
257: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %D \n",(double) real, index);
259: VecDestroy(&X);
260: VecDestroy(&Y);
261: return(0);
262: }
264: int main(int argc, char **args)
265: {
268: PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
269: test_view();
270: test_axpy_dot_max();
271: PetscFinalize();
272: return ierr;
273: }
275: /*TEST
277: test:
278: args: -explicit_is 0
280: test:
281: suffix: 2
282: args: -explicit_is 1
283: output_file: output/ex37_1.out
285: test:
286: suffix: 3
287: nsize: 2
288: args: -explicit_is 0
290: testset:
291: nsize: 2
292: args: -explicit_is 1
293: output_file: output/ex37_4.out
294: filter: grep -v -e "type: mpi" -e "type=mpi"
296: test:
297: suffix: 4
299: test:
300: requires: kokkos_kernels
301: suffix: kokkos
302: args: -vec_type kokkos
304: test:
305: requires: hip
306: suffix: hip
307: args: -vec_type hip
309: TEST*/