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