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