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