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 <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: PetscViewerSetFormat( PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL );
105: VecView( X, PETSC_VIEWER_STDOUT_WORLD );
107: VecDestroy(&X);
108: return(0);
109: }
111: #if 0
114: PetscErrorCode test_vec_ops( void )
115: {
116: Vec X, a,b;
117: Vec c,d,e,f;
118: PetscScalar val;
122: PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
124: VecCreate( PETSC_COMM_WORLD, &X );
125: VecSetSizes( X, 2, 2 );
126: VecSetType( X, VECNEST );
128: VecCreate( PETSC_COMM_WORLD, &a );
129: VecSetSizes( a, 2, 2 );
130: VecSetType( a, VECNEST );
132: VecCreate( PETSC_COMM_WORLD, &b );
133: VecSetSizes( b, 2, 2 );
134: VecSetType( b, VECNEST );
136: /* assemble X */
137: VecNestSetSubVec( X, 0, a );
138: VecNestSetSubVec( X, 1, b );
139: VecAssemblyBegin(X);
140: VecAssemblyEnd(X);
142: VecCreate( PETSC_COMM_WORLD, &c );
143: VecSetSizes( c, 3, 3 );
144: VecSetType( c, VECSEQ );
145: VecDuplicate( c, &d );
146: VecDuplicate( c, &e );
147: VecDuplicate( c, &f );
149: VecSet( c, 1.0 );
150: VecSet( d, 2.0 );
151: VecSet( e, 3.0 );
152: VecSet( f, 4.0 );
154: /* assemble a */
155: VecNestSetSubVec( a, 0, c );
156: VecNestSetSubVec( a, 1, d );
157: VecAssemblyBegin(a);
158: VecAssemblyEnd(a);
160: /* assemble b */
161: VecNestSetSubVec( b, 0, e );
162: VecNestSetSubVec( b, 1, f );
163: VecAssemblyBegin(b);
164: VecAssemblyEnd(b);
166: // PetscPrintf( PETSC_COMM_WORLD, "X \n");
167: // VecView( X, PETSC_VIEWER_STDOUT_WORLD );
169: VecDot( X,X, &val );
170: PetscPrintf( PETSC_COMM_WORLD, "X.X = %f \n", val );
172: return(0);
173: }
174: #endif
178: PetscErrorCode gen_test_vector( MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v )
179: {
180: int nproc;
181: Vec v;
182: PetscInt i;
183: PetscScalar vx;
186: MPI_Comm_size( comm, &nproc );
188: VecCreate( comm, &v );
189: VecSetSizes( v, PETSC_DECIDE, length );
190: if( nproc == 1 ) { VecSetType( v, VECSEQ ); }
191: else { VecSetType( v, VECMPI ); }
193: for( i=0; i<length; i++ ) {
194: vx = (PetscScalar)( start_value + i * stride );
195: VecSetValue( v, i, vx, INSERT_VALUES );
196: }
197: VecAssemblyBegin( v );
198: VecAssemblyEnd( v );
200: *_v = v;
202: return(0);
203: }
206: /*
207: X = ( [0,1,2,3], [10,12,14,16,18] )
208: Y = ( [4,7,10,13], [5,6,7,8,9] )
210: Y = aX + y = ( [4,8,12,16], (15,18,21,24,27] )
211: Y = aX + y = ( [4,9,14,19], (25,30,35,40,45] )
213: */
216: PetscErrorCode test_axpy_dot_max( void )
217: {
218: Vec x1,y1, x2,y2;
219: Vec tmp_buf[2];
220: Vec X, Y;
221: PetscReal real;
222: PetscScalar scalar;
223: PetscInt index;
227: PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME );
229: gen_test_vector( PETSC_COMM_WORLD, 4, 0, 1, &x1 );
230: gen_test_vector( PETSC_COMM_WORLD, 5, 10, 2, &x2 );
232: gen_test_vector( PETSC_COMM_WORLD, 4, 4, 3, &y1 );
233: gen_test_vector( PETSC_COMM_WORLD, 5, 5, 1, &y2 );
235: tmp_buf[0] = x1;
236: tmp_buf[1] = x2;
237: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&X);
238: VecAssemblyBegin(X);
239: VecAssemblyEnd(X);
240: VecDestroy(&x1);
241: VecDestroy(&x2);
244: tmp_buf[0] = y1;
245: tmp_buf[1] = y2;
246: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&Y);
247: VecAssemblyBegin(Y);
248: VecAssemblyEnd(Y);
249: VecDestroy(&y1);
250: VecDestroy(&y2);
253: PetscPrintf( PETSC_COMM_WORLD, "VecAXPY \n");
254: VecAXPY( Y, 1.0, X ); /* Y <- a X + Y */
255: VecNestGetSubVec( Y, 0, &y1 );
256: VecNestGetSubVec( Y, 1, &y2 );
257: PetscPrintf( PETSC_COMM_WORLD, "(1) y1 = \n" ); VecView( y1, PETSC_VIEWER_STDOUT_WORLD );
258: PetscPrintf( PETSC_COMM_WORLD, "(1) y2 = \n" ); VecView( y2, PETSC_VIEWER_STDOUT_WORLD );
259: VecDot( X,Y, &scalar );
261: PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar) );
264: VecAXPY( Y, 1.0, X ); /* Y <- a X + Y */
265: VecNestGetSubVec( Y, 0, &y1 );
266: VecNestGetSubVec( Y, 1, &y2 );
267: PetscPrintf( PETSC_COMM_WORLD, "(2) y1 = \n" ); VecView( y1, PETSC_VIEWER_STDOUT_WORLD );
268: PetscPrintf( PETSC_COMM_WORLD, "(2) y2 = \n" ); VecView( y2, PETSC_VIEWER_STDOUT_WORLD );
269: VecDot( X,Y, &scalar );
271: PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar) );
275: VecMax( X, &index, &real );
276: PetscPrintf( PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", real, index );
277: VecMin( X, &index, &real );
278: PetscPrintf( PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", real, index );
280: VecDestroy(&X);
281: VecDestroy(&Y);
283: return(0);
284: }
287: int main( int argc, char **args )
288: {
289: PetscInitialize( &argc, &args,(char *)0, help );
291: test_view();
292: test_axpy_dot_max();
294: PetscFinalize();
295: return 0;
296: }