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