Actual source code: ex22.c

petsc-3.3-p6 2013-02-11
  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: {

229:   PetscInitialize( &argc, &args,(char *)0, help);
230:   test_solve();
231:   test_solve_matgetvecs();
232:   PetscFinalize();
233:   return 0;
234: }