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