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