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: }