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