Actual source code: ex9.c

petsc-3.13.4 2020-08-01
Report Typos and Errors

  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: {
  8:   IS             *r,*c;
  9:   PetscInt       i,j,nr,nc;

 12:   MatNestGetSubMats(A,&nr,&nc,NULL);
 13:   PetscMalloc1(nr,&r);
 14:   PetscMalloc1(nc,&c);
 15:   MatNestGetISs(A,r,c);
 16:   for (i=0;i<nr;i++) {
 17:     for (j=0;j<nc;j++) {
 18:       Mat        sA,nA;
 19:       const char *prefix;

 21:       MatCreateSubMatrix(A,r[i],c[j],MAT_INITIAL_MATRIX,&sA);
 22:       MatDuplicate(sA,MAT_COPY_VALUES,&nA);
 23:       MatGetOptionsPrefix(sA,&prefix);
 24:       MatSetOptionsPrefix(nA,prefix);
 25:       MatNestSetSubMat(A,i,j,nA);
 26:       MatDestroy(&nA);
 27:       MatDestroy(&sA);
 28:     }
 29:   }
 30:   PetscFree(r);
 31:   PetscFree(c);
 32:   return(0);
 33: }

 35: int main(int argc, char *argv[])
 36: {
 37:    KSP            ksp;
 38:    PC             pc;
 39:    Mat            M,A,P,sA[2][2],sP[2][2];
 40:    Vec            x,b;
 41:    IS             f[2];
 42:    PetscInt       i,j,rstart,rend;
 43:    PetscBool      missA,missM;

 46:    PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 47:    MatCreateAIJ(PETSC_COMM_WORLD,10,10,PETSC_DECIDE,PETSC_DECIDE,1,NULL,0,NULL,&M);
 48:    MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
 49:    MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
 50:    MatShift(M,1.);
 51:    MatGetOwnershipRange(M,&rstart,&rend);
 52:    ISCreateStride(PetscObjectComm((PetscObject)M),7,rstart,1,&f[0]);
 53:    ISComplement(f[0],rstart,rend,&f[1]);
 54:    for (i=0;i<2;i++) {
 55:      for (j=0;j<2;j++) {
 56:        MatCreateSubMatrix(M,f[i],f[j],MAT_INITIAL_MATRIX,&sA[i][j]);
 57:        MatCreateSubMatrix(M,f[i],f[j],MAT_INITIAL_MATRIX,&sP[i][j]);
 58:      }
 59:    }
 60:    MatCreateNest(PetscObjectComm((PetscObject)M),2,f,2,f,&sA[0][0],&A);
 61:    MatCreateNest(PetscObjectComm((PetscObject)M),2,f,2,f,&sP[0][0],&P);

 63:    /* Tests MatMissingDiagonal_Nest */
 64:    MatMissingDiagonal(M,&missM,NULL);
 65:    MatMissingDiagonal(A,&missA,NULL);
 66:    if (missM != missA) {
 67:      PetscPrintf(PETSC_COMM_WORLD,"Unexpected %s != %s\n",missM ? "true": "false",missA ? "true" : "false");
 68:    }

 70:    MatDestroy(&M);

 72:    KSPCreate(PetscObjectComm((PetscObject)A),&ksp);
 73:    KSPSetOperators(ksp,A,P);
 74:    KSPGetPC(ksp,&pc);
 75:    PCSetType(pc,PCFIELDSPLIT);
 76:    KSPSetFromOptions(ksp);
 77:    MatCreateVecs(A,&x,&b);
 78:    VecSetRandom(b,NULL);
 79:    KSPSolve(ksp,b,x);
 80:    replace_submats(A);
 81:    replace_submats(P);
 82:    KSPSolve(ksp,b,x);

 84:    KSPDestroy(&ksp);
 85:    VecDestroy(&x);
 86:    VecDestroy(&b);
 87:    MatDestroy(&A);
 88:    MatDestroy(&P);
 89:    for (i=0;i<2;i++) {
 90:      ISDestroy(&f[i]);
 91:      for (j=0;j<2;j++) {
 92:        MatDestroy(&sA[i][j]);
 93:        MatDestroy(&sP[i][j]);
 94:      }
 95:    }
 96:    PetscFinalize();
 97:    return ierr;
 98: }

100: /*TEST

102:    test:
103:      nsize: 1
104:      filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
105:      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

107:    test:
108:      suffix: schur
109:      nsize: 1
110:      filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
111:      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

113: TEST*/