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