Actual source code: ex9.c
2: static char help[] = "Tests MatCreateComposite()\n\n";
4: /*T
5: Concepts: Mat^composite matrices
6: Processors: n
7: T*/
9: /*
10: Include "petscmat.h" so that we can use matrices.
11: automatically includes:
12: petscsys.h - base PETSc routines petscvec.h - vectors
13: petscmat.h - matrices
14: petscis.h - index sets petscviewer.h - viewers
15: */
16: #include <petscmat.h>
18: int main(int argc,char **args)
19: {
20: Mat *A,B; /* matrix */
21: PetscErrorCode ierr;
22: Vec x,y,v,v2,z,z2;
23: PetscReal rnorm;
24: PetscInt n = 20; /* size of the matrix */
25: PetscInt nmat = 3; /* number of matrices */
26: PetscInt i;
27: PetscRandom rctx;
28: MatCompositeType type;
29: PetscScalar scalings[5]={2,3,4,5,6};
31: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
32: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
33: PetscOptionsGetInt(NULL,NULL,"-nmat",&nmat,NULL);
35: /*
36: Create random matrices
37: */
38: PetscMalloc1(nmat+3,&A);
39: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
40: MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n/2,3,NULL,3,NULL,&A[0]);
41: for (i = 1; i < nmat+1; i++) {
42: MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,3,NULL,3,NULL,&A[i]);
43: }
44: MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n/2,n,3,NULL,3,NULL,&A[nmat+1]);
45: for (i = 0; i < nmat+2; i++) {
46: MatSetRandom(A[i],rctx);
47: }
49: MatCreateVecs(A[1],&x,&y);
50: VecDuplicate(y,&z);
51: VecDuplicate(z,&z2);
52: MatCreateVecs(A[0],&v,NULL);
53: VecDuplicate(v,&v2);
55: /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */
57: /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */
58: VecSet(x,1.0);
59: MatMult(A[1],x,z);
60: VecScale(z,scalings[1]);
61: for (i = 2; i < nmat+1; i++) {
62: MatMult(A[i],x,z2);
63: VecAXPY(z,scalings[i],z2);
64: }
66: /* Do MatMult using MatComposite and store the result in y */
67: VecSet(y,0.0);
68: MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B);
69: MatSetFromOptions(B);
70: MatCompositeSetScalings(B,&scalings[1]);
71: MatMultAdd(B,x,y,y);
73: /* Diff y and z */
74: VecAXPY(y,-1.0,z);
75: VecNorm(y,NORM_2,&rnorm);
76: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
77: PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);
78: }
80: /* Test MatCompositeMerge on ADDITIVE MatComposite */
81: MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN); /* default */
82: MatCompositeMerge(B);
83: MatMult(B,x,y);
84: MatDestroy(&B);
85: VecAXPY(y,-1.0,z);
86: VecNorm(y,NORM_2,&rnorm);
87: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
88: PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);
89: }
91: /*
92: Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings
93: */
95: /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */
96: VecSet(v,1.0);
97: MatMult(A[0],v,z);
98: VecScale(z,scalings[0]);
99: for (i = 1; i < nmat; i++) {
100: MatMult(A[i],z,y);
101: VecScale(y,scalings[i]);
102: VecCopy(y,z);
103: }
105: /* Do MatMult using MatComposite and store the result in y */
106: MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
107: MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
108: MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT);
109: MatSetFromOptions(B);
110: MatCompositeSetScalings(B,&scalings[0]);
111: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
112: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
113: MatMult(B,v,y);
114: MatDestroy(&B);
116: /* Diff y and z */
117: VecAXPY(y,-1.0,z);
118: VecNorm(y,NORM_2,&rnorm);
119: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
120: PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
121: }
123: /*
124: Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings
125: */
126: VecSet(x,1.0);
127: MatMult(A[2],x,z);
128: for (i = 3; i < nmat+1; i++) {
129: MatMult(A[i],z,y);
130: VecCopy(y,z);
131: }
132: MatMult(A[nmat+1],z,v);
134: MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B);
135: MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
136: MatSetFromOptions(B);
137: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
139: MatMult(B,x,v2);
140: MatDestroy(&B);
142: VecAXPY(v2,-1.0,v);
143: VecNorm(v2,NORM_2,&rnorm);
144: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
145: PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
146: }
148: /*
149: Test get functions
150: */
151: MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
152: MatCompositeGetNumberMat(B,&n);
153: if (nmat != n) {
154: PetscPrintf(PETSC_COMM_WORLD,"Error with GetNumberMat %D != %D\n",nmat,n);
155: }
156: MatCompositeGetMat(B,0,&A[nmat+2]);
157: if (A[0] != A[nmat+2]) {
158: PetscPrintf(PETSC_COMM_WORLD,"Error with GetMat\n");
159: }
160: MatCompositeGetType(B,&type);
161: if (type != MAT_COMPOSITE_ADDITIVE) {
162: PetscPrintf(PETSC_COMM_WORLD,"Error with GetType\n");
163: }
164: MatDestroy(&B);
166: /*
167: Free work space. All PETSc objects should be destroyed when they
168: are no longer needed.
169: */
170: VecDestroy(&x);
171: VecDestroy(&y);
172: VecDestroy(&v);
173: VecDestroy(&v2);
174: VecDestroy(&z);
175: VecDestroy(&z2);
176: PetscRandomDestroy(&rctx);
177: for (i = 0; i < nmat+2; i++) {
178: MatDestroy(&A[i]);
179: }
180: PetscFree(A);
182: PetscFinalize();
183: return ierr;
184: }
186: /*TEST
188: test:
189: nsize: 2
190: requires: double
191: args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output}
193: TEST*/