1 2 3 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 4 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 5 #include <petscdm.h> 6 7 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 8 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 9 10 PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4; 11 12 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 13 struct _PC_FieldSplitLink { 14 KSP ksp; 15 Vec x,y,z; 16 char *splitname; 17 PetscInt nfields; 18 PetscInt *fields,*fields_col; 19 VecScatter sctx; 20 IS is,is_col; 21 PC_FieldSplitLink next,previous; 22 PetscLogEvent event; 23 }; 24 25 typedef struct { 26 PCCompositeType type; 27 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 28 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 29 PetscInt bs; /* Block size for IS and Mat structures */ 30 PetscInt nsplits; /* Number of field divisions defined */ 31 Vec *x,*y,w1,w2; 32 Mat *mat; /* The diagonal block for each split */ 33 Mat *pmat; /* The preconditioning diagonal block for each split */ 34 Mat *Afield; /* The rows of the matrix associated with each split */ 35 PetscBool issetup; 36 37 /* Only used when Schur complement preconditioning is used */ 38 Mat B; /* The (0,1) block */ 39 Mat C; /* The (1,0) block */ 40 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 41 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 42 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 43 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 44 PCFieldSplitSchurFactType schurfactorization; 45 KSP kspschur; /* The solver for S */ 46 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 47 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 48 49 PC_FieldSplitLink head; 50 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 51 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 52 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 53 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 54 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 55 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 56 } PC_FieldSplit; 57 58 /* 59 Notes: 60 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 61 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 62 PC you could change this. 63 */ 64 65 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 66 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 67 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 68 { 69 switch (jac->schurpre) { 70 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 71 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 72 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 73 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 74 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 75 default: 76 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 77 } 78 } 79 80 81 #include <petscdraw.h> 82 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 83 { 84 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 85 PetscErrorCode ierr; 86 PetscBool iascii,isdraw; 87 PetscInt i,j; 88 PC_FieldSplitLink ilink = jac->head; 89 90 PetscFunctionBegin; 91 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 92 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 93 if (iascii) { 94 if (jac->bs > 0) { 95 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 96 } else { 97 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 98 } 99 if (pc->useAmat) { 100 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 101 } 102 if (jac->diag_use_amat) { 103 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr); 104 } 105 if (jac->offdiag_use_amat) { 106 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr); 107 } 108 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 109 for (i=0; i<jac->nsplits; i++) { 110 if (ilink->fields) { 111 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 112 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 113 for (j=0; j<ilink->nfields; j++) { 114 if (j > 0) { 115 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 116 } 117 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 118 } 119 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 120 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 121 } else { 122 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 123 } 124 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 125 ilink = ilink->next; 126 } 127 } 128 129 if (isdraw) { 130 PetscDraw draw; 131 PetscReal x,y,w,wd; 132 133 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 134 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 135 w = 2*PetscMin(1.0 - x,x); 136 wd = w/(jac->nsplits + 1); 137 x = x - wd*(jac->nsplits-1)/2.0; 138 for (i=0; i<jac->nsplits; i++) { 139 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 140 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 141 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 142 x += wd; 143 ilink = ilink->next; 144 } 145 } 146 PetscFunctionReturn(0); 147 } 148 149 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 150 { 151 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 152 PetscErrorCode ierr; 153 PetscBool iascii,isdraw; 154 PetscInt i,j; 155 PC_FieldSplitLink ilink = jac->head; 156 MatSchurComplementAinvType atype; 157 158 PetscFunctionBegin; 159 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 160 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 161 if (iascii) { 162 if (jac->bs > 0) { 163 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 164 } else { 165 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 166 } 167 if (pc->useAmat) { 168 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 169 } 170 switch (jac->schurpre) { 171 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 172 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr); 173 break; 174 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 175 ierr = MatSchurComplementGetAinvType(jac->schur,&atype);CHKERRQ(ierr); 176 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped "));CHKERRQ(ierr);break; 177 case PC_FIELDSPLIT_SCHUR_PRE_A11: 178 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 179 break; 180 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 181 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr); 182 break; 183 case PC_FIELDSPLIT_SCHUR_PRE_USER: 184 if (jac->schur_user) { 185 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 186 } else { 187 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 188 } 189 break; 190 default: 191 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 192 } 193 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 194 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 195 for (i=0; i<jac->nsplits; i++) { 196 if (ilink->fields) { 197 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 198 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 199 for (j=0; j<ilink->nfields; j++) { 200 if (j > 0) { 201 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 202 } 203 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 204 } 205 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 206 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 207 } else { 208 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 209 } 210 ilink = ilink->next; 211 } 212 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 214 if (jac->head) { 215 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 216 } else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 217 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 218 if (jac->head && jac->kspupper != jac->head->ksp) { 219 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 220 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 221 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);} 222 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 223 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 224 } 225 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 226 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 227 if (jac->kspschur) { 228 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 229 } else { 230 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 231 } 232 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 233 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 234 } else if (isdraw && jac->head) { 235 PetscDraw draw; 236 PetscReal x,y,w,wd,h; 237 PetscInt cnt = 2; 238 char str[32]; 239 240 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 241 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 242 if (jac->kspupper != jac->head->ksp) cnt++; 243 w = 2*PetscMin(1.0 - x,x); 244 wd = w/(cnt + 1); 245 246 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 247 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 248 y -= h; 249 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 250 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 251 } else { 252 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 253 } 254 ierr = PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 255 y -= h; 256 x = x - wd*(cnt-1)/2.0; 257 258 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 259 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 260 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 261 if (jac->kspupper != jac->head->ksp) { 262 x += wd; 263 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 264 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 265 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 266 } 267 x += wd; 268 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 269 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 270 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 271 } 272 PetscFunctionReturn(0); 273 } 274 275 /* Precondition: jac->bs is set to a meaningful value */ 276 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 277 { 278 PetscErrorCode ierr; 279 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 280 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 281 PetscBool flg,flg_col; 282 char optionname[128],splitname[8],optionname_col[128]; 283 284 PetscFunctionBegin; 285 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr); 286 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr); 287 for (i=0,flg=PETSC_TRUE;; i++) { 288 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 289 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 290 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 291 nfields = jac->bs; 292 nfields_col = jac->bs; 293 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 294 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 295 if (!flg) break; 296 else if (flg && !flg_col) { 297 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 298 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 299 } else { 300 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 301 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 302 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 303 } 304 } 305 if (i > 0) { 306 /* Makes command-line setting of splits take precedence over setting them in code. 307 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 308 create new splits, which would probably not be what the user wanted. */ 309 jac->splitdefined = PETSC_TRUE; 310 } 311 ierr = PetscFree(ifields);CHKERRQ(ierr); 312 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 317 { 318 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 319 PetscErrorCode ierr; 320 PC_FieldSplitLink ilink = jac->head; 321 PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE; 322 PetscInt i; 323 324 PetscFunctionBegin; 325 /* 326 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 327 Should probably be rewritten. 328 */ 329 if (!ilink) { 330 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr); 331 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 332 PetscInt numFields, f, i, j; 333 char **fieldNames; 334 IS *fields; 335 DM *dms; 336 DM subdm[128]; 337 PetscBool flg; 338 339 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 340 /* Allow the user to prescribe the splits */ 341 for (i = 0, flg = PETSC_TRUE;; i++) { 342 PetscInt ifields[128]; 343 IS compField; 344 char optionname[128], splitname[8]; 345 PetscInt nfields = numFields; 346 347 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 348 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 349 if (!flg) break; 350 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 351 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 352 if (nfields == 1) { 353 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 354 } else { 355 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 356 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 357 } 358 ierr = ISDestroy(&compField);CHKERRQ(ierr); 359 for (j = 0; j < nfields; ++j) { 360 f = ifields[j]; 361 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 362 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 363 } 364 } 365 if (i == 0) { 366 for (f = 0; f < numFields; ++f) { 367 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 368 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 369 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 370 } 371 } else { 372 for (j=0; j<numFields; j++) { 373 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 374 } 375 ierr = PetscFree(dms);CHKERRQ(ierr); 376 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr); 377 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 378 } 379 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 380 ierr = PetscFree(fields);CHKERRQ(ierr); 381 if (dms) { 382 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 383 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 384 const char *prefix; 385 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 386 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 387 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 388 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 389 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 390 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 391 } 392 ierr = PetscFree(dms);CHKERRQ(ierr); 393 } 394 } else { 395 if (jac->bs <= 0) { 396 if (pc->pmat) { 397 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 398 } else jac->bs = 1; 399 } 400 401 if (jac->detect) { 402 IS zerodiags,rest; 403 PetscInt nmin,nmax; 404 405 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 406 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 407 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 408 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 409 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 410 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 411 ierr = ISDestroy(&rest);CHKERRQ(ierr); 412 } else if (coupling) { 413 IS coupling,rest; 414 PetscInt nmin,nmax; 415 416 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 417 ierr = MatFindOffBlockDiagonalEntries(pc->mat,&coupling);CHKERRQ(ierr); 418 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);CHKERRQ(ierr); 419 ierr = ISSetIdentity(rest);CHKERRQ(ierr); 420 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 421 ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr); 422 ierr = ISDestroy(&coupling);CHKERRQ(ierr); 423 ierr = ISDestroy(&rest);CHKERRQ(ierr); 424 } else { 425 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 426 if (!fieldsplit_default) { 427 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 428 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 429 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 430 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 431 } 432 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 433 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 434 for (i=0; i<jac->bs; i++) { 435 char splitname[8]; 436 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 437 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 438 } 439 jac->defaultsplit = PETSC_TRUE; 440 } 441 } 442 } 443 } else if (jac->nsplits == 1) { 444 if (ilink->is) { 445 IS is2; 446 PetscInt nmin,nmax; 447 448 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 449 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 450 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 451 ierr = ISDestroy(&is2);CHKERRQ(ierr); 452 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 453 } 454 455 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 456 PetscFunctionReturn(0); 457 } 458 459 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg); 460 461 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 462 { 463 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 464 PetscErrorCode ierr; 465 PC_FieldSplitLink ilink; 466 PetscInt i,nsplit; 467 PetscBool sorted, sorted_col; 468 469 PetscFunctionBegin; 470 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 471 nsplit = jac->nsplits; 472 ilink = jac->head; 473 474 /* get the matrices for each split */ 475 if (!jac->issetup) { 476 PetscInt rstart,rend,nslots,bs; 477 478 jac->issetup = PETSC_TRUE; 479 480 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 481 if (jac->defaultsplit || !ilink->is) { 482 if (jac->bs <= 0) jac->bs = nsplit; 483 } 484 bs = jac->bs; 485 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 486 nslots = (rend - rstart)/bs; 487 for (i=0; i<nsplit; i++) { 488 if (jac->defaultsplit) { 489 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 490 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 491 } else if (!ilink->is) { 492 if (ilink->nfields > 1) { 493 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 494 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 495 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 496 for (j=0; j<nslots; j++) { 497 for (k=0; k<nfields; k++) { 498 ii[nfields*j + k] = rstart + bs*j + fields[k]; 499 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 500 } 501 } 502 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 503 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 504 ierr = ISSetBlockSize(ilink->is, nfields);CHKERRQ(ierr); 505 ierr = ISSetBlockSize(ilink->is_col, nfields);CHKERRQ(ierr); 506 } else { 507 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 508 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 509 } 510 } 511 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 512 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 513 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 514 ilink = ilink->next; 515 } 516 } 517 518 ilink = jac->head; 519 if (!jac->pmat) { 520 Vec xtmp; 521 522 ierr = MatCreateVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 523 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 524 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 525 for (i=0; i<nsplit; i++) { 526 MatNullSpace sp; 527 528 /* Check for preconditioning matrix attached to IS */ 529 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 530 if (jac->pmat[i]) { 531 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 532 if (jac->type == PC_COMPOSITE_SCHUR) { 533 jac->schur_user = jac->pmat[i]; 534 535 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 536 } 537 } else { 538 const char *prefix; 539 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 540 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr); 541 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr); 542 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr); 543 } 544 /* create work vectors for each split */ 545 ierr = MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 546 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 547 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 548 ierr = VecScatterCreateWithData(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 549 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 550 if (sp) { 551 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 552 } 553 ilink = ilink->next; 554 } 555 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 556 } else { 557 MatReuse scall; 558 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 559 for (i=0; i<nsplit; i++) { 560 ierr = MatDestroy(&jac->pmat[i]);CHKERRQ(ierr); 561 } 562 scall = MAT_INITIAL_MATRIX; 563 } else scall = MAT_REUSE_MATRIX; 564 565 for (i=0; i<nsplit; i++) { 566 Mat pmat; 567 568 /* Check for preconditioning matrix attached to IS */ 569 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 570 if (!pmat) { 571 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);CHKERRQ(ierr); 572 } 573 ilink = ilink->next; 574 } 575 } 576 if (jac->diag_use_amat) { 577 ilink = jac->head; 578 if (!jac->mat) { 579 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 580 for (i=0; i<nsplit; i++) { 581 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 582 ilink = ilink->next; 583 } 584 } else { 585 MatReuse scall; 586 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 587 for (i=0; i<nsplit; i++) { 588 ierr = MatDestroy(&jac->mat[i]);CHKERRQ(ierr); 589 } 590 scall = MAT_INITIAL_MATRIX; 591 } else scall = MAT_REUSE_MATRIX; 592 593 for (i=0; i<nsplit; i++) { 594 if (jac->mat[i]) {ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);CHKERRQ(ierr);} 595 ilink = ilink->next; 596 } 597 } 598 } else { 599 jac->mat = jac->pmat; 600 } 601 602 /* Check for null space attached to IS */ 603 ilink = jac->head; 604 for (i=0; i<nsplit; i++) { 605 MatNullSpace sp; 606 607 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 608 if (sp) { 609 ierr = MatSetNullSpace(jac->mat[i], sp);CHKERRQ(ierr); 610 } 611 ilink = ilink->next; 612 } 613 614 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 615 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 616 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 617 ilink = jac->head; 618 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 619 /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */ 620 if (!jac->Afield) { 621 ierr = PetscCalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 622 if (jac->offdiag_use_amat) { 623 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 624 } else { 625 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 626 } 627 } else { 628 MatReuse scall; 629 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 630 for (i=0; i<nsplit; i++) { 631 ierr = MatDestroy(&jac->Afield[1]);CHKERRQ(ierr); 632 } 633 scall = MAT_INITIAL_MATRIX; 634 } else scall = MAT_REUSE_MATRIX; 635 636 if (jac->offdiag_use_amat) { 637 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr); 638 } else { 639 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr); 640 } 641 } 642 } else { 643 if (!jac->Afield) { 644 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 645 for (i=0; i<nsplit; i++) { 646 if (jac->offdiag_use_amat) { 647 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 648 } else { 649 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 650 } 651 ilink = ilink->next; 652 } 653 } else { 654 MatReuse scall; 655 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 656 for (i=0; i<nsplit; i++) { 657 ierr = MatDestroy(&jac->Afield[i]);CHKERRQ(ierr); 658 } 659 scall = MAT_INITIAL_MATRIX; 660 } else scall = MAT_REUSE_MATRIX; 661 662 for (i=0; i<nsplit; i++) { 663 if (jac->offdiag_use_amat) { 664 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr); 665 } else { 666 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr); 667 } 668 ilink = ilink->next; 669 } 670 } 671 } 672 } 673 674 if (jac->type == PC_COMPOSITE_SCHUR) { 675 IS ccis; 676 PetscBool isspd; 677 PetscInt rstart,rend; 678 char lscname[256]; 679 PetscObject LSC_L; 680 681 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 682 683 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 684 if (jac->schurscale == (PetscScalar)-1.0) { 685 ierr = MatGetOption(pc->pmat,MAT_SPD,&isspd);CHKERRQ(ierr); 686 jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0; 687 } 688 689 /* When extracting off-diagonal submatrices, we take complements from this range */ 690 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 691 692 /* need to handle case when one is resetting up the preconditioner */ 693 if (jac->schur) { 694 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 695 696 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 697 ilink = jac->head; 698 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 699 if (jac->offdiag_use_amat) { 700 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 701 } else { 702 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 703 } 704 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 705 ilink = ilink->next; 706 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 707 if (jac->offdiag_use_amat) { 708 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 709 } else { 710 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 711 } 712 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 713 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 714 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 715 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 716 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 717 } 718 if (kspA != kspInner) { 719 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 720 } 721 if (kspUpper != kspA) { 722 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 723 } 724 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 725 } else { 726 const char *Dprefix; 727 char schurprefix[256], schurmatprefix[256]; 728 char schurtestoption[256]; 729 MatNullSpace sp; 730 PetscBool flg; 731 KSP kspt; 732 733 /* extract the A01 and A10 matrices */ 734 ilink = jac->head; 735 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 736 if (jac->offdiag_use_amat) { 737 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 738 } else { 739 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 740 } 741 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 742 ilink = ilink->next; 743 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 744 if (jac->offdiag_use_amat) { 745 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 746 } else { 747 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 748 } 749 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 750 751 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 752 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 753 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 754 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 755 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 756 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 757 ierr = MatSchurComplementGetKSP(jac->schur,&kspt);CHKERRQ(ierr); 758 ierr = KSPSetOptionsPrefix(kspt,schurmatprefix);CHKERRQ(ierr); 759 760 /* Note: this is not true in general */ 761 ierr = MatGetNullSpace(jac->mat[1], &sp);CHKERRQ(ierr); 762 if (sp) { 763 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 764 } 765 766 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 767 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 768 if (flg) { 769 DM dmInner; 770 KSP kspInner; 771 PC pcInner; 772 773 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 774 ierr = KSPReset(kspInner);CHKERRQ(ierr); 775 ierr = KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 776 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 777 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 778 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 779 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);CHKERRQ(ierr); 780 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 781 782 /* Set DM for new solver */ 783 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 784 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 785 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 786 787 /* Defaults to PCKSP as preconditioner */ 788 ierr = KSPGetPC(kspInner, &pcInner);CHKERRQ(ierr); 789 ierr = PCSetType(pcInner, PCKSP);CHKERRQ(ierr); 790 ierr = PCKSPSetKSP(pcInner, jac->head->ksp);CHKERRQ(ierr); 791 } else { 792 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 793 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 794 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 795 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 796 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 797 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 798 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 799 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 800 } 801 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 802 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 803 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 804 805 ierr = PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);CHKERRQ(ierr); 806 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */ 807 KSP kspInner; 808 PC pcInner; 809 810 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 811 ierr = KSPGetPC(kspInner, &pcInner);CHKERRQ(ierr); 812 ierr = PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);CHKERRQ(ierr); 813 if (flg) { 814 KSP ksp; 815 816 ierr = PCKSPGetKSP(pcInner, &ksp);CHKERRQ(ierr); 817 if (ksp == jac->head->ksp) { 818 ierr = PCSetUseAmat(pcInner, PETSC_TRUE);CHKERRQ(ierr); 819 } 820 } 821 } 822 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 823 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 824 if (flg) { 825 DM dmInner; 826 827 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 828 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 829 ierr = KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);CHKERRQ(ierr); 830 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 831 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);CHKERRQ(ierr); 832 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);CHKERRQ(ierr); 833 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 834 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 835 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 836 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 837 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 838 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 839 } else { 840 jac->kspupper = jac->head->ksp; 841 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 842 } 843 844 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 845 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 846 } 847 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 848 ierr = KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);CHKERRQ(ierr); 849 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 850 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 851 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 852 PC pcschur; 853 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 854 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 855 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 856 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 857 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 858 } 859 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 860 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 861 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 862 /* propagate DM */ 863 { 864 DM sdm; 865 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 866 if (sdm) { 867 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 868 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 869 } 870 } 871 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 872 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 873 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 874 } 875 ierr = MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 876 ierr = MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 877 878 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 879 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 880 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 881 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 882 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 883 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 884 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 885 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 886 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 887 } else { 888 /* set up the individual splits' PCs */ 889 i = 0; 890 ilink = jac->head; 891 while (ilink) { 892 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 893 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 894 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 895 i++; 896 ilink = ilink->next; 897 } 898 } 899 900 jac->suboptionsset = PETSC_TRUE; 901 PetscFunctionReturn(0); 902 } 903 904 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 905 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 906 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 907 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 908 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 909 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 910 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 911 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 912 913 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 914 { 915 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 916 PetscErrorCode ierr; 917 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 918 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 919 920 PetscFunctionBegin; 921 switch (jac->schurfactorization) { 922 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 923 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 924 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 925 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 926 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 927 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 928 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 929 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 930 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 931 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 932 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 933 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 934 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 935 ierr = VecScale(ilinkD->y,jac->schurscale);CHKERRQ(ierr); 936 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 937 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 938 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 939 break; 940 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 941 /* [A00 0; A10 S], suitable for left preconditioning */ 942 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 943 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 944 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 945 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 946 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 947 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 948 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 949 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 950 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 951 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 952 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 953 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 954 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 955 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 956 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 957 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 958 break; 959 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 960 /* [A00 A01; 0 S], suitable for right preconditioning */ 961 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 962 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 963 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 964 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 965 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 966 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 967 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 968 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 969 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 970 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 971 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 972 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 973 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 974 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 975 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 976 break; 977 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 978 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */ 979 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 980 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 981 ierr = PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 982 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 983 ierr = PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 984 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 985 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 986 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 987 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 988 989 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 990 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 991 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 992 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 993 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 994 995 if (kspUpper == kspA) { 996 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 997 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 998 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 999 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1000 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1001 } else { 1002 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1003 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1004 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1005 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 1006 ierr = PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 1007 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 1008 ierr = PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 1009 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 1010 } 1011 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1012 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1013 } 1014 PetscFunctionReturn(0); 1015 } 1016 1017 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 1018 { 1019 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1020 PetscErrorCode ierr; 1021 PC_FieldSplitLink ilink = jac->head; 1022 PetscInt cnt,bs; 1023 KSPConvergedReason reason; 1024 1025 PetscFunctionBegin; 1026 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1027 if (jac->defaultsplit) { 1028 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1029 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1030 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1031 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1032 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1033 while (ilink) { 1034 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1035 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1036 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1037 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1038 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1039 pc->failedreason = PC_SUBPC_ERROR; 1040 } 1041 ilink = ilink->next; 1042 } 1043 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1044 } else { 1045 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1046 while (ilink) { 1047 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1048 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1049 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1050 pc->failedreason = PC_SUBPC_ERROR; 1051 } 1052 ilink = ilink->next; 1053 } 1054 } 1055 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1056 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1057 /* solve on first block for first block variables */ 1058 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1059 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1060 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1061 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1062 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1063 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1064 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1065 pc->failedreason = PC_SUBPC_ERROR; 1066 } 1067 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1068 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1069 1070 /* compute the residual only onto second block variables using first block variables */ 1071 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 1072 ilink = ilink->next; 1073 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1074 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1075 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1076 1077 /* solve on second block variables */ 1078 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1079 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1080 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1081 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1082 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1083 pc->failedreason = PC_SUBPC_ERROR; 1084 } 1085 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1086 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1087 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1088 if (!jac->w1) { 1089 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1090 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1091 } 1092 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1093 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1094 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1095 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1096 pc->failedreason = PC_SUBPC_ERROR; 1097 } 1098 cnt = 1; 1099 while (ilink->next) { 1100 ilink = ilink->next; 1101 /* compute the residual only over the part of the vector needed */ 1102 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 1103 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1104 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1105 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1106 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1107 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1108 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1109 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1110 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1111 pc->failedreason = PC_SUBPC_ERROR; 1112 } 1113 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1114 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1115 } 1116 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1117 cnt -= 2; 1118 while (ilink->previous) { 1119 ilink = ilink->previous; 1120 /* compute the residual only over the part of the vector needed */ 1121 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 1122 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1123 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1124 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1125 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1126 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1127 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1128 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1129 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1130 pc->failedreason = PC_SUBPC_ERROR; 1131 } 1132 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1133 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1134 } 1135 } 1136 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1141 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1142 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1143 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1144 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1145 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1146 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1147 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1148 1149 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1150 { 1151 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1152 PetscErrorCode ierr; 1153 PC_FieldSplitLink ilink = jac->head; 1154 PetscInt bs; 1155 KSPConvergedReason reason; 1156 1157 PetscFunctionBegin; 1158 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1159 if (jac->defaultsplit) { 1160 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1161 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1162 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1163 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1164 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1165 while (ilink) { 1166 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1167 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1168 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1169 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1170 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1171 pc->failedreason = PC_SUBPC_ERROR; 1172 } 1173 ilink = ilink->next; 1174 } 1175 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1176 } else { 1177 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1178 while (ilink) { 1179 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1180 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1181 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1182 pc->failedreason = PC_SUBPC_ERROR; 1183 } 1184 ilink = ilink->next; 1185 } 1186 } 1187 } else { 1188 if (!jac->w1) { 1189 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1190 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1191 } 1192 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1193 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1194 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1195 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1196 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1197 pc->failedreason = PC_SUBPC_ERROR; 1198 } 1199 while (ilink->next) { 1200 ilink = ilink->next; 1201 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1202 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1203 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1204 } 1205 while (ilink->previous) { 1206 ilink = ilink->previous; 1207 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1208 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1209 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1210 } 1211 } else { 1212 while (ilink->next) { /* get to last entry in linked list */ 1213 ilink = ilink->next; 1214 } 1215 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1216 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1217 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1218 pc->failedreason = PC_SUBPC_ERROR; 1219 } 1220 while (ilink->previous) { 1221 ilink = ilink->previous; 1222 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1223 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1224 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1225 } 1226 } 1227 } 1228 PetscFunctionReturn(0); 1229 } 1230 1231 static PetscErrorCode PCReset_FieldSplit(PC pc) 1232 { 1233 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1234 PetscErrorCode ierr; 1235 PC_FieldSplitLink ilink = jac->head,next; 1236 1237 PetscFunctionBegin; 1238 while (ilink) { 1239 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1240 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1241 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1242 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1243 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1244 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1245 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1246 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1247 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1248 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1249 next = ilink->next; 1250 ierr = PetscFree(ilink);CHKERRQ(ierr); 1251 ilink = next; 1252 } 1253 jac->head = NULL; 1254 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1255 if (jac->mat && jac->mat != jac->pmat) { 1256 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1257 } else if (jac->mat) { 1258 jac->mat = NULL; 1259 } 1260 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1261 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1262 jac->nsplits = 0; 1263 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1264 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1265 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1266 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1267 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1268 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1269 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1270 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1271 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1272 jac->isrestrict = PETSC_FALSE; 1273 PetscFunctionReturn(0); 1274 } 1275 1276 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1277 { 1278 PetscErrorCode ierr; 1279 1280 PetscFunctionBegin; 1281 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1282 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1283 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);CHKERRQ(ierr); 1284 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1285 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1286 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1287 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1288 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1289 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1290 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1291 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1292 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);CHKERRQ(ierr); 1293 PetscFunctionReturn(0); 1294 } 1295 1296 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1297 { 1298 PetscErrorCode ierr; 1299 PetscInt bs; 1300 PetscBool flg; 1301 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1302 PCCompositeType ctype; 1303 1304 PetscFunctionBegin; 1305 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1306 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1307 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1308 if (flg) { 1309 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1310 } 1311 jac->diag_use_amat = pc->useAmat; 1312 ierr = PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",jac->diag_use_amat,&jac->diag_use_amat,NULL);CHKERRQ(ierr); 1313 jac->offdiag_use_amat = pc->useAmat; 1314 ierr = PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat","Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat",jac->offdiag_use_amat,&jac->offdiag_use_amat,NULL);CHKERRQ(ierr); 1315 ierr = PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);CHKERRQ(ierr); 1316 ierr = PCFieldSplitSetDetectSaddlePoint(pc,jac->detect);CHKERRQ(ierr); /* Sets split type and Schur PC type */ 1317 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1318 if (flg) { 1319 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1320 } 1321 /* Only setup fields once */ 1322 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1323 /* only allow user to set fields from command line if bs is already known. 1324 otherwise user can set them in PCFieldSplitSetDefaults() */ 1325 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1326 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1327 } 1328 if (jac->type == PC_COMPOSITE_SCHUR) { 1329 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1330 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1331 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_fact_type","Which off-diagonal parts of the block factorization to use","PCFieldSplitSetSchurFactType",PCFieldSplitSchurFactTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,NULL);CHKERRQ(ierr); 1332 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1333 ierr = PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);CHKERRQ(ierr); 1334 } 1335 ierr = PetscOptionsTail();CHKERRQ(ierr); 1336 PetscFunctionReturn(0); 1337 } 1338 1339 /*------------------------------------------------------------------------------------*/ 1340 1341 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1342 { 1343 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1344 PetscErrorCode ierr; 1345 PC_FieldSplitLink ilink,next = jac->head; 1346 char prefix[128]; 1347 PetscInt i; 1348 1349 PetscFunctionBegin; 1350 if (jac->splitdefined) { 1351 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1352 PetscFunctionReturn(0); 1353 } 1354 for (i=0; i<n; i++) { 1355 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1356 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1357 } 1358 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1359 if (splitname) { 1360 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1361 } else { 1362 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1363 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1364 } 1365 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1366 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1367 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1368 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1369 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1370 1371 ilink->nfields = n; 1372 ilink->next = NULL; 1373 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1374 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1375 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1376 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1377 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1378 1379 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1380 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1381 1382 if (!next) { 1383 jac->head = ilink; 1384 ilink->previous = NULL; 1385 } else { 1386 while (next->next) { 1387 next = next->next; 1388 } 1389 next->next = ilink; 1390 ilink->previous = next; 1391 } 1392 jac->nsplits++; 1393 PetscFunctionReturn(0); 1394 } 1395 1396 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1397 { 1398 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBegin; 1402 *subksp = NULL; 1403 if (n) *n = 0; 1404 if (jac->type == PC_COMPOSITE_SCHUR) { 1405 PetscInt nn; 1406 1407 if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1408 if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits); 1409 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1410 ierr = PetscMalloc1(nn,subksp);CHKERRQ(ierr); 1411 (*subksp)[0] = jac->head->ksp; 1412 (*subksp)[1] = jac->kspschur; 1413 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1414 if (n) *n = nn; 1415 } 1416 PetscFunctionReturn(0); 1417 } 1418 1419 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1420 { 1421 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1422 PetscErrorCode ierr; 1423 1424 PetscFunctionBegin; 1425 if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1426 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1427 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1428 1429 (*subksp)[1] = jac->kspschur; 1430 if (n) *n = jac->nsplits; 1431 PetscFunctionReturn(0); 1432 } 1433 1434 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1435 { 1436 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1437 PetscErrorCode ierr; 1438 PetscInt cnt = 0; 1439 PC_FieldSplitLink ilink = jac->head; 1440 1441 PetscFunctionBegin; 1442 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1443 while (ilink) { 1444 (*subksp)[cnt++] = ilink->ksp; 1445 ilink = ilink->next; 1446 } 1447 if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number of splits in linked list %D does not match number in object %D",cnt,jac->nsplits); 1448 if (n) *n = jac->nsplits; 1449 PetscFunctionReturn(0); 1450 } 1451 1452 /*@C 1453 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1454 1455 Input Parameters: 1456 + pc - the preconditioner context 1457 + is - the index set that defines the indices to which the fieldsplit is to be restricted 1458 1459 Level: advanced 1460 1461 @*/ 1462 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1463 { 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1468 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1469 ierr = PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));CHKERRQ(ierr); 1470 PetscFunctionReturn(0); 1471 } 1472 1473 1474 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1475 { 1476 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1477 PetscErrorCode ierr; 1478 PC_FieldSplitLink ilink = jac->head, next; 1479 PetscInt localsize,size,sizez,i; 1480 const PetscInt *ind, *indz; 1481 PetscInt *indc, *indcz; 1482 PetscBool flg; 1483 1484 PetscFunctionBegin; 1485 ierr = ISGetLocalSize(isy,&localsize);CHKERRQ(ierr); 1486 ierr = MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));CHKERRQ(ierr); 1487 size -= localsize; 1488 while(ilink) { 1489 IS isrl,isr; 1490 PC subpc; 1491 ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1492 ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr); 1493 ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr); 1494 ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr); 1495 ierr = PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1496 ierr = ISRestoreIndices(isrl,&ind);CHKERRQ(ierr); 1497 ierr = ISDestroy(&isrl);CHKERRQ(ierr); 1498 for (i=0; i<localsize; i++) *(indc+i) += size; 1499 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);CHKERRQ(ierr); 1500 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1501 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1502 ilink->is = isr; 1503 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1504 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1505 ilink->is_col = isr; 1506 ierr = ISDestroy(&isr);CHKERRQ(ierr); 1507 ierr = KSPGetPC(ilink->ksp, &subpc);CHKERRQ(ierr); 1508 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1509 if(flg) { 1510 IS iszl,isz; 1511 MPI_Comm comm; 1512 ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr); 1513 comm = PetscObjectComm((PetscObject)ilink->is); 1514 ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1515 ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 1516 sizez -= localsize; 1517 ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr); 1518 ierr = PetscMalloc1(localsize,&indcz);CHKERRQ(ierr); 1519 ierr = ISGetIndices(iszl,&indz);CHKERRQ(ierr); 1520 ierr = PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1521 ierr = ISRestoreIndices(iszl,&indz);CHKERRQ(ierr); 1522 ierr = ISDestroy(&iszl);CHKERRQ(ierr); 1523 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1524 ierr = ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);CHKERRQ(ierr); 1525 ierr = PCFieldSplitRestrictIS(subpc,isz);CHKERRQ(ierr); 1526 ierr = ISDestroy(&isz);CHKERRQ(ierr); 1527 } 1528 next = ilink->next; 1529 ilink = next; 1530 } 1531 jac->isrestrict = PETSC_TRUE; 1532 PetscFunctionReturn(0); 1533 } 1534 1535 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1536 { 1537 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1538 PetscErrorCode ierr; 1539 PC_FieldSplitLink ilink, next = jac->head; 1540 char prefix[128]; 1541 1542 PetscFunctionBegin; 1543 if (jac->splitdefined) { 1544 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1545 PetscFunctionReturn(0); 1546 } 1547 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1548 if (splitname) { 1549 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1550 } else { 1551 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1552 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1553 } 1554 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1555 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1556 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1557 ilink->is = is; 1558 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1559 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1560 ilink->is_col = is; 1561 ilink->next = NULL; 1562 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1563 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1564 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1565 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1566 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1567 1568 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1569 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1570 1571 if (!next) { 1572 jac->head = ilink; 1573 ilink->previous = NULL; 1574 } else { 1575 while (next->next) { 1576 next = next->next; 1577 } 1578 next->next = ilink; 1579 ilink->previous = next; 1580 } 1581 jac->nsplits++; 1582 PetscFunctionReturn(0); 1583 } 1584 1585 /*@C 1586 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1587 1588 Logically Collective on PC 1589 1590 Input Parameters: 1591 + pc - the preconditioner context 1592 . splitname - name of this split, if NULL the number of the split is used 1593 . n - the number of fields in this split 1594 - fields - the fields in this split 1595 1596 Level: intermediate 1597 1598 Notes: 1599 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1600 1601 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1602 size is three then one can define a field as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean 1603 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1604 where the numbered entries indicate what is in the field. 1605 1606 This function is called once per split (it creates a new split each time). Solve options 1607 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1608 1609 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1610 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1611 available when this routine is called. 1612 1613 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1614 1615 @*/ 1616 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1617 { 1618 PetscErrorCode ierr; 1619 1620 PetscFunctionBegin; 1621 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1622 PetscValidCharPointer(splitname,2); 1623 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1624 PetscValidIntPointer(fields,3); 1625 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1626 PetscFunctionReturn(0); 1627 } 1628 1629 /*@ 1630 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1631 1632 Logically Collective on PC 1633 1634 Input Parameters: 1635 + pc - the preconditioner object 1636 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1637 1638 Options Database: 1639 . -pc_fieldsplit_diag_use_amat 1640 1641 Level: intermediate 1642 1643 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1644 1645 @*/ 1646 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1647 { 1648 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1649 PetscBool isfs; 1650 PetscErrorCode ierr; 1651 1652 PetscFunctionBegin; 1653 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1654 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1655 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1656 jac->diag_use_amat = flg; 1657 PetscFunctionReturn(0); 1658 } 1659 1660 /*@ 1661 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1662 1663 Logically Collective on PC 1664 1665 Input Parameters: 1666 . pc - the preconditioner object 1667 1668 Output Parameters: 1669 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1670 1671 1672 Level: intermediate 1673 1674 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1675 1676 @*/ 1677 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1678 { 1679 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1680 PetscBool isfs; 1681 PetscErrorCode ierr; 1682 1683 PetscFunctionBegin; 1684 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1685 PetscValidPointer(flg,2); 1686 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1687 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1688 *flg = jac->diag_use_amat; 1689 PetscFunctionReturn(0); 1690 } 1691 1692 /*@ 1693 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1694 1695 Logically Collective on PC 1696 1697 Input Parameters: 1698 + pc - the preconditioner object 1699 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1700 1701 Options Database: 1702 . -pc_fieldsplit_off_diag_use_amat 1703 1704 Level: intermediate 1705 1706 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1707 1708 @*/ 1709 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1710 { 1711 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1712 PetscBool isfs; 1713 PetscErrorCode ierr; 1714 1715 PetscFunctionBegin; 1716 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1717 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1718 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1719 jac->offdiag_use_amat = flg; 1720 PetscFunctionReturn(0); 1721 } 1722 1723 /*@ 1724 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1725 1726 Logically Collective on PC 1727 1728 Input Parameters: 1729 . pc - the preconditioner object 1730 1731 Output Parameters: 1732 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1733 1734 1735 Level: intermediate 1736 1737 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1738 1739 @*/ 1740 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1741 { 1742 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1743 PetscBool isfs; 1744 PetscErrorCode ierr; 1745 1746 PetscFunctionBegin; 1747 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1748 PetscValidPointer(flg,2); 1749 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1750 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1751 *flg = jac->offdiag_use_amat; 1752 PetscFunctionReturn(0); 1753 } 1754 1755 1756 1757 /*@C 1758 PCFieldSplitSetIS - Sets the exact elements for field 1759 1760 Logically Collective on PC 1761 1762 Input Parameters: 1763 + pc - the preconditioner context 1764 . splitname - name of this split, if NULL the number of the split is used 1765 - is - the index set that defines the vector elements in this field 1766 1767 1768 Notes: 1769 Use PCFieldSplitSetFields(), for fields defined by strided types. 1770 1771 This function is called once per split (it creates a new split each time). Solve options 1772 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1773 1774 Level: intermediate 1775 1776 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1777 1778 @*/ 1779 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1780 { 1781 PetscErrorCode ierr; 1782 1783 PetscFunctionBegin; 1784 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1785 if (splitname) PetscValidCharPointer(splitname,2); 1786 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1787 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1788 PetscFunctionReturn(0); 1789 } 1790 1791 /*@C 1792 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1793 1794 Logically Collective on PC 1795 1796 Input Parameters: 1797 + pc - the preconditioner context 1798 - splitname - name of this split 1799 1800 Output Parameter: 1801 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1802 1803 Level: intermediate 1804 1805 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1806 1807 @*/ 1808 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1809 { 1810 PetscErrorCode ierr; 1811 1812 PetscFunctionBegin; 1813 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1814 PetscValidCharPointer(splitname,2); 1815 PetscValidPointer(is,3); 1816 { 1817 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1818 PC_FieldSplitLink ilink = jac->head; 1819 PetscBool found; 1820 1821 *is = NULL; 1822 while (ilink) { 1823 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1824 if (found) { 1825 *is = ilink->is; 1826 break; 1827 } 1828 ilink = ilink->next; 1829 } 1830 } 1831 PetscFunctionReturn(0); 1832 } 1833 1834 /*@ 1835 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1836 fieldsplit preconditioner. If not set the matrix block size is used. 1837 1838 Logically Collective on PC 1839 1840 Input Parameters: 1841 + pc - the preconditioner context 1842 - bs - the block size 1843 1844 Level: intermediate 1845 1846 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1847 1848 @*/ 1849 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1850 { 1851 PetscErrorCode ierr; 1852 1853 PetscFunctionBegin; 1854 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1855 PetscValidLogicalCollectiveInt(pc,bs,2); 1856 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1857 PetscFunctionReturn(0); 1858 } 1859 1860 /*@C 1861 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1862 1863 Collective on KSP 1864 1865 Input Parameter: 1866 . pc - the preconditioner context 1867 1868 Output Parameters: 1869 + n - the number of splits 1870 - subksp - the array of KSP contexts 1871 1872 Note: 1873 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 1874 (not the KSP just the array that contains them). 1875 1876 You must call PCSetUp() before calling PCFieldSplitGetSubKSP(). 1877 1878 If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the 1879 Schur complement and the KSP object used to iterate over the Schur complement. 1880 To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP() 1881 1882 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1883 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 1884 KSP array must be. 1885 1886 1887 Level: advanced 1888 1889 .seealso: PCFIELDSPLIT 1890 @*/ 1891 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1892 { 1893 PetscErrorCode ierr; 1894 1895 PetscFunctionBegin; 1896 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1897 if (n) PetscValidIntPointer(n,2); 1898 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 /*@C 1903 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT 1904 1905 Collective on KSP 1906 1907 Input Parameter: 1908 . pc - the preconditioner context 1909 1910 Output Parameters: 1911 + n - the number of splits 1912 - subksp - the array of KSP contexts 1913 1914 Note: 1915 After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 1916 (not the KSP just the array that contains them). 1917 1918 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP(). 1919 1920 If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order) 1921 - the KSP used for the (1,1) block 1922 - the KSP used for the Schur complement (not the one used for the interior Schur solver) 1923 - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 1924 1925 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP(). 1926 1927 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1928 You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 1929 KSP array must be. 1930 1931 Level: advanced 1932 1933 .seealso: PCFIELDSPLIT 1934 @*/ 1935 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1936 { 1937 PetscErrorCode ierr; 1938 1939 PetscFunctionBegin; 1940 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1941 if (n) PetscValidIntPointer(n,2); 1942 ierr = PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1943 PetscFunctionReturn(0); 1944 } 1945 1946 /*@ 1947 PCFieldSplitSetSchurPre - Indicates what operator is used to construct the preconditioner for the Schur complement. 1948 A11 matrix. Otherwise no preconditioner is used. 1949 1950 Collective on PC 1951 1952 Input Parameters: 1953 + pc - the preconditioner context 1954 . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER 1955 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 1956 - userpre - matrix to use for preconditioning, or NULL 1957 1958 Options Database: 1959 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 1960 1961 Notes: 1962 $ If ptype is 1963 $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 1964 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 1965 $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 1966 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 1967 $ preconditioner 1968 $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 1969 $ to this function). 1970 $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1971 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 1972 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 1973 $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive) 1974 $ useful mostly as a test that the Schur complement approach can work for your problem 1975 1976 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1977 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1978 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1979 1980 Level: intermediate 1981 1982 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1983 MatSchurComplementSetAinvType(), PCLSC 1984 1985 @*/ 1986 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1987 { 1988 PetscErrorCode ierr; 1989 1990 PetscFunctionBegin; 1991 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1992 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1993 PetscFunctionReturn(0); 1994 } 1995 1996 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1997 1998 /*@ 1999 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2000 preconditioned. See PCFieldSplitSetSchurPre() for details. 2001 2002 Logically Collective on PC 2003 2004 Input Parameters: 2005 . pc - the preconditioner context 2006 2007 Output Parameters: 2008 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 2009 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 2010 2011 Level: intermediate 2012 2013 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 2014 2015 @*/ 2016 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2017 { 2018 PetscErrorCode ierr; 2019 2020 PetscFunctionBegin; 2021 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2022 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 2023 PetscFunctionReturn(0); 2024 } 2025 2026 /*@ 2027 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 2028 2029 Not collective 2030 2031 Input Parameter: 2032 . pc - the preconditioner context 2033 2034 Output Parameter: 2035 . S - the Schur complement matrix 2036 2037 Notes: 2038 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 2039 2040 Level: advanced 2041 2042 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 2043 2044 @*/ 2045 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 2046 { 2047 PetscErrorCode ierr; 2048 const char* t; 2049 PetscBool isfs; 2050 PC_FieldSplit *jac; 2051 2052 PetscFunctionBegin; 2053 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2054 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 2055 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2056 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2057 jac = (PC_FieldSplit*)pc->data; 2058 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2059 if (S) *S = jac->schur; 2060 PetscFunctionReturn(0); 2061 } 2062 2063 /*@ 2064 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 2065 2066 Not collective 2067 2068 Input Parameters: 2069 + pc - the preconditioner context 2070 . S - the Schur complement matrix 2071 2072 Level: advanced 2073 2074 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 2075 2076 @*/ 2077 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2078 { 2079 PetscErrorCode ierr; 2080 const char* t; 2081 PetscBool isfs; 2082 PC_FieldSplit *jac; 2083 2084 PetscFunctionBegin; 2085 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2086 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 2087 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2088 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2089 jac = (PC_FieldSplit*)pc->data; 2090 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2091 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2092 PetscFunctionReturn(0); 2093 } 2094 2095 2096 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2097 { 2098 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2099 PetscErrorCode ierr; 2100 2101 PetscFunctionBegin; 2102 jac->schurpre = ptype; 2103 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2104 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 2105 jac->schur_user = pre; 2106 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 2107 } 2108 PetscFunctionReturn(0); 2109 } 2110 2111 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2112 { 2113 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2114 2115 PetscFunctionBegin; 2116 *ptype = jac->schurpre; 2117 *pre = jac->schur_user; 2118 PetscFunctionReturn(0); 2119 } 2120 2121 /*@ 2122 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner 2123 2124 Collective on PC 2125 2126 Input Parameters: 2127 + pc - the preconditioner context 2128 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2129 2130 Options Database: 2131 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 2132 2133 2134 Level: intermediate 2135 2136 Notes: 2137 The FULL factorization is 2138 2139 $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U 2140 $ (C E) (C*Ainv 1) (0 S) (0 1 ) 2141 2142 where S = E - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D, 2143 and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale(). 2144 2145 $ If A and S are solved exactly 2146 $ *) FULL factorization is a direct solver. 2147 $ *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations. 2148 $ *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 2149 2150 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2151 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2152 2153 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. 2154 2155 Note that a flexible method like KSPFGMRES or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S). 2156 2157 References: 2158 + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2159 - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2160 2161 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale() 2162 @*/ 2163 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2164 { 2165 PetscErrorCode ierr; 2166 2167 PetscFunctionBegin; 2168 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2169 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 2170 PetscFunctionReturn(0); 2171 } 2172 2173 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2174 { 2175 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2176 2177 PetscFunctionBegin; 2178 jac->schurfactorization = ftype; 2179 PetscFunctionReturn(0); 2180 } 2181 2182 /*@ 2183 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG. 2184 2185 Collective on PC 2186 2187 Input Parameters: 2188 + pc - the preconditioner context 2189 - scale - scaling factor for the Schur complement 2190 2191 Options Database: 2192 . -pc_fieldsplit_schur_scale - default is -1.0 2193 2194 Level: intermediate 2195 2196 .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale() 2197 @*/ 2198 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale) 2199 { 2200 PetscErrorCode ierr; 2201 2202 PetscFunctionBegin; 2203 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2204 PetscValidLogicalCollectiveScalar(pc,scale,2); 2205 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));CHKERRQ(ierr); 2206 PetscFunctionReturn(0); 2207 } 2208 2209 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale) 2210 { 2211 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2212 2213 PetscFunctionBegin; 2214 jac->schurscale = scale; 2215 PetscFunctionReturn(0); 2216 } 2217 2218 /*@C 2219 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2220 2221 Collective on KSP 2222 2223 Input Parameter: 2224 . pc - the preconditioner context 2225 2226 Output Parameters: 2227 + A00 - the (0,0) block 2228 . A01 - the (0,1) block 2229 . A10 - the (1,0) block 2230 - A11 - the (1,1) block 2231 2232 Level: advanced 2233 2234 .seealso: PCFIELDSPLIT 2235 @*/ 2236 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2237 { 2238 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2239 2240 PetscFunctionBegin; 2241 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2242 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2243 if (A00) *A00 = jac->pmat[0]; 2244 if (A01) *A01 = jac->B; 2245 if (A10) *A10 = jac->C; 2246 if (A11) *A11 = jac->pmat[1]; 2247 PetscFunctionReturn(0); 2248 } 2249 2250 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2251 { 2252 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2253 PetscErrorCode ierr; 2254 2255 PetscFunctionBegin; 2256 jac->type = type; 2257 if (type == PC_COMPOSITE_SCHUR) { 2258 pc->ops->apply = PCApply_FieldSplit_Schur; 2259 pc->ops->view = PCView_FieldSplit_Schur; 2260 2261 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 2262 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 2263 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 2264 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 2265 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);CHKERRQ(ierr); 2266 2267 } else { 2268 pc->ops->apply = PCApply_FieldSplit; 2269 pc->ops->view = PCView_FieldSplit; 2270 2271 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2272 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 2273 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 2274 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 2275 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);CHKERRQ(ierr); 2276 } 2277 PetscFunctionReturn(0); 2278 } 2279 2280 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2281 { 2282 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2283 2284 PetscFunctionBegin; 2285 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2286 if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 2287 jac->bs = bs; 2288 PetscFunctionReturn(0); 2289 } 2290 2291 /*@ 2292 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2293 2294 Collective on PC 2295 2296 Input Parameter: 2297 . pc - the preconditioner context 2298 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2299 2300 Options Database Key: 2301 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2302 2303 Level: Intermediate 2304 2305 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2306 2307 .seealso: PCCompositeSetType() 2308 2309 @*/ 2310 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2311 { 2312 PetscErrorCode ierr; 2313 2314 PetscFunctionBegin; 2315 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2316 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2317 PetscFunctionReturn(0); 2318 } 2319 2320 /*@ 2321 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2322 2323 Not collective 2324 2325 Input Parameter: 2326 . pc - the preconditioner context 2327 2328 Output Parameter: 2329 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2330 2331 Level: Intermediate 2332 2333 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2334 .seealso: PCCompositeSetType() 2335 @*/ 2336 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2337 { 2338 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2339 2340 PetscFunctionBegin; 2341 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2342 PetscValidIntPointer(type,2); 2343 *type = jac->type; 2344 PetscFunctionReturn(0); 2345 } 2346 2347 /*@ 2348 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2349 2350 Logically Collective 2351 2352 Input Parameters: 2353 + pc - the preconditioner context 2354 - flg - boolean indicating whether to use field splits defined by the DM 2355 2356 Options Database Key: 2357 . -pc_fieldsplit_dm_splits 2358 2359 Level: Intermediate 2360 2361 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2362 2363 .seealso: PCFieldSplitGetDMSplits() 2364 2365 @*/ 2366 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2367 { 2368 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2369 PetscBool isfs; 2370 PetscErrorCode ierr; 2371 2372 PetscFunctionBegin; 2373 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2374 PetscValidLogicalCollectiveBool(pc,flg,2); 2375 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2376 if (isfs) { 2377 jac->dm_splits = flg; 2378 } 2379 PetscFunctionReturn(0); 2380 } 2381 2382 2383 /*@ 2384 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2385 2386 Logically Collective 2387 2388 Input Parameter: 2389 . pc - the preconditioner context 2390 2391 Output Parameter: 2392 . flg - boolean indicating whether to use field splits defined by the DM 2393 2394 Level: Intermediate 2395 2396 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2397 2398 .seealso: PCFieldSplitSetDMSplits() 2399 2400 @*/ 2401 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2402 { 2403 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2404 PetscBool isfs; 2405 PetscErrorCode ierr; 2406 2407 PetscFunctionBegin; 2408 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2409 PetscValidPointer(flg,2); 2410 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2411 if (isfs) { 2412 if(flg) *flg = jac->dm_splits; 2413 } 2414 PetscFunctionReturn(0); 2415 } 2416 2417 /*@ 2418 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2419 2420 Logically Collective 2421 2422 Input Parameter: 2423 . pc - the preconditioner context 2424 2425 Output Parameter: 2426 . flg - boolean indicating whether to detect fields or not 2427 2428 Level: Intermediate 2429 2430 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint() 2431 2432 @*/ 2433 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg) 2434 { 2435 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2436 2437 PetscFunctionBegin; 2438 *flg = jac->detect; 2439 PetscFunctionReturn(0); 2440 } 2441 2442 /*@ 2443 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2444 2445 Logically Collective 2446 2447 Notes: 2448 Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()). 2449 2450 Input Parameter: 2451 . pc - the preconditioner context 2452 2453 Output Parameter: 2454 . flg - boolean indicating whether to detect fields or not 2455 2456 Options Database Key: 2457 . -pc_fieldsplit_detect_saddle_point 2458 2459 Level: Intermediate 2460 2461 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre() 2462 2463 @*/ 2464 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg) 2465 { 2466 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2467 PetscErrorCode ierr; 2468 2469 PetscFunctionBegin; 2470 jac->detect = flg; 2471 if (jac->detect) { 2472 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 2473 ierr = PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);CHKERRQ(ierr); 2474 } 2475 PetscFunctionReturn(0); 2476 } 2477 2478 /* -------------------------------------------------------------------------------------*/ 2479 /*MC 2480 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2481 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2482 2483 To set options on the solvers for each block append -fieldsplit_ to all the PC 2484 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2485 2486 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2487 and set the options directly on the resulting KSP object 2488 2489 Level: intermediate 2490 2491 Options Database Keys: 2492 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2493 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2494 been supplied explicitly by -pc_fieldsplit_%d_fields 2495 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2496 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2497 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 2498 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 2499 2500 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2501 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2502 2503 Notes: 2504 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2505 to define a field by an arbitrary collection of entries. 2506 2507 If no fields are set the default is used. The fields are defined by entries strided by bs, 2508 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2509 if this is not called the block size defaults to the blocksize of the second matrix passed 2510 to KSPSetOperators()/PCSetOperators(). 2511 2512 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2513 $ ( A10 A11 ) 2514 $ the preconditioner using full factorization is 2515 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2516 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2517 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2518 $ S = A11 - A10 ksp(A00) A01 2519 which is usually dense and not stored explicitly. The action of ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 2520 in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0, 2521 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2522 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 2523 2524 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2525 diag gives 2526 $ ( inv(A00) 0 ) 2527 $ ( 0 -ksp(S) ) 2528 note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip 2529 can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of 2530 $ ( A00 0 ) 2531 $ ( A10 S ) 2532 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2533 $ ( A00 A01 ) 2534 $ ( 0 S ) 2535 where again the inverses of A00 and S are applied using KSPs. 2536 2537 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2538 is used automatically for a second block. 2539 2540 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2541 Generally it should be used with the AIJ format. 2542 2543 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2544 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2545 inside a smoother resulting in "Distributive Smoothers". 2546 2547 Concepts: physics based preconditioners, block preconditioners 2548 2549 There is a nice discussion of block preconditioners in 2550 2551 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2552 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2553 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2554 2555 The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the 2556 residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables. 2557 2558 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2559 PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2560 MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), 2561 PCFieldSplitSetDetectSaddlePoint() 2562 M*/ 2563 2564 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2565 { 2566 PetscErrorCode ierr; 2567 PC_FieldSplit *jac; 2568 2569 PetscFunctionBegin; 2570 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2571 2572 jac->bs = -1; 2573 jac->nsplits = 0; 2574 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2575 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2576 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2577 jac->schurscale = -1.0; 2578 jac->dm_splits = PETSC_TRUE; 2579 jac->detect = PETSC_FALSE; 2580 2581 pc->data = (void*)jac; 2582 2583 pc->ops->apply = PCApply_FieldSplit; 2584 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2585 pc->ops->setup = PCSetUp_FieldSplit; 2586 pc->ops->reset = PCReset_FieldSplit; 2587 pc->ops->destroy = PCDestroy_FieldSplit; 2588 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2589 pc->ops->view = PCView_FieldSplit; 2590 pc->ops->applyrichardson = 0; 2591 2592 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);CHKERRQ(ierr); 2593 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2594 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2595 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2596 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2597 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2598 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);CHKERRQ(ierr); 2599 PetscFunctionReturn(0); 2600 } 2601