1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 3 #include <petscdm.h> 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",NULL}; 6 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",NULL}; 7 8 PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4; 9 10 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 11 struct _PC_FieldSplitLink { 12 KSP ksp; 13 Vec x,y,z; 14 char *splitname; 15 PetscInt nfields; 16 PetscInt *fields,*fields_col; 17 VecScatter sctx; 18 IS is,is_col; 19 PC_FieldSplitLink next,previous; 20 PetscLogEvent event; 21 22 /* Used only when setting coordinates with PCSetCoordinates */ 23 PetscInt dim; 24 PetscInt ndofs; 25 PetscReal *coords; 26 }; 27 28 typedef struct { 29 PCCompositeType type; 30 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 31 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 32 PetscInt bs; /* Block size for IS and Mat structures */ 33 PetscInt nsplits; /* Number of field divisions defined */ 34 Vec *x,*y,w1,w2; 35 Mat *mat; /* The diagonal block for each split */ 36 Mat *pmat; /* The preconditioning diagonal block for each split */ 37 Mat *Afield; /* The rows of the matrix associated with each split */ 38 PetscBool issetup; 39 40 /* Only used when Schur complement preconditioning is used */ 41 Mat B; /* The (0,1) block */ 42 Mat C; /* The (1,0) block */ 43 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 44 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 45 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 46 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 47 PCFieldSplitSchurFactType schurfactorization; 48 KSP kspschur; /* The solver for S */ 49 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 50 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 51 52 /* Only used when Golub-Kahan bidiagonalization preconditioning is used */ 53 Mat H; /* The modified matrix H = A00 + nu*A01*A01' */ 54 PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */ 55 PetscInt gkbdelay; /* The delay window for the stopping criterion */ 56 PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */ 57 PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */ 58 PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */ 59 PetscViewer gkbviewer; /* Viewer context for gkbmonitor */ 60 Vec u,v,d,Hu; /* Work vectors for the GKB algorithm */ 61 PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */ 62 63 PC_FieldSplitLink head; 64 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 65 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 66 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 67 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 68 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 69 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 70 PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */ 71 } PC_FieldSplit; 72 73 /* 74 Notes: 75 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 76 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 77 PC you could change this. 78 */ 79 80 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 81 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 82 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 83 { 84 switch (jac->schurpre) { 85 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 86 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 87 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 88 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 89 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 90 default: 91 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 92 } 93 } 94 95 #include <petscdraw.h> 96 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 97 { 98 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 99 PetscBool iascii,isdraw; 100 PetscInt i,j; 101 PC_FieldSplitLink ilink = jac->head; 102 103 PetscFunctionBegin; 104 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 105 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 106 if (iascii) { 107 if (jac->bs > 0) { 108 CHKERRQ(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs)); 109 } else { 110 CHKERRQ(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits)); 111 } 112 if (pc->useAmat) { 113 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n")); 114 } 115 if (jac->diag_use_amat) { 116 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n")); 117 } 118 if (jac->offdiag_use_amat) { 119 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n")); 120 } 121 CHKERRQ(PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n")); 122 for (i=0; i<jac->nsplits; i++) { 123 if (ilink->fields) { 124 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i)); 125 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 126 for (j=0; j<ilink->nfields; j++) { 127 if (j > 0) { 128 CHKERRQ(PetscViewerASCIIPrintf(viewer,",")); 129 } 130 CHKERRQ(PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j])); 131 } 132 CHKERRQ(PetscViewerASCIIPrintf(viewer,"\n")); 133 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 134 } else { 135 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i)); 136 } 137 CHKERRQ(KSPView(ilink->ksp,viewer)); 138 ilink = ilink->next; 139 } 140 } 141 142 if (isdraw) { 143 PetscDraw draw; 144 PetscReal x,y,w,wd; 145 146 CHKERRQ(PetscViewerDrawGetDraw(viewer,0,&draw)); 147 CHKERRQ(PetscDrawGetCurrentPoint(draw,&x,&y)); 148 w = 2*PetscMin(1.0 - x,x); 149 wd = w/(jac->nsplits + 1); 150 x = x - wd*(jac->nsplits-1)/2.0; 151 for (i=0; i<jac->nsplits; i++) { 152 CHKERRQ(PetscDrawPushCurrentPoint(draw,x,y)); 153 CHKERRQ(KSPView(ilink->ksp,viewer)); 154 CHKERRQ(PetscDrawPopCurrentPoint(draw)); 155 x += wd; 156 ilink = ilink->next; 157 } 158 } 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 163 { 164 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 165 PetscBool iascii,isdraw; 166 PetscInt i,j; 167 PC_FieldSplitLink ilink = jac->head; 168 MatSchurComplementAinvType atype; 169 170 PetscFunctionBegin; 171 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 172 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 173 if (iascii) { 174 if (jac->bs > 0) { 175 CHKERRQ(PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization])); 176 } else { 177 CHKERRQ(PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization])); 178 } 179 if (pc->useAmat) { 180 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n")); 181 } 182 switch (jac->schurpre) { 183 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 184 CHKERRQ(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n")); 185 break; 186 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 187 CHKERRQ(MatSchurComplementGetAinvType(jac->schur,&atype)); 188 CHKERRQ(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped ")));break; 189 case PC_FIELDSPLIT_SCHUR_PRE_A11: 190 CHKERRQ(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n")); 191 break; 192 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 193 CHKERRQ(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n")); 194 break; 195 case PC_FIELDSPLIT_SCHUR_PRE_USER: 196 if (jac->schur_user) { 197 CHKERRQ(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n")); 198 } else { 199 CHKERRQ(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n")); 200 } 201 break; 202 default: 203 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 204 } 205 CHKERRQ(PetscViewerASCIIPrintf(viewer," Split info:\n")); 206 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 207 for (i=0; i<jac->nsplits; i++) { 208 if (ilink->fields) { 209 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i)); 210 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 211 for (j=0; j<ilink->nfields; j++) { 212 if (j > 0) { 213 CHKERRQ(PetscViewerASCIIPrintf(viewer,",")); 214 } 215 CHKERRQ(PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j])); 216 } 217 CHKERRQ(PetscViewerASCIIPrintf(viewer,"\n")); 218 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 219 } else { 220 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i)); 221 } 222 ilink = ilink->next; 223 } 224 CHKERRQ(PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n")); 225 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 226 if (jac->head) { 227 CHKERRQ(KSPView(jac->head->ksp,viewer)); 228 } else CHKERRQ(PetscViewerASCIIPrintf(viewer," not yet available\n")); 229 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 230 if (jac->head && jac->kspupper != jac->head->ksp) { 231 CHKERRQ(PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n")); 232 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 233 if (jac->kspupper) CHKERRQ(KSPView(jac->kspupper,viewer)); 234 else CHKERRQ(PetscViewerASCIIPrintf(viewer," not yet available\n")); 235 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 236 } 237 CHKERRQ(PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n")); 238 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 239 if (jac->kspschur) { 240 CHKERRQ(KSPView(jac->kspschur,viewer)); 241 } else { 242 CHKERRQ(PetscViewerASCIIPrintf(viewer," not yet available\n")); 243 } 244 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 245 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 246 } else if (isdraw && jac->head) { 247 PetscDraw draw; 248 PetscReal x,y,w,wd,h; 249 PetscInt cnt = 2; 250 char str[32]; 251 252 CHKERRQ(PetscViewerDrawGetDraw(viewer,0,&draw)); 253 CHKERRQ(PetscDrawGetCurrentPoint(draw,&x,&y)); 254 if (jac->kspupper != jac->head->ksp) cnt++; 255 w = 2*PetscMin(1.0 - x,x); 256 wd = w/(cnt + 1); 257 258 CHKERRQ(PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization])); 259 CHKERRQ(PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h)); 260 y -= h; 261 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 262 CHKERRQ(PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11])); 263 } else { 264 CHKERRQ(PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre])); 265 } 266 CHKERRQ(PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h)); 267 y -= h; 268 x = x - wd*(cnt-1)/2.0; 269 270 CHKERRQ(PetscDrawPushCurrentPoint(draw,x,y)); 271 CHKERRQ(KSPView(jac->head->ksp,viewer)); 272 CHKERRQ(PetscDrawPopCurrentPoint(draw)); 273 if (jac->kspupper != jac->head->ksp) { 274 x += wd; 275 CHKERRQ(PetscDrawPushCurrentPoint(draw,x,y)); 276 CHKERRQ(KSPView(jac->kspupper,viewer)); 277 CHKERRQ(PetscDrawPopCurrentPoint(draw)); 278 } 279 x += wd; 280 CHKERRQ(PetscDrawPushCurrentPoint(draw,x,y)); 281 CHKERRQ(KSPView(jac->kspschur,viewer)); 282 CHKERRQ(PetscDrawPopCurrentPoint(draw)); 283 } 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer) 288 { 289 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 290 PetscBool iascii,isdraw; 291 PetscInt i,j; 292 PC_FieldSplitLink ilink = jac->head; 293 294 PetscFunctionBegin; 295 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 296 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 297 if (iascii) { 298 if (jac->bs > 0) { 299 CHKERRQ(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs)); 300 } else { 301 CHKERRQ(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits)); 302 } 303 if (pc->useAmat) { 304 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n")); 305 } 306 if (jac->diag_use_amat) { 307 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n")); 308 } 309 if (jac->offdiag_use_amat) { 310 CHKERRQ(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n")); 311 } 312 313 CHKERRQ(PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%D, maximum iterations=%D\n",jac->gkbtol,jac->gkbdelay,jac->gkbmaxit)); 314 CHKERRQ(PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n")); 315 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 316 317 if (ilink->fields) { 318 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0)); 319 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 320 for (j=0; j<ilink->nfields; j++) { 321 if (j > 0) { 322 CHKERRQ(PetscViewerASCIIPrintf(viewer,",")); 323 } 324 CHKERRQ(PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j])); 325 } 326 CHKERRQ(PetscViewerASCIIPrintf(viewer,"\n")); 327 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 328 } else { 329 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0)); 330 } 331 CHKERRQ(KSPView(ilink->ksp,viewer)); 332 333 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 334 } 335 336 if (isdraw) { 337 PetscDraw draw; 338 PetscReal x,y,w,wd; 339 340 CHKERRQ(PetscViewerDrawGetDraw(viewer,0,&draw)); 341 CHKERRQ(PetscDrawGetCurrentPoint(draw,&x,&y)); 342 w = 2*PetscMin(1.0 - x,x); 343 wd = w/(jac->nsplits + 1); 344 x = x - wd*(jac->nsplits-1)/2.0; 345 for (i=0; i<jac->nsplits; i++) { 346 CHKERRQ(PetscDrawPushCurrentPoint(draw,x,y)); 347 CHKERRQ(KSPView(ilink->ksp,viewer)); 348 CHKERRQ(PetscDrawPopCurrentPoint(draw)); 349 x += wd; 350 ilink = ilink->next; 351 } 352 } 353 PetscFunctionReturn(0); 354 } 355 356 /* Precondition: jac->bs is set to a meaningful value */ 357 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 358 { 359 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 360 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 361 PetscBool flg,flg_col; 362 char optionname[128],splitname[8],optionname_col[128]; 363 364 PetscFunctionBegin; 365 CHKERRQ(PetscMalloc1(jac->bs,&ifields)); 366 CHKERRQ(PetscMalloc1(jac->bs,&ifields_col)); 367 for (i=0,flg=PETSC_TRUE;; i++) { 368 CHKERRQ(PetscSNPrintf(splitname,sizeof(splitname),"%D",i)); 369 CHKERRQ(PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i)); 370 CHKERRQ(PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i)); 371 nfields = jac->bs; 372 nfields_col = jac->bs; 373 CHKERRQ(PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg)); 374 CHKERRQ(PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col)); 375 if (!flg) break; 376 else if (flg && !flg_col) { 377 PetscCheck(nfields,PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 378 CHKERRQ(PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields)); 379 } else { 380 PetscCheckFalse(!nfields || !nfields_col,PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 381 PetscCheckFalse(nfields != nfields_col,PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 382 CHKERRQ(PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col)); 383 } 384 } 385 if (i > 0) { 386 /* Makes command-line setting of splits take precedence over setting them in code. 387 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 388 create new splits, which would probably not be what the user wanted. */ 389 jac->splitdefined = PETSC_TRUE; 390 } 391 CHKERRQ(PetscFree(ifields)); 392 CHKERRQ(PetscFree(ifields_col)); 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 397 { 398 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 399 PC_FieldSplitLink ilink = jac->head; 400 PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE; 401 PetscInt i; 402 403 PetscFunctionBegin; 404 /* 405 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 406 Should probably be rewritten. 407 */ 408 if (!ilink) { 409 CHKERRQ(PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL)); 410 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 411 PetscInt numFields, f, i, j; 412 char **fieldNames; 413 IS *fields; 414 DM *dms; 415 DM subdm[128]; 416 PetscBool flg; 417 418 CHKERRQ(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms)); 419 /* Allow the user to prescribe the splits */ 420 for (i = 0, flg = PETSC_TRUE;; i++) { 421 PetscInt ifields[128]; 422 IS compField; 423 char optionname[128], splitname[8]; 424 PetscInt nfields = numFields; 425 426 CHKERRQ(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i)); 427 CHKERRQ(PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg)); 428 if (!flg) break; 429 PetscCheckFalse(numFields > 128,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 430 CHKERRQ(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i])); 431 if (nfields == 1) { 432 CHKERRQ(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField)); 433 } else { 434 CHKERRQ(PetscSNPrintf(splitname, sizeof(splitname), "%D", i)); 435 CHKERRQ(PCFieldSplitSetIS(pc, splitname, compField)); 436 } 437 CHKERRQ(ISDestroy(&compField)); 438 for (j = 0; j < nfields; ++j) { 439 f = ifields[j]; 440 CHKERRQ(PetscFree(fieldNames[f])); 441 CHKERRQ(ISDestroy(&fields[f])); 442 } 443 } 444 if (i == 0) { 445 for (f = 0; f < numFields; ++f) { 446 CHKERRQ(PCFieldSplitSetIS(pc, fieldNames[f], fields[f])); 447 CHKERRQ(PetscFree(fieldNames[f])); 448 CHKERRQ(ISDestroy(&fields[f])); 449 } 450 } else { 451 for (j=0; j<numFields; j++) { 452 CHKERRQ(DMDestroy(dms+j)); 453 } 454 CHKERRQ(PetscFree(dms)); 455 CHKERRQ(PetscMalloc1(i, &dms)); 456 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 457 } 458 CHKERRQ(PetscFree(fieldNames)); 459 CHKERRQ(PetscFree(fields)); 460 if (dms) { 461 CHKERRQ(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n")); 462 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 463 const char *prefix; 464 CHKERRQ(PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix)); 465 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix)); 466 CHKERRQ(KSPSetDM(ilink->ksp, dms[i])); 467 CHKERRQ(KSPSetDMActive(ilink->ksp, PETSC_FALSE)); 468 { 469 PetscErrorCode (*func)(KSP,Mat,Mat,void*); 470 void *ctx; 471 472 CHKERRQ(DMKSPGetComputeOperators(pc->dm, &func, &ctx)); 473 CHKERRQ(DMKSPSetComputeOperators(dms[i], func, ctx)); 474 } 475 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0)); 476 CHKERRQ(DMDestroy(&dms[i])); 477 } 478 CHKERRQ(PetscFree(dms)); 479 } 480 } else { 481 if (jac->bs <= 0) { 482 if (pc->pmat) { 483 CHKERRQ(MatGetBlockSize(pc->pmat,&jac->bs)); 484 } else jac->bs = 1; 485 } 486 487 if (jac->detect) { 488 IS zerodiags,rest; 489 PetscInt nmin,nmax; 490 491 CHKERRQ(MatGetOwnershipRange(pc->mat,&nmin,&nmax)); 492 if (jac->diag_use_amat) { 493 CHKERRQ(MatFindZeroDiagonals(pc->mat,&zerodiags)); 494 } else { 495 CHKERRQ(MatFindZeroDiagonals(pc->pmat,&zerodiags)); 496 } 497 CHKERRQ(ISComplement(zerodiags,nmin,nmax,&rest)); 498 CHKERRQ(PCFieldSplitSetIS(pc,"0",rest)); 499 CHKERRQ(PCFieldSplitSetIS(pc,"1",zerodiags)); 500 CHKERRQ(ISDestroy(&zerodiags)); 501 CHKERRQ(ISDestroy(&rest)); 502 } else if (coupling) { 503 IS coupling,rest; 504 PetscInt nmin,nmax; 505 506 CHKERRQ(MatGetOwnershipRange(pc->mat,&nmin,&nmax)); 507 if (jac->offdiag_use_amat) { 508 CHKERRQ(MatFindOffBlockDiagonalEntries(pc->mat,&coupling)); 509 } else { 510 CHKERRQ(MatFindOffBlockDiagonalEntries(pc->pmat,&coupling)); 511 } 512 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest)); 513 CHKERRQ(ISSetIdentity(rest)); 514 CHKERRQ(PCFieldSplitSetIS(pc,"0",rest)); 515 CHKERRQ(PCFieldSplitSetIS(pc,"1",coupling)); 516 CHKERRQ(ISDestroy(&coupling)); 517 CHKERRQ(ISDestroy(&rest)); 518 } else { 519 CHKERRQ(PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL)); 520 if (!fieldsplit_default) { 521 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 522 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 523 CHKERRQ(PCFieldSplitSetRuntimeSplits_Private(pc)); 524 if (jac->splitdefined) CHKERRQ(PetscInfo(pc,"Splits defined using the options database\n")); 525 } 526 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 527 Mat M = pc->pmat; 528 PetscBool isnest; 529 530 CHKERRQ(PetscInfo(pc,"Using default splitting of fields\n")); 531 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest)); 532 if (!isnest) { 533 M = pc->mat; 534 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest)); 535 } 536 if (isnest) { 537 IS *fields; 538 PetscInt nf; 539 540 CHKERRQ(MatNestGetSize(M,&nf,NULL)); 541 CHKERRQ(PetscMalloc1(nf,&fields)); 542 CHKERRQ(MatNestGetISs(M,fields,NULL)); 543 for (i=0;i<nf;i++) { 544 CHKERRQ(PCFieldSplitSetIS(pc,NULL,fields[i])); 545 } 546 CHKERRQ(PetscFree(fields)); 547 } else { 548 for (i=0; i<jac->bs; i++) { 549 char splitname[8]; 550 CHKERRQ(PetscSNPrintf(splitname,sizeof(splitname),"%D",i)); 551 CHKERRQ(PCFieldSplitSetFields(pc,splitname,1,&i,&i)); 552 } 553 jac->defaultsplit = PETSC_TRUE; 554 } 555 } 556 } 557 } 558 } else if (jac->nsplits == 1) { 559 if (ilink->is) { 560 IS is2; 561 PetscInt nmin,nmax; 562 563 CHKERRQ(MatGetOwnershipRange(pc->mat,&nmin,&nmax)); 564 CHKERRQ(ISComplement(ilink->is,nmin,nmax,&is2)); 565 CHKERRQ(PCFieldSplitSetIS(pc,"1",is2)); 566 CHKERRQ(ISDestroy(&is2)); 567 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 568 } 569 570 PetscCheckFalse(jac->nsplits < 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 571 PetscFunctionReturn(0); 572 } 573 574 static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu) 575 { 576 Mat BT,T; 577 PetscReal nrmT,nrmB; 578 579 PetscFunctionBegin; 580 CHKERRQ(MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T)); /* Test if augmented matrix is symmetric */ 581 CHKERRQ(MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN)); 582 CHKERRQ(MatNorm(T,NORM_1,&nrmT)); 583 CHKERRQ(MatNorm(B,NORM_1,&nrmB)); 584 if (nrmB > 0) { 585 if (nrmT/nrmB >= PETSC_SMALL) { 586 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable."); 587 } 588 } 589 /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */ 590 /* setting N := 1/nu*I in [Ar13]. */ 591 CHKERRQ(MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT)); 592 CHKERRQ(MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H)); /* H = A01*A01' */ 593 CHKERRQ(MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN)); /* H = A00 + nu*A01*A01' */ 594 595 CHKERRQ(MatDestroy(&BT)); 596 CHKERRQ(MatDestroy(&T)); 597 PetscFunctionReturn(0); 598 } 599 600 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg); 601 602 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 603 { 604 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 605 PC_FieldSplitLink ilink; 606 PetscInt i,nsplit; 607 PetscBool sorted, sorted_col; 608 609 PetscFunctionBegin; 610 pc->failedreason = PC_NOERROR; 611 CHKERRQ(PCFieldSplitSetDefaults(pc)); 612 nsplit = jac->nsplits; 613 ilink = jac->head; 614 615 /* get the matrices for each split */ 616 if (!jac->issetup) { 617 PetscInt rstart,rend,nslots,bs; 618 619 jac->issetup = PETSC_TRUE; 620 621 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 622 if (jac->defaultsplit || !ilink->is) { 623 if (jac->bs <= 0) jac->bs = nsplit; 624 } 625 bs = jac->bs; 626 CHKERRQ(MatGetOwnershipRange(pc->pmat,&rstart,&rend)); 627 nslots = (rend - rstart)/bs; 628 for (i=0; i<nsplit; i++) { 629 if (jac->defaultsplit) { 630 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is)); 631 CHKERRQ(ISDuplicate(ilink->is,&ilink->is_col)); 632 } else if (!ilink->is) { 633 if (ilink->nfields > 1) { 634 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 635 CHKERRQ(PetscMalloc1(ilink->nfields*nslots,&ii)); 636 CHKERRQ(PetscMalloc1(ilink->nfields*nslots,&jj)); 637 for (j=0; j<nslots; j++) { 638 for (k=0; k<nfields; k++) { 639 ii[nfields*j + k] = rstart + bs*j + fields[k]; 640 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 641 } 642 } 643 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is)); 644 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col)); 645 CHKERRQ(ISSetBlockSize(ilink->is, nfields)); 646 CHKERRQ(ISSetBlockSize(ilink->is_col, nfields)); 647 } else { 648 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is)); 649 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col)); 650 } 651 } 652 CHKERRQ(ISSorted(ilink->is,&sorted)); 653 if (ilink->is_col) CHKERRQ(ISSorted(ilink->is_col,&sorted_col)); 654 PetscCheckFalse(!sorted || !sorted_col,PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 655 ilink = ilink->next; 656 } 657 } 658 659 ilink = jac->head; 660 if (!jac->pmat) { 661 Vec xtmp; 662 663 CHKERRQ(MatCreateVecs(pc->pmat,&xtmp,NULL)); 664 CHKERRQ(PetscMalloc1(nsplit,&jac->pmat)); 665 CHKERRQ(PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y)); 666 for (i=0; i<nsplit; i++) { 667 MatNullSpace sp; 668 669 /* Check for preconditioning matrix attached to IS */ 670 CHKERRQ(PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i])); 671 if (jac->pmat[i]) { 672 CHKERRQ(PetscObjectReference((PetscObject) jac->pmat[i])); 673 if (jac->type == PC_COMPOSITE_SCHUR) { 674 jac->schur_user = jac->pmat[i]; 675 676 CHKERRQ(PetscObjectReference((PetscObject) jac->schur_user)); 677 } 678 } else { 679 const char *prefix; 680 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i])); 681 CHKERRQ(KSPGetOptionsPrefix(ilink->ksp,&prefix)); 682 CHKERRQ(MatSetOptionsPrefix(jac->pmat[i],prefix)); 683 CHKERRQ(MatViewFromOptions(jac->pmat[i],NULL,"-mat_view")); 684 } 685 /* create work vectors for each split */ 686 CHKERRQ(MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i])); 687 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 688 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 689 CHKERRQ(VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx)); 690 CHKERRQ(PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp)); 691 if (sp) { 692 CHKERRQ(MatSetNearNullSpace(jac->pmat[i], sp)); 693 } 694 ilink = ilink->next; 695 } 696 CHKERRQ(VecDestroy(&xtmp)); 697 } else { 698 MatReuse scall; 699 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 700 for (i=0; i<nsplit; i++) { 701 CHKERRQ(MatDestroy(&jac->pmat[i])); 702 } 703 scall = MAT_INITIAL_MATRIX; 704 } else scall = MAT_REUSE_MATRIX; 705 706 for (i=0; i<nsplit; i++) { 707 Mat pmat; 708 709 /* Check for preconditioning matrix attached to IS */ 710 CHKERRQ(PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat)); 711 if (!pmat) { 712 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i])); 713 } 714 ilink = ilink->next; 715 } 716 } 717 if (jac->diag_use_amat) { 718 ilink = jac->head; 719 if (!jac->mat) { 720 CHKERRQ(PetscMalloc1(nsplit,&jac->mat)); 721 for (i=0; i<nsplit; i++) { 722 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i])); 723 ilink = ilink->next; 724 } 725 } else { 726 MatReuse scall; 727 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 728 for (i=0; i<nsplit; i++) { 729 CHKERRQ(MatDestroy(&jac->mat[i])); 730 } 731 scall = MAT_INITIAL_MATRIX; 732 } else scall = MAT_REUSE_MATRIX; 733 734 for (i=0; i<nsplit; i++) { 735 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i])); 736 ilink = ilink->next; 737 } 738 } 739 } else { 740 jac->mat = jac->pmat; 741 } 742 743 /* Check for null space attached to IS */ 744 ilink = jac->head; 745 for (i=0; i<nsplit; i++) { 746 MatNullSpace sp; 747 748 CHKERRQ(PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp)); 749 if (sp) { 750 CHKERRQ(MatSetNullSpace(jac->mat[i], sp)); 751 } 752 ilink = ilink->next; 753 } 754 755 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) { 756 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 757 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 758 ilink = jac->head; 759 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 760 /* 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 */ 761 if (!jac->Afield) { 762 CHKERRQ(PetscCalloc1(nsplit,&jac->Afield)); 763 if (jac->offdiag_use_amat) { 764 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1])); 765 } else { 766 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1])); 767 } 768 } else { 769 MatReuse scall; 770 771 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 772 CHKERRQ(MatDestroy(&jac->Afield[1])); 773 scall = MAT_INITIAL_MATRIX; 774 } else scall = MAT_REUSE_MATRIX; 775 776 if (jac->offdiag_use_amat) { 777 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1])); 778 } else { 779 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1])); 780 } 781 } 782 } else { 783 if (!jac->Afield) { 784 CHKERRQ(PetscMalloc1(nsplit,&jac->Afield)); 785 for (i=0; i<nsplit; i++) { 786 if (jac->offdiag_use_amat) { 787 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i])); 788 } else { 789 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i])); 790 } 791 ilink = ilink->next; 792 } 793 } else { 794 MatReuse scall; 795 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 796 for (i=0; i<nsplit; i++) { 797 CHKERRQ(MatDestroy(&jac->Afield[i])); 798 } 799 scall = MAT_INITIAL_MATRIX; 800 } else scall = MAT_REUSE_MATRIX; 801 802 for (i=0; i<nsplit; i++) { 803 if (jac->offdiag_use_amat) { 804 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i])); 805 } else { 806 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i])); 807 } 808 ilink = ilink->next; 809 } 810 } 811 } 812 } 813 814 if (jac->type == PC_COMPOSITE_SCHUR) { 815 IS ccis; 816 PetscBool isspd; 817 PetscInt rstart,rend; 818 char lscname[256]; 819 PetscObject LSC_L; 820 821 PetscCheckFalse(nsplit != 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 822 823 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 824 if (jac->schurscale == (PetscScalar)-1.0) { 825 CHKERRQ(MatGetOption(pc->pmat,MAT_SPD,&isspd)); 826 jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0; 827 } 828 829 /* When extracting off-diagonal submatrices, we take complements from this range */ 830 CHKERRQ(MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend)); 831 832 if (jac->schur) { 833 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 834 MatReuse scall; 835 836 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 837 scall = MAT_INITIAL_MATRIX; 838 CHKERRQ(MatDestroy(&jac->B)); 839 CHKERRQ(MatDestroy(&jac->C)); 840 } else scall = MAT_REUSE_MATRIX; 841 842 CHKERRQ(MatSchurComplementGetKSP(jac->schur, &kspInner)); 843 ilink = jac->head; 844 CHKERRQ(ISComplement(ilink->is_col,rstart,rend,&ccis)); 845 if (jac->offdiag_use_amat) { 846 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B)); 847 } else { 848 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B)); 849 } 850 CHKERRQ(ISDestroy(&ccis)); 851 ilink = ilink->next; 852 CHKERRQ(ISComplement(ilink->is_col,rstart,rend,&ccis)); 853 if (jac->offdiag_use_amat) { 854 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C)); 855 } else { 856 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C)); 857 } 858 CHKERRQ(ISDestroy(&ccis)); 859 CHKERRQ(MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1])); 860 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 861 CHKERRQ(MatDestroy(&jac->schurp)); 862 CHKERRQ(MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp)); 863 } 864 if (kspA != kspInner) { 865 CHKERRQ(KSPSetOperators(kspA,jac->mat[0],jac->pmat[0])); 866 } 867 if (kspUpper != kspA) { 868 CHKERRQ(KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0])); 869 } 870 CHKERRQ(KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac))); 871 } else { 872 const char *Dprefix; 873 char schurprefix[256], schurmatprefix[256]; 874 char schurtestoption[256]; 875 MatNullSpace sp; 876 PetscBool flg; 877 KSP kspt; 878 879 /* extract the A01 and A10 matrices */ 880 ilink = jac->head; 881 CHKERRQ(ISComplement(ilink->is_col,rstart,rend,&ccis)); 882 if (jac->offdiag_use_amat) { 883 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 884 } else { 885 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 886 } 887 CHKERRQ(ISDestroy(&ccis)); 888 ilink = ilink->next; 889 CHKERRQ(ISComplement(ilink->is_col,rstart,rend,&ccis)); 890 if (jac->offdiag_use_amat) { 891 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 892 } else { 893 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 894 } 895 CHKERRQ(ISDestroy(&ccis)); 896 897 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 898 CHKERRQ(MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur)); 899 CHKERRQ(MatSetType(jac->schur,MATSCHURCOMPLEMENT)); 900 CHKERRQ(MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1])); 901 CHKERRQ(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 902 CHKERRQ(MatSetOptionsPrefix(jac->schur,schurmatprefix)); 903 CHKERRQ(MatSchurComplementGetKSP(jac->schur,&kspt)); 904 CHKERRQ(KSPSetOptionsPrefix(kspt,schurmatprefix)); 905 906 /* Note: this is not true in general */ 907 CHKERRQ(MatGetNullSpace(jac->mat[1], &sp)); 908 if (sp) { 909 CHKERRQ(MatSetNullSpace(jac->schur, sp)); 910 } 911 912 CHKERRQ(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname)); 913 CHKERRQ(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg)); 914 if (flg) { 915 DM dmInner; 916 KSP kspInner; 917 PC pcInner; 918 919 CHKERRQ(MatSchurComplementGetKSP(jac->schur, &kspInner)); 920 CHKERRQ(KSPReset(kspInner)); 921 CHKERRQ(KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0])); 922 CHKERRQ(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 923 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 924 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2)); 925 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2)); 926 CHKERRQ(KSPSetOptionsPrefix(kspInner, schurprefix)); 927 928 /* Set DM for new solver */ 929 CHKERRQ(KSPGetDM(jac->head->ksp, &dmInner)); 930 CHKERRQ(KSPSetDM(kspInner, dmInner)); 931 CHKERRQ(KSPSetDMActive(kspInner, PETSC_FALSE)); 932 933 /* Defaults to PCKSP as preconditioner */ 934 CHKERRQ(KSPGetPC(kspInner, &pcInner)); 935 CHKERRQ(PCSetType(pcInner, PCKSP)); 936 CHKERRQ(PCKSPSetKSP(pcInner, jac->head->ksp)); 937 } else { 938 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 939 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 940 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 941 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 942 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 943 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 944 CHKERRQ(KSPSetType(jac->head->ksp,KSPGMRES)); 945 CHKERRQ(MatSchurComplementSetKSP(jac->schur,jac->head->ksp)); 946 } 947 CHKERRQ(KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0])); 948 CHKERRQ(KSPSetFromOptions(jac->head->ksp)); 949 CHKERRQ(MatSetFromOptions(jac->schur)); 950 951 CHKERRQ(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg)); 952 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */ 953 KSP kspInner; 954 PC pcInner; 955 956 CHKERRQ(MatSchurComplementGetKSP(jac->schur, &kspInner)); 957 CHKERRQ(KSPGetPC(kspInner, &pcInner)); 958 CHKERRQ(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg)); 959 if (flg) { 960 KSP ksp; 961 962 CHKERRQ(PCKSPGetKSP(pcInner, &ksp)); 963 if (ksp == jac->head->ksp) { 964 CHKERRQ(PCSetUseAmat(pcInner, PETSC_TRUE)); 965 } 966 } 967 } 968 CHKERRQ(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname)); 969 CHKERRQ(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg)); 970 if (flg) { 971 DM dmInner; 972 973 CHKERRQ(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 974 CHKERRQ(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper)); 975 CHKERRQ(KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure)); 976 CHKERRQ(KSPSetOptionsPrefix(jac->kspupper, schurprefix)); 977 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1)); 978 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1)); 979 CHKERRQ(KSPGetDM(jac->head->ksp, &dmInner)); 980 CHKERRQ(KSPSetDM(jac->kspupper, dmInner)); 981 CHKERRQ(KSPSetDMActive(jac->kspupper, PETSC_FALSE)); 982 CHKERRQ(KSPSetFromOptions(jac->kspupper)); 983 CHKERRQ(KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0])); 984 CHKERRQ(VecDuplicate(jac->head->x, &jac->head->z)); 985 } else { 986 jac->kspupper = jac->head->ksp; 987 CHKERRQ(PetscObjectReference((PetscObject) jac->head->ksp)); 988 } 989 990 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 991 CHKERRQ(MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp)); 992 } 993 CHKERRQ(KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur)); 994 CHKERRQ(KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure)); 995 CHKERRQ(PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur)); 996 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1)); 997 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 998 PC pcschur; 999 CHKERRQ(KSPGetPC(jac->kspschur,&pcschur)); 1000 CHKERRQ(PCSetType(pcschur,PCNONE)); 1001 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 1002 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 1003 CHKERRQ(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 1004 } 1005 CHKERRQ(KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac))); 1006 CHKERRQ(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix)); 1007 CHKERRQ(KSPSetOptionsPrefix(jac->kspschur, Dprefix)); 1008 /* propagate DM */ 1009 { 1010 DM sdm; 1011 CHKERRQ(KSPGetDM(jac->head->next->ksp, &sdm)); 1012 if (sdm) { 1013 CHKERRQ(KSPSetDM(jac->kspschur, sdm)); 1014 CHKERRQ(KSPSetDMActive(jac->kspschur, PETSC_FALSE)); 1015 } 1016 } 1017 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1018 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 1019 CHKERRQ(KSPSetFromOptions(jac->kspschur)); 1020 } 1021 CHKERRQ(MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY)); 1022 CHKERRQ(MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY)); 1023 1024 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 1025 CHKERRQ(PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname)); 1026 CHKERRQ(PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L)); 1027 if (!LSC_L) CHKERRQ(PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L)); 1028 if (LSC_L) CHKERRQ(PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L)); 1029 CHKERRQ(PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname)); 1030 CHKERRQ(PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L)); 1031 if (!LSC_L) CHKERRQ(PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L)); 1032 if (LSC_L) CHKERRQ(PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L)); 1033 } else if (jac->type == PC_COMPOSITE_GKB) { 1034 IS ccis; 1035 PetscInt rstart,rend; 1036 1037 PetscCheckFalse(nsplit != 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields"); 1038 1039 ilink = jac->head; 1040 1041 /* When extracting off-diagonal submatrices, we take complements from this range */ 1042 CHKERRQ(MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend)); 1043 1044 CHKERRQ(ISComplement(ilink->is_col,rstart,rend,&ccis)); 1045 if (jac->offdiag_use_amat) { 1046 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 1047 } else { 1048 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 1049 } 1050 CHKERRQ(ISDestroy(&ccis)); 1051 /* Create work vectors for GKB algorithm */ 1052 CHKERRQ(VecDuplicate(ilink->x,&jac->u)); 1053 CHKERRQ(VecDuplicate(ilink->x,&jac->Hu)); 1054 CHKERRQ(VecDuplicate(ilink->x,&jac->w2)); 1055 ilink = ilink->next; 1056 CHKERRQ(ISComplement(ilink->is_col,rstart,rend,&ccis)); 1057 if (jac->offdiag_use_amat) { 1058 CHKERRQ(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 1059 } else { 1060 CHKERRQ(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 1061 } 1062 CHKERRQ(ISDestroy(&ccis)); 1063 /* Create work vectors for GKB algorithm */ 1064 CHKERRQ(VecDuplicate(ilink->x,&jac->v)); 1065 CHKERRQ(VecDuplicate(ilink->x,&jac->d)); 1066 CHKERRQ(VecDuplicate(ilink->x,&jac->w1)); 1067 CHKERRQ(MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu)); 1068 CHKERRQ(PetscCalloc1(jac->gkbdelay,&jac->vecz)); 1069 1070 ilink = jac->head; 1071 CHKERRQ(KSPSetOperators(ilink->ksp,jac->H,jac->H)); 1072 if (!jac->suboptionsset) CHKERRQ(KSPSetFromOptions(ilink->ksp)); 1073 /* Create gkb_monitor context */ 1074 if (jac->gkbmonitor) { 1075 PetscInt tablevel; 1076 CHKERRQ(PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer)); 1077 CHKERRQ(PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII)); 1078 CHKERRQ(PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel)); 1079 CHKERRQ(PetscViewerASCIISetTab(jac->gkbviewer,tablevel)); 1080 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1)); 1081 } 1082 } else { 1083 /* set up the individual splits' PCs */ 1084 i = 0; 1085 ilink = jac->head; 1086 while (ilink) { 1087 CHKERRQ(KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i])); 1088 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1089 if (!jac->suboptionsset) CHKERRQ(KSPSetFromOptions(ilink->ksp)); 1090 i++; 1091 ilink = ilink->next; 1092 } 1093 } 1094 1095 /* Set coordinates to the sub PC objects whenever these are set */ 1096 if (jac->coordinates_set) { 1097 PC pc_coords; 1098 if (jac->type == PC_COMPOSITE_SCHUR) { 1099 // Head is first block. 1100 CHKERRQ(KSPGetPC(jac->head->ksp, &pc_coords)); 1101 CHKERRQ(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords)); 1102 // Second one is Schur block, but its KSP object is in kspschur. 1103 CHKERRQ(KSPGetPC(jac->kspschur, &pc_coords)); 1104 CHKERRQ(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords)); 1105 } else if (jac->type == PC_COMPOSITE_GKB) { 1106 CHKERRQ(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner")); 1107 } else { 1108 ilink = jac->head; 1109 while (ilink) { 1110 CHKERRQ(KSPGetPC(ilink->ksp, &pc_coords)); 1111 CHKERRQ(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords)); 1112 ilink = ilink->next; 1113 } 1114 } 1115 } 1116 1117 jac->suboptionsset = PETSC_TRUE; 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 1122 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 1123 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 1124 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 1125 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 1126 KSPCheckSolve(ilink->ksp,pc,ilink->y) || \ 1127 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 1128 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1129 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 1130 1131 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 1132 { 1133 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1134 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1135 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1136 1137 PetscFunctionBegin; 1138 switch (jac->schurfactorization) { 1139 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1140 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1141 CHKERRQ(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1142 CHKERRQ(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1143 CHKERRQ(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1144 CHKERRQ(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1145 CHKERRQ(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1146 CHKERRQ(KSPCheckSolve(kspA,pc,ilinkA->y)); 1147 CHKERRQ(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1148 CHKERRQ(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1149 CHKERRQ(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1150 CHKERRQ(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1151 CHKERRQ(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1152 CHKERRQ(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1153 CHKERRQ(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1154 CHKERRQ(VecScale(ilinkD->y,jac->schurscale)); 1155 CHKERRQ(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1156 CHKERRQ(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1157 CHKERRQ(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1158 break; 1159 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1160 /* [A00 0; A10 S], suitable for left preconditioning */ 1161 CHKERRQ(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1162 CHKERRQ(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1163 CHKERRQ(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1164 CHKERRQ(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1165 CHKERRQ(KSPCheckSolve(kspA,pc,ilinkA->y)); 1166 CHKERRQ(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1167 CHKERRQ(MatMult(jac->C,ilinkA->y,ilinkD->x)); 1168 CHKERRQ(VecScale(ilinkD->x,-1.)); 1169 CHKERRQ(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1170 CHKERRQ(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1171 CHKERRQ(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1172 CHKERRQ(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1173 CHKERRQ(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1174 CHKERRQ(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1175 CHKERRQ(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1176 CHKERRQ(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1177 CHKERRQ(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1178 CHKERRQ(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1179 break; 1180 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1181 /* [A00 A01; 0 S], suitable for right preconditioning */ 1182 CHKERRQ(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1183 CHKERRQ(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1184 CHKERRQ(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1185 CHKERRQ(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1186 CHKERRQ(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1187 CHKERRQ(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); CHKERRQ(MatMult(jac->B,ilinkD->y,ilinkA->x)); 1188 CHKERRQ(VecScale(ilinkA->x,-1.)); 1189 CHKERRQ(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD)); 1190 CHKERRQ(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1191 CHKERRQ(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD)); 1192 CHKERRQ(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1193 CHKERRQ(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1194 CHKERRQ(KSPCheckSolve(kspA,pc,ilinkA->y)); 1195 CHKERRQ(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1196 CHKERRQ(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1197 CHKERRQ(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1198 CHKERRQ(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1199 break; 1200 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1201 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */ 1202 CHKERRQ(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1203 CHKERRQ(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1204 CHKERRQ(PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL)); 1205 CHKERRQ(KSPSolve(kspLower,ilinkA->x,ilinkA->y)); 1206 CHKERRQ(KSPCheckSolve(kspLower,pc,ilinkA->y)); 1207 CHKERRQ(PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL)); 1208 CHKERRQ(MatMult(jac->C,ilinkA->y,ilinkD->x)); 1209 CHKERRQ(VecScale(ilinkD->x,-1.0)); 1210 CHKERRQ(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1211 CHKERRQ(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1212 1213 CHKERRQ(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1214 CHKERRQ(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1215 CHKERRQ(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1216 CHKERRQ(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1217 CHKERRQ(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1218 1219 if (kspUpper == kspA) { 1220 CHKERRQ(MatMult(jac->B,ilinkD->y,ilinkA->y)); 1221 CHKERRQ(VecAXPY(ilinkA->x,-1.0,ilinkA->y)); 1222 CHKERRQ(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1223 CHKERRQ(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1224 CHKERRQ(KSPCheckSolve(kspA,pc,ilinkA->y)); 1225 CHKERRQ(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1226 } else { 1227 CHKERRQ(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1228 CHKERRQ(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1229 CHKERRQ(KSPCheckSolve(kspA,pc,ilinkA->y)); 1230 CHKERRQ(MatMult(jac->B,ilinkD->y,ilinkA->x)); 1231 CHKERRQ(PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL)); 1232 CHKERRQ(KSPSolve(kspUpper,ilinkA->x,ilinkA->z)); 1233 CHKERRQ(KSPCheckSolve(kspUpper,pc,ilinkA->z)); 1234 CHKERRQ(PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL)); 1235 CHKERRQ(VecAXPY(ilinkA->y,-1.0,ilinkA->z)); 1236 } 1237 CHKERRQ(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1238 CHKERRQ(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1239 CHKERRQ(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1240 } 1241 PetscFunctionReturn(0); 1242 } 1243 1244 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 1245 { 1246 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1247 PC_FieldSplitLink ilink = jac->head; 1248 PetscInt cnt,bs; 1249 1250 PetscFunctionBegin; 1251 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1252 if (jac->defaultsplit) { 1253 CHKERRQ(VecGetBlockSize(x,&bs)); 1254 PetscCheckFalse(jac->bs > 0 && bs != jac->bs,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1255 CHKERRQ(VecGetBlockSize(y,&bs)); 1256 PetscCheckFalse(jac->bs > 0 && bs != jac->bs,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1257 CHKERRQ(VecStrideGatherAll(x,jac->x,INSERT_VALUES)); 1258 while (ilink) { 1259 CHKERRQ(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1260 CHKERRQ(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1261 CHKERRQ(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1262 CHKERRQ(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1263 ilink = ilink->next; 1264 } 1265 CHKERRQ(VecStrideScatterAll(jac->y,y,INSERT_VALUES)); 1266 } else { 1267 CHKERRQ(VecSet(y,0.0)); 1268 while (ilink) { 1269 CHKERRQ(FieldSplitSplitSolveAdd(ilink,x,y)); 1270 ilink = ilink->next; 1271 } 1272 } 1273 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1274 CHKERRQ(VecSet(y,0.0)); 1275 /* solve on first block for first block variables */ 1276 CHKERRQ(VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD)); 1277 CHKERRQ(VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD)); 1278 CHKERRQ(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1279 CHKERRQ(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1280 CHKERRQ(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1281 CHKERRQ(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1282 CHKERRQ(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1283 CHKERRQ(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1284 1285 /* compute the residual only onto second block variables using first block variables */ 1286 CHKERRQ(MatMult(jac->Afield[1],ilink->y,ilink->next->x)); 1287 ilink = ilink->next; 1288 CHKERRQ(VecScale(ilink->x,-1.0)); 1289 CHKERRQ(VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1290 CHKERRQ(VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1291 1292 /* solve on second block variables */ 1293 CHKERRQ(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1294 CHKERRQ(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1295 CHKERRQ(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1296 CHKERRQ(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1297 CHKERRQ(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1298 CHKERRQ(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1299 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1300 if (!jac->w1) { 1301 CHKERRQ(VecDuplicate(x,&jac->w1)); 1302 CHKERRQ(VecDuplicate(x,&jac->w2)); 1303 } 1304 CHKERRQ(VecSet(y,0.0)); 1305 CHKERRQ(FieldSplitSplitSolveAdd(ilink,x,y)); 1306 cnt = 1; 1307 while (ilink->next) { 1308 ilink = ilink->next; 1309 /* compute the residual only over the part of the vector needed */ 1310 CHKERRQ(MatMult(jac->Afield[cnt++],y,ilink->x)); 1311 CHKERRQ(VecScale(ilink->x,-1.0)); 1312 CHKERRQ(VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1313 CHKERRQ(VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1314 CHKERRQ(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1315 CHKERRQ(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1316 CHKERRQ(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1317 CHKERRQ(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1318 CHKERRQ(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1319 CHKERRQ(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1320 } 1321 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1322 cnt -= 2; 1323 while (ilink->previous) { 1324 ilink = ilink->previous; 1325 /* compute the residual only over the part of the vector needed */ 1326 CHKERRQ(MatMult(jac->Afield[cnt--],y,ilink->x)); 1327 CHKERRQ(VecScale(ilink->x,-1.0)); 1328 CHKERRQ(VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1329 CHKERRQ(VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1330 CHKERRQ(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1331 CHKERRQ(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1332 CHKERRQ(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1333 CHKERRQ(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1334 CHKERRQ(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1335 CHKERRQ(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1336 } 1337 } 1338 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1339 PetscFunctionReturn(0); 1340 } 1341 1342 static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y) 1343 { 1344 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1345 PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next; 1346 KSP ksp = ilinkA->ksp; 1347 Vec u,v,Hu,d,work1,work2; 1348 PetscScalar alpha,z,nrmz2,*vecz; 1349 PetscReal lowbnd,nu,beta; 1350 PetscInt j,iterGKB; 1351 1352 PetscFunctionBegin; 1353 CHKERRQ(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1354 CHKERRQ(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1355 CHKERRQ(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1356 CHKERRQ(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1357 1358 u = jac->u; 1359 v = jac->v; 1360 Hu = jac->Hu; 1361 d = jac->d; 1362 work1 = jac->w1; 1363 work2 = jac->w2; 1364 vecz = jac->vecz; 1365 1366 /* Change RHS to comply with matrix regularization H = A + nu*B*B' */ 1367 /* Add q = q + nu*B*b */ 1368 if (jac->gkbnu) { 1369 nu = jac->gkbnu; 1370 CHKERRQ(VecScale(ilinkD->x,jac->gkbnu)); 1371 CHKERRQ(MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x)); /* q = q + nu*B*b */ 1372 } else { 1373 /* Situation when no augmented Lagrangian is used. Then we set inner */ 1374 /* matrix N = I in [Ar13], and thus nu = 1. */ 1375 nu = 1; 1376 } 1377 1378 /* Transform rhs from [q,tilde{b}] to [0,b] */ 1379 CHKERRQ(PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL)); 1380 CHKERRQ(KSPSolve(ksp,ilinkA->x,ilinkA->y)); 1381 CHKERRQ(KSPCheckSolve(ksp,pc,ilinkA->y)); 1382 CHKERRQ(PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL)); 1383 CHKERRQ(MatMultHermitianTranspose(jac->B,ilinkA->y,work1)); 1384 CHKERRQ(VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x)); /* c = b - B'*x */ 1385 1386 /* First step of algorithm */ 1387 CHKERRQ(VecNorm(work1,NORM_2,&beta)); /* beta = sqrt(nu*c'*c)*/ 1388 KSPCheckDot(ksp,beta); 1389 beta = PetscSqrtReal(nu)*beta; 1390 CHKERRQ(VecAXPBY(v,nu/beta,0.0,work1)); /* v = nu/beta *c */ 1391 CHKERRQ(MatMult(jac->B,v,work2)); /* u = H^{-1}*B*v */ 1392 CHKERRQ(PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL)); 1393 CHKERRQ(KSPSolve(ksp,work2,u)); 1394 CHKERRQ(KSPCheckSolve(ksp,pc,u)); 1395 CHKERRQ(PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL)); 1396 CHKERRQ(MatMult(jac->H,u,Hu)); /* alpha = u'*H*u */ 1397 CHKERRQ(VecDot(Hu,u,&alpha)); 1398 KSPCheckDot(ksp,alpha); 1399 PetscCheckFalse(PetscRealPart(alpha) <= 0.0,PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite"); 1400 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1401 CHKERRQ(VecScale(u,1.0/alpha)); 1402 CHKERRQ(VecAXPBY(d,1.0/alpha,0.0,v)); /* v = nu/beta *c */ 1403 1404 z = beta/alpha; 1405 vecz[1] = z; 1406 1407 /* Computation of first iterate x(1) and p(1) */ 1408 CHKERRQ(VecAXPY(ilinkA->y,z,u)); 1409 CHKERRQ(VecCopy(d,ilinkD->y)); 1410 CHKERRQ(VecScale(ilinkD->y,-z)); 1411 1412 iterGKB = 1; lowbnd = 2*jac->gkbtol; 1413 if (jac->gkbmonitor) { 1414 CHKERRQ(PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd)); 1415 } 1416 1417 while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) { 1418 iterGKB += 1; 1419 CHKERRQ(MatMultHermitianTranspose(jac->B,u,work1)); /* v <- nu*(B'*u-alpha/nu*v) */ 1420 CHKERRQ(VecAXPBY(v,nu,-alpha,work1)); 1421 CHKERRQ(VecNorm(v,NORM_2,&beta)); /* beta = sqrt(nu)*v'*v */ 1422 beta = beta/PetscSqrtReal(nu); 1423 CHKERRQ(VecScale(v,1.0/beta)); 1424 CHKERRQ(MatMult(jac->B,v,work2)); /* u <- H^{-1}*(B*v-beta*H*u) */ 1425 CHKERRQ(MatMult(jac->H,u,Hu)); 1426 CHKERRQ(VecAXPY(work2,-beta,Hu)); 1427 CHKERRQ(PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL)); 1428 CHKERRQ(KSPSolve(ksp,work2,u)); 1429 CHKERRQ(KSPCheckSolve(ksp,pc,u)); 1430 CHKERRQ(PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL)); 1431 CHKERRQ(MatMult(jac->H,u,Hu)); /* alpha = u'*H*u */ 1432 CHKERRQ(VecDot(Hu,u,&alpha)); 1433 KSPCheckDot(ksp,alpha); 1434 PetscCheckFalse(PetscRealPart(alpha) <= 0.0,PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite"); 1435 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1436 CHKERRQ(VecScale(u,1.0/alpha)); 1437 1438 z = -beta/alpha*z; /* z <- beta/alpha*z */ 1439 vecz[0] = z; 1440 1441 /* Computation of new iterate x(i+1) and p(i+1) */ 1442 CHKERRQ(VecAXPBY(d,1.0/alpha,-beta/alpha,v)); /* d = (v-beta*d)/alpha */ 1443 CHKERRQ(VecAXPY(ilinkA->y,z,u)); /* r = r + z*u */ 1444 CHKERRQ(VecAXPY(ilinkD->y,-z,d)); /* p = p - z*d */ 1445 CHKERRQ(MatMult(jac->H,ilinkA->y,Hu)); /* ||u||_H = u'*H*u */ 1446 CHKERRQ(VecDot(Hu,ilinkA->y,&nrmz2)); 1447 1448 /* Compute Lower Bound estimate */ 1449 if (iterGKB > jac->gkbdelay) { 1450 lowbnd = 0.0; 1451 for (j=0; j<jac->gkbdelay; j++) { 1452 lowbnd += PetscAbsScalar(vecz[j]*vecz[j]); 1453 } 1454 lowbnd = PetscSqrtReal(lowbnd/PetscAbsScalar(nrmz2)); 1455 } 1456 1457 for (j=0; j<jac->gkbdelay-1; j++) { 1458 vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2]; 1459 } 1460 if (jac->gkbmonitor) { 1461 CHKERRQ(PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd)); 1462 } 1463 } 1464 1465 CHKERRQ(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1466 CHKERRQ(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1467 CHKERRQ(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1468 CHKERRQ(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1469 1470 PetscFunctionReturn(0); 1471 } 1472 1473 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1474 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1475 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1476 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1477 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1478 KSPCheckSolve(ilink->ksp,pc,ilink->x) || \ 1479 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1480 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1481 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1482 1483 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1484 { 1485 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1486 PC_FieldSplitLink ilink = jac->head; 1487 PetscInt bs; 1488 1489 PetscFunctionBegin; 1490 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1491 if (jac->defaultsplit) { 1492 CHKERRQ(VecGetBlockSize(x,&bs)); 1493 PetscCheckFalse(jac->bs > 0 && bs != jac->bs,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1494 CHKERRQ(VecGetBlockSize(y,&bs)); 1495 PetscCheckFalse(jac->bs > 0 && bs != jac->bs,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1496 CHKERRQ(VecStrideGatherAll(x,jac->x,INSERT_VALUES)); 1497 while (ilink) { 1498 CHKERRQ(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1499 CHKERRQ(KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y)); 1500 CHKERRQ(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1501 CHKERRQ(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1502 ilink = ilink->next; 1503 } 1504 CHKERRQ(VecStrideScatterAll(jac->y,y,INSERT_VALUES)); 1505 } else { 1506 CHKERRQ(VecSet(y,0.0)); 1507 while (ilink) { 1508 CHKERRQ(FieldSplitSplitSolveAddTranspose(ilink,x,y)); 1509 ilink = ilink->next; 1510 } 1511 } 1512 } else { 1513 if (!jac->w1) { 1514 CHKERRQ(VecDuplicate(x,&jac->w1)); 1515 CHKERRQ(VecDuplicate(x,&jac->w2)); 1516 } 1517 CHKERRQ(VecSet(y,0.0)); 1518 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1519 CHKERRQ(FieldSplitSplitSolveAddTranspose(ilink,x,y)); 1520 while (ilink->next) { 1521 ilink = ilink->next; 1522 CHKERRQ(MatMultTranspose(pc->mat,y,jac->w1)); 1523 CHKERRQ(VecWAXPY(jac->w2,-1.0,jac->w1,x)); 1524 CHKERRQ(FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y)); 1525 } 1526 while (ilink->previous) { 1527 ilink = ilink->previous; 1528 CHKERRQ(MatMultTranspose(pc->mat,y,jac->w1)); 1529 CHKERRQ(VecWAXPY(jac->w2,-1.0,jac->w1,x)); 1530 CHKERRQ(FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y)); 1531 } 1532 } else { 1533 while (ilink->next) { /* get to last entry in linked list */ 1534 ilink = ilink->next; 1535 } 1536 CHKERRQ(FieldSplitSplitSolveAddTranspose(ilink,x,y)); 1537 while (ilink->previous) { 1538 ilink = ilink->previous; 1539 CHKERRQ(MatMultTranspose(pc->mat,y,jac->w1)); 1540 CHKERRQ(VecWAXPY(jac->w2,-1.0,jac->w1,x)); 1541 CHKERRQ(FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y)); 1542 } 1543 } 1544 } 1545 PetscFunctionReturn(0); 1546 } 1547 1548 static PetscErrorCode PCReset_FieldSplit(PC pc) 1549 { 1550 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1551 PC_FieldSplitLink ilink = jac->head,next; 1552 1553 PetscFunctionBegin; 1554 while (ilink) { 1555 CHKERRQ(KSPDestroy(&ilink->ksp)); 1556 CHKERRQ(VecDestroy(&ilink->x)); 1557 CHKERRQ(VecDestroy(&ilink->y)); 1558 CHKERRQ(VecDestroy(&ilink->z)); 1559 CHKERRQ(VecScatterDestroy(&ilink->sctx)); 1560 CHKERRQ(ISDestroy(&ilink->is)); 1561 CHKERRQ(ISDestroy(&ilink->is_col)); 1562 CHKERRQ(PetscFree(ilink->splitname)); 1563 CHKERRQ(PetscFree(ilink->fields)); 1564 CHKERRQ(PetscFree(ilink->fields_col)); 1565 next = ilink->next; 1566 CHKERRQ(PetscFree(ilink)); 1567 ilink = next; 1568 } 1569 jac->head = NULL; 1570 CHKERRQ(PetscFree2(jac->x,jac->y)); 1571 if (jac->mat && jac->mat != jac->pmat) { 1572 CHKERRQ(MatDestroyMatrices(jac->nsplits,&jac->mat)); 1573 } else if (jac->mat) { 1574 jac->mat = NULL; 1575 } 1576 if (jac->pmat) CHKERRQ(MatDestroyMatrices(jac->nsplits,&jac->pmat)); 1577 if (jac->Afield) CHKERRQ(MatDestroyMatrices(jac->nsplits,&jac->Afield)); 1578 jac->nsplits = 0; 1579 CHKERRQ(VecDestroy(&jac->w1)); 1580 CHKERRQ(VecDestroy(&jac->w2)); 1581 CHKERRQ(MatDestroy(&jac->schur)); 1582 CHKERRQ(MatDestroy(&jac->schurp)); 1583 CHKERRQ(MatDestroy(&jac->schur_user)); 1584 CHKERRQ(KSPDestroy(&jac->kspschur)); 1585 CHKERRQ(KSPDestroy(&jac->kspupper)); 1586 CHKERRQ(MatDestroy(&jac->B)); 1587 CHKERRQ(MatDestroy(&jac->C)); 1588 CHKERRQ(MatDestroy(&jac->H)); 1589 CHKERRQ(VecDestroy(&jac->u)); 1590 CHKERRQ(VecDestroy(&jac->v)); 1591 CHKERRQ(VecDestroy(&jac->Hu)); 1592 CHKERRQ(VecDestroy(&jac->d)); 1593 CHKERRQ(PetscFree(jac->vecz)); 1594 CHKERRQ(PetscViewerDestroy(&jac->gkbviewer)); 1595 jac->isrestrict = PETSC_FALSE; 1596 PetscFunctionReturn(0); 1597 } 1598 1599 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1600 { 1601 PetscFunctionBegin; 1602 CHKERRQ(PCReset_FieldSplit(pc)); 1603 CHKERRQ(PetscFree(pc->data)); 1604 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL)); 1605 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 1606 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL)); 1607 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL)); 1608 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL)); 1609 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL)); 1610 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL)); 1611 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL)); 1612 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL)); 1613 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL)); 1614 PetscFunctionReturn(0); 1615 } 1616 1617 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1618 { 1619 PetscInt bs; 1620 PetscBool flg; 1621 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1622 PCCompositeType ctype; 1623 1624 PetscFunctionBegin; 1625 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"FieldSplit options")); 1626 CHKERRQ(PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL)); 1627 CHKERRQ(PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg)); 1628 if (flg) { 1629 CHKERRQ(PCFieldSplitSetBlockSize(pc,bs)); 1630 } 1631 jac->diag_use_amat = pc->useAmat; 1632 CHKERRQ(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)); 1633 jac->offdiag_use_amat = pc->useAmat; 1634 CHKERRQ(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)); 1635 CHKERRQ(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL)); 1636 CHKERRQ(PCFieldSplitSetDetectSaddlePoint(pc,jac->detect)); /* Sets split type and Schur PC type */ 1637 CHKERRQ(PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg)); 1638 if (flg) { 1639 CHKERRQ(PCFieldSplitSetType(pc,ctype)); 1640 } 1641 /* Only setup fields once */ 1642 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1643 /* only allow user to set fields from command line if bs is already known. 1644 otherwise user can set them in PCFieldSplitSetDefaults() */ 1645 CHKERRQ(PCFieldSplitSetRuntimeSplits_Private(pc)); 1646 if (jac->splitdefined) CHKERRQ(PetscInfo(pc,"Splits defined using the options database\n")); 1647 } 1648 if (jac->type == PC_COMPOSITE_SCHUR) { 1649 CHKERRQ(PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg)); 1650 if (flg) CHKERRQ(PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n")); 1651 CHKERRQ(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)); 1652 CHKERRQ(PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL)); 1653 CHKERRQ(PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL)); 1654 } else if (jac->type == PC_COMPOSITE_GKB) { 1655 CHKERRQ(PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL)); 1656 CHKERRQ(PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL)); 1657 CHKERRQ(PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL)); 1658 PetscCheckFalse(jac->gkbnu < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu); 1659 CHKERRQ(PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL)); 1660 CHKERRQ(PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL)); 1661 } 1662 CHKERRQ(PetscOptionsTail()); 1663 PetscFunctionReturn(0); 1664 } 1665 1666 /*------------------------------------------------------------------------------------*/ 1667 1668 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1669 { 1670 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1671 PC_FieldSplitLink ilink,next = jac->head; 1672 char prefix[128]; 1673 PetscInt i; 1674 1675 PetscFunctionBegin; 1676 if (jac->splitdefined) { 1677 CHKERRQ(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1678 PetscFunctionReturn(0); 1679 } 1680 for (i=0; i<n; i++) { 1681 PetscCheckFalse(fields[i] >= jac->bs,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1682 PetscCheckFalse(fields[i] < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1683 } 1684 CHKERRQ(PetscNew(&ilink)); 1685 if (splitname) { 1686 CHKERRQ(PetscStrallocpy(splitname,&ilink->splitname)); 1687 } else { 1688 CHKERRQ(PetscMalloc1(3,&ilink->splitname)); 1689 CHKERRQ(PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits)); 1690 } 1691 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1692 CHKERRQ(PetscMalloc1(n,&ilink->fields)); 1693 CHKERRQ(PetscArraycpy(ilink->fields,fields,n)); 1694 CHKERRQ(PetscMalloc1(n,&ilink->fields_col)); 1695 CHKERRQ(PetscArraycpy(ilink->fields_col,fields_col,n)); 1696 1697 ilink->nfields = n; 1698 ilink->next = NULL; 1699 CHKERRQ(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1700 CHKERRQ(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1701 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1702 CHKERRQ(KSPSetType(ilink->ksp,KSPPREONLY)); 1703 CHKERRQ(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1704 1705 CHKERRQ(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1706 CHKERRQ(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1707 1708 if (!next) { 1709 jac->head = ilink; 1710 ilink->previous = NULL; 1711 } else { 1712 while (next->next) { 1713 next = next->next; 1714 } 1715 next->next = ilink; 1716 ilink->previous = next; 1717 } 1718 jac->nsplits++; 1719 PetscFunctionReturn(0); 1720 } 1721 1722 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1723 { 1724 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1725 1726 PetscFunctionBegin; 1727 *subksp = NULL; 1728 if (n) *n = 0; 1729 if (jac->type == PC_COMPOSITE_SCHUR) { 1730 PetscInt nn; 1731 1732 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1733 PetscCheckFalse(jac->nsplits != 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits); 1734 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1735 CHKERRQ(PetscMalloc1(nn,subksp)); 1736 (*subksp)[0] = jac->head->ksp; 1737 (*subksp)[1] = jac->kspschur; 1738 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1739 if (n) *n = nn; 1740 } 1741 PetscFunctionReturn(0); 1742 } 1743 1744 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1745 { 1746 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1747 1748 PetscFunctionBegin; 1749 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1750 CHKERRQ(PetscMalloc1(jac->nsplits,subksp)); 1751 CHKERRQ(MatSchurComplementGetKSP(jac->schur,*subksp)); 1752 1753 (*subksp)[1] = jac->kspschur; 1754 if (n) *n = jac->nsplits; 1755 PetscFunctionReturn(0); 1756 } 1757 1758 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1759 { 1760 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1761 PetscInt cnt = 0; 1762 PC_FieldSplitLink ilink = jac->head; 1763 1764 PetscFunctionBegin; 1765 CHKERRQ(PetscMalloc1(jac->nsplits,subksp)); 1766 while (ilink) { 1767 (*subksp)[cnt++] = ilink->ksp; 1768 ilink = ilink->next; 1769 } 1770 PetscCheckFalse(cnt != jac->nsplits,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); 1771 if (n) *n = jac->nsplits; 1772 PetscFunctionReturn(0); 1773 } 1774 1775 /*@C 1776 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1777 1778 Input Parameters: 1779 + pc - the preconditioner context 1780 - is - the index set that defines the indices to which the fieldsplit is to be restricted 1781 1782 Level: advanced 1783 1784 @*/ 1785 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1786 { 1787 PetscFunctionBegin; 1788 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1789 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1790 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy))); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1795 { 1796 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1797 PC_FieldSplitLink ilink = jac->head, next; 1798 PetscInt localsize,size,sizez,i; 1799 const PetscInt *ind, *indz; 1800 PetscInt *indc, *indcz; 1801 PetscBool flg; 1802 1803 PetscFunctionBegin; 1804 CHKERRQ(ISGetLocalSize(isy,&localsize)); 1805 CHKERRMPI(MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy))); 1806 size -= localsize; 1807 while (ilink) { 1808 IS isrl,isr; 1809 PC subpc; 1810 CHKERRQ(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl)); 1811 CHKERRQ(ISGetLocalSize(isrl,&localsize)); 1812 CHKERRQ(PetscMalloc1(localsize,&indc)); 1813 CHKERRQ(ISGetIndices(isrl,&ind)); 1814 CHKERRQ(PetscArraycpy(indc,ind,localsize)); 1815 CHKERRQ(ISRestoreIndices(isrl,&ind)); 1816 CHKERRQ(ISDestroy(&isrl)); 1817 for (i=0; i<localsize; i++) *(indc+i) += size; 1818 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr)); 1819 CHKERRQ(PetscObjectReference((PetscObject)isr)); 1820 CHKERRQ(ISDestroy(&ilink->is)); 1821 ilink->is = isr; 1822 CHKERRQ(PetscObjectReference((PetscObject)isr)); 1823 CHKERRQ(ISDestroy(&ilink->is_col)); 1824 ilink->is_col = isr; 1825 CHKERRQ(ISDestroy(&isr)); 1826 CHKERRQ(KSPGetPC(ilink->ksp, &subpc)); 1827 CHKERRQ(PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg)); 1828 if (flg) { 1829 IS iszl,isz; 1830 MPI_Comm comm; 1831 CHKERRQ(ISGetLocalSize(ilink->is,&localsize)); 1832 comm = PetscObjectComm((PetscObject)ilink->is); 1833 CHKERRQ(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl)); 1834 CHKERRMPI(MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm)); 1835 sizez -= localsize; 1836 CHKERRQ(ISGetLocalSize(iszl,&localsize)); 1837 CHKERRQ(PetscMalloc1(localsize,&indcz)); 1838 CHKERRQ(ISGetIndices(iszl,&indz)); 1839 CHKERRQ(PetscArraycpy(indcz,indz,localsize)); 1840 CHKERRQ(ISRestoreIndices(iszl,&indz)); 1841 CHKERRQ(ISDestroy(&iszl)); 1842 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1843 CHKERRQ(ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz)); 1844 CHKERRQ(PCFieldSplitRestrictIS(subpc,isz)); 1845 CHKERRQ(ISDestroy(&isz)); 1846 } 1847 next = ilink->next; 1848 ilink = next; 1849 } 1850 jac->isrestrict = PETSC_TRUE; 1851 PetscFunctionReturn(0); 1852 } 1853 1854 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1855 { 1856 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1857 PC_FieldSplitLink ilink, next = jac->head; 1858 char prefix[128]; 1859 1860 PetscFunctionBegin; 1861 if (jac->splitdefined) { 1862 CHKERRQ(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1863 PetscFunctionReturn(0); 1864 } 1865 CHKERRQ(PetscNew(&ilink)); 1866 if (splitname) { 1867 CHKERRQ(PetscStrallocpy(splitname,&ilink->splitname)); 1868 } else { 1869 CHKERRQ(PetscMalloc1(8,&ilink->splitname)); 1870 CHKERRQ(PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits)); 1871 } 1872 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1873 CHKERRQ(PetscObjectReference((PetscObject)is)); 1874 CHKERRQ(ISDestroy(&ilink->is)); 1875 ilink->is = is; 1876 CHKERRQ(PetscObjectReference((PetscObject)is)); 1877 CHKERRQ(ISDestroy(&ilink->is_col)); 1878 ilink->is_col = is; 1879 ilink->next = NULL; 1880 CHKERRQ(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1881 CHKERRQ(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1882 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1883 CHKERRQ(KSPSetType(ilink->ksp,KSPPREONLY)); 1884 CHKERRQ(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1885 1886 CHKERRQ(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1887 CHKERRQ(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1888 1889 if (!next) { 1890 jac->head = ilink; 1891 ilink->previous = NULL; 1892 } else { 1893 while (next->next) { 1894 next = next->next; 1895 } 1896 next->next = ilink; 1897 ilink->previous = next; 1898 } 1899 jac->nsplits++; 1900 PetscFunctionReturn(0); 1901 } 1902 1903 /*@C 1904 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1905 1906 Logically Collective on PC 1907 1908 Input Parameters: 1909 + pc - the preconditioner context 1910 . splitname - name of this split, if NULL the number of the split is used 1911 . n - the number of fields in this split 1912 - fields - the fields in this split 1913 1914 Level: intermediate 1915 1916 Notes: 1917 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1918 1919 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1920 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 1921 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1922 where the numbered entries indicate what is in the field. 1923 1924 This function is called once per split (it creates a new split each time). Solve options 1925 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1926 1927 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1928 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1929 available when this routine is called. 1930 1931 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1932 1933 @*/ 1934 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1935 { 1936 PetscFunctionBegin; 1937 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1938 PetscValidCharPointer(splitname,2); 1939 PetscCheckFalse(n < 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1940 PetscValidIntPointer(fields,4); 1941 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col))); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 /*@ 1946 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1947 1948 Logically Collective on PC 1949 1950 Input Parameters: 1951 + pc - the preconditioner object 1952 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1953 1954 Options Database: 1955 . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks 1956 1957 Level: intermediate 1958 1959 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1960 1961 @*/ 1962 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1963 { 1964 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1965 PetscBool isfs; 1966 1967 PetscFunctionBegin; 1968 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1969 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 1970 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1971 jac->diag_use_amat = flg; 1972 PetscFunctionReturn(0); 1973 } 1974 1975 /*@ 1976 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1977 1978 Logically Collective on PC 1979 1980 Input Parameters: 1981 . pc - the preconditioner object 1982 1983 Output Parameters: 1984 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1985 1986 Level: intermediate 1987 1988 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1989 1990 @*/ 1991 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1992 { 1993 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1994 PetscBool isfs; 1995 1996 PetscFunctionBegin; 1997 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1998 PetscValidBoolPointer(flg,2); 1999 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2000 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2001 *flg = jac->diag_use_amat; 2002 PetscFunctionReturn(0); 2003 } 2004 2005 /*@ 2006 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2007 2008 Logically Collective on PC 2009 2010 Input Parameters: 2011 + pc - the preconditioner object 2012 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2013 2014 Options Database: 2015 . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks 2016 2017 Level: intermediate 2018 2019 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 2020 2021 @*/ 2022 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 2023 { 2024 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2025 PetscBool isfs; 2026 2027 PetscFunctionBegin; 2028 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2029 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2030 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2031 jac->offdiag_use_amat = flg; 2032 PetscFunctionReturn(0); 2033 } 2034 2035 /*@ 2036 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2037 2038 Logically Collective on PC 2039 2040 Input Parameters: 2041 . pc - the preconditioner object 2042 2043 Output Parameters: 2044 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2045 2046 Level: intermediate 2047 2048 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 2049 2050 @*/ 2051 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 2052 { 2053 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2054 PetscBool isfs; 2055 2056 PetscFunctionBegin; 2057 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2058 PetscValidBoolPointer(flg,2); 2059 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2060 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2061 *flg = jac->offdiag_use_amat; 2062 PetscFunctionReturn(0); 2063 } 2064 2065 /*@C 2066 PCFieldSplitSetIS - Sets the exact elements for field 2067 2068 Logically Collective on PC 2069 2070 Input Parameters: 2071 + pc - the preconditioner context 2072 . splitname - name of this split, if NULL the number of the split is used 2073 - is - the index set that defines the vector elements in this field 2074 2075 Notes: 2076 Use PCFieldSplitSetFields(), for fields defined by strided types. 2077 2078 This function is called once per split (it creates a new split each time). Solve options 2079 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2080 2081 Level: intermediate 2082 2083 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 2084 2085 @*/ 2086 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 2087 { 2088 PetscFunctionBegin; 2089 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2090 if (splitname) PetscValidCharPointer(splitname,2); 2091 PetscValidHeaderSpecific(is,IS_CLASSID,3); 2092 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is))); 2093 PetscFunctionReturn(0); 2094 } 2095 2096 /*@C 2097 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 2098 2099 Logically Collective on PC 2100 2101 Input Parameters: 2102 + pc - the preconditioner context 2103 - splitname - name of this split 2104 2105 Output Parameter: 2106 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 2107 2108 Level: intermediate 2109 2110 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 2111 2112 @*/ 2113 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 2114 { 2115 PetscFunctionBegin; 2116 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2117 PetscValidCharPointer(splitname,2); 2118 PetscValidPointer(is,3); 2119 { 2120 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2121 PC_FieldSplitLink ilink = jac->head; 2122 PetscBool found; 2123 2124 *is = NULL; 2125 while (ilink) { 2126 CHKERRQ(PetscStrcmp(ilink->splitname, splitname, &found)); 2127 if (found) { 2128 *is = ilink->is; 2129 break; 2130 } 2131 ilink = ilink->next; 2132 } 2133 } 2134 PetscFunctionReturn(0); 2135 } 2136 2137 /*@C 2138 PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS 2139 2140 Logically Collective on PC 2141 2142 Input Parameters: 2143 + pc - the preconditioner context 2144 - index - index of this split 2145 2146 Output Parameter: 2147 - is - the index set that defines the vector elements in this field 2148 2149 Level: intermediate 2150 2151 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS() 2152 2153 @*/ 2154 PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is) 2155 { 2156 PetscFunctionBegin; 2157 PetscCheckFalse(index < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index); 2158 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2159 PetscValidPointer(is,3); 2160 { 2161 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2162 PC_FieldSplitLink ilink = jac->head; 2163 PetscInt i = 0; 2164 PetscCheckFalse(index >= jac->nsplits,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",index,jac->nsplits); 2165 2166 while (i < index) { 2167 ilink = ilink->next; 2168 ++i; 2169 } 2170 CHKERRQ(PCFieldSplitGetIS(pc,ilink->splitname,is)); 2171 } 2172 PetscFunctionReturn(0); 2173 } 2174 2175 /*@ 2176 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2177 fieldsplit preconditioner. If not set the matrix block size is used. 2178 2179 Logically Collective on PC 2180 2181 Input Parameters: 2182 + pc - the preconditioner context 2183 - bs - the block size 2184 2185 Level: intermediate 2186 2187 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 2188 2189 @*/ 2190 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 2191 { 2192 PetscFunctionBegin; 2193 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2194 PetscValidLogicalCollectiveInt(pc,bs,2); 2195 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs))); 2196 PetscFunctionReturn(0); 2197 } 2198 2199 /*@C 2200 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 2201 2202 Collective on KSP 2203 2204 Input Parameter: 2205 . pc - the preconditioner context 2206 2207 Output Parameters: 2208 + n - the number of splits 2209 - subksp - the array of KSP contexts 2210 2211 Note: 2212 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2213 (not the KSP just the array that contains them). 2214 2215 You must call PCSetUp() before calling PCFieldSplitGetSubKSP(). 2216 2217 If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the 2218 Schur complement and the KSP object used to iterate over the Schur complement. 2219 To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP(). 2220 2221 If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the 2222 inner linear system defined by the matrix H in each loop. 2223 2224 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2225 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2226 KSP array must be. 2227 2228 Level: advanced 2229 2230 .seealso: PCFIELDSPLIT 2231 @*/ 2232 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2233 { 2234 PetscFunctionBegin; 2235 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2236 if (n) PetscValidIntPointer(n,2); 2237 CHKERRQ(PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp))); 2238 PetscFunctionReturn(0); 2239 } 2240 2241 /*@C 2242 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT 2243 2244 Collective on KSP 2245 2246 Input Parameter: 2247 . pc - the preconditioner context 2248 2249 Output Parameters: 2250 + n - the number of splits 2251 - subksp - the array of KSP contexts 2252 2253 Note: 2254 After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2255 (not the KSP just the array that contains them). 2256 2257 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP(). 2258 2259 If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order) 2260 - the KSP used for the (1,1) block 2261 - the KSP used for the Schur complement (not the one used for the interior Schur solver) 2262 - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2263 2264 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP(). 2265 2266 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2267 You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2268 KSP array must be. 2269 2270 Level: advanced 2271 2272 .seealso: PCFIELDSPLIT 2273 @*/ 2274 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2275 { 2276 PetscFunctionBegin; 2277 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2278 if (n) PetscValidIntPointer(n,2); 2279 CHKERRQ(PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp))); 2280 PetscFunctionReturn(0); 2281 } 2282 2283 /*@ 2284 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement. 2285 The default is the A11 matrix. 2286 2287 Collective on PC 2288 2289 Input Parameters: 2290 + pc - the preconditioner context 2291 . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER 2292 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 2293 - userpre - matrix to use for preconditioning, or NULL 2294 2295 Options Database: 2296 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 2297 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2298 2299 Notes: 2300 $ If ptype is 2301 $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 2302 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 2303 $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 2304 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 2305 $ preconditioner 2306 $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2307 $ to this function). 2308 $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 2309 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 2310 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 2311 $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive) 2312 $ useful mostly as a test that the Schur complement approach can work for your problem 2313 2314 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 2315 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 2316 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 2317 2318 Level: intermediate 2319 2320 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 2321 MatSchurComplementSetAinvType(), PCLSC 2322 2323 @*/ 2324 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2325 { 2326 PetscFunctionBegin; 2327 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2328 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre))); 2329 PetscFunctionReturn(0); 2330 } 2331 2332 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 2333 2334 /*@ 2335 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2336 preconditioned. See PCFieldSplitSetSchurPre() for details. 2337 2338 Logically Collective on PC 2339 2340 Input Parameter: 2341 . pc - the preconditioner context 2342 2343 Output Parameters: 2344 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 2345 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 2346 2347 Level: intermediate 2348 2349 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 2350 2351 @*/ 2352 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2353 { 2354 PetscFunctionBegin; 2355 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2356 CHKERRQ(PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre))); 2357 PetscFunctionReturn(0); 2358 } 2359 2360 /*@ 2361 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 2362 2363 Not collective 2364 2365 Input Parameter: 2366 . pc - the preconditioner context 2367 2368 Output Parameter: 2369 . S - the Schur complement matrix 2370 2371 Notes: 2372 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 2373 2374 Level: advanced 2375 2376 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 2377 2378 @*/ 2379 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 2380 { 2381 const char* t; 2382 PetscBool isfs; 2383 PC_FieldSplit *jac; 2384 2385 PetscFunctionBegin; 2386 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2387 CHKERRQ(PetscObjectGetType((PetscObject)pc,&t)); 2388 CHKERRQ(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2389 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2390 jac = (PC_FieldSplit*)pc->data; 2391 PetscCheckFalse(jac->type != PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2392 if (S) *S = jac->schur; 2393 PetscFunctionReturn(0); 2394 } 2395 2396 /*@ 2397 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 2398 2399 Not collective 2400 2401 Input Parameters: 2402 + pc - the preconditioner context 2403 - S - the Schur complement matrix 2404 2405 Level: advanced 2406 2407 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 2408 2409 @*/ 2410 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2411 { 2412 const char* t; 2413 PetscBool isfs; 2414 PC_FieldSplit *jac; 2415 2416 PetscFunctionBegin; 2417 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2418 CHKERRQ(PetscObjectGetType((PetscObject)pc,&t)); 2419 CHKERRQ(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2420 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2421 jac = (PC_FieldSplit*)pc->data; 2422 PetscCheckFalse(jac->type != PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2423 PetscCheckFalse(!S || *S != jac->schur,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2424 PetscFunctionReturn(0); 2425 } 2426 2427 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2428 { 2429 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2430 2431 PetscFunctionBegin; 2432 jac->schurpre = ptype; 2433 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2434 CHKERRQ(MatDestroy(&jac->schur_user)); 2435 jac->schur_user = pre; 2436 CHKERRQ(PetscObjectReference((PetscObject)jac->schur_user)); 2437 } 2438 PetscFunctionReturn(0); 2439 } 2440 2441 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2442 { 2443 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2444 2445 PetscFunctionBegin; 2446 *ptype = jac->schurpre; 2447 *pre = jac->schur_user; 2448 PetscFunctionReturn(0); 2449 } 2450 2451 /*@ 2452 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner 2453 2454 Collective on PC 2455 2456 Input Parameters: 2457 + pc - the preconditioner context 2458 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2459 2460 Options Database: 2461 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full 2462 2463 Level: intermediate 2464 2465 Notes: 2466 The FULL factorization is 2467 2468 .vb 2469 (A B) = (1 0) (A 0) (1 Ainv*B) = L D U 2470 (C E) (C*Ainv 1) (0 S) (0 1) 2471 .vb 2472 where S = E - 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, 2473 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). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale(). 2474 2475 If A and S are solved exactly 2476 .vb 2477 *) FULL factorization is a direct solver. 2478 *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations. 2479 *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 2480 .ve 2481 2482 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2483 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2484 2485 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. 2486 2487 Note that a flexible method like KSPFGMRES or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S). 2488 2489 References: 2490 + * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2491 - * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2492 2493 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale() 2494 @*/ 2495 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2496 { 2497 PetscFunctionBegin; 2498 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2499 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype))); 2500 PetscFunctionReturn(0); 2501 } 2502 2503 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2504 { 2505 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2506 2507 PetscFunctionBegin; 2508 jac->schurfactorization = ftype; 2509 PetscFunctionReturn(0); 2510 } 2511 2512 /*@ 2513 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG. 2514 2515 Collective on PC 2516 2517 Input Parameters: 2518 + pc - the preconditioner context 2519 - scale - scaling factor for the Schur complement 2520 2521 Options Database: 2522 . -pc_fieldsplit_schur_scale - default is -1.0 2523 2524 Level: intermediate 2525 2526 .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale() 2527 @*/ 2528 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale) 2529 { 2530 PetscFunctionBegin; 2531 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2532 PetscValidLogicalCollectiveScalar(pc,scale,2); 2533 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale))); 2534 PetscFunctionReturn(0); 2535 } 2536 2537 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale) 2538 { 2539 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2540 2541 PetscFunctionBegin; 2542 jac->schurscale = scale; 2543 PetscFunctionReturn(0); 2544 } 2545 2546 /*@C 2547 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2548 2549 Collective on KSP 2550 2551 Input Parameter: 2552 . pc - the preconditioner context 2553 2554 Output Parameters: 2555 + A00 - the (0,0) block 2556 . A01 - the (0,1) block 2557 . A10 - the (1,0) block 2558 - A11 - the (1,1) block 2559 2560 Level: advanced 2561 2562 .seealso: PCFIELDSPLIT 2563 @*/ 2564 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2565 { 2566 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2567 2568 PetscFunctionBegin; 2569 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2570 PetscCheckFalse(jac->type != PC_COMPOSITE_SCHUR,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2571 if (A00) *A00 = jac->pmat[0]; 2572 if (A01) *A01 = jac->B; 2573 if (A10) *A10 = jac->C; 2574 if (A11) *A11 = jac->pmat[1]; 2575 PetscFunctionReturn(0); 2576 } 2577 2578 /*@ 2579 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner. 2580 2581 Collective on PC 2582 2583 Notes: 2584 The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. 2585 It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than 2586 this estimate, the stopping criterion is satisfactory in practical cases [A13]. 2587 2588 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2589 2590 Input Parameters: 2591 + pc - the preconditioner context 2592 - tolerance - the solver tolerance 2593 2594 Options Database: 2595 . -pc_fieldsplit_gkb_tol - default is 1e-5 2596 2597 Level: intermediate 2598 2599 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit() 2600 @*/ 2601 PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance) 2602 { 2603 PetscFunctionBegin; 2604 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2605 PetscValidLogicalCollectiveReal(pc,tolerance,2); 2606 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance))); 2607 PetscFunctionReturn(0); 2608 } 2609 2610 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance) 2611 { 2612 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2613 2614 PetscFunctionBegin; 2615 jac->gkbtol = tolerance; 2616 PetscFunctionReturn(0); 2617 } 2618 2619 /*@ 2620 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization 2621 preconditioner. 2622 2623 Collective on PC 2624 2625 Input Parameters: 2626 + pc - the preconditioner context 2627 - maxit - the maximum number of iterations 2628 2629 Options Database: 2630 . -pc_fieldsplit_gkb_maxit - default is 100 2631 2632 Level: intermediate 2633 2634 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu() 2635 @*/ 2636 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit) 2637 { 2638 PetscFunctionBegin; 2639 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2640 PetscValidLogicalCollectiveInt(pc,maxit,2); 2641 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit))); 2642 PetscFunctionReturn(0); 2643 } 2644 2645 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit) 2646 { 2647 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2648 2649 PetscFunctionBegin; 2650 jac->gkbmaxit = maxit; 2651 PetscFunctionReturn(0); 2652 } 2653 2654 /*@ 2655 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization 2656 preconditioner. 2657 2658 Collective on PC 2659 2660 Notes: 2661 The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error ||u-u^k||_H 2662 is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs 2663 at least (delay + 1) iterations to stop. For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to 2664 2665 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2666 2667 Input Parameters: 2668 + pc - the preconditioner context 2669 - delay - the delay window in the lower bound estimate 2670 2671 Options Database: 2672 . -pc_fieldsplit_gkb_delay - default is 5 2673 2674 Level: intermediate 2675 2676 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit() 2677 @*/ 2678 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay) 2679 { 2680 PetscFunctionBegin; 2681 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2682 PetscValidLogicalCollectiveInt(pc,delay,2); 2683 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay))); 2684 PetscFunctionReturn(0); 2685 } 2686 2687 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay) 2688 { 2689 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2690 2691 PetscFunctionBegin; 2692 jac->gkbdelay = delay; 2693 PetscFunctionReturn(0); 2694 } 2695 2696 /*@ 2697 PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the Golub-Kahan bidiagonalization preconditioner. 2698 2699 Collective on PC 2700 2701 Notes: 2702 This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by chosing nu sufficiently big. However, 2703 if nu is chosen too big, the matrix H might be badly conditioned and the solution of the linear system Hx = b in the inner loop gets difficult. It is therefore 2704 necessary to find a good balance in between the convergence of the inner and outer loop. 2705 2706 For nu = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in [Ar13] is then chosen as identity. 2707 2708 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2709 2710 Input Parameters: 2711 + pc - the preconditioner context 2712 - nu - the shift parameter 2713 2714 Options Database: 2715 . -pc_fieldsplit_gkb_nu - default is 1 2716 2717 Level: intermediate 2718 2719 .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit() 2720 @*/ 2721 PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu) 2722 { 2723 PetscFunctionBegin; 2724 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2725 PetscValidLogicalCollectiveReal(pc,nu,2); 2726 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu))); 2727 PetscFunctionReturn(0); 2728 } 2729 2730 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu) 2731 { 2732 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2733 2734 PetscFunctionBegin; 2735 jac->gkbnu = nu; 2736 PetscFunctionReturn(0); 2737 } 2738 2739 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2740 { 2741 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2742 2743 PetscFunctionBegin; 2744 jac->type = type; 2745 2746 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0)); 2747 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0)); 2748 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0)); 2749 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0)); 2750 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0)); 2751 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0)); 2752 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0)); 2753 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0)); 2754 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0)); 2755 2756 if (type == PC_COMPOSITE_SCHUR) { 2757 pc->ops->apply = PCApply_FieldSplit_Schur; 2758 pc->ops->view = PCView_FieldSplit_Schur; 2759 2760 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur)); 2761 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit)); 2762 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit)); 2763 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit)); 2764 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit)); 2765 } else if (type == PC_COMPOSITE_GKB) { 2766 pc->ops->apply = PCApply_FieldSplit_GKB; 2767 pc->ops->view = PCView_FieldSplit_GKB; 2768 2769 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2770 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit)); 2771 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit)); 2772 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit)); 2773 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit)); 2774 } else { 2775 pc->ops->apply = PCApply_FieldSplit; 2776 pc->ops->view = PCView_FieldSplit; 2777 2778 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2779 } 2780 PetscFunctionReturn(0); 2781 } 2782 2783 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2784 { 2785 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2786 2787 PetscFunctionBegin; 2788 PetscCheckFalse(bs < 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2789 PetscCheckFalse(jac->bs > 0 && jac->bs != bs,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 2790 jac->bs = bs; 2791 PetscFunctionReturn(0); 2792 } 2793 2794 static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[]) 2795 { 2796 PC_FieldSplit * jac = (PC_FieldSplit*)pc->data; 2797 PC_FieldSplitLink ilink_current = jac->head; 2798 PetscInt ii; 2799 IS is_owned; 2800 2801 PetscFunctionBegin; 2802 jac->coordinates_set = PETSC_TRUE; // Internal flag 2803 CHKERRQ(MatGetOwnershipIS(pc->mat,&is_owned,PETSC_NULL)); 2804 2805 ii=0; 2806 while (ilink_current) { 2807 // For each IS, embed it to get local coords indces 2808 IS is_coords; 2809 PetscInt ndofs_block; 2810 const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block 2811 2812 // Setting drop to true for safety. It should make no difference. 2813 CHKERRQ(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords)); 2814 CHKERRQ(ISGetLocalSize(is_coords, &ndofs_block)); 2815 CHKERRQ(ISGetIndices(is_coords, &block_dofs_enumeration)); 2816 2817 // Allocate coordinates vector and set it directly 2818 CHKERRQ(PetscMalloc1(ndofs_block * dim, &(ilink_current->coords))); 2819 for (PetscInt dof=0;dof<ndofs_block;++dof) { 2820 for (PetscInt d=0;d<dim;++d) { 2821 (ilink_current->coords)[dim*dof + d] = coords[dim * block_dofs_enumeration[dof] + d]; 2822 } 2823 } 2824 ilink_current->dim = dim; 2825 ilink_current->ndofs = ndofs_block; 2826 CHKERRQ(ISRestoreIndices(is_coords, &block_dofs_enumeration)); 2827 CHKERRQ(ISDestroy(&is_coords)); 2828 ilink_current = ilink_current->next; 2829 ++ii; 2830 } 2831 CHKERRQ(ISDestroy(&is_owned)); 2832 PetscFunctionReturn(0); 2833 } 2834 2835 /*@ 2836 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2837 2838 Collective on PC 2839 2840 Input Parameters: 2841 + pc - the preconditioner context 2842 - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2843 2844 Options Database Key: 2845 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2846 2847 Level: Intermediate 2848 2849 .seealso: PCCompositeSetType() 2850 2851 @*/ 2852 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2853 { 2854 PetscFunctionBegin; 2855 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2856 CHKERRQ(PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type))); 2857 PetscFunctionReturn(0); 2858 } 2859 2860 /*@ 2861 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2862 2863 Not collective 2864 2865 Input Parameter: 2866 . pc - the preconditioner context 2867 2868 Output Parameter: 2869 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2870 2871 Level: Intermediate 2872 2873 .seealso: PCCompositeSetType() 2874 @*/ 2875 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2876 { 2877 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2878 2879 PetscFunctionBegin; 2880 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2881 PetscValidIntPointer(type,2); 2882 *type = jac->type; 2883 PetscFunctionReturn(0); 2884 } 2885 2886 /*@ 2887 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2888 2889 Logically Collective 2890 2891 Input Parameters: 2892 + pc - the preconditioner context 2893 - flg - boolean indicating whether to use field splits defined by the DM 2894 2895 Options Database Key: 2896 . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the DM 2897 2898 Level: Intermediate 2899 2900 .seealso: PCFieldSplitGetDMSplits() 2901 2902 @*/ 2903 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2904 { 2905 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2906 PetscBool isfs; 2907 2908 PetscFunctionBegin; 2909 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2910 PetscValidLogicalCollectiveBool(pc,flg,2); 2911 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2912 if (isfs) { 2913 jac->dm_splits = flg; 2914 } 2915 PetscFunctionReturn(0); 2916 } 2917 2918 /*@ 2919 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2920 2921 Logically Collective 2922 2923 Input Parameter: 2924 . pc - the preconditioner context 2925 2926 Output Parameter: 2927 . flg - boolean indicating whether to use field splits defined by the DM 2928 2929 Level: Intermediate 2930 2931 .seealso: PCFieldSplitSetDMSplits() 2932 2933 @*/ 2934 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2935 { 2936 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2937 PetscBool isfs; 2938 2939 PetscFunctionBegin; 2940 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2941 PetscValidBoolPointer(flg,2); 2942 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2943 if (isfs) { 2944 if (flg) *flg = jac->dm_splits; 2945 } 2946 PetscFunctionReturn(0); 2947 } 2948 2949 /*@ 2950 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2951 2952 Logically Collective 2953 2954 Input Parameter: 2955 . pc - the preconditioner context 2956 2957 Output Parameter: 2958 . flg - boolean indicating whether to detect fields or not 2959 2960 Level: Intermediate 2961 2962 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint() 2963 2964 @*/ 2965 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg) 2966 { 2967 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2968 2969 PetscFunctionBegin; 2970 *flg = jac->detect; 2971 PetscFunctionReturn(0); 2972 } 2973 2974 /*@ 2975 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2976 2977 Logically Collective 2978 2979 Input Parameter: 2980 . pc - the preconditioner context 2981 2982 Output Parameter: 2983 . flg - boolean indicating whether to detect fields or not 2984 2985 Options Database Key: 2986 . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point 2987 2988 Notes: 2989 Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()). 2990 2991 Level: Intermediate 2992 2993 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre() 2994 2995 @*/ 2996 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg) 2997 { 2998 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2999 3000 PetscFunctionBegin; 3001 jac->detect = flg; 3002 if (jac->detect) { 3003 CHKERRQ(PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR)); 3004 CHKERRQ(PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL)); 3005 } 3006 PetscFunctionReturn(0); 3007 } 3008 3009 /* -------------------------------------------------------------------------------------*/ 3010 /*MC 3011 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 3012 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 3013 3014 To set options on the solvers for each block append -fieldsplit_ to all the PC 3015 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 3016 3017 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 3018 and set the options directly on the resulting KSP object 3019 3020 Level: intermediate 3021 3022 Options Database Keys: 3023 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 3024 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 3025 been supplied explicitly by -pc_fieldsplit_%d_fields 3026 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 3027 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting 3028 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 3029 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using -pc_fieldsplit_type schur; see PCFieldSplitSetSchurFactType() 3030 - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 3031 3032 Options prefixes for inner solvers when using the Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ . 3033 For all other solvers they are -fieldsplit_%d_ for the dth field; use -fieldsplit_ for all fields. 3034 The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is -fieldsplit_0_ 3035 3036 Notes: 3037 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 3038 to define a field by an arbitrary collection of entries. 3039 3040 If no fields are set the default is used. The fields are defined by entries strided by bs, 3041 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 3042 if this is not called the block size defaults to the blocksize of the second matrix passed 3043 to KSPSetOperators()/PCSetOperators(). 3044 3045 $ For the Schur complement preconditioner if J = [ A00 A01] 3046 $ [ A10 A11] 3047 $ the preconditioner using full factorization is 3048 $ [ I -ksp(A00) A01] [ inv(A00) 0 ] [ I 0] 3049 $ [ 0 I ] [ 0 ksp(S)] [ -A10 ksp(A00) I] 3050 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 3051 $ S = A11 - A10 ksp(A00) A01 3052 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 3053 in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0, 3054 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 3055 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 3056 3057 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 3058 diag gives 3059 $ [ inv(A00) 0 ] 3060 $ [ 0 -ksp(S)] 3061 note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip 3062 can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of 3063 $ [ A00 0] 3064 $ [ A10 S] 3065 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 3066 $ [ A00 A01] 3067 $ [ 0 S ] 3068 where again the inverses of A00 and S are applied using KSPs. 3069 3070 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 3071 is used automatically for a second block. 3072 3073 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 3074 Generally it should be used with the AIJ format. 3075 3076 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 3077 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 3078 inside a smoother resulting in "Distributive Smoothers". 3079 3080 References: 3081 . * - A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 3082 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 3083 http://chess.cs.umd.edu/~elman/papers/tax.pdf 3084 3085 The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the 3086 residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables. 3087 3088 The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape 3089 $ [ A00 A01] 3090 $ [ A01' 0 ] 3091 with A00 positive semi-definite. The implementation follows [Ar13]. Therein, we choose N := 1/nu * I and the (1,1)-block of the matrix is modified to H = A00 + nu*A01*A01'. 3092 A linear system Hx = b has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix -fieldsplit_0_. 3093 3094 . * - [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 3095 3096 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCLSC, 3097 PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(), 3098 PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), PCFieldSplitSetSchurFactType(), 3099 MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), PCFieldSplitSetDetectSaddlePoint() 3100 M*/ 3101 3102 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 3103 { 3104 PC_FieldSplit *jac; 3105 3106 PetscFunctionBegin; 3107 CHKERRQ(PetscNewLog(pc,&jac)); 3108 3109 jac->bs = -1; 3110 jac->nsplits = 0; 3111 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 3112 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 3113 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 3114 jac->schurscale = -1.0; 3115 jac->dm_splits = PETSC_TRUE; 3116 jac->detect = PETSC_FALSE; 3117 jac->gkbtol = 1e-5; 3118 jac->gkbdelay = 5; 3119 jac->gkbnu = 1; 3120 jac->gkbmaxit = 100; 3121 jac->gkbmonitor = PETSC_FALSE; 3122 jac->coordinates_set = PETSC_FALSE; 3123 3124 pc->data = (void*)jac; 3125 3126 pc->ops->apply = PCApply_FieldSplit; 3127 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 3128 pc->ops->setup = PCSetUp_FieldSplit; 3129 pc->ops->reset = PCReset_FieldSplit; 3130 pc->ops->destroy = PCDestroy_FieldSplit; 3131 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 3132 pc->ops->view = PCView_FieldSplit; 3133 pc->ops->applyrichardson = NULL; 3134 3135 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit)); 3136 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 3137 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit)); 3138 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit)); 3139 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit)); 3140 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit)); 3141 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit)); 3142 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_FieldSplit)); 3143 PetscFunctionReturn(0); 3144 } 3145