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: }