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