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