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

 20: int main(int argc,char **args)
 21: {
 22:   Mat                   A[3],B;                /* matrix */
 23:   PetscViewer           fd;               /* viewer */
 24:   char                  file[PETSC_MAX_PATH_LEN];     /* input file name */
 25:   PetscErrorCode        ierr;
 26:   PetscTruth            flg;
 27:   Vec                   x,y,z,work;
 28:   PetscReal             rnorm;

 30:   PetscInitialize(&argc,&args,(char *)0,help);


 33:   /* 
 34:      Determine files from which we read the two linear systems
 35:      (matrix and right-hand-side vector).
 36:   */
 37:   PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);
 38:   if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");

 40:     /* 
 41:        Open binary file.  Note that we use FILE_MODE_READ to indicate
 42:        reading from this file.
 43:     */
 44:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 46:     /*
 47:        Load the matrix; then destroy the viewer.
 48:     */
 49:     MatLoad(fd,MATAIJ,&A[0]);
 50:     PetscViewerDestroy(fd);

 52:     MatDuplicate(A[0],MAT_COPY_VALUES,&A[1]);
 53:     MatDuplicate(A[0],MAT_COPY_VALUES,&A[2]);
 54:     MatShift(A[1],1.0);
 55:     MatShift(A[1],2.0);

 57:     MatGetVecs(A[0],&x,&y);
 58:     VecDuplicate(y,&work);
 59:     VecDuplicate(y,&z);
 60: 
 61:     VecSet(x,1.0);
 62:     MatMult(A[0],x,z);
 63:     MatMultAdd(A[1],x,z,z);
 64:     MatMultAdd(A[2],x,z,z);

 66:     MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 67:     MatMult(B,x,y);
 68:     MatDestroy(B);
 69:     VecAXPY(y,-1.0,y);
 70:     VecNorm(y,NORM_2,&rnorm);
 71:     if (rnorm > 1.e-10) {
 72:       PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %G\n",rnorm);
 73:     }

 75:     MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 76:     MatCompositeMerge(B);
 77:     MatMult(B,x,y);
 78:     MatDestroy(B);
 79:     VecAXPY(y,-1.0,y);
 80:     VecNorm(y,NORM_2,&rnorm);
 81:     if (rnorm > 1.e-10) {
 82:       PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %G\n",rnorm);
 83:     }

 85:     VecSet(x,1.0);
 86:     MatMult(A[0],x,z);
 87:     MatMult(A[1],z,work);
 88:     MatMult(A[2],work,z);

 90:     MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 91:     MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
 92:     MatMult(B,x,y);
 93:     MatDestroy(B);
 94:     VecAXPY(y,-1.0,y);
 95:     VecNorm(y,NORM_2,&rnorm);
 96:     if (rnorm > 1.e-10) {
 97:       PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %G\n",rnorm);
 98:     }

100:     MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
101:     MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
102:     MatCompositeMerge(B);
103:     MatMult(B,x,y);
104:     MatDestroy(B);
105:     VecAXPY(y,-1.0,y);
106:     VecNorm(y,NORM_2,&rnorm);
107:     if (rnorm > 1.e-10) {
108:       PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative after merge %G\n",rnorm);
109:     }

111:     /* 
112:        Free work space.  All PETSc objects should be destroyed when they
113:        are no longer needed.
114:     */
115:     VecDestroy(x);
116:     VecDestroy(y);
117:     VecDestroy(work);
118:     VecDestroy(z);
119:     MatDestroy(A[0]);
120:     MatDestroy(A[1]);
121:     MatDestroy(A[2]);

123:   PetscFinalize();
124:   return 0;
125: }