Actual source code: ex9.c

petsc-master 2019-08-17
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;
 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;

 30:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 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:   MatCreateVecs(A[0],&v,NULL);
 51:   VecDuplicate(v,&v2);

 53:   VecSet(x,1.0);
 54:   VecSet(y,0.0);
 55:   MatMult(A[1],x,z);
 56:   for (i = 2; i < nmat+1; i++) {
 57:     MatMultAdd(A[i],x,z,z);
 58:   }

 60:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B);
 61:   MatMultAdd(B,x,y,y);
 62:   VecAXPY(y,-1.0,z);
 63:   VecNorm(y,NORM_2,&rnorm);
 64:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
 65:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);
 66:   }

 68:   MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN); /* default */
 69:   MatCompositeMerge(B);
 70:   MatMult(B,x,y);
 71:   MatDestroy(&B);
 72:   VecAXPY(y,-1.0,z);
 73:   VecNorm(y,NORM_2,&rnorm);
 74:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
 75:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);
 76:   }

 78:   /*
 79:      Test n x n/2 multiplicative composite
 80:   */
 81:   VecSet(v,1.0);
 82:   MatMult(A[0],v,z);
 83:   for (i = 1; i < nmat; i++) {
 84:     MatMult(A[i],z,y);
 85:     VecCopy(y,z);
 86:   }

 88:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
 89:   MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
 90:   MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT);
 91:   MatSetFromOptions(B);
 92:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 93:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
 94:   MatMult(B,v,y);
 95:   MatDestroy(&B);
 96:   VecAXPY(y,-1.0,z);
 97:   VecNorm(y,NORM_2,&rnorm);
 98:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
 99:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
100:   }

102:   /*
103:      Test n/2 x n multiplicative composite
104:   */
105:   VecSet(x,1.0);
106:   MatMult(A[2],x,z);
107:   for (i = 3; i < nmat+1; i++) {
108:     MatMult(A[i],z,y);
109:     VecCopy(y,z);
110:   }
111:   MatMult(A[nmat+1],z,v);

113:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B);
114:   MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
115:   MatSetFromOptions(B);
116:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
117:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
118:   MatMult(B,x,v2);
119:   MatDestroy(&B);
120:   VecAXPY(v2,-1.0,v);
121:   VecNorm(v2,NORM_2,&rnorm);
122:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
123:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
124:   }

126:   /*
127:      Test get functions
128:   */
129:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
130:   MatCompositeGetNumberMat(B,&n);
131:   if (nmat != n) {
132:     PetscPrintf(PETSC_COMM_WORLD,"Error with GetNumberMat %D != %D\n",nmat,n);
133:   }
134:   MatCompositeGetMat(B,0,&A[nmat+2]);
135:   if (A[0] != A[nmat+2]) {
136:     PetscPrintf(PETSC_COMM_WORLD,"Error with GetMat\n");
137:   }
138:   MatCompositeGetType(B,&type);
139:   if (type != MAT_COMPOSITE_ADDITIVE) {
140:     PetscPrintf(PETSC_COMM_WORLD,"Error with GetType\n");
141:   }
142:   MatDestroy(&B);

144:   /*
145:      Free work space.  All PETSc objects should be destroyed when they
146:      are no longer needed.
147:   */
148:   VecDestroy(&x);
149:   VecDestroy(&y);
150:   VecDestroy(&v);
151:   VecDestroy(&v2);
152:   VecDestroy(&z);
153:   PetscRandomDestroy(&rctx);
154:   for (i = 0; i < nmat+2; i++) {
155:     MatDestroy(&A[i]);
156:   }
157:   PetscFree(A);

159:   PetscFinalize();
160:   return ierr;
161: }

163: /*TEST

165:    test:
166:       nsize: 2
167:       requires: double
168:       args: -mat_composite_merge {{0 1}shared output}

170: TEST*/