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