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*/