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