Actual source code: ex37.c
petsc-3.3-p6 2013-02-11
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 <stdio.h>
11: #include <stdlib.h>
12: #include <math.h>
14: #include <petsc.h>
15: #include <petscvec.h>
16: #include <petscmat.h>
17: #include <petscksp.h>
21: static PetscErrorCode GetISs(Vec vecs[],IS is[])
22: {
24: PetscInt rstart[2],rend[2];
27: VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
28: VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
29: ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
30: ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
31: return(0);
32: }
37: PetscErrorCode test_view( void )
38: {
39: Vec X, a,b;
40: Vec c,d,e,f;
41: Vec tmp_buf[2];
42: IS tmp_is[2];
43: PetscInt index;
44: PetscReal val;
45: PetscInt list[]={0,1,2};
46: PetscScalar vals[]={0.720032,0.061794,0.0100223};
48: PetscBool explcit = PETSC_FALSE;
51: PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME );
53: VecCreate( PETSC_COMM_WORLD, &c );
54: VecSetSizes( c, PETSC_DECIDE, 3 );
55: VecSetFromOptions( c );
56: VecDuplicate( c, &d );
57: VecDuplicate( c, &e );
58: VecDuplicate( c, &f );
60: VecSet( c, 1.0 );
61: VecSet( d, 2.0 );
62: VecSet( e, 3.0 );
63: VecSetValues(f,3,list,vals,INSERT_VALUES);
64: VecAssemblyBegin(f);
65: VecAssemblyEnd(f);
66: VecScale( f, 10.0 );
68: tmp_buf[0] = e;
69: tmp_buf[1] = f;
70: PetscOptionsGetBool(0,"-explicit_is",&explcit,0);
71: GetISs(tmp_buf,tmp_is);
72: VecCreateNest(PETSC_COMM_WORLD,2,explcit?tmp_is:PETSC_NULL,tmp_buf,&b);
73: VecDestroy(&e);
74: VecDestroy(&f);
75: ISDestroy(&tmp_is[0]);
76: ISDestroy(&tmp_is[1]);
78: tmp_buf[0] = c;
79: tmp_buf[1] = d;
80: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&a);
81: VecDestroy(&c); VecDestroy(&d);
83: tmp_buf[0] = a;
84: tmp_buf[1] = b;
85: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&X);
86: VecDestroy(&a);
88: VecAssemblyBegin(X);
89: VecAssemblyEnd(X);
91: VecMax( b, &index, &val );
92: PetscPrintf( PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n", val, index );
94: VecMin( b, &index, &val );
95: PetscPrintf( PETSC_COMM_WORLD, "(min-b) = %f : index = %d \n", val, index );
97: VecDestroy(&b);
99: VecMax( X, &index, &val );
100: PetscPrintf( PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", val, index );
101: VecMin( X, &index, &val );
102: PetscPrintf( PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", val, index );
104: VecView( X, PETSC_VIEWER_STDOUT_WORLD );
106: VecDestroy(&X);
107: return(0);
108: }
110: #if 0
113: PetscErrorCode test_vec_ops( void )
114: {
115: Vec X, a,b;
116: Vec c,d,e,f;
117: PetscScalar val;
121: PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
123: VecCreate( PETSC_COMM_WORLD, &X );
124: VecSetSizes( X, 2, 2 );
125: VecSetType( X, VECNEST );
127: VecCreate( PETSC_COMM_WORLD, &a );
128: VecSetSizes( a, 2, 2 );
129: VecSetType( a, VECNEST );
131: VecCreate( PETSC_COMM_WORLD, &b );
132: VecSetSizes( b, 2, 2 );
133: VecSetType( b, VECNEST );
135: /* assemble X */
136: VecNestSetSubVec( X, 0, a );
137: VecNestSetSubVec( X, 1, b );
138: VecAssemblyBegin(X);
139: VecAssemblyEnd(X);
141: VecCreate( PETSC_COMM_WORLD, &c );
142: VecSetSizes( c, 3, 3 );
143: VecSetType( c, VECSEQ );
144: VecDuplicate( c, &d );
145: VecDuplicate( c, &e );
146: VecDuplicate( c, &f );
148: VecSet( c, 1.0 );
149: VecSet( d, 2.0 );
150: VecSet( e, 3.0 );
151: VecSet( f, 4.0 );
153: /* assemble a */
154: VecNestSetSubVec( a, 0, c );
155: VecNestSetSubVec( a, 1, d );
156: VecAssemblyBegin(a);
157: VecAssemblyEnd(a);
159: /* assemble b */
160: VecNestSetSubVec( b, 0, e );
161: VecNestSetSubVec( b, 1, f );
162: VecAssemblyBegin(b);
163: VecAssemblyEnd(b);
165: // PetscPrintf( PETSC_COMM_WORLD, "X \n");
166: // VecView( X, PETSC_VIEWER_STDOUT_WORLD );
168: VecDot( X,X, &val );
169: PetscPrintf( PETSC_COMM_WORLD, "X.X = %f \n", val );
171: return(0);
172: }
173: #endif
177: PetscErrorCode gen_test_vector( MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v )
178: {
179: int nproc;
180: Vec v;
181: PetscInt i;
182: PetscScalar vx;
185: MPI_Comm_size( comm, &nproc );
187: VecCreate( comm, &v );
188: VecSetSizes( v, PETSC_DECIDE, length );
189: if( nproc == 1 ) { VecSetType( v, VECSEQ ); }
190: else { VecSetType( v, VECMPI ); }
192: for( i=0; i<length; i++ ) {
193: vx = (PetscScalar)( start_value + i * stride );
194: VecSetValue( v, i, vx, INSERT_VALUES );
195: }
196: VecAssemblyBegin( v );
197: VecAssemblyEnd( v );
199: *_v = v;
201: return(0);
202: }
205: /*
206: X = ( [0,1,2,3], [10,12,14,16,18] )
207: Y = ( [4,7,10,13], [5,6,7,8,9] )
209: Y = aX + y = ( [4,8,12,16], (15,18,21,24,27] )
210: Y = aX + y = ( [4,9,14,19], (25,30,35,40,45] )
212: */
215: PetscErrorCode test_axpy_dot_max( void )
216: {
217: Vec x1,y1, x2,y2;
218: Vec tmp_buf[2];
219: Vec X, Y;
220: PetscReal real;
221: PetscScalar scalar, scalar2;
222: PetscInt index;
226: PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME );
228: gen_test_vector( PETSC_COMM_WORLD, 4, 0, 1, &x1 );
229: gen_test_vector( PETSC_COMM_WORLD, 5, 10, 2, &x2 );
231: gen_test_vector( PETSC_COMM_WORLD, 4, 4, 3, &y1 );
232: gen_test_vector( PETSC_COMM_WORLD, 5, 5, 1, &y2 );
234: tmp_buf[0] = x1;
235: tmp_buf[1] = x2;
236: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&X);
237: VecAssemblyBegin(X);
238: VecAssemblyEnd(X);
239: VecDestroy(&x1);
240: VecDestroy(&x2);
243: tmp_buf[0] = y1;
244: tmp_buf[1] = y2;
245: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&Y);
246: VecAssemblyBegin(Y);
247: VecAssemblyEnd(Y);
248: VecDestroy(&y1);
249: VecDestroy(&y2);
252: PetscPrintf( PETSC_COMM_WORLD, "VecAXPY \n");
253: VecAXPY( Y, 1.0, X ); /* Y <- a X + Y */
254: VecNestGetSubVec( Y, 0, &y1 );
255: VecNestGetSubVec( Y, 1, &y2 );
256: PetscPrintf( PETSC_COMM_WORLD, "(1) y1 = \n" ); VecView( y1, PETSC_VIEWER_STDOUT_WORLD );
257: PetscPrintf( PETSC_COMM_WORLD, "(1) y2 = \n" ); VecView( y2, PETSC_VIEWER_STDOUT_WORLD );
258: VecDot( X,Y, &scalar );
260: PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar) );
262: VecDotNorm2( X,Y, &scalar, &scalar2 );
263: PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf + %lfi\n", PetscRealPart(scalar), PetscImaginaryPart(scalar), PetscRealPart(scalar2), PetscImaginaryPart(scalar2) );
266: VecAXPY( Y, 1.0, X ); /* Y <- a X + Y */
267: VecNestGetSubVec( Y, 0, &y1 );
268: VecNestGetSubVec( Y, 1, &y2 );
269: PetscPrintf( PETSC_COMM_WORLD, "(2) y1 = \n" ); VecView( y1, PETSC_VIEWER_STDOUT_WORLD );
270: PetscPrintf( PETSC_COMM_WORLD, "(2) y2 = \n" ); VecView( y2, PETSC_VIEWER_STDOUT_WORLD );
271: VecDot( X,Y, &scalar );
273: PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar) );
274: VecDotNorm2( X,Y, &scalar, &scalar2 );
275: PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf + %lfi\n", PetscRealPart(scalar), PetscImaginaryPart(scalar), PetscRealPart(scalar2), PetscImaginaryPart(scalar2) );
278: VecMax( X, &index, &real );
279: PetscPrintf( PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", real, index );
280: VecMin( X, &index, &real );
281: PetscPrintf( PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", real, index );
283: VecDestroy(&X);
284: VecDestroy(&Y);
286: return(0);
287: }
290: int main( int argc, char **args )
291: {
292: PetscInitialize( &argc, &args,(char *)0, help );
294: test_view();
295: test_axpy_dot_max();
297: PetscFinalize();
298: return 0;
299: }