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 /* 583 Do we really want to attach the A11-block's nullspace to S? Note that this may generate 584 "false positives", when the A11's nullspace isn't S's: Stokes or A = [1, 1; 1, 0]. 585 Using a false nullspace may prevent KSP(S) from converging, since it might force an inconsistent rhs. 586 */ 587 /* if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);} */ 588 /* set tabbing, options prefix and DM of KSP inside the MatSchurComplement */ 589 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 590 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 591 { 592 PC pcinner; 593 ierr = KSPGetPC(ksp, &pcinner); CHKERRQ(ierr); 594 ierr = PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);CHKERRQ(ierr); 595 } 596 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",jac->head->splitname);CHKERRQ(ierr); 597 ierr = KSPSetOptionsPrefix(ksp,schurprefix);CHKERRQ(ierr); 598 /* 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. */ 599 ierr = KSPSetDM(ksp,jac->head->dm); CHKERRQ(ierr); 600 ierr = KSPSetDMActive(ksp, PETSC_FALSE); CHKERRQ(ierr); 601 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 602 ierr = KSPSetUp(ksp); CHKERRQ(ierr); 603 /* Need to call this everytime because new matrix is being created */ 604 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 605 ierr = MatSetUp(jac->schur);CHKERRQ(ierr); 606 /* 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. */ 607 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 608 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 609 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 610 { 611 PC pcschur; 612 ierr = KSPGetPC(jac->kspschur, &pcschur); CHKERRQ(ierr); 613 ierr = PetscObjectIncrementTabLevel((PetscObject)pcschur,(PetscObject)pc,1);CHKERRQ(ierr); 614 } 615 /* 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. */ 616 ierr = KSPSetDM(jac->kspschur,ilink->dm); CHKERRQ(ierr); 617 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE); CHKERRQ(ierr); 618 619 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 620 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 621 PC pcschur; 622 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 623 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 624 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 625 } 626 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 627 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 628 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 629 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 630 ierr = KSPSetFromOptions(jac->kspschur); CHKERRQ(ierr); 631 ierr = KSPSetUp(jac->kspschur); CHKERRQ(ierr); 632 } 633 634 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 635 ierr = PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 636 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 637 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 638 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 639 ierr = PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 640 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 641 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 642 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 643 } 644 /* set up the individual PCs */ 645 i = 0; 646 ilink = jac->head; 647 while (ilink) { 648 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag); CHKERRQ(ierr); 649 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 650 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp); CHKERRQ(ierr);} 651 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 652 i++; 653 ilink = ilink->next; 654 } 655 656 jac->suboptionsset = PETSC_TRUE; 657 PetscFunctionReturn(0); 658 } 659 660 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 661 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 662 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 663 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 664 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 665 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 666 667 #undef __FUNCT__ 668 #define __FUNCT__ "PCApply_FieldSplit_Schur" 669 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 670 { 671 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 672 PetscErrorCode ierr; 673 KSP ksp; 674 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 675 676 PetscFunctionBegin; 677 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 678 679 switch (jac->schurfactorization) { 680 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 681 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 682 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 683 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 684 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 685 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 686 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 687 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 688 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 689 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 690 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 691 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 692 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 693 break; 694 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 695 /* [A00 0; A10 S], suitable for left preconditioning */ 696 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 697 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 698 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 699 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 700 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 701 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 702 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 703 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 704 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 705 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 706 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 707 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 708 break; 709 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 710 /* [A00 A01; 0 S], suitable for right preconditioning */ 711 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 712 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 713 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 714 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 715 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 716 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 717 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 718 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 719 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 720 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 721 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 722 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 723 break; 724 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 725 /* [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 */ 726 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 727 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 728 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 729 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 730 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 731 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 732 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 733 734 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 735 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 736 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 737 738 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 739 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 740 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 741 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 742 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 743 } 744 PetscFunctionReturn(0); 745 } 746 747 #undef __FUNCT__ 748 #define __FUNCT__ "PCApply_FieldSplit" 749 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 750 { 751 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 752 PetscErrorCode ierr; 753 PC_FieldSplitLink ilink = jac->head; 754 PetscInt cnt,bs; 755 756 PetscFunctionBegin; 757 CHKMEMQ; 758 if (jac->type == PC_COMPOSITE_ADDITIVE) { 759 if (jac->defaultsplit) { 760 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 761 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); 762 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 763 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); 764 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 765 while (ilink) { 766 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 767 ilink = ilink->next; 768 } 769 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 770 } else { 771 ierr = VecSet(y,0.0);CHKERRQ(ierr); 772 while (ilink) { 773 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 774 ilink = ilink->next; 775 } 776 } 777 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 778 if (!jac->w1) { 779 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 780 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 781 } 782 ierr = VecSet(y,0.0);CHKERRQ(ierr); 783 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 784 cnt = 1; 785 while (ilink->next) { 786 ilink = ilink->next; 787 /* compute the residual only over the part of the vector needed */ 788 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 789 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 790 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 791 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 792 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 793 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 794 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 795 } 796 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 797 cnt -= 2; 798 while (ilink->previous) { 799 ilink = ilink->previous; 800 /* compute the residual only over the part of the vector needed */ 801 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 802 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 803 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 804 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 805 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 806 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 807 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 808 } 809 } 810 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 811 CHKMEMQ; 812 PetscFunctionReturn(0); 813 } 814 815 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 816 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 817 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 818 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 819 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 820 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 824 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 825 { 826 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 827 PetscErrorCode ierr; 828 PC_FieldSplitLink ilink = jac->head; 829 PetscInt bs; 830 831 PetscFunctionBegin; 832 CHKMEMQ; 833 if (jac->type == PC_COMPOSITE_ADDITIVE) { 834 if (jac->defaultsplit) { 835 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 836 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); 837 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 838 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); 839 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 840 while (ilink) { 841 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 842 ilink = ilink->next; 843 } 844 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 845 } else { 846 ierr = VecSet(y,0.0);CHKERRQ(ierr); 847 while (ilink) { 848 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 849 ilink = ilink->next; 850 } 851 } 852 } else { 853 if (!jac->w1) { 854 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 855 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 856 } 857 ierr = VecSet(y,0.0);CHKERRQ(ierr); 858 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 859 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 860 while (ilink->next) { 861 ilink = ilink->next; 862 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 863 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 864 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 865 } 866 while (ilink->previous) { 867 ilink = ilink->previous; 868 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 869 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 870 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 871 } 872 } else { 873 while (ilink->next) { /* get to last entry in linked list */ 874 ilink = ilink->next; 875 } 876 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 877 while (ilink->previous) { 878 ilink = ilink->previous; 879 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 880 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 881 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 882 } 883 } 884 } 885 CHKMEMQ; 886 PetscFunctionReturn(0); 887 } 888 889 #undef __FUNCT__ 890 #define __FUNCT__ "PCReset_FieldSplit" 891 static PetscErrorCode PCReset_FieldSplit(PC pc) 892 { 893 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 894 PetscErrorCode ierr; 895 PC_FieldSplitLink ilink = jac->head,next; 896 897 PetscFunctionBegin; 898 while (ilink) { 899 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 900 ierr = DMDestroy(&ilink->dm);CHKERRQ(ierr); 901 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 902 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 903 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 904 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 905 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 906 next = ilink->next; 907 ilink = next; 908 } 909 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 910 if (jac->mat && jac->mat != jac->pmat) { 911 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 912 } else if (jac->mat) { 913 jac->mat = PETSC_NULL; 914 } 915 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 916 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 917 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 918 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 919 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 920 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 921 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 922 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 923 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 924 jac->reset = PETSC_TRUE; 925 PetscFunctionReturn(0); 926 } 927 928 #undef __FUNCT__ 929 #define __FUNCT__ "PCDestroy_FieldSplit" 930 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 931 { 932 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 933 PetscErrorCode ierr; 934 PC_FieldSplitLink ilink = jac->head,next; 935 936 PetscFunctionBegin; 937 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 938 while (ilink) { 939 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 940 next = ilink->next; 941 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 942 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 943 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 944 ierr = PetscFree(ilink);CHKERRQ(ierr); 945 ilink = next; 946 } 947 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 948 ierr = PetscFree(pc->data);CHKERRQ(ierr); 949 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 950 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 951 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 952 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 953 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 954 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",PETSC_NULL);CHKERRQ(ierr); 955 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",PETSC_NULL);CHKERRQ(ierr); 956 PetscFunctionReturn(0); 957 } 958 959 #undef __FUNCT__ 960 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 961 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 962 { 963 PetscErrorCode ierr; 964 PetscInt bs; 965 PetscBool flg,stokes = PETSC_FALSE; 966 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 967 PCCompositeType ctype; 968 DM ddm; 969 char ddm_name[1025]; 970 971 PetscFunctionBegin; 972 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 973 if(pc->dm) { 974 /* Allow the user to request a decomposition DM by name */ 975 ierr = PetscStrncpy(ddm_name, "", 1024); CHKERRQ(ierr); 976 ierr = PetscOptionsString("-pc_fieldsplit_decomposition", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); 977 if(flg) { 978 ierr = DMCreateFieldDecompositionDM(pc->dm, ddm_name, &ddm); CHKERRQ(ierr); 979 if(!ddm) { 980 SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown field decomposition name %s", ddm_name); 981 } 982 ierr = PetscInfo(pc,"Using field decomposition DM defined using options database\n");CHKERRQ(ierr); 983 ierr = PCSetDM(pc,ddm); CHKERRQ(ierr); 984 } 985 } 986 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 987 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 988 if (flg) { 989 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 990 } 991 992 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 993 if (stokes) { 994 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 995 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 996 } 997 998 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 999 if (flg) { 1000 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1001 } 1002 /* Only setup fields once */ 1003 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1004 /* only allow user to set fields from command line if bs is already known. 1005 otherwise user can set them in PCFieldSplitSetDefaults() */ 1006 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1007 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1008 } 1009 if (jac->type == PC_COMPOSITE_SCHUR) { 1010 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1011 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1012 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); 1013 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); 1014 } 1015 ierr = PetscOptionsTail();CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 /*------------------------------------------------------------------------------------*/ 1020 1021 EXTERN_C_BEGIN 1022 #undef __FUNCT__ 1023 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1024 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1025 { 1026 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1027 PetscErrorCode ierr; 1028 PC_FieldSplitLink ilink,next = jac->head; 1029 char prefix[128]; 1030 PetscInt i; 1031 1032 PetscFunctionBegin; 1033 if (jac->splitdefined) { 1034 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1035 PetscFunctionReturn(0); 1036 } 1037 for (i=0; i<n; i++) { 1038 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1039 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1040 } 1041 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1042 if (splitname) { 1043 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1044 } else { 1045 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1046 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1047 } 1048 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 1049 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1050 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 1051 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1052 ilink->nfields = n; 1053 ilink->next = PETSC_NULL; 1054 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1055 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1056 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1057 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1058 1059 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1060 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1061 1062 if (!next) { 1063 jac->head = ilink; 1064 ilink->previous = PETSC_NULL; 1065 } else { 1066 while (next->next) { 1067 next = next->next; 1068 } 1069 next->next = ilink; 1070 ilink->previous = next; 1071 } 1072 jac->nsplits++; 1073 PetscFunctionReturn(0); 1074 } 1075 EXTERN_C_END 1076 1077 EXTERN_C_BEGIN 1078 #undef __FUNCT__ 1079 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1080 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1081 { 1082 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1083 PetscErrorCode ierr; 1084 1085 PetscFunctionBegin; 1086 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1087 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1088 (*subksp)[1] = jac->kspschur; 1089 if (n) *n = jac->nsplits; 1090 PetscFunctionReturn(0); 1091 } 1092 EXTERN_C_END 1093 1094 EXTERN_C_BEGIN 1095 #undef __FUNCT__ 1096 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1097 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1098 { 1099 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1100 PetscErrorCode ierr; 1101 PetscInt cnt = 0; 1102 PC_FieldSplitLink ilink = jac->head; 1103 1104 PetscFunctionBegin; 1105 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1106 while (ilink) { 1107 (*subksp)[cnt++] = ilink->ksp; 1108 ilink = ilink->next; 1109 } 1110 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); 1111 if (n) *n = jac->nsplits; 1112 PetscFunctionReturn(0); 1113 } 1114 EXTERN_C_END 1115 1116 EXTERN_C_BEGIN 1117 #undef __FUNCT__ 1118 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1119 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1120 { 1121 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1122 PetscErrorCode ierr; 1123 PC_FieldSplitLink ilink, next = jac->head; 1124 char prefix[128]; 1125 1126 PetscFunctionBegin; 1127 if (jac->splitdefined) { 1128 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1132 if (splitname) { 1133 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1134 } else { 1135 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1136 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1137 } 1138 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1139 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1140 ilink->is = is; 1141 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1142 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1143 ilink->is_col = is; 1144 ilink->next = PETSC_NULL; 1145 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1146 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1147 { 1148 PC ilinkpc; 1149 ierr = KSPGetPC(ilink->ksp, &ilinkpc); CHKERRQ(ierr); 1150 ierr = PetscObjectIncrementTabLevel((PetscObject)ilinkpc,(PetscObject)pc,1);CHKERRQ(ierr); 1151 } 1152 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1153 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1154 1155 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1156 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1157 1158 if (!next) { 1159 jac->head = ilink; 1160 ilink->previous = PETSC_NULL; 1161 } else { 1162 while (next->next) { 1163 next = next->next; 1164 } 1165 next->next = ilink; 1166 ilink->previous = next; 1167 } 1168 jac->nsplits++; 1169 1170 PetscFunctionReturn(0); 1171 } 1172 EXTERN_C_END 1173 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "PCFieldSplitSetFields" 1176 /*@ 1177 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1178 1179 Logically Collective on PC 1180 1181 Input Parameters: 1182 + pc - the preconditioner context 1183 . splitname - name of this split, if PETSC_NULL the number of the split is used 1184 . n - the number of fields in this split 1185 - fields - the fields in this split 1186 1187 Level: intermediate 1188 1189 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1190 1191 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1192 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 1193 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1194 where the numbered entries indicate what is in the field. 1195 1196 This function is called once per split (it creates a new split each time). Solve options 1197 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1198 1199 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1200 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1201 available when this routine is called. 1202 1203 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1204 1205 @*/ 1206 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1207 { 1208 PetscErrorCode ierr; 1209 1210 PetscFunctionBegin; 1211 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1212 PetscValidCharPointer(splitname,2); 1213 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1214 PetscValidIntPointer(fields,3); 1215 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1216 PetscFunctionReturn(0); 1217 } 1218 1219 #undef __FUNCT__ 1220 #define __FUNCT__ "PCFieldSplitSetIS" 1221 /*@ 1222 PCFieldSplitSetIS - Sets the exact elements for field 1223 1224 Logically Collective on PC 1225 1226 Input Parameters: 1227 + pc - the preconditioner context 1228 . splitname - name of this split, if PETSC_NULL the number of the split is used 1229 - is - the index set that defines the vector elements in this field 1230 1231 1232 Notes: 1233 Use PCFieldSplitSetFields(), for fields defined by strided types. 1234 1235 This function is called once per split (it creates a new split each time). Solve options 1236 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1237 1238 Level: intermediate 1239 1240 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1241 1242 @*/ 1243 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1244 { 1245 PetscErrorCode ierr; 1246 1247 PetscFunctionBegin; 1248 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1249 if (splitname) PetscValidCharPointer(splitname,2); 1250 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1251 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1252 PetscFunctionReturn(0); 1253 } 1254 1255 #undef __FUNCT__ 1256 #define __FUNCT__ "PCFieldSplitGetIS" 1257 /*@ 1258 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1259 1260 Logically Collective on PC 1261 1262 Input Parameters: 1263 + pc - the preconditioner context 1264 - splitname - name of this split 1265 1266 Output Parameter: 1267 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1268 1269 Level: intermediate 1270 1271 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1272 1273 @*/ 1274 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1275 { 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1280 PetscValidCharPointer(splitname,2); 1281 PetscValidPointer(is,3); 1282 { 1283 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1284 PC_FieldSplitLink ilink = jac->head; 1285 PetscBool found; 1286 1287 *is = PETSC_NULL; 1288 while(ilink) { 1289 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1290 if (found) { 1291 *is = ilink->is; 1292 break; 1293 } 1294 ilink = ilink->next; 1295 } 1296 } 1297 PetscFunctionReturn(0); 1298 } 1299 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1302 /*@ 1303 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1304 fieldsplit preconditioner. If not set the matrix block size is used. 1305 1306 Logically Collective on PC 1307 1308 Input Parameters: 1309 + pc - the preconditioner context 1310 - bs - the block size 1311 1312 Level: intermediate 1313 1314 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1315 1316 @*/ 1317 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1318 { 1319 PetscErrorCode ierr; 1320 1321 PetscFunctionBegin; 1322 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1323 PetscValidLogicalCollectiveInt(pc,bs,2); 1324 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1325 PetscFunctionReturn(0); 1326 } 1327 1328 #undef __FUNCT__ 1329 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1330 /*@C 1331 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1332 1333 Collective on KSP 1334 1335 Input Parameter: 1336 . pc - the preconditioner context 1337 1338 Output Parameters: 1339 + n - the number of splits 1340 - pc - the array of KSP contexts 1341 1342 Note: 1343 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1344 (not the KSP just the array that contains them). 1345 1346 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1347 1348 Level: advanced 1349 1350 .seealso: PCFIELDSPLIT 1351 @*/ 1352 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1353 { 1354 PetscErrorCode ierr; 1355 1356 PetscFunctionBegin; 1357 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1358 if (n) PetscValidIntPointer(n,2); 1359 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1360 PetscFunctionReturn(0); 1361 } 1362 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1365 /*@ 1366 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1367 A11 matrix. Otherwise no preconditioner is used. 1368 1369 Collective on PC 1370 1371 Input Parameters: 1372 + pc - the preconditioner context 1373 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1374 - userpre - matrix to use for preconditioning, or PETSC_NULL 1375 1376 Options Database: 1377 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1378 1379 Notes: 1380 $ If ptype is 1381 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1382 $ to this function). 1383 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1384 $ matrix associated with the Schur complement (i.e. A11) 1385 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1386 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1387 $ preconditioner 1388 1389 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1390 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1391 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1392 1393 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1394 the name since it sets a proceedure to use. 1395 1396 Level: intermediate 1397 1398 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1399 1400 @*/ 1401 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1402 { 1403 PetscErrorCode ierr; 1404 1405 PetscFunctionBegin; 1406 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1407 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1408 PetscFunctionReturn(0); 1409 } 1410 1411 EXTERN_C_BEGIN 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1414 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1415 { 1416 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1417 PetscErrorCode ierr; 1418 1419 PetscFunctionBegin; 1420 jac->schurpre = ptype; 1421 if (pre) { 1422 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1423 jac->schur_user = pre; 1424 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1425 } 1426 PetscFunctionReturn(0); 1427 } 1428 EXTERN_C_END 1429 1430 #undef __FUNCT__ 1431 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1432 /*@ 1433 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1434 1435 Collective on PC 1436 1437 Input Parameters: 1438 + pc - the preconditioner context 1439 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1440 1441 Options Database: 1442 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1443 1444 1445 Level: intermediate 1446 1447 Notes: 1448 The FULL factorization is 1449 1450 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1451 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1452 1453 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, 1454 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). 1455 1456 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 1457 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 1458 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, 1459 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1460 1461 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 1462 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). 1463 1464 References: 1465 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1466 1467 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1468 1469 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1470 @*/ 1471 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1472 { 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1477 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 #undef __FUNCT__ 1482 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1483 PETSC_EXTERN_C PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1484 { 1485 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1486 1487 PetscFunctionBegin; 1488 jac->schurfactorization = ftype; 1489 PetscFunctionReturn(0); 1490 } 1491 1492 #undef __FUNCT__ 1493 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1494 /*@C 1495 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1496 1497 Collective on KSP 1498 1499 Input Parameter: 1500 . pc - the preconditioner context 1501 1502 Output Parameters: 1503 + A00 - the (0,0) block 1504 . A01 - the (0,1) block 1505 . A10 - the (1,0) block 1506 - A11 - the (1,1) block 1507 1508 Level: advanced 1509 1510 .seealso: PCFIELDSPLIT 1511 @*/ 1512 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1513 { 1514 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1515 1516 PetscFunctionBegin; 1517 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1518 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1519 if (A00) *A00 = jac->pmat[0]; 1520 if (A01) *A01 = jac->B; 1521 if (A10) *A10 = jac->C; 1522 if (A11) *A11 = jac->pmat[1]; 1523 PetscFunctionReturn(0); 1524 } 1525 1526 EXTERN_C_BEGIN 1527 #undef __FUNCT__ 1528 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1529 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1530 { 1531 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1532 PetscErrorCode ierr; 1533 1534 PetscFunctionBegin; 1535 jac->type = type; 1536 if (type == PC_COMPOSITE_SCHUR) { 1537 pc->ops->apply = PCApply_FieldSplit_Schur; 1538 pc->ops->view = PCView_FieldSplit_Schur; 1539 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1540 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1541 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","PCFieldSplitSetSchurFactType_FieldSplit",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1542 1543 } else { 1544 pc->ops->apply = PCApply_FieldSplit; 1545 pc->ops->view = PCView_FieldSplit; 1546 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1547 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1548 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",0);CHKERRQ(ierr); 1549 } 1550 PetscFunctionReturn(0); 1551 } 1552 EXTERN_C_END 1553 1554 EXTERN_C_BEGIN 1555 #undef __FUNCT__ 1556 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1557 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1558 { 1559 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1560 1561 PetscFunctionBegin; 1562 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1563 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); 1564 jac->bs = bs; 1565 PetscFunctionReturn(0); 1566 } 1567 EXTERN_C_END 1568 1569 #undef __FUNCT__ 1570 #define __FUNCT__ "PCFieldSplitSetType" 1571 /*@ 1572 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1573 1574 Collective on PC 1575 1576 Input Parameter: 1577 . pc - the preconditioner context 1578 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1579 1580 Options Database Key: 1581 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1582 1583 Level: Intermediate 1584 1585 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1586 1587 .seealso: PCCompositeSetType() 1588 1589 @*/ 1590 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1591 { 1592 PetscErrorCode ierr; 1593 1594 PetscFunctionBegin; 1595 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1596 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1597 PetscFunctionReturn(0); 1598 } 1599 1600 #undef __FUNCT__ 1601 #define __FUNCT__ "PCFieldSplitGetType" 1602 /*@ 1603 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1604 1605 Not collective 1606 1607 Input Parameter: 1608 . pc - the preconditioner context 1609 1610 Output Parameter: 1611 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1612 1613 Level: Intermediate 1614 1615 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1616 .seealso: PCCompositeSetType() 1617 @*/ 1618 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1619 { 1620 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1621 1622 PetscFunctionBegin; 1623 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1624 PetscValidIntPointer(type,2); 1625 *type = jac->type; 1626 PetscFunctionReturn(0); 1627 } 1628 1629 /* -------------------------------------------------------------------------------------*/ 1630 /*MC 1631 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1632 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1633 1634 To set options on the solvers for each block append -fieldsplit_ to all the PC 1635 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1636 1637 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1638 and set the options directly on the resulting KSP object 1639 1640 Level: intermediate 1641 1642 Options Database Keys: 1643 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1644 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1645 been supplied explicitly by -pc_fieldsplit_%d_fields 1646 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1647 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1648 . -pc_fieldsplit_schur_precondition <self,user,diag> - default is diag 1649 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1650 1651 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1652 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1653 1654 Notes: 1655 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1656 to define a field by an arbitrary collection of entries. 1657 1658 If no fields are set the default is used. The fields are defined by entries strided by bs, 1659 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1660 if this is not called the block size defaults to the blocksize of the second matrix passed 1661 to KSPSetOperators()/PCSetOperators(). 1662 1663 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1664 $ ( A10 A11 ) 1665 $ the preconditioner using full factorization is 1666 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1667 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1668 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1669 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1670 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1671 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1672 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1673 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1674 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1675 diag gives 1676 $ ( inv(A00) 0 ) 1677 $ ( 0 -ksp(S) ) 1678 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 1679 $ ( A00 0 ) 1680 $ ( A10 S ) 1681 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1682 $ ( A00 A01 ) 1683 $ ( 0 S ) 1684 where again the inverses of A00 and S are applied using KSPs. 1685 1686 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1687 is used automatically for a second block. 1688 1689 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1690 Generally it should be used with the AIJ format. 1691 1692 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1693 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1694 inside a smoother resulting in "Distributive Smoothers". 1695 1696 Concepts: physics based preconditioners, block preconditioners 1697 1698 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1699 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1700 M*/ 1701 1702 EXTERN_C_BEGIN 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "PCCreate_FieldSplit" 1705 PetscErrorCode PCCreate_FieldSplit(PC pc) 1706 { 1707 PetscErrorCode ierr; 1708 PC_FieldSplit *jac; 1709 1710 PetscFunctionBegin; 1711 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1712 jac->bs = -1; 1713 jac->nsplits = 0; 1714 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1715 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1716 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1717 1718 pc->data = (void*)jac; 1719 1720 pc->ops->apply = PCApply_FieldSplit; 1721 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1722 pc->ops->setup = PCSetUp_FieldSplit; 1723 pc->ops->reset = PCReset_FieldSplit; 1724 pc->ops->destroy = PCDestroy_FieldSplit; 1725 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1726 pc->ops->view = PCView_FieldSplit; 1727 pc->ops->applyrichardson = 0; 1728 1729 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1730 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1732 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1734 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1735 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1736 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1738 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1739 PetscFunctionReturn(0); 1740 } 1741 EXTERN_C_END 1742 1743 1744 1745