Actual source code: ex9.c

petsc-master 2016-05-30
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>

 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 */
 26:   PetscBool      flg;
 27:   Vec            x,y,z,work;
 28:   PetscReal      rnorm;

 30:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 31:   /*
 32:      Determine files from which we read the two linear systems
 33:      (matrix and right-hand-side vector).
 34:   */
 35:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 36:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

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

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

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

 56:   MatCreateVecs(A[0],&x,&y);
 57:   VecDuplicate(y,&work);
 58:   VecDuplicate(y,&z);

 60:   VecSet(x,1.0);
 61:   MatMult(A[0],x,z);
 62:   MatMultAdd(A[1],x,z,z);
 63:   MatMultAdd(A[2],x,z,z);

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

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

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

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

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

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

122:   PetscFinalize();
123:   return ierr;
124: }