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

  5:  #include <petscmat.h>

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

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

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

 28:   MatCreate(PETSC_COMM_WORLD,&C);
 29:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);
 30:   MatSetFromOptions(C);
 31:   MatSetUp(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);
 75:   CHKMEMQ;
 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(C,alpha);
 82:   VecAssemblyBegin(u);
 83:   VecAssemblyEnd(u);

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

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

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

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

113:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");
114:   MatMultTransposeAdd(C,x,u,s);
115:   VecAXPY(y,one,u);
116:   VecAXPY(y,negone,s);
117:   VecNorm(y,NORM_2,&norm);
118:   if (norm > tol) {
119:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
120:   }

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

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

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

146:   PetscFinalize();
147:   return ierr;
148: }