Actual source code: ex9.c
2: static char help[] = "Tests repeated setups and solves of PCFIELDSPLIT.\n\n";
3: #include <petscksp.h>
5: static PetscErrorCode replace_submats(Mat A)
6: {
7: IS *r,*c;
8: PetscInt i,j,nr,nc;
11: MatNestGetSubMats(A,&nr,&nc,NULL);
12: PetscMalloc1(nr,&r);
13: PetscMalloc1(nc,&c);
14: MatNestGetISs(A,r,c);
15: for (i=0;i<nr;i++) {
16: for (j=0;j<nc;j++) {
17: Mat sA,nA;
18: const char *prefix;
20: MatCreateSubMatrix(A,r[i],c[j],MAT_INITIAL_MATRIX,&sA);
21: MatDuplicate(sA,MAT_COPY_VALUES,&nA);
22: MatGetOptionsPrefix(sA,&prefix);
23: MatSetOptionsPrefix(nA,prefix);
24: MatNestSetSubMat(A,i,j,nA);
25: MatDestroy(&nA);
26: MatDestroy(&sA);
27: }
28: }
29: PetscFree(r);
30: PetscFree(c);
31: return 0;
32: }
34: int main(int argc, char *argv[])
35: {
36: KSP ksp;
37: PC pc;
38: Mat M,A,P,sA[2][2],sP[2][2];
39: Vec x,b;
40: IS f[2];
41: PetscInt i,j,rstart,rend;
42: PetscBool missA,missM;
44: PetscInitialize(&argc,&argv,(char*)0,help);
45: MatCreateAIJ(PETSC_COMM_WORLD,10,10,PETSC_DECIDE,PETSC_DECIDE,1,NULL,0,NULL,&M);
46: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
47: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
48: MatShift(M,1.);
49: MatGetOwnershipRange(M,&rstart,&rend);
50: ISCreateStride(PetscObjectComm((PetscObject)M),7,rstart,1,&f[0]);
51: ISComplement(f[0],rstart,rend,&f[1]);
52: for (i=0;i<2;i++) {
53: for (j=0;j<2;j++) {
54: MatCreateSubMatrix(M,f[i],f[j],MAT_INITIAL_MATRIX,&sA[i][j]);
55: MatCreateSubMatrix(M,f[i],f[j],MAT_INITIAL_MATRIX,&sP[i][j]);
56: }
57: }
58: MatCreateNest(PetscObjectComm((PetscObject)M),2,f,2,f,&sA[0][0],&A);
59: MatCreateNest(PetscObjectComm((PetscObject)M),2,f,2,f,&sP[0][0],&P);
61: /* Tests MatMissingDiagonal_Nest */
62: MatMissingDiagonal(M,&missM,NULL);
63: MatMissingDiagonal(A,&missA,NULL);
64: if (missM != missA) {
65: PetscPrintf(PETSC_COMM_WORLD,"Unexpected %s != %s\n",missM ? "true": "false",missA ? "true" : "false");
66: }
68: MatDestroy(&M);
70: KSPCreate(PetscObjectComm((PetscObject)A),&ksp);
71: KSPSetOperators(ksp,A,P);
72: KSPGetPC(ksp,&pc);
73: PCSetType(pc,PCFIELDSPLIT);
74: KSPSetFromOptions(ksp);
75: MatCreateVecs(A,&x,&b);
76: VecSetRandom(b,NULL);
77: KSPSolve(ksp,b,x);
78: replace_submats(A);
79: replace_submats(P);
80: KSPSolve(ksp,b,x);
82: KSPDestroy(&ksp);
83: VecDestroy(&x);
84: VecDestroy(&b);
85: MatDestroy(&A);
86: MatDestroy(&P);
87: for (i=0;i<2;i++) {
88: ISDestroy(&f[i]);
89: for (j=0;j<2;j++) {
90: MatDestroy(&sA[i][j]);
91: MatDestroy(&sP[i][j]);
92: }
93: }
94: PetscFinalize();
95: return 0;
96: }
98: /*TEST
100: test:
101: nsize: 1
102: filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
103: args: -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_type {{additive multiplicative}} -ksp_converged_reason -ksp_error_if_not_converged
105: test:
106: suffix: schur
107: nsize: 1
108: filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
109: args: -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_type schur -pc_fieldsplit_schur_scale 1.0 -pc_fieldsplit_schur_fact_type {{diag lower upper full}} -ksp_converged_reason -ksp_error_if_not_converged
111: TEST*/