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

      2: static char help[] = "Tests MatMult(), MatMultAdd(), MatMultTranspose().\n\
      3: Also MatMultTransposeAdd(), MatScale(), MatGetDiagonal(), and MatDiagonalScale().\n\n";

      5: #include <petscmat.h>

      9: int main(int argc,char **args)
     10: {
     11:   Mat            C;
     12:   Vec            s,u,w,x,y,z;
     14:   PetscInt       i,j,m = 8,n,rstart,rend,vstart,vend;
     15:   PetscScalar    one = 1.0,negone = -1.0,v,alpha=0.1;
     16:   PetscReal      norm;
     17:   PetscBool      flg;

     19:   PetscInitialize(&argc,&args,(char*)0,help);
     20:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
     21:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
     22:   n    = m;
     23:   PetscOptionsHasName(NULL,"-rectA",&flg);
     24:   if (flg) n += 2;
     25:   PetscOptionsHasName(NULL,"-rectB",&flg);
     26:   if (flg) n -= 2;

     28:   /* ---------- Assemble matrix and vectors ----------- */

     30:   MatCreate(PETSC_COMM_WORLD,&C);
     31:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);
     32:   MatSetFromOptions(C);
     33:   MatSetUp(C);
     34:   MatGetOwnershipRange(C,&rstart,&rend);
     35:   VecCreate(PETSC_COMM_WORLD,&x);
     36:   VecSetSizes(x,PETSC_DECIDE,m);
     37:   VecSetFromOptions(x);
     38:   VecDuplicate(x,&z);
     39:   VecDuplicate(x,&w);
     40:   VecCreate(PETSC_COMM_WORLD,&y);
     41:   VecSetSizes(y,PETSC_DECIDE,n);
     42:   VecSetFromOptions(y);
     43:   VecDuplicate(y,&u);
     44:   VecDuplicate(y,&s);
     45:   VecGetOwnershipRange(y,&vstart,&vend);

     47:   /* Assembly */
     48:   for (i=rstart; i<rend; i++) {
     49:     v    = 100*(i+1);
     50:     VecSetValues(z,1,&i,&v,INSERT_VALUES);
     51:     for (j=0; j<n; j++) {
     52:       v    = 10*(i+1)+j+1;
     53:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
     54:     }
     55:   }

     57:   /* Flush off proc Vec values and do more assembly */
     58:   VecAssemblyBegin(z);
     59:   for (i=vstart; i<vend; i++) {
     60:     v    = one*((PetscReal)i);
     61:     VecSetValues(y,1,&i,&v,INSERT_VALUES);
     62:     v    = 100.0*i;
     63:     VecSetValues(u,1,&i,&v,INSERT_VALUES);
     64:   }

     66:   /* Flush off proc Mat values and do more assembly */
     67:   MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);
     68:   for (i=rstart; i<rend; i++) {
     69:     for (j=0; j<n; j++) {
     70:       v    = 10*(i+1)+j+1;
     71:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
     72:     }
     73:   }
     74:   /* Try overlap Coomunication with the next stage XXXSetValues */
     75:   VecAssemblyEnd(z);
     76:   MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);
     77:   CHKMEMQ;
     78:   /* The Assembly for the second Stage */
     79:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
     80:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
     81:   VecAssemblyBegin(y);
     82:   VecAssemblyEnd(y);
     83:   MatScale(C,alpha);
     84:   VecAssemblyBegin(u);
     85:   VecAssemblyEnd(u);

     87:   /* ------------ Test MatMult(), MatMultAdd()  ---------- */

     89:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n");
     90:   CHKMEMQ;
     91:   MatMult(C,y,x);
     92:   CHKMEMQ;
     93:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
     94:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n");
     95:   MatMultAdd(C,y,z,w);
     96:   VecAXPY(x,one,z);
     97:   VecAXPY(x,negone,w);
     98:   VecNorm(x,NORM_2,&norm);
     99:   if (norm > 1.e-8) {
    100:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %G\n",norm);
    101:   }

    103:   /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */

    105:   for (i=rstart; i<rend; i++) {
    106:     v    = one*((PetscReal)i);
    107:     VecSetValues(x,1,&i,&v,INSERT_VALUES);
    108:   }
    109:   VecAssemblyBegin(x);
    110:   VecAssemblyEnd(x);
    111:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n");
    112:   MatMultTranspose(C,x,y);
    113:   VecView(y,PETSC_VIEWER_STDOUT_WORLD);

    115:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");
    116:   MatMultTransposeAdd(C,x,u,s);
    117:   VecAXPY(y,one,u);
    118:   VecAXPY(y,negone,s);
    119:   VecNorm(y,NORM_2,&norm);
    120:   if (norm > 1.e-8) {
    121:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %G\n",norm);
    122:   }

    124:   /* -------------------- Test MatGetDiagonal() ------------------ */

    126:   PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n");
    127:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
    128:   VecSet(x,one);
    129:   MatGetDiagonal(C,x);
    130:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    131:   for (i=vstart; i<vend; i++) {
    132:     v    = one*((PetscReal)(i+1));
    133:     VecSetValues(y,1,&i,&v,INSERT_VALUES);
    134:   }

    136:   /* -------------------- Test () MatDiagonalScale ------------------ */
    137:   PetscOptionsHasName(NULL,"-test_diagonalscale",&flg);
    138:   if (flg) {
    139:     MatDiagonalScale(C,x,y);
    140:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
    141:   }
    142:   /* Free data structures */
    143:   VecDestroy(&u); VecDestroy(&s);
    144:   VecDestroy(&w); VecDestroy(&x);
    145:   VecDestroy(&y); VecDestroy(&z);
    146:   MatDestroy(&C);

    148:   PetscFinalize();
    149:   return 0;
    150: }