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)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 289 ierr = PetscOptionsGetIntArray(((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)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 328 ierr = PetscOptionsGetBool(((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)->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 = coupling; 429 jac->head->next->is = rest; 430 } else { 431 ierr = PCFieldSplitSetIS(pc,"0",coupling);CHKERRQ(ierr); 432 ierr = PCFieldSplitSetIS(pc,"1",rest);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)->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(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 /* Check for null space attached to IS */ 566 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 567 if (sp) { 568 ierr = MatSetNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 569 } 570 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 571 if (sp) { 572 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 573 } 574 ilink = ilink->next; 575 } 576 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 577 } else { 578 for (i=0; i<nsplit; i++) { 579 Mat pmat; 580 581 /* Check for preconditioning matrix attached to IS */ 582 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 583 if (!pmat) { 584 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 585 } 586 ilink = ilink->next; 587 } 588 } 589 if (jac->diag_use_amat) { 590 ilink = jac->head; 591 if (!jac->mat) { 592 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 593 for (i=0; i<nsplit; i++) { 594 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 595 ilink = ilink->next; 596 } 597 } else { 598 for (i=0; i<nsplit; i++) { 599 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 600 ilink = ilink->next; 601 } 602 } 603 } else { 604 jac->mat = jac->pmat; 605 } 606 607 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 608 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 609 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 610 ilink = jac->head; 611 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 612 /* 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 */ 613 if (!jac->Afield) { 614 ierr = PetscCalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 615 ierr = MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 616 } else { 617 ierr = MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 618 } 619 } else { 620 if (!jac->Afield) { 621 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 622 for (i=0; i<nsplit; i++) { 623 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 624 ilink = ilink->next; 625 } 626 } else { 627 for (i=0; i<nsplit; i++) { 628 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 629 ilink = ilink->next; 630 } 631 } 632 } 633 } 634 635 if (jac->type == PC_COMPOSITE_SCHUR) { 636 IS ccis; 637 PetscInt rstart,rend; 638 char lscname[256]; 639 PetscObject LSC_L; 640 641 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 642 643 /* When extracting off-diagonal submatrices, we take complements from this range */ 644 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 645 646 /* need to handle case when one is resetting up the preconditioner */ 647 if (jac->schur) { 648 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 649 650 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 651 ilink = jac->head; 652 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 653 if (jac->offdiag_use_amat) { 654 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 655 } else { 656 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 657 } 658 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 659 ilink = ilink->next; 660 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 661 if (jac->offdiag_use_amat) { 662 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 663 } else { 664 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 665 } 666 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 667 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 668 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 669 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 670 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 671 } 672 if (kspA != kspInner) { 673 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 674 } 675 if (kspUpper != kspA) { 676 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 677 } 678 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 679 } else { 680 const char *Dprefix; 681 char schurprefix[256], schurmatprefix[256]; 682 char schurtestoption[256]; 683 MatNullSpace sp; 684 PetscBool flg; 685 686 /* extract the A01 and A10 matrices */ 687 ilink = jac->head; 688 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 689 if (jac->offdiag_use_amat) { 690 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 691 } else { 692 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 693 } 694 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 695 ilink = ilink->next; 696 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 697 if (jac->offdiag_use_amat) { 698 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 699 } else { 700 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 701 } 702 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 703 704 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 705 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 706 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 707 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 708 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 709 /* 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? */ 710 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 711 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 712 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 713 if (sp) { 714 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 715 } 716 717 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 718 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 719 if (flg) { 720 DM dmInner; 721 KSP kspInner; 722 723 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 724 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 725 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 726 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 727 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 728 729 /* Set DM for new solver */ 730 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 731 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 732 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 733 } else { 734 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 735 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 736 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 737 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 738 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 739 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 740 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 741 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 742 } 743 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 744 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 745 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 746 747 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 748 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 749 if (flg) { 750 DM dmInner; 751 752 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 753 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 754 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 755 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 756 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 757 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 758 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 759 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 760 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 761 } else { 762 jac->kspupper = jac->head->ksp; 763 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 764 } 765 766 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 767 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 768 } 769 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 770 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 771 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 772 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 773 PC pcschur; 774 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 775 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 776 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 777 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 778 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 779 } 780 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 781 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 782 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 783 /* propogate DM */ 784 { 785 DM sdm; 786 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 787 if (sdm) { 788 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 789 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 790 } 791 } 792 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 793 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 794 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 795 } 796 797 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 798 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 799 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 800 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 801 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 802 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 803 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 804 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 805 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 806 } else { 807 /* set up the individual splits' PCs */ 808 i = 0; 809 ilink = jac->head; 810 while (ilink) { 811 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 812 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 813 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 814 i++; 815 ilink = ilink->next; 816 } 817 } 818 819 jac->suboptionsset = PETSC_TRUE; 820 PetscFunctionReturn(0); 821 } 822 823 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 824 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 825 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 826 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 827 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 828 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 829 830 #undef __FUNCT__ 831 #define __FUNCT__ "PCApply_FieldSplit_Schur" 832 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 833 { 834 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 835 PetscErrorCode ierr; 836 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 837 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 838 839 PetscFunctionBegin; 840 switch (jac->schurfactorization) { 841 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 842 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 843 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 844 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 845 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 846 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 847 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 848 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 849 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 850 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 851 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 852 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 853 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 854 break; 855 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 856 /* [A00 0; A10 S], suitable for left preconditioning */ 857 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 858 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 859 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 860 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 861 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 862 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 863 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 864 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 865 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 866 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 867 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 868 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 869 break; 870 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 871 /* [A00 A01; 0 S], suitable for right preconditioning */ 872 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 873 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 874 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 875 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 876 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 877 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 878 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 879 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 880 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 881 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 882 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 883 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 884 break; 885 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 886 /* [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 */ 887 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 888 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 889 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 890 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 891 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 892 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 893 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 894 895 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 896 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 897 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 898 899 if (kspUpper == kspA) { 900 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 901 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 902 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 903 } else { 904 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 905 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 906 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 907 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 908 } 909 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 910 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 911 } 912 PetscFunctionReturn(0); 913 } 914 915 #undef __FUNCT__ 916 #define __FUNCT__ "PCApply_FieldSplit" 917 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 918 { 919 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 920 PetscErrorCode ierr; 921 PC_FieldSplitLink ilink = jac->head; 922 PetscInt cnt,bs; 923 924 PetscFunctionBegin; 925 if (jac->type == PC_COMPOSITE_ADDITIVE) { 926 if (jac->defaultsplit) { 927 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 928 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); 929 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 930 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); 931 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 932 while (ilink) { 933 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 934 ilink = ilink->next; 935 } 936 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 937 } else { 938 ierr = VecSet(y,0.0);CHKERRQ(ierr); 939 while (ilink) { 940 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 941 ilink = ilink->next; 942 } 943 } 944 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 945 ierr = VecSet(y,0.0);CHKERRQ(ierr); 946 /* solve on first block for first block variables */ 947 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 948 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 949 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 950 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 951 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 952 953 /* compute the residual only onto second block variables using first block variables */ 954 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 955 ilink = ilink->next; 956 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 957 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 958 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 959 960 /* solve on second block variables */ 961 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 962 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 963 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 964 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 965 if (!jac->w1) { 966 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 967 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 968 } 969 ierr = VecSet(y,0.0);CHKERRQ(ierr); 970 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 971 cnt = 1; 972 while (ilink->next) { 973 ilink = ilink->next; 974 /* compute the residual only over the part of the vector needed */ 975 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 976 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 977 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 978 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 979 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 980 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 981 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 982 } 983 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 984 cnt -= 2; 985 while (ilink->previous) { 986 ilink = ilink->previous; 987 /* compute the residual only over the part of the vector needed */ 988 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 989 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 990 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 991 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 992 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 993 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 994 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 995 } 996 } 997 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 998 PetscFunctionReturn(0); 999 } 1000 1001 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1002 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1003 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1004 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1005 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1006 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1007 1008 #undef __FUNCT__ 1009 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 1010 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1011 { 1012 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1013 PetscErrorCode ierr; 1014 PC_FieldSplitLink ilink = jac->head; 1015 PetscInt bs; 1016 1017 PetscFunctionBegin; 1018 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1019 if (jac->defaultsplit) { 1020 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1021 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); 1022 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1023 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); 1024 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1025 while (ilink) { 1026 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1027 ilink = ilink->next; 1028 } 1029 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1030 } else { 1031 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1032 while (ilink) { 1033 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1034 ilink = ilink->next; 1035 } 1036 } 1037 } else { 1038 if (!jac->w1) { 1039 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1040 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1041 } 1042 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1043 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1044 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1045 while (ilink->next) { 1046 ilink = ilink->next; 1047 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1048 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1049 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1050 } 1051 while (ilink->previous) { 1052 ilink = ilink->previous; 1053 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1054 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1055 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1056 } 1057 } else { 1058 while (ilink->next) { /* get to last entry in linked list */ 1059 ilink = ilink->next; 1060 } 1061 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1062 while (ilink->previous) { 1063 ilink = ilink->previous; 1064 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1065 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1066 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1067 } 1068 } 1069 } 1070 PetscFunctionReturn(0); 1071 } 1072 1073 #undef __FUNCT__ 1074 #define __FUNCT__ "PCReset_FieldSplit" 1075 static PetscErrorCode PCReset_FieldSplit(PC pc) 1076 { 1077 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1078 PetscErrorCode ierr; 1079 PC_FieldSplitLink ilink = jac->head,next; 1080 1081 PetscFunctionBegin; 1082 while (ilink) { 1083 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1084 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1085 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1086 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1087 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1088 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1089 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1090 next = ilink->next; 1091 ilink = next; 1092 } 1093 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1094 if (jac->mat && jac->mat != jac->pmat) { 1095 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1096 } else if (jac->mat) { 1097 jac->mat = NULL; 1098 } 1099 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1100 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1101 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1102 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1103 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1104 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1105 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1106 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1107 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1108 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1109 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1110 jac->reset = PETSC_TRUE; 1111 PetscFunctionReturn(0); 1112 } 1113 1114 #undef __FUNCT__ 1115 #define __FUNCT__ "PCDestroy_FieldSplit" 1116 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1117 { 1118 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1119 PetscErrorCode ierr; 1120 PC_FieldSplitLink ilink = jac->head,next; 1121 1122 PetscFunctionBegin; 1123 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1124 while (ilink) { 1125 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1126 next = ilink->next; 1127 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1128 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1129 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1130 ierr = PetscFree(ilink);CHKERRQ(ierr); 1131 ilink = next; 1132 } 1133 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1134 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1135 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1136 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1137 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1138 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1139 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1140 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1141 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1142 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1148 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptions *PetscOptionsObject,PC pc) 1149 { 1150 PetscErrorCode ierr; 1151 PetscInt bs; 1152 PetscBool flg,stokes = PETSC_FALSE; 1153 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1154 PCCompositeType ctype; 1155 1156 PetscFunctionBegin; 1157 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1158 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1159 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1160 if (flg) { 1161 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1162 } 1163 jac->diag_use_amat = pc->useAmat; 1164 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); 1165 jac->offdiag_use_amat = pc->useAmat; 1166 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); 1167 /* FIXME: No programmatic equivalent to the following. */ 1168 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1169 if (stokes) { 1170 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1171 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1172 } 1173 1174 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1175 if (flg) { 1176 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1177 } 1178 /* Only setup fields once */ 1179 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1180 /* only allow user to set fields from command line if bs is already known. 1181 otherwise user can set them in PCFieldSplitSetDefaults() */ 1182 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1183 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1184 } 1185 if (jac->type == PC_COMPOSITE_SCHUR) { 1186 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1187 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1188 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); 1189 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1190 } 1191 ierr = PetscOptionsTail();CHKERRQ(ierr); 1192 PetscFunctionReturn(0); 1193 } 1194 1195 /*------------------------------------------------------------------------------------*/ 1196 1197 #undef __FUNCT__ 1198 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1199 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1200 { 1201 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1202 PetscErrorCode ierr; 1203 PC_FieldSplitLink ilink,next = jac->head; 1204 char prefix[128]; 1205 PetscInt i; 1206 1207 PetscFunctionBegin; 1208 if (jac->splitdefined) { 1209 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1210 PetscFunctionReturn(0); 1211 } 1212 for (i=0; i<n; i++) { 1213 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1214 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1215 } 1216 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1217 if (splitname) { 1218 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1219 } else { 1220 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1221 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1222 } 1223 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1224 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1225 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1226 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1227 1228 ilink->nfields = n; 1229 ilink->next = NULL; 1230 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1231 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1232 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1233 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1234 1235 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1236 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1237 1238 if (!next) { 1239 jac->head = ilink; 1240 ilink->previous = NULL; 1241 } else { 1242 while (next->next) { 1243 next = next->next; 1244 } 1245 next->next = ilink; 1246 ilink->previous = next; 1247 } 1248 jac->nsplits++; 1249 PetscFunctionReturn(0); 1250 } 1251 1252 #undef __FUNCT__ 1253 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1254 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1255 { 1256 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1257 PetscErrorCode ierr; 1258 1259 PetscFunctionBegin; 1260 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1261 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1262 1263 (*subksp)[1] = jac->kspschur; 1264 if (n) *n = jac->nsplits; 1265 PetscFunctionReturn(0); 1266 } 1267 1268 #undef __FUNCT__ 1269 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1270 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1271 { 1272 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1273 PetscErrorCode ierr; 1274 PetscInt cnt = 0; 1275 PC_FieldSplitLink ilink = jac->head; 1276 1277 PetscFunctionBegin; 1278 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1279 while (ilink) { 1280 (*subksp)[cnt++] = ilink->ksp; 1281 ilink = ilink->next; 1282 } 1283 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); 1284 if (n) *n = jac->nsplits; 1285 PetscFunctionReturn(0); 1286 } 1287 1288 #undef __FUNCT__ 1289 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1290 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1291 { 1292 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1293 PetscErrorCode ierr; 1294 PC_FieldSplitLink ilink, next = jac->head; 1295 char prefix[128]; 1296 1297 PetscFunctionBegin; 1298 if (jac->splitdefined) { 1299 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1300 PetscFunctionReturn(0); 1301 } 1302 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1303 if (splitname) { 1304 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1305 } else { 1306 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1307 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1308 } 1309 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1310 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1311 ilink->is = is; 1312 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1313 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1314 ilink->is_col = is; 1315 ilink->next = NULL; 1316 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1317 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1318 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1319 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1320 1321 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1322 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1323 1324 if (!next) { 1325 jac->head = ilink; 1326 ilink->previous = NULL; 1327 } else { 1328 while (next->next) { 1329 next = next->next; 1330 } 1331 next->next = ilink; 1332 ilink->previous = next; 1333 } 1334 jac->nsplits++; 1335 PetscFunctionReturn(0); 1336 } 1337 1338 #undef __FUNCT__ 1339 #define __FUNCT__ "PCFieldSplitSetFields" 1340 /*@ 1341 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1342 1343 Logically Collective on PC 1344 1345 Input Parameters: 1346 + pc - the preconditioner context 1347 . splitname - name of this split, if NULL the number of the split is used 1348 . n - the number of fields in this split 1349 - fields - the fields in this split 1350 1351 Level: intermediate 1352 1353 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1354 1355 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1356 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 1357 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1358 where the numbered entries indicate what is in the field. 1359 1360 This function is called once per split (it creates a new split each time). Solve options 1361 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1362 1363 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1364 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1365 available when this routine is called. 1366 1367 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1368 1369 @*/ 1370 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1371 { 1372 PetscErrorCode ierr; 1373 1374 PetscFunctionBegin; 1375 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1376 PetscValidCharPointer(splitname,2); 1377 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1378 PetscValidIntPointer(fields,3); 1379 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1380 PetscFunctionReturn(0); 1381 } 1382 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "PCFieldSplitSetDiagUseAmat" 1385 /*@ 1386 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1387 1388 Logically Collective on PC 1389 1390 Input Parameters: 1391 + pc - the preconditioner object 1392 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1393 1394 Options Database: 1395 . -pc_fieldsplit_diag_use_amat 1396 1397 Level: intermediate 1398 1399 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1400 1401 @*/ 1402 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1403 { 1404 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1405 PetscBool isfs; 1406 PetscErrorCode ierr; 1407 1408 PetscFunctionBegin; 1409 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1410 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1411 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1412 jac->diag_use_amat = flg; 1413 PetscFunctionReturn(0); 1414 } 1415 1416 #undef __FUNCT__ 1417 #define __FUNCT__ "PCFieldSplitGetDiagUseAmat" 1418 /*@ 1419 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1420 1421 Logically Collective on PC 1422 1423 Input Parameters: 1424 . pc - the preconditioner object 1425 1426 Output Parameters: 1427 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1428 1429 1430 Level: intermediate 1431 1432 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1433 1434 @*/ 1435 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1436 { 1437 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1438 PetscBool isfs; 1439 PetscErrorCode ierr; 1440 1441 PetscFunctionBegin; 1442 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1443 PetscValidPointer(flg,2); 1444 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1445 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1446 *flg = jac->diag_use_amat; 1447 PetscFunctionReturn(0); 1448 } 1449 1450 #undef __FUNCT__ 1451 #define __FUNCT__ "PCFieldSplitSetOffDiagUseAmat" 1452 /*@ 1453 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1454 1455 Logically Collective on PC 1456 1457 Input Parameters: 1458 + pc - the preconditioner object 1459 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1460 1461 Options Database: 1462 . -pc_fieldsplit_off_diag_use_amat 1463 1464 Level: intermediate 1465 1466 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1467 1468 @*/ 1469 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1470 { 1471 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1472 PetscBool isfs; 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1477 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1478 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1479 jac->offdiag_use_amat = flg; 1480 PetscFunctionReturn(0); 1481 } 1482 1483 #undef __FUNCT__ 1484 #define __FUNCT__ "PCFieldSplitGetOffDiagUseAmat" 1485 /*@ 1486 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1487 1488 Logically Collective on PC 1489 1490 Input Parameters: 1491 . pc - the preconditioner object 1492 1493 Output Parameters: 1494 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1495 1496 1497 Level: intermediate 1498 1499 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1500 1501 @*/ 1502 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1503 { 1504 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1505 PetscBool isfs; 1506 PetscErrorCode ierr; 1507 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1510 PetscValidPointer(flg,2); 1511 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1512 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1513 *flg = jac->offdiag_use_amat; 1514 PetscFunctionReturn(0); 1515 } 1516 1517 1518 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "PCFieldSplitSetIS" 1521 /*@C 1522 PCFieldSplitSetIS - Sets the exact elements for field 1523 1524 Logically Collective on PC 1525 1526 Input Parameters: 1527 + pc - the preconditioner context 1528 . splitname - name of this split, if NULL the number of the split is used 1529 - is - the index set that defines the vector elements in this field 1530 1531 1532 Notes: 1533 Use PCFieldSplitSetFields(), for fields defined by strided types. 1534 1535 This function is called once per split (it creates a new split each time). Solve options 1536 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1537 1538 Level: intermediate 1539 1540 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1541 1542 @*/ 1543 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1544 { 1545 PetscErrorCode ierr; 1546 1547 PetscFunctionBegin; 1548 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1549 if (splitname) PetscValidCharPointer(splitname,2); 1550 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1551 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1552 PetscFunctionReturn(0); 1553 } 1554 1555 #undef __FUNCT__ 1556 #define __FUNCT__ "PCFieldSplitGetIS" 1557 /*@ 1558 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1559 1560 Logically Collective on PC 1561 1562 Input Parameters: 1563 + pc - the preconditioner context 1564 - splitname - name of this split 1565 1566 Output Parameter: 1567 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1568 1569 Level: intermediate 1570 1571 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1572 1573 @*/ 1574 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1575 { 1576 PetscErrorCode ierr; 1577 1578 PetscFunctionBegin; 1579 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1580 PetscValidCharPointer(splitname,2); 1581 PetscValidPointer(is,3); 1582 { 1583 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1584 PC_FieldSplitLink ilink = jac->head; 1585 PetscBool found; 1586 1587 *is = NULL; 1588 while (ilink) { 1589 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1590 if (found) { 1591 *is = ilink->is; 1592 break; 1593 } 1594 ilink = ilink->next; 1595 } 1596 } 1597 PetscFunctionReturn(0); 1598 } 1599 1600 #undef __FUNCT__ 1601 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1602 /*@ 1603 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1604 fieldsplit preconditioner. If not set the matrix block size is used. 1605 1606 Logically Collective on PC 1607 1608 Input Parameters: 1609 + pc - the preconditioner context 1610 - bs - the block size 1611 1612 Level: intermediate 1613 1614 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1615 1616 @*/ 1617 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1618 { 1619 PetscErrorCode ierr; 1620 1621 PetscFunctionBegin; 1622 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1623 PetscValidLogicalCollectiveInt(pc,bs,2); 1624 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1625 PetscFunctionReturn(0); 1626 } 1627 1628 #undef __FUNCT__ 1629 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1630 /*@C 1631 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1632 1633 Collective on KSP 1634 1635 Input Parameter: 1636 . pc - the preconditioner context 1637 1638 Output Parameters: 1639 + n - the number of splits 1640 - pc - the array of KSP contexts 1641 1642 Note: 1643 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1644 (not the KSP just the array that contains them). 1645 1646 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1647 1648 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1649 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1650 KSP array must be. 1651 1652 1653 Level: advanced 1654 1655 .seealso: PCFIELDSPLIT 1656 @*/ 1657 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1658 { 1659 PetscErrorCode ierr; 1660 1661 PetscFunctionBegin; 1662 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1663 if (n) PetscValidIntPointer(n,2); 1664 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 #undef __FUNCT__ 1669 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1670 /*@ 1671 PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1672 A11 matrix. Otherwise no preconditioner is used. 1673 1674 Collective on PC 1675 1676 Input Parameters: 1677 + pc - the preconditioner context 1678 . 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 1679 - userpre - matrix to use for preconditioning, or NULL 1680 1681 Options Database: 1682 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11 1683 1684 Notes: 1685 $ If ptype is 1686 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1687 $ to this function). 1688 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the preconditioner 1689 $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix 1690 $ full then the preconditioner uses the exact Schur complement (this is expensive) 1691 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1692 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1693 $ preconditioner 1694 $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1695 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 1696 $ lumped before extracting the diagonal: -fieldsplit_1_mat_schur_complement_ainv_type lump; diag is the default. 1697 1698 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1699 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1700 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1701 1702 Level: intermediate 1703 1704 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1705 MatSchurComplementSetAinvType(), PCLSC 1706 1707 @*/ 1708 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1709 { 1710 PetscErrorCode ierr; 1711 1712 PetscFunctionBegin; 1713 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1714 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1715 PetscFunctionReturn(0); 1716 } 1717 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1718 1719 #undef __FUNCT__ 1720 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1721 /*@ 1722 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1723 preconditioned. See PCFieldSplitSetSchurPre() for details. 1724 1725 Logically Collective on PC 1726 1727 Input Parameters: 1728 . pc - the preconditioner context 1729 1730 Output Parameters: 1731 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1732 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1733 1734 Level: intermediate 1735 1736 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1737 1738 @*/ 1739 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1740 { 1741 PetscErrorCode ierr; 1742 1743 PetscFunctionBegin; 1744 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1745 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1746 PetscFunctionReturn(0); 1747 } 1748 1749 #undef __FUNCT__ 1750 #define __FUNCT__ "PCFieldSplitSchurGetS" 1751 /*@ 1752 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1753 1754 Not collective 1755 1756 Input Parameter: 1757 . pc - the preconditioner context 1758 1759 Output Parameter: 1760 . S - the Schur complement matrix 1761 1762 Notes: 1763 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1764 1765 Level: advanced 1766 1767 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1768 1769 @*/ 1770 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1771 { 1772 PetscErrorCode ierr; 1773 const char* t; 1774 PetscBool isfs; 1775 PC_FieldSplit *jac; 1776 1777 PetscFunctionBegin; 1778 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1779 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1780 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1781 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1782 jac = (PC_FieldSplit*)pc->data; 1783 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1784 if (S) *S = jac->schur; 1785 PetscFunctionReturn(0); 1786 } 1787 1788 #undef __FUNCT__ 1789 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1790 /*@ 1791 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1792 1793 Not collective 1794 1795 Input Parameters: 1796 + pc - the preconditioner context 1797 . S - the Schur complement matrix 1798 1799 Level: advanced 1800 1801 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 1802 1803 @*/ 1804 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1805 { 1806 PetscErrorCode ierr; 1807 const char* t; 1808 PetscBool isfs; 1809 PC_FieldSplit *jac; 1810 1811 PetscFunctionBegin; 1812 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1813 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1814 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1815 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1816 jac = (PC_FieldSplit*)pc->data; 1817 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1818 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1819 PetscFunctionReturn(0); 1820 } 1821 1822 1823 #undef __FUNCT__ 1824 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 1825 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1826 { 1827 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1828 PetscErrorCode ierr; 1829 1830 PetscFunctionBegin; 1831 jac->schurpre = ptype; 1832 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 1833 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1834 jac->schur_user = pre; 1835 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1836 } 1837 PetscFunctionReturn(0); 1838 } 1839 1840 #undef __FUNCT__ 1841 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 1842 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1843 { 1844 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1845 1846 PetscFunctionBegin; 1847 *ptype = jac->schurpre; 1848 *pre = jac->schur_user; 1849 PetscFunctionReturn(0); 1850 } 1851 1852 #undef __FUNCT__ 1853 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1854 /*@ 1855 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1856 1857 Collective on PC 1858 1859 Input Parameters: 1860 + pc - the preconditioner context 1861 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1862 1863 Options Database: 1864 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1865 1866 1867 Level: intermediate 1868 1869 Notes: 1870 The FULL factorization is 1871 1872 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1873 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1874 1875 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, 1876 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). 1877 1878 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 1879 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 1880 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, 1881 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1882 1883 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 1884 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). 1885 1886 References: 1887 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1888 1889 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1890 1891 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1892 @*/ 1893 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1894 { 1895 PetscErrorCode ierr; 1896 1897 PetscFunctionBegin; 1898 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1899 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1900 PetscFunctionReturn(0); 1901 } 1902 1903 #undef __FUNCT__ 1904 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1905 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1906 { 1907 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1908 1909 PetscFunctionBegin; 1910 jac->schurfactorization = ftype; 1911 PetscFunctionReturn(0); 1912 } 1913 1914 #undef __FUNCT__ 1915 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1916 /*@C 1917 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1918 1919 Collective on KSP 1920 1921 Input Parameter: 1922 . pc - the preconditioner context 1923 1924 Output Parameters: 1925 + A00 - the (0,0) block 1926 . A01 - the (0,1) block 1927 . A10 - the (1,0) block 1928 - A11 - the (1,1) block 1929 1930 Level: advanced 1931 1932 .seealso: PCFIELDSPLIT 1933 @*/ 1934 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1935 { 1936 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1937 1938 PetscFunctionBegin; 1939 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1940 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1941 if (A00) *A00 = jac->pmat[0]; 1942 if (A01) *A01 = jac->B; 1943 if (A10) *A10 = jac->C; 1944 if (A11) *A11 = jac->pmat[1]; 1945 PetscFunctionReturn(0); 1946 } 1947 1948 #undef __FUNCT__ 1949 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1950 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1951 { 1952 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1953 PetscErrorCode ierr; 1954 1955 PetscFunctionBegin; 1956 jac->type = type; 1957 if (type == PC_COMPOSITE_SCHUR) { 1958 pc->ops->apply = PCApply_FieldSplit_Schur; 1959 pc->ops->view = PCView_FieldSplit_Schur; 1960 1961 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1962 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 1963 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 1964 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1965 1966 } else { 1967 pc->ops->apply = PCApply_FieldSplit; 1968 pc->ops->view = PCView_FieldSplit; 1969 1970 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1971 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 1972 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 1973 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1974 } 1975 PetscFunctionReturn(0); 1976 } 1977 1978 #undef __FUNCT__ 1979 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1980 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1981 { 1982 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1983 1984 PetscFunctionBegin; 1985 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1986 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); 1987 jac->bs = bs; 1988 PetscFunctionReturn(0); 1989 } 1990 1991 #undef __FUNCT__ 1992 #define __FUNCT__ "PCFieldSplitSetType" 1993 /*@ 1994 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1995 1996 Collective on PC 1997 1998 Input Parameter: 1999 . pc - the preconditioner context 2000 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2001 2002 Options Database Key: 2003 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2004 2005 Level: Intermediate 2006 2007 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2008 2009 .seealso: PCCompositeSetType() 2010 2011 @*/ 2012 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2013 { 2014 PetscErrorCode ierr; 2015 2016 PetscFunctionBegin; 2017 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2018 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2019 PetscFunctionReturn(0); 2020 } 2021 2022 #undef __FUNCT__ 2023 #define __FUNCT__ "PCFieldSplitGetType" 2024 /*@ 2025 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2026 2027 Not collective 2028 2029 Input Parameter: 2030 . pc - the preconditioner context 2031 2032 Output Parameter: 2033 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2034 2035 Level: Intermediate 2036 2037 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2038 .seealso: PCCompositeSetType() 2039 @*/ 2040 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2041 { 2042 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2043 2044 PetscFunctionBegin; 2045 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2046 PetscValidIntPointer(type,2); 2047 *type = jac->type; 2048 PetscFunctionReturn(0); 2049 } 2050 2051 #undef __FUNCT__ 2052 #define __FUNCT__ "PCFieldSplitSetDMSplits" 2053 /*@ 2054 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2055 2056 Logically Collective 2057 2058 Input Parameters: 2059 + pc - the preconditioner context 2060 - flg - boolean indicating whether to use field splits defined by the DM 2061 2062 Options Database Key: 2063 . -pc_fieldsplit_dm_splits 2064 2065 Level: Intermediate 2066 2067 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2068 2069 .seealso: PCFieldSplitGetDMSplits() 2070 2071 @*/ 2072 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2073 { 2074 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2075 PetscBool isfs; 2076 PetscErrorCode ierr; 2077 2078 PetscFunctionBegin; 2079 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2080 PetscValidLogicalCollectiveBool(pc,flg,2); 2081 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2082 if (isfs) { 2083 jac->dm_splits = flg; 2084 } 2085 PetscFunctionReturn(0); 2086 } 2087 2088 2089 #undef __FUNCT__ 2090 #define __FUNCT__ "PCFieldSplitGetDMSplits" 2091 /*@ 2092 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2093 2094 Logically Collective 2095 2096 Input Parameter: 2097 . pc - the preconditioner context 2098 2099 Output Parameter: 2100 . flg - boolean indicating whether to use field splits defined by the DM 2101 2102 Level: Intermediate 2103 2104 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2105 2106 .seealso: PCFieldSplitSetDMSplits() 2107 2108 @*/ 2109 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2110 { 2111 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2112 PetscBool isfs; 2113 PetscErrorCode ierr; 2114 2115 PetscFunctionBegin; 2116 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2117 PetscValidPointer(flg,2); 2118 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2119 if (isfs) { 2120 if(flg) *flg = jac->dm_splits; 2121 } 2122 PetscFunctionReturn(0); 2123 } 2124 2125 /* -------------------------------------------------------------------------------------*/ 2126 /*MC 2127 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2128 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2129 2130 To set options on the solvers for each block append -fieldsplit_ to all the PC 2131 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2132 2133 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2134 and set the options directly on the resulting KSP object 2135 2136 Level: intermediate 2137 2138 Options Database Keys: 2139 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2140 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2141 been supplied explicitly by -pc_fieldsplit_%d_fields 2142 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2143 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2144 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11 2145 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2146 2147 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2148 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2149 2150 Notes: 2151 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2152 to define a field by an arbitrary collection of entries. 2153 2154 If no fields are set the default is used. The fields are defined by entries strided by bs, 2155 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2156 if this is not called the block size defaults to the blocksize of the second matrix passed 2157 to KSPSetOperators()/PCSetOperators(). 2158 2159 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2160 $ ( A10 A11 ) 2161 $ the preconditioner using full factorization is 2162 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2163 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2164 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2165 $ S = A11 - A10 ksp(A00) A01 2166 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 2167 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2168 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2169 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this 2170 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. 2171 When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled -- 2172 Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners 2173 (e.g., -fieldsplit_1_pc_type asm). Optionally, A00 can be lumped before extracting the diagonal: 2174 -fieldsplit_1_mat_schur_complement_ainv_type lump; diag is the default. 2175 2176 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2177 diag gives 2178 $ ( inv(A00) 0 ) 2179 $ ( 0 -ksp(S) ) 2180 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 2181 $ ( A00 0 ) 2182 $ ( A10 S ) 2183 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2184 $ ( A00 A01 ) 2185 $ ( 0 S ) 2186 where again the inverses of A00 and S are applied using KSPs. 2187 2188 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2189 is used automatically for a second block. 2190 2191 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2192 Generally it should be used with the AIJ format. 2193 2194 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2195 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2196 inside a smoother resulting in "Distributive Smoothers". 2197 2198 Concepts: physics based preconditioners, block preconditioners 2199 2200 There is a nice discussion of block preconditioners in 2201 2202 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2203 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2204 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2205 2206 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2207 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2208 MatSchurComplementSetAinvType() 2209 M*/ 2210 2211 #undef __FUNCT__ 2212 #define __FUNCT__ "PCCreate_FieldSplit" 2213 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2214 { 2215 PetscErrorCode ierr; 2216 PC_FieldSplit *jac; 2217 2218 PetscFunctionBegin; 2219 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2220 2221 jac->bs = -1; 2222 jac->nsplits = 0; 2223 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2224 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2225 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2226 jac->dm_splits = PETSC_TRUE; 2227 2228 pc->data = (void*)jac; 2229 2230 pc->ops->apply = PCApply_FieldSplit; 2231 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2232 pc->ops->setup = PCSetUp_FieldSplit; 2233 pc->ops->reset = PCReset_FieldSplit; 2234 pc->ops->destroy = PCDestroy_FieldSplit; 2235 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2236 pc->ops->view = PCView_FieldSplit; 2237 pc->ops->applyrichardson = 0; 2238 2239 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2240 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2241 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2242 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2243 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2244 PetscFunctionReturn(0); 2245 } 2246 2247 2248 2249