Actual source code: ex37.c
petsc-3.12.2 2019-11-22
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: }
26: PetscErrorCode test_view(void)
27: {
28: Vec X, a,b;
29: Vec c,d,e,f;
30: Vec tmp_buf[2];
31: IS tmp_is[2];
32: PetscInt index;
33: PetscReal val;
34: PetscInt list[]={0,1,2};
35: PetscScalar vals[]={0.720032,0.061794,0.0100223};
37: PetscBool explcit = PETSC_FALSE;
40: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
42: VecCreate(PETSC_COMM_WORLD, &c);
43: VecSetSizes(c, PETSC_DECIDE, 3);
44: VecSetFromOptions(c);
45: VecDuplicate(c, &d);
46: VecDuplicate(c, &e);
47: VecDuplicate(c, &f);
49: VecSet(c, 1.0);
50: VecSet(d, 2.0);
51: VecSet(e, 3.0);
52: VecSetValues(f,3,list,vals,INSERT_VALUES);
53: VecAssemblyBegin(f);
54: VecAssemblyEnd(f);
55: VecScale(f, 10.0);
57: tmp_buf[0] = e;
58: tmp_buf[1] = f;
59: PetscOptionsGetBool(NULL,0,"-explicit_is",&explcit,0);
60: GetISs(tmp_buf,tmp_is);
61: VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
62: VecDestroy(&e);
63: VecDestroy(&f);
64: ISDestroy(&tmp_is[0]);
65: ISDestroy(&tmp_is[1]);
67: tmp_buf[0] = c;
68: tmp_buf[1] = d;
69: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
70: VecDestroy(&c); VecDestroy(&d);
72: tmp_buf[0] = a;
73: tmp_buf[1] = b;
74: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
75: VecDestroy(&a);
77: VecAssemblyBegin(X);
78: VecAssemblyEnd(X);
80: VecMax(b, &index, &val);
81: PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n",(double) val, index);
83: VecMin(b, &index, &val);
84: PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %d \n",(double) val, index);
86: VecDestroy(&b);
88: VecMax(X, &index, &val);
89: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n",(double) val, index);
90: VecMin(X, &index, &val);
91: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n",(double) val, index);
93: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
95: VecDestroy(&X);
96: return(0);
97: }
99: #if 0
100: PetscErrorCode test_vec_ops(void)
101: {
102: Vec X, a,b;
103: Vec c,d,e,f;
104: PetscScalar val;
108: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
110: VecCreate(PETSC_COMM_WORLD, &X);
111: VecSetSizes(X, 2, 2);
112: VecSetType(X, VECNEST);
114: VecCreate(PETSC_COMM_WORLD, &a);
115: VecSetSizes(a, 2, 2);
116: VecSetType(a, VECNEST);
118: VecCreate(PETSC_COMM_WORLD, &b);
119: VecSetSizes(b, 2, 2);
120: VecSetType(b, VECNEST);
122: /* assemble X */
123: VecNestSetSubVec(X, 0, a);
124: VecNestSetSubVec(X, 1, b);
125: VecAssemblyBegin(X);
126: VecAssemblyEnd(X);
128: VecCreate(PETSC_COMM_WORLD, &c);
129: VecSetSizes(c, 3, 3);
130: VecSetType(c, VECSEQ);
131: VecDuplicate(c, &d);
132: VecDuplicate(c, &e);
133: VecDuplicate(c, &f);
135: VecSet(c, 1.0);
136: VecSet(d, 2.0);
137: VecSet(e, 3.0);
138: VecSet(f, 4.0);
140: /* assemble a */
141: VecNestSetSubVec(a, 0, c);
142: VecNestSetSubVec(a, 1, d);
143: VecAssemblyBegin(a);
144: VecAssemblyEnd(a);
146: /* assemble b */
147: VecNestSetSubVec(b, 0, e);
148: VecNestSetSubVec(b, 1, f);
149: VecAssemblyBegin(b);
150: VecAssemblyEnd(b);
152: VecDot(X,X, &val);
153: PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
154: return(0);
155: }
156: #endif
158: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
159: {
160: int size;
161: Vec v;
162: PetscInt i;
163: PetscScalar vx;
166: MPI_Comm_size(comm, &size);
168: VecCreate(comm, &v);
169: VecSetSizes(v, PETSC_DECIDE, length);
170: if (size == 1) { VecSetType(v, VECSEQ); }
171: else { VecSetType(v, VECMPI); }
173: for (i=0; i<length; i++) {
174: vx = (PetscScalar)(start_value + i * stride);
175: VecSetValue(v, i, vx, INSERT_VALUES);
176: }
177: VecAssemblyBegin(v);
178: VecAssemblyEnd(v);
180: *_v = v;
181: return(0);
182: }
185: /*
186: X = ([0,1,2,3], [10,12,14,16,18])
187: Y = ([4,7,10,13], [5,6,7,8,9])
189: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
190: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
192: */
193: PetscErrorCode test_axpy_dot_max(void)
194: {
195: Vec x1,y1, x2,y2;
196: Vec tmp_buf[2];
197: Vec X, Y;
198: PetscReal real,real2;
199: PetscScalar scalar;
200: PetscInt index;
204: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
206: gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
207: gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);
209: gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
210: gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);
212: tmp_buf[0] = x1;
213: tmp_buf[1] = x2;
214: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
215: VecAssemblyBegin(X);
216: VecAssemblyEnd(X);
217: VecDestroy(&x1);
218: VecDestroy(&x2);
221: tmp_buf[0] = y1;
222: tmp_buf[1] = y2;
223: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
224: VecAssemblyBegin(Y);
225: VecAssemblyEnd(Y);
226: VecDestroy(&y1);
227: VecDestroy(&y2);
229: PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
230: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
231: VecNestGetSubVec(Y, 0, &y1);
232: VecNestGetSubVec(Y, 1, &y2);
233: PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
234: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
235: PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
236: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
237: VecDot(X,Y, &scalar);
239: 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);
245: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
246: VecNestGetSubVec(Y, 0, &y1);
247: VecNestGetSubVec(Y, 1, &y2);
248: PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
249: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
250: PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
251: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
252: VecDot(X,Y, &scalar);
254: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
255: VecDotNorm2(X,Y, &scalar, &real2);
256: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
259: VecMax(X, &index, &real);
260: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n",(double) real, index);
261: VecMin(X, &index, &real);
262: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n",(double) real, index);
264: VecDestroy(&X);
265: VecDestroy(&Y);
266: return(0);
267: }
270: int main(int argc, char **args)
271: {
274: PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
275: test_view();
276: test_axpy_dot_max();
277: PetscFinalize();
278: return ierr;
279: }
282: /*TEST
284: test:
285: args: -explicit_is 0
287: test:
288: suffix: 2
289: args: -explicit_is 1
290: output_file: output/ex37_1.out
292: test:
293: suffix: 3
294: nsize: 2
295: args: -explicit_is 0
297: test:
298: suffix: 4
299: nsize: 2
300: args: -explicit_is 1
302: TEST*/