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