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 = MatDestroy(&jac->B);CHKERRQ(ierr); 950 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 951 jac->reset = PETSC_TRUE; 952 PetscFunctionReturn(0); 953 } 954 955 #undef __FUNCT__ 956 #define __FUNCT__ "PCDestroy_FieldSplit" 957 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 958 { 959 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 960 PetscErrorCode ierr; 961 PC_FieldSplitLink ilink = jac->head,next; 962 963 PetscFunctionBegin; 964 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 965 while (ilink) { 966 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 967 next = ilink->next; 968 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 969 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 970 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 971 ierr = PetscFree(ilink);CHKERRQ(ierr); 972 ilink = next; 973 } 974 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 975 ierr = PetscFree(pc->data);CHKERRQ(ierr); 976 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 977 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 978 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 979 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 980 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 981 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",PETSC_NULL);CHKERRQ(ierr); 982 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",PETSC_NULL);CHKERRQ(ierr); 983 PetscFunctionReturn(0); 984 } 985 986 #undef __FUNCT__ 987 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 988 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 989 { 990 PetscErrorCode ierr; 991 PetscInt bs; 992 PetscBool flg,stokes = PETSC_FALSE; 993 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 994 PCCompositeType ctype; 995 DM ddm; 996 char ddm_name[1025]; 997 998 PetscFunctionBegin; 999 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1000 if (pc->dm) { 1001 /* Allow the user to request a decomposition DM by name */ 1002 ierr = PetscStrncpy(ddm_name, "", 1024); CHKERRQ(ierr); 1003 ierr = PetscOptionsString("-pc_fieldsplit_decomposition", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); 1004 if (flg) { 1005 ierr = DMCreateFieldDecompositionDM(pc->dm, ddm_name, &ddm); CHKERRQ(ierr); 1006 if (!ddm) { 1007 SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown field decomposition name %s", ddm_name); 1008 } 1009 ierr = PetscInfo(pc,"Using field decomposition DM defined using options database\n");CHKERRQ(ierr); 1010 ierr = PCSetDM(pc,ddm); CHKERRQ(ierr); 1011 } 1012 } 1013 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 1014 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1015 if (flg) { 1016 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1017 } 1018 1019 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 1020 if (stokes) { 1021 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1022 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1023 } 1024 1025 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1026 if (flg) { 1027 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1028 } 1029 /* Only setup fields once */ 1030 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1031 /* only allow user to set fields from command line if bs is already known. 1032 otherwise user can set them in PCFieldSplitSetDefaults() */ 1033 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1034 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1035 } 1036 if (jac->type == PC_COMPOSITE_SCHUR) { 1037 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1038 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1039 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); 1040 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); 1041 } 1042 ierr = PetscOptionsTail();CHKERRQ(ierr); 1043 PetscFunctionReturn(0); 1044 } 1045 1046 /*------------------------------------------------------------------------------------*/ 1047 1048 EXTERN_C_BEGIN 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1051 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1052 { 1053 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1054 PetscErrorCode ierr; 1055 PC_FieldSplitLink ilink,next = jac->head; 1056 char prefix[128]; 1057 PetscInt i; 1058 1059 PetscFunctionBegin; 1060 if (jac->splitdefined) { 1061 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 for (i=0; i<n; i++) { 1065 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1066 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1067 } 1068 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1069 if (splitname) { 1070 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1071 } else { 1072 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1073 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1074 } 1075 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 1076 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1077 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 1078 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1079 ilink->nfields = n; 1080 ilink->next = PETSC_NULL; 1081 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1082 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1083 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1084 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1085 1086 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1087 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1088 1089 if (!next) { 1090 jac->head = ilink; 1091 ilink->previous = PETSC_NULL; 1092 } else { 1093 while (next->next) { 1094 next = next->next; 1095 } 1096 next->next = ilink; 1097 ilink->previous = next; 1098 } 1099 jac->nsplits++; 1100 PetscFunctionReturn(0); 1101 } 1102 EXTERN_C_END 1103 1104 EXTERN_C_BEGIN 1105 #undef __FUNCT__ 1106 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1107 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1108 { 1109 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1110 PetscErrorCode ierr; 1111 1112 PetscFunctionBegin; 1113 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1114 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1115 (*subksp)[1] = jac->kspschur; 1116 if (n) *n = jac->nsplits; 1117 PetscFunctionReturn(0); 1118 } 1119 EXTERN_C_END 1120 1121 EXTERN_C_BEGIN 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1124 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1125 { 1126 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1127 PetscErrorCode ierr; 1128 PetscInt cnt = 0; 1129 PC_FieldSplitLink ilink = jac->head; 1130 1131 PetscFunctionBegin; 1132 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1133 while (ilink) { 1134 (*subksp)[cnt++] = ilink->ksp; 1135 ilink = ilink->next; 1136 } 1137 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); 1138 if (n) *n = jac->nsplits; 1139 PetscFunctionReturn(0); 1140 } 1141 EXTERN_C_END 1142 1143 EXTERN_C_BEGIN 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1146 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1147 { 1148 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1149 PetscErrorCode ierr; 1150 PC_FieldSplitLink ilink, next = jac->head; 1151 char prefix[128]; 1152 1153 PetscFunctionBegin; 1154 if (jac->splitdefined) { 1155 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1156 PetscFunctionReturn(0); 1157 } 1158 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1159 if (splitname) { 1160 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1161 } else { 1162 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1163 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1164 } 1165 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1166 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1167 ilink->is = is; 1168 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1169 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1170 ilink->is_col = is; 1171 ilink->next = PETSC_NULL; 1172 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1173 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1174 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1175 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1176 1177 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1178 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1179 1180 if (!next) { 1181 jac->head = ilink; 1182 ilink->previous = PETSC_NULL; 1183 } else { 1184 while (next->next) { 1185 next = next->next; 1186 } 1187 next->next = ilink; 1188 ilink->previous = next; 1189 } 1190 jac->nsplits++; 1191 1192 PetscFunctionReturn(0); 1193 } 1194 EXTERN_C_END 1195 1196 #undef __FUNCT__ 1197 #define __FUNCT__ "PCFieldSplitSetFields" 1198 /*@ 1199 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1200 1201 Logically Collective on PC 1202 1203 Input Parameters: 1204 + pc - the preconditioner context 1205 . splitname - name of this split, if PETSC_NULL the number of the split is used 1206 . n - the number of fields in this split 1207 - fields - the fields in this split 1208 1209 Level: intermediate 1210 1211 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1212 1213 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1214 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 1215 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1216 where the numbered entries indicate what is in the field. 1217 1218 This function is called once per split (it creates a new split each time). Solve options 1219 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1220 1221 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1222 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1223 available when this routine is called. 1224 1225 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1226 1227 @*/ 1228 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1229 { 1230 PetscErrorCode ierr; 1231 1232 PetscFunctionBegin; 1233 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1234 PetscValidCharPointer(splitname,2); 1235 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1236 PetscValidIntPointer(fields,3); 1237 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1238 PetscFunctionReturn(0); 1239 } 1240 1241 #undef __FUNCT__ 1242 #define __FUNCT__ "PCFieldSplitSetIS" 1243 /*@ 1244 PCFieldSplitSetIS - Sets the exact elements for field 1245 1246 Logically Collective on PC 1247 1248 Input Parameters: 1249 + pc - the preconditioner context 1250 . splitname - name of this split, if PETSC_NULL the number of the split is used 1251 - is - the index set that defines the vector elements in this field 1252 1253 1254 Notes: 1255 Use PCFieldSplitSetFields(), for fields defined by strided types. 1256 1257 This function is called once per split (it creates a new split each time). Solve options 1258 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1259 1260 Level: intermediate 1261 1262 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1263 1264 @*/ 1265 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1266 { 1267 PetscErrorCode ierr; 1268 1269 PetscFunctionBegin; 1270 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1271 if (splitname) PetscValidCharPointer(splitname,2); 1272 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1273 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1274 PetscFunctionReturn(0); 1275 } 1276 1277 #undef __FUNCT__ 1278 #define __FUNCT__ "PCFieldSplitGetIS" 1279 /*@ 1280 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1281 1282 Logically Collective on PC 1283 1284 Input Parameters: 1285 + pc - the preconditioner context 1286 - splitname - name of this split 1287 1288 Output Parameter: 1289 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1290 1291 Level: intermediate 1292 1293 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1294 1295 @*/ 1296 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1297 { 1298 PetscErrorCode ierr; 1299 1300 PetscFunctionBegin; 1301 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1302 PetscValidCharPointer(splitname,2); 1303 PetscValidPointer(is,3); 1304 { 1305 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1306 PC_FieldSplitLink ilink = jac->head; 1307 PetscBool found; 1308 1309 *is = PETSC_NULL; 1310 while(ilink) { 1311 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1312 if (found) { 1313 *is = ilink->is; 1314 break; 1315 } 1316 ilink = ilink->next; 1317 } 1318 } 1319 PetscFunctionReturn(0); 1320 } 1321 1322 #undef __FUNCT__ 1323 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1324 /*@ 1325 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1326 fieldsplit preconditioner. If not set the matrix block size is used. 1327 1328 Logically Collective on PC 1329 1330 Input Parameters: 1331 + pc - the preconditioner context 1332 - bs - the block size 1333 1334 Level: intermediate 1335 1336 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1337 1338 @*/ 1339 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1340 { 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1345 PetscValidLogicalCollectiveInt(pc,bs,2); 1346 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1347 PetscFunctionReturn(0); 1348 } 1349 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1352 /*@C 1353 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1354 1355 Collective on KSP 1356 1357 Input Parameter: 1358 . pc - the preconditioner context 1359 1360 Output Parameters: 1361 + n - the number of splits 1362 - pc - the array of KSP contexts 1363 1364 Note: 1365 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1366 (not the KSP just the array that contains them). 1367 1368 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1369 1370 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1371 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_OBJECT,ierr) to determine how large the 1372 KSP array must be. 1373 1374 1375 Level: advanced 1376 1377 .seealso: PCFIELDSPLIT 1378 @*/ 1379 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1380 { 1381 PetscErrorCode ierr; 1382 1383 PetscFunctionBegin; 1384 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1385 if (n) PetscValidIntPointer(n,2); 1386 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1387 PetscFunctionReturn(0); 1388 } 1389 1390 #undef __FUNCT__ 1391 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1392 /*@ 1393 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1394 A11 matrix. Otherwise no preconditioner is used. 1395 1396 Collective on PC 1397 1398 Input Parameters: 1399 + pc - the preconditioner context 1400 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1401 - userpre - matrix to use for preconditioning, or PETSC_NULL 1402 1403 Options Database: 1404 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1405 1406 Notes: 1407 $ If ptype is 1408 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1409 $ to this function). 1410 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1411 $ matrix associated with the Schur complement (i.e. A11) 1412 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1413 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1414 $ preconditioner 1415 1416 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1417 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1418 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1419 1420 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1421 the name since it sets a proceedure to use. 1422 1423 Level: intermediate 1424 1425 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1426 1427 @*/ 1428 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1429 { 1430 PetscErrorCode ierr; 1431 1432 PetscFunctionBegin; 1433 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1434 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1435 PetscFunctionReturn(0); 1436 } 1437 1438 EXTERN_C_BEGIN 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1441 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1442 { 1443 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1444 PetscErrorCode ierr; 1445 1446 PetscFunctionBegin; 1447 jac->schurpre = ptype; 1448 if (pre) { 1449 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1450 jac->schur_user = pre; 1451 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1452 } 1453 PetscFunctionReturn(0); 1454 } 1455 EXTERN_C_END 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1459 /*@ 1460 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1461 1462 Collective on PC 1463 1464 Input Parameters: 1465 + pc - the preconditioner context 1466 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1467 1468 Options Database: 1469 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1470 1471 1472 Level: intermediate 1473 1474 Notes: 1475 The FULL factorization is 1476 1477 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1478 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1479 1480 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, 1481 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). 1482 1483 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 1484 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 1485 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, 1486 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1487 1488 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 1489 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). 1490 1491 References: 1492 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1493 1494 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1495 1496 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1497 @*/ 1498 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1499 { 1500 PetscErrorCode ierr; 1501 1502 PetscFunctionBegin; 1503 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1504 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1505 PetscFunctionReturn(0); 1506 } 1507 1508 #undef __FUNCT__ 1509 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1510 PETSC_EXTERN_C PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1511 { 1512 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1513 1514 PetscFunctionBegin; 1515 jac->schurfactorization = ftype; 1516 PetscFunctionReturn(0); 1517 } 1518 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1521 /*@C 1522 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1523 1524 Collective on KSP 1525 1526 Input Parameter: 1527 . pc - the preconditioner context 1528 1529 Output Parameters: 1530 + A00 - the (0,0) block 1531 . A01 - the (0,1) block 1532 . A10 - the (1,0) block 1533 - A11 - the (1,1) block 1534 1535 Level: advanced 1536 1537 .seealso: PCFIELDSPLIT 1538 @*/ 1539 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1540 { 1541 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1542 1543 PetscFunctionBegin; 1544 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1545 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1546 if (A00) *A00 = jac->pmat[0]; 1547 if (A01) *A01 = jac->B; 1548 if (A10) *A10 = jac->C; 1549 if (A11) *A11 = jac->pmat[1]; 1550 PetscFunctionReturn(0); 1551 } 1552 1553 EXTERN_C_BEGIN 1554 #undef __FUNCT__ 1555 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1556 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1557 { 1558 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1559 PetscErrorCode ierr; 1560 1561 PetscFunctionBegin; 1562 jac->type = type; 1563 if (type == PC_COMPOSITE_SCHUR) { 1564 pc->ops->apply = PCApply_FieldSplit_Schur; 1565 pc->ops->view = PCView_FieldSplit_Schur; 1566 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1567 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1568 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","PCFieldSplitSetSchurFactType_FieldSplit",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1569 1570 } else { 1571 pc->ops->apply = PCApply_FieldSplit; 1572 pc->ops->view = PCView_FieldSplit; 1573 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1574 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1575 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",0);CHKERRQ(ierr); 1576 } 1577 PetscFunctionReturn(0); 1578 } 1579 EXTERN_C_END 1580 1581 EXTERN_C_BEGIN 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1584 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1585 { 1586 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1587 1588 PetscFunctionBegin; 1589 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1590 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); 1591 jac->bs = bs; 1592 PetscFunctionReturn(0); 1593 } 1594 EXTERN_C_END 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "PCFieldSplitSetType" 1598 /*@ 1599 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1600 1601 Collective on PC 1602 1603 Input Parameter: 1604 . pc - the preconditioner context 1605 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1606 1607 Options Database Key: 1608 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1609 1610 Level: Intermediate 1611 1612 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1613 1614 .seealso: PCCompositeSetType() 1615 1616 @*/ 1617 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1618 { 1619 PetscErrorCode ierr; 1620 1621 PetscFunctionBegin; 1622 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1623 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1624 PetscFunctionReturn(0); 1625 } 1626 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "PCFieldSplitGetType" 1629 /*@ 1630 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1631 1632 Not collective 1633 1634 Input Parameter: 1635 . pc - the preconditioner context 1636 1637 Output Parameter: 1638 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1639 1640 Level: Intermediate 1641 1642 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1643 .seealso: PCCompositeSetType() 1644 @*/ 1645 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1646 { 1647 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1648 1649 PetscFunctionBegin; 1650 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1651 PetscValidIntPointer(type,2); 1652 *type = jac->type; 1653 PetscFunctionReturn(0); 1654 } 1655 1656 /* -------------------------------------------------------------------------------------*/ 1657 /*MC 1658 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1659 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1660 1661 To set options on the solvers for each block append -fieldsplit_ to all the PC 1662 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1663 1664 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1665 and set the options directly on the resulting KSP object 1666 1667 Level: intermediate 1668 1669 Options Database Keys: 1670 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1671 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1672 been supplied explicitly by -pc_fieldsplit_%d_fields 1673 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1674 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1675 . -pc_fieldsplit_schur_precondition <self,user,diag> - default is diag 1676 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1677 1678 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1679 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1680 1681 Notes: 1682 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1683 to define a field by an arbitrary collection of entries. 1684 1685 If no fields are set the default is used. The fields are defined by entries strided by bs, 1686 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1687 if this is not called the block size defaults to the blocksize of the second matrix passed 1688 to KSPSetOperators()/PCSetOperators(). 1689 1690 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1691 $ ( A10 A11 ) 1692 $ the preconditioner using full factorization is 1693 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1694 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1695 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1696 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1697 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1698 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1699 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1700 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1701 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1702 diag gives 1703 $ ( inv(A00) 0 ) 1704 $ ( 0 -ksp(S) ) 1705 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 1706 $ ( A00 0 ) 1707 $ ( A10 S ) 1708 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1709 $ ( A00 A01 ) 1710 $ ( 0 S ) 1711 where again the inverses of A00 and S are applied using KSPs. 1712 1713 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1714 is used automatically for a second block. 1715 1716 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1717 Generally it should be used with the AIJ format. 1718 1719 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1720 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1721 inside a smoother resulting in "Distributive Smoothers". 1722 1723 Concepts: physics based preconditioners, block preconditioners 1724 1725 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1726 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1727 M*/ 1728 1729 EXTERN_C_BEGIN 1730 #undef __FUNCT__ 1731 #define __FUNCT__ "PCCreate_FieldSplit" 1732 PetscErrorCode PCCreate_FieldSplit(PC pc) 1733 { 1734 PetscErrorCode ierr; 1735 PC_FieldSplit *jac; 1736 1737 PetscFunctionBegin; 1738 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1739 jac->bs = -1; 1740 jac->nsplits = 0; 1741 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1742 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1743 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1744 1745 pc->data = (void*)jac; 1746 1747 pc->ops->apply = PCApply_FieldSplit; 1748 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1749 pc->ops->setup = PCSetUp_FieldSplit; 1750 pc->ops->reset = PCReset_FieldSplit; 1751 pc->ops->destroy = PCDestroy_FieldSplit; 1752 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1753 pc->ops->view = PCView_FieldSplit; 1754 pc->ops->applyrichardson = 0; 1755 1756 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1757 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1758 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1759 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1760 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1761 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1762 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1763 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1764 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1765 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1766 PetscFunctionReturn(0); 1767 } 1768 EXTERN_C_END 1769 1770 1771 1772