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[3],B;                       /* matrix */
 21:   PetscViewer    fd;                           /* viewer */
 22:   char           file[PETSC_MAX_PATH_LEN];     /* input file name */
 24:   PetscBool      flg;
 25:   Vec            x,y,z,work;
 26:   PetscReal      rnorm;

 28:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 29:   /*
 30:      Determine files from which we read the matrix
 31:   */
 32:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 33:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

 35:   /*
 36:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 37:      reading from this file.
 38:   */
 39:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 41:   /*
 42:      Load the matrix; then destroy the viewer.
 43:   */
 44:   MatCreate(PETSC_COMM_WORLD,&A[0]);
 45:   MatLoad(A[0],fd);
 46:   PetscViewerDestroy(&fd);

 48:   MatDuplicate(A[0],MAT_COPY_VALUES,&A[1]);
 49:   MatDuplicate(A[0],MAT_COPY_VALUES,&A[2]);
 50:   MatShift(A[1],1.0);
 51:   MatShift(A[1],2.0);

 53:   MatCreateVecs(A[0],&x,&y);
 54:   VecDuplicate(y,&work);
 55:   VecDuplicate(y,&z);

 57:   VecSet(x,1.0);
 58:   MatMult(A[0],x,z);
 59:   MatMultAdd(A[1],x,z,z);
 60:   MatMultAdd(A[2],x,z,z);

 62:   MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 63:   MatMult(B,x,y);
 64:   MatDestroy(&B);
 65:   VecAXPY(y,-1.0,z);
 66:   VecNorm(y,NORM_2,&rnorm);
 67:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
 68:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);
 69:   }

 71:   MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 72:   MatCompositeMerge(B);
 73:   MatMult(B,x,y);
 74:   MatDestroy(&B);
 75:   VecAXPY(y,-1.0,z);
 76:   VecNorm(y,NORM_2,&rnorm);
 77:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
 78:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);
 79:   }

 81:   VecSet(x,1.0);
 82:   MatMult(A[0],x,z);
 83:   MatMult(A[1],z,work);
 84:   MatMult(A[2],work,z);

 86:   MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 87:   MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
 88:   MatMult(B,x,y);
 89:   MatDestroy(&B);
 90:   VecAXPY(y,-1.0,z);
 91:   VecNorm(y,NORM_2,&rnorm);
 92:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
 93:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
 94:   }

 96:   MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
 97:   MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
 98:   MatCompositeMerge(B);
 99:   MatMult(B,x,y);
100:   MatDestroy(&B);
101:   VecAXPY(y,-1.0,z);
102:   VecNorm(y,NORM_2,&rnorm);
103:   if (rnorm > 1000000.0*PETSC_MACHINE_EPSILON) {
104:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative after merge %g\n",(double)rnorm);
105:   }

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

119:   PetscFinalize();
120:   return ierr;
121: }

123: /*TEST

125:    test:
126:       nsize: 2
127:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
128:       args: -f ${DATAFILESPATH}/matrices/medium -viewer_binary_skip_info

130: TEST*/