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