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: Vec x,y,v,v2,z,z2;
22: PetscReal rnorm;
23: PetscInt n = 20; /* size of the matrix */
24: PetscInt nmat = 3; /* number of matrices */
25: PetscInt i;
26: PetscRandom rctx;
27: MatCompositeType type;
28: PetscScalar scalings[5]={2,3,4,5,6};
30: PetscInitialize(&argc,&args,(char*)0,help);
31: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
32: PetscOptionsGetInt(NULL,NULL,"-nmat",&nmat,NULL);
34: /*
35: Create random matrices
36: */
37: PetscMalloc1(nmat+3,&A);
38: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
39: MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n/2,3,NULL,3,NULL,&A[0]);
40: for (i = 1; i < nmat+1; i++) {
41: MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,3,NULL,3,NULL,&A[i]);
42: }
43: MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n/2,n,3,NULL,3,NULL,&A[nmat+1]);
44: for (i = 0; i < nmat+2; i++) {
45: MatSetRandom(A[i],rctx);
46: }
48: MatCreateVecs(A[1],&x,&y);
49: VecDuplicate(y,&z);
50: VecDuplicate(z,&z2);
51: MatCreateVecs(A[0],&v,NULL);
52: VecDuplicate(v,&v2);
54: /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */
56: /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */
57: VecSet(x,1.0);
58: MatMult(A[1],x,z);
59: VecScale(z,scalings[1]);
60: for (i = 2; i < nmat+1; i++) {
61: MatMult(A[i],x,z2);
62: VecAXPY(z,scalings[i],z2);
63: }
65: /* Do MatMult using MatComposite and store the result in y */
66: VecSet(y,0.0);
67: MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B);
68: MatSetFromOptions(B);
69: MatCompositeSetScalings(B,&scalings[1]);
70: MatMultAdd(B,x,y,y);
72: /* Diff y and z */
73: VecAXPY(y,-1.0,z);
74: VecNorm(y,NORM_2,&rnorm);
75: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
76: PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);
77: }
79: /* Test MatCompositeMerge on ADDITIVE MatComposite */
80: MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN); /* default */
81: MatCompositeMerge(B);
82: MatMult(B,x,y);
83: MatDestroy(&B);
84: VecAXPY(y,-1.0,z);
85: VecNorm(y,NORM_2,&rnorm);
86: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
87: PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);
88: }
90: /*
91: Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings
92: */
94: /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */
95: VecSet(v,1.0);
96: MatMult(A[0],v,z);
97: VecScale(z,scalings[0]);
98: for (i = 1; i < nmat; i++) {
99: MatMult(A[i],z,y);
100: VecScale(y,scalings[i]);
101: VecCopy(y,z);
102: }
104: /* Do MatMult using MatComposite and store the result in y */
105: MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
106: MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
107: MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT);
108: MatSetFromOptions(B);
109: MatCompositeSetScalings(B,&scalings[0]);
110: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
112: MatMult(B,v,y);
113: MatDestroy(&B);
115: /* Diff y and z */
116: VecAXPY(y,-1.0,z);
117: VecNorm(y,NORM_2,&rnorm);
118: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
119: PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
120: }
122: /*
123: Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings
124: */
125: VecSet(x,1.0);
126: MatMult(A[2],x,z);
127: for (i = 3; i < nmat+1; i++) {
128: MatMult(A[i],z,y);
129: VecCopy(y,z);
130: }
131: MatMult(A[nmat+1],z,v);
133: MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B);
134: MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
135: MatSetFromOptions(B);
136: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
137: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
138: MatMult(B,x,v2);
139: MatDestroy(&B);
141: VecAXPY(v2,-1.0,v);
142: VecNorm(v2,NORM_2,&rnorm);
143: if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
144: PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
145: }
147: /*
148: Test get functions
149: */
150: MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
151: MatCompositeGetNumberMat(B,&n);
152: if (nmat != n) {
153: PetscPrintf(PETSC_COMM_WORLD,"Error with GetNumberMat %" PetscInt_FMT " != %" PetscInt_FMT "\n",nmat,n);
154: }
155: MatCompositeGetMat(B,0,&A[nmat+2]);
156: if (A[0] != A[nmat+2]) {
157: PetscPrintf(PETSC_COMM_WORLD,"Error with GetMat\n");
158: }
159: MatCompositeGetType(B,&type);
160: if (type != MAT_COMPOSITE_ADDITIVE) {
161: PetscPrintf(PETSC_COMM_WORLD,"Error with GetType\n");
162: }
163: MatDestroy(&B);
165: /*
166: Free work space. All PETSc objects should be destroyed when they
167: are no longer needed.
168: */
169: VecDestroy(&x);
170: VecDestroy(&y);
171: VecDestroy(&v);
172: VecDestroy(&v2);
173: VecDestroy(&z);
174: VecDestroy(&z2);
175: PetscRandomDestroy(&rctx);
176: for (i = 0; i < nmat+2; i++) {
177: MatDestroy(&A[i]);
178: }
179: PetscFree(A);
181: PetscFinalize();
182: return 0;
183: }
185: /*TEST
187: test:
188: nsize: 2
189: requires: double
190: args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output}
192: TEST*/