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