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