Actual source code: ex37.c
petsc-3.4.2 2013-07-02

      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>

     14: static PetscErrorCode GetISs(Vec vecs[],IS is[])
     15: {
     17:   PetscInt       rstart[2],rend[2];

     20:   VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
     21:   VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
     22:   ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
     23:   ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
     24:   return(0);
     25: }

     30: PetscErrorCode test_view(void)
     31: {
     32:   Vec            X, a,b;
     33:   Vec            c,d,e,f;
     34:   Vec            tmp_buf[2];
     35:   IS             tmp_is[2];
     36:   PetscInt       index;
     37:   PetscReal      val;
     38:   PetscInt       list[]={0,1,2};
     39:   PetscScalar    vals[]={0.720032,0.061794,0.0100223};
     41:   PetscBool      explcit = PETSC_FALSE;

     44:   PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);

     46:   VecCreate(PETSC_COMM_WORLD, &c);
     47:   VecSetSizes(c, PETSC_DECIDE, 3);
     48:   VecSetFromOptions(c);
     49:   VecDuplicate(c, &d);
     50:   VecDuplicate(c, &e);
     51:   VecDuplicate(c, &f);

     53:   VecSet(c, 1.0);
     54:   VecSet(d, 2.0);
     55:   VecSet(e, 3.0);
     56:   VecSetValues(f,3,list,vals,INSERT_VALUES);
     57:   VecAssemblyBegin(f);
     58:   VecAssemblyEnd(f);
     59:   VecScale(f, 10.0);

     61:   tmp_buf[0] = e;
     62:   tmp_buf[1] = f;
     63:   PetscOptionsGetBool(0,"-explicit_is",&explcit,0);
     64:   GetISs(tmp_buf,tmp_is);
     65:   VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
     66:   VecDestroy(&e);
     67:   VecDestroy(&f);
     68:   ISDestroy(&tmp_is[0]);
     69:   ISDestroy(&tmp_is[1]);

     71:   tmp_buf[0] = c;
     72:   tmp_buf[1] = d;
     73:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
     74:   VecDestroy(&c);   VecDestroy(&d);

     76:   tmp_buf[0] = a;
     77:   tmp_buf[1] = b;
     78:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
     79:   VecDestroy(&a);

     81:   VecAssemblyBegin(X);
     82:   VecAssemblyEnd(X);

     84:   VecMax(b, &index, &val);
     85:   PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n", val, index);

     87:   VecMin(b, &index, &val);
     88:   PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %d \n", val, index);

     90:   VecDestroy(&b);

     92:   VecMax(X, &index, &val);
     93:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", val, index);
     94:   VecMin(X, &index, &val);
     95:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", val, index);

     97:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);

     99:   VecDestroy(&X);
    100:   return(0);
    101: }

    103: #if 0
    106: PetscErrorCode test_vec_ops(void)
    107: {
    108:   Vec            X, a,b;
    109:   Vec            c,d,e,f;
    110:   PetscScalar    val;

    114:   PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);

    116:   VecCreate(PETSC_COMM_WORLD, &X);
    117:   VecSetSizes(X, 2, 2);
    118:   VecSetType(X, VECNEST);

    120:   VecCreate(PETSC_COMM_WORLD, &a);
    121:   VecSetSizes(a, 2, 2);
    122:   VecSetType(a, VECNEST);

    124:   VecCreate(PETSC_COMM_WORLD, &b);
    125:   VecSetSizes(b, 2, 2);
    126:   VecSetType(b, VECNEST);

    128:   /* assemble X */
    129:   VecNestSetSubVec(X, 0, a);
    130:   VecNestSetSubVec(X, 1, b);
    131:   VecAssemblyBegin(X);
    132:   VecAssemblyEnd(X);

    134:   VecCreate(PETSC_COMM_WORLD, &c);
    135:   VecSetSizes(c, 3, 3);
    136:   VecSetType(c, VECSEQ);
    137:   VecDuplicate(c, &d);
    138:   VecDuplicate(c, &e);
    139:   VecDuplicate(c, &f);

    141:   VecSet(c, 1.0);
    142:   VecSet(d, 2.0);
    143:   VecSet(e, 3.0);
    144:   VecSet(f, 4.0);

    146:   /* assemble a */
    147:   VecNestSetSubVec(a, 0, c);
    148:   VecNestSetSubVec(a, 1, d);
    149:   VecAssemblyBegin(a);
    150:   VecAssemblyEnd(a);

    152:   /* assemble b */
    153:   VecNestSetSubVec(b, 0, e);
    154:   VecNestSetSubVec(b, 1, f);
    155:   VecAssemblyBegin(b);
    156:   VecAssemblyEnd(b);

    158:   /*PetscPrintf(PETSC_COMM_WORLD, "X \n");*/
    159:   /*VecView(X, PETSC_VIEWER_STDOUT_WORLD);*/

    161:   VecDot(X,X, &val);
    162:   PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n", val);
    163:   return(0);
    164: }
    165: #endif

    169: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
    170: {
    171:   int            nproc;
    172:   Vec            v;
    173:   PetscInt       i;
    174:   PetscScalar    vx;

    177:   MPI_Comm_size(comm, &nproc);

    179:   VecCreate(comm, &v);
    180:   VecSetSizes(v, PETSC_DECIDE, length);
    181:   if (nproc == 1) { VecSetType(v, VECSEQ); }
    182:   else { VecSetType(v, VECMPI); }

    184:   for (i=0; i<length; i++) {
    185:     vx   = (PetscScalar)(start_value + i * stride);
    186:     VecSetValue(v, i, vx, INSERT_VALUES);
    187:   }
    188:   VecAssemblyBegin(v);
    189:   VecAssemblyEnd(v);

    191:   *_v = v;
    192:   return(0);
    193: }

    196: /*
    197: X = ([0,1,2,3], [10,12,14,16,18])
    198: Y = ([4,7,10,13], [5,6,7,8,9])

    200: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
    201: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])

    203: */
    206: PetscErrorCode test_axpy_dot_max(void)
    207: {
    208:   Vec            x1,y1, x2,y2;
    209:   Vec            tmp_buf[2];
    210:   Vec            X, Y;
    211:   PetscReal      real,real2;
    212:   PetscScalar    scalar;
    213:   PetscInt       index;

    217:   PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);

    219:   gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
    220:   gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);

    222:   gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
    223:   gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);

    225:   tmp_buf[0] = x1;
    226:   tmp_buf[1] = x2;
    227:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
    228:   VecAssemblyBegin(X);
    229:   VecAssemblyEnd(X);
    230:   VecDestroy(&x1);
    231:   VecDestroy(&x2);

    234:   tmp_buf[0] = y1;
    235:   tmp_buf[1] = y2;
    236:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
    237:   VecAssemblyBegin(Y);
    238:   VecAssemblyEnd(Y);
    239:   VecDestroy(&y1);
    240:   VecDestroy(&y2);

    243:   PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
    244:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
    245:   VecNestGetSubVec(Y, 0, &y1);
    246:   VecNestGetSubVec(Y, 1, &y2);
    247:   PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
    248:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
    249:   PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
    250:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
    251:   VecDot(X,Y, &scalar);

    253:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar));

    255:   VecDotNorm2(X,Y, &scalar, &real2);
    256:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi     norm2(Y) = %lf\n", PetscRealPart(scalar), PetscImaginaryPart(scalar), real2);

    259:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
    260:   VecNestGetSubVec(Y, 0, &y1);
    261:   VecNestGetSubVec(Y, 1, &y2);
    262:   PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
    263:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
    264:   PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
    265:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
    266:   VecDot(X,Y, &scalar);

    268:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar));
    269:   VecDotNorm2(X,Y, &scalar, &real2);
    270:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi     norm2(Y) = %lf\n", PetscRealPart(scalar), PetscImaginaryPart(scalar), real2);

    273:   VecMax(X, &index, &real);
    274:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", real, index);
    275:   VecMin(X, &index, &real);
    276:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", real, index);

    278:   VecDestroy(&X);
    279:   VecDestroy(&Y);
    280:   return(0);
    281: }

    284: int main(int argc, char **args)
    285: {
    286:   PetscInitialize(&argc, &args,(char*)0, help);

    288:   test_view();
    289:   test_axpy_dot_max();

    291:   PetscFinalize();
    292:   return 0;
    293: }