Actual source code: ex22.c
1: static const char help[] = "Test MatNest solving a linear system\n\n";
3: #include <petscksp.h>
7: PetscErrorCode test_solve( void )
8: {
9: Mat A11, A12,A21,A22, A, tmp[2][2];
10: KSP ksp;
11: PC pc;
12: Vec b,x , f,h, diag, x1,x2;
13: Vec tmp_x[2],*_tmp_x;
14: int n, np, i,j;
18: PetscPrintf( PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME );
20: n = 3;
21: np = 2;
22: /* Create matrices */
23: /* A11 */
24: VecCreate( PETSC_COMM_WORLD, &diag );
25: VecSetSizes( diag, PETSC_DECIDE, n );
26: VecSetFromOptions(diag);
28: VecSet( diag, (1.0/10.0) ); /* so inverse = diag(10) */
30: /* As a test, create a diagonal matrix for A11 */
31: MatCreate( PETSC_COMM_WORLD, &A11 );
32: MatSetSizes( A11, PETSC_DECIDE, PETSC_DECIDE, n, n );
33: MatSetType( A11, MATAIJ );
34: MatSeqAIJSetPreallocation( A11, n, PETSC_NULL );
35: MatMPIAIJSetPreallocation( A11, np, PETSC_NULL,np, PETSC_NULL );
36: MatDiagonalSet( A11, diag, INSERT_VALUES );
38: VecDestroy(& diag );
40: /* A12 */
41: MatCreate( PETSC_COMM_WORLD, &A12 );
42: MatSetSizes( A12, PETSC_DECIDE, PETSC_DECIDE, n, np );
43: MatSetType( A12, MATAIJ );
44: MatSeqAIJSetPreallocation( A12, np, PETSC_NULL );
45: MatMPIAIJSetPreallocation( A12, np, PETSC_NULL,np, PETSC_NULL );
47: for( i=0; i<n; i++ ) {
48: for( j=0; j<np; j++ ) {
49: MatSetValue( A12, i,j, (double)(i+j*n), INSERT_VALUES );
50: }
51: }
52: MatSetValue( A12, 2,1, (double)(4), INSERT_VALUES );
53: MatAssemblyBegin( A12, MAT_FINAL_ASSEMBLY);
54: MatAssemblyEnd( A12, MAT_FINAL_ASSEMBLY);
56: /* A21 */
57: MatTranspose( A12, MAT_INITIAL_MATRIX, &A21 );
59: A22 = PETSC_NULL;
61: /* Create block matrix */
62: tmp[0][0] = A11;
63: tmp[0][1] = A12;
64: tmp[1][0] = A21;
65: tmp[1][1] = A22;
66: MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,2,PETSC_NULL,&tmp[0][0],&A);
67: MatNestSetVecType(A,VECNEST);
68: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
71: /* Create vectors */
72: MatGetVecs( A12, &h, &f );
74: VecSet( f, 1.0 );
75: VecSet( h, 0.0 );
77: /* Create block vector */
78: tmp_x[0] = f;
79: tmp_x[1] = h;
80: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_x,&b);
81: VecAssemblyBegin(b);
82: VecAssemblyEnd(b);
83: VecDuplicate( b, &x );
85: KSPCreate( PETSC_COMM_WORLD, &ksp );
86: KSPSetOperators( ksp, A, A, SAME_NONZERO_PATTERN );
87: KSPSetType( ksp, "gmres" );
88: KSPGetPC( ksp, &pc );
89: PCSetType( pc, "none" );
90: KSPSetFromOptions( ksp );
92: KSPSolve( ksp, b, x );
94: VecNestGetSubVecs(x,PETSC_NULL,&_tmp_x);
95: x1 = _tmp_x[0];
96: x2 = _tmp_x[1];
98: PetscPrintf( PETSC_COMM_WORLD, "x1 \n");
99: PetscViewerSetFormat( PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL );
100: VecView( x1, PETSC_VIEWER_STDOUT_WORLD );
101: PetscPrintf( PETSC_COMM_WORLD, "x2 \n");
102: PetscViewerSetFormat( PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL );
103: VecView( x2, PETSC_VIEWER_STDOUT_WORLD );
105: KSPDestroy(& ksp );
106: VecDestroy(& x );
107: VecDestroy(& b );
108: MatDestroy(& A11 );
109: MatDestroy(& A12 );
110: MatDestroy(& A21 );
111: VecDestroy(& f );
112: VecDestroy(& h );
114: MatDestroy(& A );
116: return(0);
117: }
122: PetscErrorCode test_solve_matgetvecs( void )
123: {
124: Mat A11, A12,A21, A;
125: KSP ksp;
126: PC pc;
127: Vec b,x , f,h, diag, x1,x2;
128: int n, np, i,j;
129: Mat tmp[2][2];
130: Vec *tmp_x;
134: PetscPrintf( PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME );
136: n = 3;
137: np = 2;
138: /* Create matrices */
139: /* A11 */
140: VecCreate( PETSC_COMM_WORLD, &diag );
141: VecSetSizes( diag, PETSC_DECIDE, n );
142: VecSetFromOptions(diag);
144: VecSet( diag, (1.0/10.0) ); /* so inverse = diag(10) */
146: /* As a test, create a diagonal matrix for A11 */
147: MatCreate( PETSC_COMM_WORLD, &A11 );
148: MatSetSizes( A11, PETSC_DECIDE, PETSC_DECIDE, n, n );
149: MatSetType( A11, MATAIJ );
150: MatSeqAIJSetPreallocation( A11, n, PETSC_NULL );
151: MatMPIAIJSetPreallocation( A11, np, PETSC_NULL,np, PETSC_NULL );
152: MatDiagonalSet( A11, diag, INSERT_VALUES );
154: VecDestroy(& diag );
156: /* A12 */
157: MatCreate( PETSC_COMM_WORLD, &A12 );
158: MatSetSizes( A12, PETSC_DECIDE, PETSC_DECIDE, n, np );
159: MatSetType( A12, MATAIJ );
160: MatSeqAIJSetPreallocation( A12, np, PETSC_NULL );
161: MatMPIAIJSetPreallocation( A12, np, PETSC_NULL,np, PETSC_NULL );
163: for( i=0; i<n; i++ ) {
164: for( j=0; j<np; j++ ) {
165: MatSetValue( A12, i,j, (double)(i+j*n), INSERT_VALUES );
166: }
167: }
168: MatSetValue( A12, 2,1, (double)(4), INSERT_VALUES );
169: MatAssemblyBegin( A12, MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd( A12, MAT_FINAL_ASSEMBLY);
172: /* A21 */
173: MatTranspose( A12, MAT_INITIAL_MATRIX, &A21 );
175: /* Create block matrix */
176: tmp[0][0] = A11;
177: tmp[0][1] = A12;
178: tmp[1][0] = A21;
179: tmp[1][1] = PETSC_NULL;
180: MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,2,PETSC_NULL,&tmp[0][0],&A);
181: MatNestSetVecType(A,VECNEST);
182: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
183: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
185: /* Create vectors */
186: MatGetVecs( A, &b, &x );
187: VecNestGetSubVecs(b,PETSC_NULL,&tmp_x);
188: f = tmp_x[0];
189: h = tmp_x[1];
191: VecSet( f, 1.0 );
192: VecSet( h, 0.0 );
194: KSPCreate( PETSC_COMM_WORLD, &ksp );
195: KSPSetOperators( ksp, A, A, SAME_NONZERO_PATTERN );
196: KSPGetPC( ksp, &pc );
197: PCSetType( pc, PCNONE );
198: KSPSetFromOptions( ksp );
200: KSPSolve( ksp, b, x );
201: VecNestGetSubVecs(x,PETSC_NULL,&tmp_x);
202: x1 = tmp_x[0];
203: x2 = tmp_x[1];
205: PetscPrintf( PETSC_COMM_WORLD, "x1 \n");
206: PetscViewerSetFormat( PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL );
207: VecView( x1, PETSC_VIEWER_STDOUT_WORLD );
208: PetscPrintf( PETSC_COMM_WORLD, "x2 \n");
209: PetscViewerSetFormat( PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL );
210: VecView( x2, PETSC_VIEWER_STDOUT_WORLD );
212: KSPDestroy(& ksp );
213: VecDestroy(& x );
214: VecDestroy(& b );
215: MatDestroy(& A11 );
216: MatDestroy(& A12 );
217: MatDestroy(& A21 );
219: MatDestroy(& A );
220: return(0);
221: }
225: int main( int argc, char **args )
226: {
227: PetscInitialize( &argc, &args,(char *)0, help);
229: test_solve();
230: test_solve_matgetvecs();
232: PetscFinalize();
233: return 0;
234: }