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