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 /* 553 Set from options only the A00 split. The other split's solver won't be used with Schur. 554 Should it be destroyed? Should KSPCreate() be moved here from PCFieldSplitSetIS() and invoked 555 only when necessary? 556 */ 557 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 558 559 /* need to handle case when one is resetting up the preconditioner */ 560 if (jac->schur) { 561 ilink = jac->head; 562 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 563 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 564 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 565 ilink = ilink->next; 566 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 567 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 568 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 569 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr); 570 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 571 572 } else { 573 KSP ksp; 574 char schurprefix[256]; 575 MatNullSpace sp; 576 577 /* extract the A01 and A10 matrices */ 578 ilink = jac->head; 579 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 580 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 581 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 582 ilink = ilink->next; 583 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 584 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 585 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 586 /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */ 587 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr); 588 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 589 /* 590 Do we really want to attach the A11-block's nullspace to S? Note that this may generate 591 "false positives", when the A11's nullspace isn't S's: Stokes or A = [1, 1; 1, 0]. 592 Using a false nullspace may prevent KSP(S) from converging, since it might force an inconsistent rhs. 593 */ 594 /* if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);} */ 595 /* set tabbing, options prefix and DM of KSP inside the MatSchurComplement */ 596 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 597 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 598 { 599 PC pcinner; 600 ierr = KSPGetPC(ksp, &pcinner); CHKERRQ(ierr); 601 ierr = PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);CHKERRQ(ierr); 602 } 603 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",jac->head->splitname);CHKERRQ(ierr); 604 ierr = KSPSetOptionsPrefix(ksp,schurprefix);CHKERRQ(ierr); 605 /* 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. */ 606 ierr = KSPSetDM(ksp,jac->head->dm); CHKERRQ(ierr); 607 ierr = KSPSetDMActive(ksp, PETSC_FALSE); CHKERRQ(ierr); 608 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 609 /* Need to call this everytime because new matrix is being created */ 610 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 611 ierr = MatSetUp(jac->schur);CHKERRQ(ierr); 612 /* Create and set up the KSP for the Schur complement; forward the second split's DM and set up tabbing, including for the contained PC. */ 613 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 614 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 615 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 616 { 617 PC pcschur; 618 ierr = KSPGetPC(jac->kspschur, &pcschur); CHKERRQ(ierr); 619 ierr = PetscObjectIncrementTabLevel((PetscObject)pcschur,(PetscObject)pc,1);CHKERRQ(ierr); 620 } 621 /* 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. */ 622 ierr = KSPSetDM(jac->kspschur,ilink->dm); CHKERRQ(ierr); 623 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE); CHKERRQ(ierr); 624 625 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 626 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 627 PC pcschur; 628 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 629 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 630 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 631 } 632 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 633 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 634 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 635 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 636 ierr = KSPSetFromOptions(jac->kspschur); CHKERRQ(ierr); 637 } 638 639 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 640 ierr = PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 641 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 642 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 643 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 644 ierr = PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 645 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 646 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 647 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 648 } 649 else { 650 /* set up the individual splits' PCs */ 651 i = 0; 652 ilink = jac->head; 653 while (ilink) { 654 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag); CHKERRQ(ierr); 655 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 656 if (!jac->suboptionsset) { 657 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 658 } 659 i++; 660 ilink = ilink->next; 661 } 662 } 663 664 jac->suboptionsset = PETSC_TRUE; 665 PetscFunctionReturn(0); 666 } 667 668 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 669 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 670 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 671 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 672 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 673 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "PCApply_FieldSplit_Schur" 677 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 678 { 679 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 680 PetscErrorCode ierr; 681 KSP ksp; 682 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 683 684 PetscFunctionBegin; 685 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 686 687 switch (jac->schurfactorization) { 688 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 689 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 690 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 691 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 692 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 693 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 694 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 695 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 696 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 697 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 698 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 699 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 700 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 701 break; 702 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 703 /* [A00 0; A10 S], suitable for left preconditioning */ 704 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 705 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 706 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 707 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 708 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 709 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 710 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 711 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 712 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 713 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 714 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 715 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 716 break; 717 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 718 /* [A00 A01; 0 S], suitable for right preconditioning */ 719 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 720 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 721 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 722 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 723 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 724 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 725 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 726 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 727 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 728 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 729 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 730 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 731 break; 732 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 733 /* [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 */ 734 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 735 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 736 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 737 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 738 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 739 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 740 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 741 742 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 743 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 744 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 745 746 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 747 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 748 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 749 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 750 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 751 } 752 PetscFunctionReturn(0); 753 } 754 755 #undef __FUNCT__ 756 #define __FUNCT__ "PCApply_FieldSplit" 757 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 758 { 759 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 760 PetscErrorCode ierr; 761 PC_FieldSplitLink ilink = jac->head; 762 PetscInt cnt,bs; 763 764 PetscFunctionBegin; 765 CHKMEMQ; 766 if (jac->type == PC_COMPOSITE_ADDITIVE) { 767 if (jac->defaultsplit) { 768 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 769 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); 770 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 771 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); 772 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 773 while (ilink) { 774 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 775 ilink = ilink->next; 776 } 777 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 778 } else { 779 ierr = VecSet(y,0.0);CHKERRQ(ierr); 780 while (ilink) { 781 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 782 ilink = ilink->next; 783 } 784 } 785 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 786 if (!jac->w1) { 787 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 788 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 789 } 790 ierr = VecSet(y,0.0);CHKERRQ(ierr); 791 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 792 cnt = 1; 793 while (ilink->next) { 794 ilink = ilink->next; 795 /* compute the residual only over the part of the vector needed */ 796 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 797 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 798 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 799 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 800 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 801 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 802 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 803 } 804 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 805 cnt -= 2; 806 while (ilink->previous) { 807 ilink = ilink->previous; 808 /* compute the residual only over the part of the vector needed */ 809 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 810 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 811 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 812 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 813 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 814 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 815 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 816 } 817 } 818 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 819 CHKMEMQ; 820 PetscFunctionReturn(0); 821 } 822 823 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 824 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 825 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 826 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 827 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 828 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 829 830 #undef __FUNCT__ 831 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 832 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 833 { 834 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 835 PetscErrorCode ierr; 836 PC_FieldSplitLink ilink = jac->head; 837 PetscInt bs; 838 839 PetscFunctionBegin; 840 CHKMEMQ; 841 if (jac->type == PC_COMPOSITE_ADDITIVE) { 842 if (jac->defaultsplit) { 843 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 844 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); 845 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 846 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); 847 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 848 while (ilink) { 849 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 850 ilink = ilink->next; 851 } 852 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 853 } else { 854 ierr = VecSet(y,0.0);CHKERRQ(ierr); 855 while (ilink) { 856 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 857 ilink = ilink->next; 858 } 859 } 860 } else { 861 if (!jac->w1) { 862 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 863 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 864 } 865 ierr = VecSet(y,0.0);CHKERRQ(ierr); 866 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 867 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 868 while (ilink->next) { 869 ilink = ilink->next; 870 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 871 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 872 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 873 } 874 while (ilink->previous) { 875 ilink = ilink->previous; 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 } else { 881 while (ilink->next) { /* get to last entry in linked list */ 882 ilink = ilink->next; 883 } 884 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 885 while (ilink->previous) { 886 ilink = ilink->previous; 887 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 888 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 889 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 890 } 891 } 892 } 893 CHKMEMQ; 894 PetscFunctionReturn(0); 895 } 896 897 #undef __FUNCT__ 898 #define __FUNCT__ "PCReset_FieldSplit" 899 static PetscErrorCode PCReset_FieldSplit(PC pc) 900 { 901 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 902 PetscErrorCode ierr; 903 PC_FieldSplitLink ilink = jac->head,next; 904 905 PetscFunctionBegin; 906 while (ilink) { 907 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 908 ierr = DMDestroy(&ilink->dm);CHKERRQ(ierr); 909 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 910 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 911 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 912 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 913 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 914 next = ilink->next; 915 ilink = next; 916 } 917 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 918 if (jac->mat && jac->mat != jac->pmat) { 919 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 920 } else if (jac->mat) { 921 jac->mat = PETSC_NULL; 922 } 923 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 924 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 925 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 926 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 927 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 928 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 929 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 930 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 931 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 932 jac->reset = PETSC_TRUE; 933 PetscFunctionReturn(0); 934 } 935 936 #undef __FUNCT__ 937 #define __FUNCT__ "PCDestroy_FieldSplit" 938 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 939 { 940 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 941 PetscErrorCode ierr; 942 PC_FieldSplitLink ilink = jac->head,next; 943 944 PetscFunctionBegin; 945 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 946 while (ilink) { 947 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 948 next = ilink->next; 949 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 950 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 951 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 952 ierr = PetscFree(ilink);CHKERRQ(ierr); 953 ilink = next; 954 } 955 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 956 ierr = PetscFree(pc->data);CHKERRQ(ierr); 957 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 958 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 959 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 960 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 961 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 962 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",PETSC_NULL);CHKERRQ(ierr); 963 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",PETSC_NULL);CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 967 #undef __FUNCT__ 968 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 969 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 970 { 971 PetscErrorCode ierr; 972 PetscInt bs; 973 PetscBool flg,stokes = PETSC_FALSE; 974 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 975 PCCompositeType ctype; 976 DM ddm; 977 char ddm_name[1025]; 978 979 PetscFunctionBegin; 980 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 981 if(pc->dm) { 982 /* Allow the user to request a decomposition DM by name */ 983 ierr = PetscStrncpy(ddm_name, "", 1024); CHKERRQ(ierr); 984 ierr = PetscOptionsString("-pc_fieldsplit_decomposition", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); 985 if(flg) { 986 ierr = DMCreateFieldDecompositionDM(pc->dm, ddm_name, &ddm); CHKERRQ(ierr); 987 if(!ddm) { 988 SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown field decomposition name %s", ddm_name); 989 } 990 ierr = PetscInfo(pc,"Using field decomposition DM defined using options database\n");CHKERRQ(ierr); 991 ierr = PCSetDM(pc,ddm); CHKERRQ(ierr); 992 } 993 } 994 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 995 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 996 if (flg) { 997 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 998 } 999 1000 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 1001 if (stokes) { 1002 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1003 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1004 } 1005 1006 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1007 if (flg) { 1008 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1009 } 1010 /* Only setup fields once */ 1011 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1012 /* only allow user to set fields from command line if bs is already known. 1013 otherwise user can set them in PCFieldSplitSetDefaults() */ 1014 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1015 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1016 } 1017 if (jac->type == PC_COMPOSITE_SCHUR) { 1018 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1019 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1020 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); 1021 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); 1022 } 1023 ierr = PetscOptionsTail();CHKERRQ(ierr); 1024 PetscFunctionReturn(0); 1025 } 1026 1027 /*------------------------------------------------------------------------------------*/ 1028 1029 EXTERN_C_BEGIN 1030 #undef __FUNCT__ 1031 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1032 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1033 { 1034 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1035 PetscErrorCode ierr; 1036 PC_FieldSplitLink ilink,next = jac->head; 1037 char prefix[128]; 1038 PetscInt i; 1039 1040 PetscFunctionBegin; 1041 if (jac->splitdefined) { 1042 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1043 PetscFunctionReturn(0); 1044 } 1045 for (i=0; i<n; i++) { 1046 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1047 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1048 } 1049 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1050 if (splitname) { 1051 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1052 } else { 1053 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1054 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1055 } 1056 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 1057 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1058 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 1059 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1060 ilink->nfields = n; 1061 ilink->next = PETSC_NULL; 1062 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1063 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1064 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1065 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1066 1067 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1068 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1069 1070 if (!next) { 1071 jac->head = ilink; 1072 ilink->previous = PETSC_NULL; 1073 } else { 1074 while (next->next) { 1075 next = next->next; 1076 } 1077 next->next = ilink; 1078 ilink->previous = next; 1079 } 1080 jac->nsplits++; 1081 PetscFunctionReturn(0); 1082 } 1083 EXTERN_C_END 1084 1085 EXTERN_C_BEGIN 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1088 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1089 { 1090 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1091 PetscErrorCode ierr; 1092 1093 PetscFunctionBegin; 1094 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1095 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1096 (*subksp)[1] = jac->kspschur; 1097 if (n) *n = jac->nsplits; 1098 PetscFunctionReturn(0); 1099 } 1100 EXTERN_C_END 1101 1102 EXTERN_C_BEGIN 1103 #undef __FUNCT__ 1104 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1105 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1106 { 1107 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1108 PetscErrorCode ierr; 1109 PetscInt cnt = 0; 1110 PC_FieldSplitLink ilink = jac->head; 1111 1112 PetscFunctionBegin; 1113 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1114 while (ilink) { 1115 (*subksp)[cnt++] = ilink->ksp; 1116 ilink = ilink->next; 1117 } 1118 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); 1119 if (n) *n = jac->nsplits; 1120 PetscFunctionReturn(0); 1121 } 1122 EXTERN_C_END 1123 1124 EXTERN_C_BEGIN 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1127 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1128 { 1129 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1130 PetscErrorCode ierr; 1131 PC_FieldSplitLink ilink, next = jac->head; 1132 char prefix[128]; 1133 1134 PetscFunctionBegin; 1135 if (jac->splitdefined) { 1136 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1140 if (splitname) { 1141 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1142 } else { 1143 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1144 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1145 } 1146 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1147 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1148 ilink->is = is; 1149 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1150 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1151 ilink->is_col = is; 1152 ilink->next = PETSC_NULL; 1153 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1154 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1155 { 1156 PC ilinkpc; 1157 ierr = KSPGetPC(ilink->ksp, &ilinkpc); CHKERRQ(ierr); 1158 ierr = PetscObjectIncrementTabLevel((PetscObject)ilinkpc,(PetscObject)pc,1);CHKERRQ(ierr); 1159 } 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