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