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