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