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