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