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, tol = PETSC_SQRT_MACHINE_EPSILON;
 17:   PetscBool      flg;

 19:   PetscInitialize(&argc,&args,(char*)0,help);
 20:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
 21:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 22:   n    = m;
 23:   PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
 24:   if (flg) n += 2;
 25:   PetscOptionsHasName(NULL,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 > tol) {
100:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)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 > tol) {
121:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)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,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: }