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: }