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