Actual source code: ex9.c

petsc-3.12.5 2020-03-29
Report Typos and Errors

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