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