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