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