1: 
  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:   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,&C);
 31:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);
 32:   MatSetFromOptions(C);
 33:   MatGetOwnershipRange(C,&rstart,&rend);
 34:   VecCreate(PETSC_COMM_WORLD,&x);
 35:   VecSetSizes(x,PETSC_DECIDE,m);
 36:   VecSetFromOptions(x);
 37:   VecDuplicate(x,&z);
 38:   VecDuplicate(x,&w);
 39:   VecCreate(PETSC_COMM_WORLD,&y);
 40:   VecSetSizes(y,PETSC_DECIDE,n);
 41:   VecSetFromOptions(y);
 42:   VecDuplicate(y,&u);
 43:   VecDuplicate(y,&s);
 44:   VecGetOwnershipRange(y,&vstart,&vend);

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

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

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

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

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

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

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

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

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

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

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

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