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