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*/