Actual source code: ex5.c

  1: /*$Id: ex5.c,v 1.27 2001/08/07 21:30:08 bsmith Exp $*/
  2: 
  3: static char help[] = "Tests MatMult(), MatMultAdd(), MatMultTranspose().\n\
  4: Also MatMultTransposeAdd(), MatScale(), MatGetDiagonal(), and MatDiagonalScale().\n\n";

 6:  #include petscmat.h

 10: int main(int argc,char **args)
 11: {
 12:   Mat          C;
 13:   Vec          s,u,w,x,y,z;
 14:   int          ierr,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:   PetscTruth   flg;

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

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

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

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

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

 64:   /* Flush off proc Mat values and do more assembly */
 65:   MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);
 66:   for (i=rstart; i<rend; i++) {
 67:     for (j=0; j<n; j++) {
 68:       v=10*(i+1)+j+1;
 69:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
 70:     }
 71:   }
 72:   /* Try overlap Coomunication with the next stage XXXSetValues */
 73:   VecAssemblyEnd(z);
 74:   MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);

 76:   /* The Assembly for the second Stage */
 77:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 78:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 79:   VecAssemblyBegin(y);
 80:   VecAssemblyEnd(y);
 81:   MatScale(&alpha,C);
 82:   VecAssemblyBegin(u);
 83:   VecAssemblyEnd(u);

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

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

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

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

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

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

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

132:   /* -------------------- Test () MatDiagonalScale ------------------ */
133:   PetscOptionsHasName(PETSC_NULL,"-test_diagonalscale",&flg);
134:   if (flg) {
135:     MatDiagonalScale(C,x,y);
136:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
137:   }
138:   /* Free data structures */
139:   VecDestroy(u); VecDestroy(s);
140:   VecDestroy(w); VecDestroy(x);
141:   VecDestroy(y); VecDestroy(z);
142:   MatDestroy(C);

144:   PetscFinalize();
145:   return 0;
146: }