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