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