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