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://chess.cs.umd.edu/~elman/papers/tax.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 132 PetscFunctionBegin; 133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 134 if (iascii) { 135 if (jac->bs > 0) { 136 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 137 } else { 138 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 139 } 140 if (jac->realdiagonal) { 141 ierr = PetscViewerASCIIPrintf(viewer," using actual matrix for blocks rather than preconditioner matrix\n");CHKERRQ(ierr); 142 } 143 switch(jac->schurpre) { 144 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 145 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);break; 146 case PC_FIELDSPLIT_SCHUR_PRE_DIAG: 147 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the block diagonal part of A11\n");CHKERRQ(ierr);break; 148 case PC_FIELDSPLIT_SCHUR_PRE_USER: 149 if (jac->schur_user) { 150 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 151 } else { 152 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the block diagonal part of A11\n");CHKERRQ(ierr); 153 } 154 break; 155 default: 156 SETERRQ1(((PetscObject) pc)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 157 } 158 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 159 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 160 for (i=0; i<jac->nsplits; i++) { 161 if (ilink->fields) { 162 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 163 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 164 for (j=0; j<ilink->nfields; j++) { 165 if (j > 0) { 166 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 167 } 168 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 169 } 170 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 171 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 172 } else { 173 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 174 } 175 ilink = ilink->next; 176 } 177 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block \n");CHKERRQ(ierr); 178 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 179 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 180 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 181 if (jac->kspupper != jac->head->ksp) { 182 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 183 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 184 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 185 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 186 } 187 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 188 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 189 if (jac->kspschur) { 190 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 191 } else { 192 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 193 } 194 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 195 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 196 } else { 197 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 198 } 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 204 /* Precondition: jac->bs is set to a meaningful value */ 205 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 206 { 207 PetscErrorCode ierr; 208 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 209 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 210 PetscBool flg,flg_col; 211 char optionname[128],splitname[8],optionname_col[128]; 212 213 PetscFunctionBegin; 214 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 215 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields_col);CHKERRQ(ierr); 216 for (i=0,flg=PETSC_TRUE; ; i++) { 217 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 218 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 219 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 220 nfields = jac->bs; 221 nfields_col = jac->bs; 222 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 223 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 224 if (!flg) break; 225 else if (flg && !flg_col) { 226 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 227 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 228 } 229 else { 230 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 231 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 232 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 233 } 234 } 235 if (i > 0) { 236 /* Makes command-line setting of splits take precedence over setting them in code. 237 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 238 create new splits, which would probably not be what the user wanted. */ 239 jac->splitdefined = PETSC_TRUE; 240 } 241 ierr = PetscFree(ifields);CHKERRQ(ierr); 242 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 #undef __FUNCT__ 247 #define __FUNCT__ "PCFieldSplitSetDefaults" 248 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 249 { 250 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 251 PetscErrorCode ierr; 252 PC_FieldSplitLink ilink = jac->head; 253 PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE; 254 PetscInt i; 255 256 PetscFunctionBegin; 257 /* 258 Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset. 259 Should probably be rewritten. 260 */ 261 if (!ilink || jac->reset) { 262 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 263 if (pc->dm && !stokes) { 264 PetscInt numFields, f, i, j; 265 char **fieldNames; 266 IS *fields; 267 DM *dms; 268 DM subdm[128]; 269 PetscBool flg; 270 271 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 272 /* Allow the user to prescribe the splits */ 273 for (i = 0, flg = PETSC_TRUE; ; i++) { 274 PetscInt ifields[128]; 275 IS compField; 276 char optionname[128], splitname[8]; 277 PetscInt nfields = numFields; 278 279 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 280 ierr = PetscOptionsGetIntArray(((PetscObject) pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 281 if (!flg) break; 282 if (numFields > 128) SETERRQ1(((PetscObject) pc)->comm,PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 283 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 284 if (nfields == 1) { 285 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 286 /* ierr = PetscPrintf(((PetscObject) pc)->comm, "%s Field Indices:", fieldNames[ifields[0]]);CHKERRQ(ierr); 287 ierr = ISView(compField, PETSC_NULL);CHKERRQ(ierr); */ 288 } else { 289 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 290 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 291 /* ierr = PetscPrintf(((PetscObject) pc)->comm, "%s Field Indices:", splitname);CHKERRQ(ierr); 292 ierr = ISView(compField, PETSC_NULL);CHKERRQ(ierr); */ 293 } 294 ierr = ISDestroy(&compField);CHKERRQ(ierr); 295 for (j = 0; j < nfields; ++j) { 296 f = ifields[j]; 297 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 298 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 299 } 300 } 301 if (i == 0) { 302 for (f = 0; f < numFields; ++f) { 303 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 304 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 305 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 306 } 307 } else { 308 ierr = PetscMalloc(i * sizeof(DM), &dms);CHKERRQ(ierr); 309 for (j = 0; j < i; ++j) { 310 dms[j] = subdm[j]; 311 } 312 } 313 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 314 ierr = PetscFree(fields);CHKERRQ(ierr); 315 if (dms) { 316 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 317 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 318 const char *prefix; 319 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix); CHKERRQ(ierr); 320 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix); CHKERRQ(ierr); 321 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 322 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 323 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0); CHKERRQ(ierr); 324 ierr = DMDestroy(&dms[i]); CHKERRQ(ierr); 325 } 326 ierr = PetscFree(dms);CHKERRQ(ierr); 327 } 328 } else { 329 if (jac->bs <= 0) { 330 if (pc->pmat) { 331 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 332 } else { 333 jac->bs = 1; 334 } 335 } 336 337 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,PETSC_NULL);CHKERRQ(ierr); 338 if (stokes) { 339 IS zerodiags,rest; 340 PetscInt nmin,nmax; 341 342 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 343 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 344 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 345 if (jac->reset) { 346 jac->head->is = rest; 347 jac->head->next->is = zerodiags; 348 } 349 else { 350 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 351 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 352 } 353 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 354 ierr = ISDestroy(&rest);CHKERRQ(ierr); 355 } else { 356 if (jac->reset) 357 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 358 if (!fieldsplit_default) { 359 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 360 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 361 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 362 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 363 } 364 if (fieldsplit_default || !jac->splitdefined) { 365 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 366 for (i=0; i<jac->bs; i++) { 367 char splitname[8]; 368 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 369 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 370 } 371 jac->defaultsplit = PETSC_TRUE; 372 } 373 } 374 } 375 } else if (jac->nsplits == 1) { 376 if (ilink->is) { 377 IS is2; 378 PetscInt nmin,nmax; 379 380 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 381 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 382 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 383 ierr = ISDestroy(&is2);CHKERRQ(ierr); 384 } else SETERRQ(((PetscObject) pc)->comm,PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 385 } 386 387 388 if (jac->nsplits < 2) SETERRQ1(((PetscObject) pc)->comm,PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 389 PetscFunctionReturn(0); 390 } 391 392 PetscErrorCode PetscOptionsFindPairPrefix_Private(const char pre[], const char name[], char *value[], PetscBool *flg); 393 394 #undef __FUNCT__ 395 #define __FUNCT__ "PCSetUp_FieldSplit" 396 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 397 { 398 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 399 PetscErrorCode ierr; 400 PC_FieldSplitLink ilink; 401 PetscInt i,nsplit; 402 MatStructure flag = pc->flag; 403 PetscBool sorted, sorted_col; 404 405 PetscFunctionBegin; 406 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 407 nsplit = jac->nsplits; 408 ilink = jac->head; 409 410 /* get the matrices for each split */ 411 if (!jac->issetup) { 412 PetscInt rstart,rend,nslots,bs; 413 414 jac->issetup = PETSC_TRUE; 415 416 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 417 if (jac->defaultsplit || !ilink->is) { 418 if (jac->bs <= 0) jac->bs = nsplit; 419 } 420 bs = jac->bs; 421 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 422 nslots = (rend - rstart)/bs; 423 for (i=0; i<nsplit; i++) { 424 if (jac->defaultsplit) { 425 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 426 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 427 } else if (!ilink->is) { 428 if (ilink->nfields > 1) { 429 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 430 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 431 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&jj);CHKERRQ(ierr); 432 for (j=0; j<nslots; j++) { 433 for (k=0; k<nfields; k++) { 434 ii[nfields*j + k] = rstart + bs*j + fields[k]; 435 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 436 } 437 } 438 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 439 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 440 } else { 441 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 442 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 443 } 444 } 445 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 446 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 447 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 448 ilink = ilink->next; 449 } 450 } 451 452 ilink = jac->head; 453 if (!jac->pmat) { 454 Vec xtmp; 455 456 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 457 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 458 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 459 for (i=0; i<nsplit; i++) { 460 MatNullSpace sp; 461 462 /* Check for preconditioning matrix attached to IS */ 463 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject *) &jac->pmat[i]);CHKERRQ(ierr); 464 if (jac->pmat[i]) { 465 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 466 if (jac->type == PC_COMPOSITE_SCHUR) { 467 jac->schur_user = jac->pmat[i]; 468 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 469 } 470 } else { 471 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 472 } 473 /* create work vectors for each split */ 474 ierr = MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 475 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = PETSC_NULL; 476 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 477 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 478 /* Check for null space attached to IS */ 479 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject *) &sp);CHKERRQ(ierr); 480 if (sp) { 481 ierr = MatSetNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 482 } 483 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject *) &sp);CHKERRQ(ierr); 484 if (sp) { 485 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 486 } 487 ilink = ilink->next; 488 } 489 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 490 } else { 491 for (i=0; i<nsplit; i++) { 492 Mat pmat; 493 494 /* Check for preconditioning matrix attached to IS */ 495 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject *) &pmat);CHKERRQ(ierr); 496 if (!pmat) { 497 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 498 } 499 ilink = ilink->next; 500 } 501 } 502 if (jac->realdiagonal) { 503 ilink = jac->head; 504 if (!jac->mat) { 505 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 506 for (i=0; i<nsplit; i++) { 507 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 508 ilink = ilink->next; 509 } 510 } else { 511 for (i=0; i<nsplit; i++) { 512 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 513 ilink = ilink->next; 514 } 515 } 516 } else { 517 jac->mat = jac->pmat; 518 } 519 520 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 521 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 522 ilink = jac->head; 523 if (!jac->Afield) { 524 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 525 for (i=0; i<nsplit; i++) { 526 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 527 ilink = ilink->next; 528 } 529 } else { 530 for (i=0; i<nsplit; i++) { 531 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 532 ilink = ilink->next; 533 } 534 } 535 } 536 537 if (jac->type == PC_COMPOSITE_SCHUR) { 538 IS ccis; 539 PetscInt rstart,rend; 540 char lscname[256]; 541 PetscObject LSC_L; 542 if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 543 544 /* When extracting off-diagonal submatrices, we take complements from this range */ 545 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 546 547 /* need to handle case when one is resetting up the preconditioner */ 548 if (jac->schur) { 549 KSP kspA = jac->head->ksp, kspInner = PETSC_NULL, kspUpper = jac->kspupper; 550 551 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 552 ilink = jac->head; 553 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 554 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 555 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 556 ilink = ilink->next; 557 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 558 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 559 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 560 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr); 561 if (kspA != kspInner) { 562 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 563 } 564 if (kspUpper != kspA) { 565 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 566 } 567 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 568 } else { 569 KSP ksp; 570 const char *Dprefix; 571 char schurprefix[256]; 572 char schurtestoption[256]; 573 MatNullSpace sp; 574 PetscBool flg; 575 576 /* extract the A01 and A10 matrices */ 577 ilink = jac->head; 578 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 579 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 580 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 581 ilink = ilink->next; 582 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 583 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 584 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 585 586 /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] to define Schur complement */ 587 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 588 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 589 ierr = MatSchurComplementGetKSP(jac->schur, &ksp);CHKERRQ(ierr); 590 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 591 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 592 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject) pc, 2);CHKERRQ(ierr); 593 ierr = KSPSetOptionsPrefix(ksp, schurprefix);CHKERRQ(ierr); 594 ierr = MatSchurComplementSet(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 595 596 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 597 if (sp) { 598 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 599 } 600 601 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 602 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, PETSC_NULL, &flg);CHKERRQ(ierr); 603 if (flg) { 604 DM dmInner; 605 606 /* Set DM for new solver */ 607 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 608 ierr = KSPSetDM(ksp, dmInner);CHKERRQ(ierr); 609 ierr = KSPSetDMActive(ksp, PETSC_FALSE);CHKERRQ(ierr); 610 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 611 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 612 } else { 613 ierr = MatSchurComplementSetKSP(jac->schur, jac->head->ksp);CHKERRQ(ierr); 614 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 615 } 616 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 617 618 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 619 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, PETSC_NULL, &flg);CHKERRQ(ierr); 620 if (flg) { 621 DM dmInner; 622 623 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 624 ierr = KSPCreate(((PetscObject) pc)->comm, &jac->kspupper);CHKERRQ(ierr); 625 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 626 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 627 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 628 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 629 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 630 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 631 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 632 } else { 633 jac->kspupper = jac->head->ksp; 634 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 635 } 636 637 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur); CHKERRQ(ierr); 638 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 639 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1); CHKERRQ(ierr); 640 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 641 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 642 PC pcschur; 643 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 644 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 645 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 646 } 647 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix); CHKERRQ(ierr); 648 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix); CHKERRQ(ierr); 649 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 650 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 651 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 652 } 653 654 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 655 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 656 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 657 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 658 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 659 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 660 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 661 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 662 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 663 } else { 664 /* set up the individual splits' PCs */ 665 i = 0; 666 ilink = jac->head; 667 while (ilink) { 668 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 669 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 670 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 671 i++; 672 ilink = ilink->next; 673 } 674 } 675 676 jac->suboptionsset = PETSC_TRUE; 677 PetscFunctionReturn(0); 678 } 679 680 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 681 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 682 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 683 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 684 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 685 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 686 687 #undef __FUNCT__ 688 #define __FUNCT__ "PCApply_FieldSplit_Schur" 689 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 690 { 691 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 692 PetscErrorCode ierr; 693 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 694 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 695 696 PetscFunctionBegin; 697 switch (jac->schurfactorization) { 698 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 699 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 700 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 701 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 702 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 703 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 704 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 705 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 706 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 707 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 708 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 709 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 710 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 711 break; 712 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 713 /* [A00 0; A10 S], suitable for left preconditioning */ 714 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 715 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 716 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 717 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 718 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 719 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 720 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 721 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 722 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 723 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 724 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 725 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 726 break; 727 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 728 /* [A00 A01; 0 S], suitable for right preconditioning */ 729 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 730 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 731 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 732 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 733 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 734 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 735 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 736 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 737 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 738 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 739 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 740 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 741 break; 742 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 743 /* [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 */ 744 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 745 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 746 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 747 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 748 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 749 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 750 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 751 752 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 753 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 754 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 755 756 if (kspUpper == kspA) { 757 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 758 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 759 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 760 } else { 761 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 762 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 763 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 764 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 765 } 766 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 767 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 768 } 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNCT__ 773 #define __FUNCT__ "PCApply_FieldSplit" 774 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 775 { 776 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 777 PetscErrorCode ierr; 778 PC_FieldSplitLink ilink = jac->head; 779 PetscInt cnt,bs; 780 781 PetscFunctionBegin; 782 x->map->bs = jac->bs; 783 y->map->bs = jac->bs; 784 CHKMEMQ; 785 if (jac->type == PC_COMPOSITE_ADDITIVE) { 786 if (jac->defaultsplit) { 787 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 788 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); 789 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 790 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); 791 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 792 while (ilink) { 793 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 794 ilink = ilink->next; 795 } 796 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 797 } else { 798 ierr = VecSet(y,0.0);CHKERRQ(ierr); 799 while (ilink) { 800 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 801 ilink = ilink->next; 802 } 803 } 804 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 805 if (!jac->w1) { 806 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 807 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 808 } 809 ierr = VecSet(y,0.0);CHKERRQ(ierr); 810 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 811 cnt = 1; 812 while (ilink->next) { 813 ilink = ilink->next; 814 /* compute the residual only over the part of the vector needed */ 815 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 816 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 817 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 818 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 819 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 820 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 821 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 822 } 823 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 824 cnt -= 2; 825 while (ilink->previous) { 826 ilink = ilink->previous; 827 /* compute the residual only over the part of the vector needed */ 828 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 829 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 830 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 831 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 832 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 833 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 834 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 835 } 836 } 837 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 838 CHKMEMQ; 839 PetscFunctionReturn(0); 840 } 841 842 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 843 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 844 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 845 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 846 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 847 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 848 849 #undef __FUNCT__ 850 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 851 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 852 { 853 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 854 PetscErrorCode ierr; 855 PC_FieldSplitLink ilink = jac->head; 856 PetscInt bs; 857 858 PetscFunctionBegin; 859 CHKMEMQ; 860 if (jac->type == PC_COMPOSITE_ADDITIVE) { 861 if (jac->defaultsplit) { 862 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 863 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); 864 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 865 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); 866 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 867 while (ilink) { 868 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 869 ilink = ilink->next; 870 } 871 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 872 } else { 873 ierr = VecSet(y,0.0);CHKERRQ(ierr); 874 while (ilink) { 875 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 876 ilink = ilink->next; 877 } 878 } 879 } else { 880 if (!jac->w1) { 881 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 882 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 883 } 884 ierr = VecSet(y,0.0);CHKERRQ(ierr); 885 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 886 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 887 while (ilink->next) { 888 ilink = ilink->next; 889 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 890 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 891 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 892 } 893 while (ilink->previous) { 894 ilink = ilink->previous; 895 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 896 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 897 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 898 } 899 } else { 900 while (ilink->next) { /* get to last entry in linked list */ 901 ilink = ilink->next; 902 } 903 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 904 while (ilink->previous) { 905 ilink = ilink->previous; 906 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 907 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 908 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 909 } 910 } 911 } 912 CHKMEMQ; 913 PetscFunctionReturn(0); 914 } 915 916 #undef __FUNCT__ 917 #define __FUNCT__ "PCReset_FieldSplit" 918 static PetscErrorCode PCReset_FieldSplit(PC pc) 919 { 920 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 921 PetscErrorCode ierr; 922 PC_FieldSplitLink ilink = jac->head,next; 923 924 PetscFunctionBegin; 925 while (ilink) { 926 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 927 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 928 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 929 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 930 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 931 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 932 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 933 next = ilink->next; 934 ilink = next; 935 } 936 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 937 if (jac->mat && jac->mat != jac->pmat) { 938 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 939 } else if (jac->mat) { 940 jac->mat = PETSC_NULL; 941 } 942 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 943 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 944 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 945 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 946 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 947 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 948 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 949 ierr = KSPDestroy(&jac->kspupper);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