Actual source code: ex9.c
petsc-3.6.1 2015-08-06
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);
33: /*
34: Determine files from which we read the two linear systems
35: (matrix and right-hand-side vector).
36: */
37: PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
38: if (!flg) SETERRQ(PETSC_COMM_WORLD,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: MatCreate(PETSC_COMM_WORLD,&A[0]);
50: MatLoad(A[0],fd);
51: PetscViewerDestroy(&fd);
53: MatDuplicate(A[0],MAT_COPY_VALUES,&A[1]);
54: MatDuplicate(A[0],MAT_COPY_VALUES,&A[2]);
55: MatShift(A[1],1.0);
56: MatShift(A[1],2.0);
58: MatCreateVecs(A[0],&x,&y);
59: VecDuplicate(y,&work);
60: VecDuplicate(y,&z);
62: VecSet(x,1.0);
63: MatMult(A[0],x,z);
64: MatMultAdd(A[1],x,z,z);
65: MatMultAdd(A[2],x,z,z);
67: MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
68: MatMult(B,x,y);
69: MatDestroy(&B);
70: VecAXPY(y,-1.0,z);
71: VecNorm(y,NORM_2,&rnorm);
72: if (rnorm > 1.e-10) {
73: PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);
74: }
76: MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
77: MatCompositeMerge(B);
78: MatMult(B,x,y);
79: MatDestroy(&B);
80: VecAXPY(y,-1.0,z);
81: VecNorm(y,NORM_2,&rnorm);
82: if (rnorm > 1.e-10) {
83: PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);
84: }
86: VecSet(x,1.0);
87: MatMult(A[0],x,z);
88: MatMult(A[1],z,work);
89: MatMult(A[2],work,z);
91: MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
92: MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
93: MatMult(B,x,y);
94: MatDestroy(&B);
95: VecAXPY(y,-1.0,z);
96: VecNorm(y,NORM_2,&rnorm);
97: if (rnorm > 1.e-10) {
98: PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
99: }
101: MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);
102: MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
103: MatCompositeMerge(B);
104: MatMult(B,x,y);
105: MatDestroy(&B);
106: VecAXPY(y,-1.0,z);
107: VecNorm(y,NORM_2,&rnorm);
108: if (rnorm > 1.e-10) {
109: PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative after merge %g\n",(double)rnorm);
110: }
112: /*
113: Free work space. All PETSc objects should be destroyed when they
114: are no longer needed.
115: */
116: VecDestroy(&x);
117: VecDestroy(&y);
118: VecDestroy(&work);
119: VecDestroy(&z);
120: MatDestroy(&A[0]);
121: MatDestroy(&A[1]);
122: MatDestroy(&A[2]);
124: PetscFinalize();
125: return 0;
126: }