Commits

Matt Knepley  committed 296d98f Merge

Merge branch 'knepley/feature-schurm-explicit-op'

* knepley/feature-schurm-explicit-op:
KSP: Fix for name change in another branch
KSP: It looks like someone forgot this declaration
Mat: Improve conversion of MATDENSE to MATAIJ - Previously, we were creating a fully dense AIJ matrix
Section: Small doc fix
PC: Allow FieldSplit to explicitly compute S for preconditioning - Yes I know it is generally slow
KSP: Added MatSchurComplementComputeExplicitOperator - This could be optimized to compute a sparse S for preconditioning (Saad has some papers here)

  • Participants
  • Parent commits e90cffd, 79bd946

Comments (0)

Files changed (7)

File include/petsc-private/kspimpl.h

 
 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;
 
-PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*);
+PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
 
 #endif

File include/petscksp.h

 PETSC_EXTERN PetscErrorCode MatSchurComplementSet(Mat,Mat,Mat,Mat,Mat,Mat);
 PETSC_EXTERN PetscErrorCode MatSchurComplementUpdate(Mat,Mat,Mat,Mat,Mat,Mat,MatStructure);
 PETSC_EXTERN PetscErrorCode MatSchurComplementGetSubmatrices(Mat,Mat*,Mat*,Mat*,Mat*,Mat*);
+PETSC_EXTERN PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat,Mat*);
 PETSC_EXTERN PetscErrorCode MatGetSchurComplement(Mat,IS,IS,IS,IS,MatReuse,Mat *,MatReuse,Mat *);
 
 PETSC_EXTERN PetscErrorCode KSPSetDM(KSP,DM);

File include/petscpc.h

 
 .seealso: PCFieldSplitSchurPrecondition()
 E*/
-typedef enum {PC_FIELDSPLIT_SCHUR_PRE_SELF,PC_FIELDSPLIT_SCHUR_PRE_A11,PC_FIELDSPLIT_SCHUR_PRE_USER} PCFieldSplitSchurPreType;
+typedef enum {PC_FIELDSPLIT_SCHUR_PRE_SELF,PC_FIELDSPLIT_SCHUR_PRE_A11,PC_FIELDSPLIT_SCHUR_PRE_USER,PC_FIELDSPLIT_SCHUR_PRE_FULL} PCFieldSplitSchurPreType;
 PETSC_EXTERN const char *const PCFieldSplitSchurPreTypes[];
 
 /*E

File src/ksp/ksp/utils/schurm.c

   PetscFunctionReturn(0);
 }
 
+#include <petsc-private/kspimpl.h>
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSchurComplementComputeExplicitOperator"
+/*@
+  MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
+
+  Collective on Mat
+
+  Input Parameter:
+. M - the matrix obtained with MatCreateSchurComplement()
+
+  Output Parameter:
+. S - the Schur complement matrix
+
+  Note: This can be expensive, so it is mainly for testing
+
+  Level: advanced
+
+.seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
+@*/
+PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
+{
+  Mat            B, C, D;
+  KSP            ksp;
+  PC             pc;
+  PetscBool      isLU, isILU;
+  PetscReal      fill = 2.0;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
+  ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
+  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
+  if (isLU || isILU) {
+    Mat       fact, Bd, AinvB, AinvBd;
+    PetscReal eps = 1.0e-10;
+
+    /* This can be sped up for banded LU */
+    ierr = KSPSetUp(ksp);CHKERRQ(ierr);
+    ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
+    ierr = MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
+    ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
+    ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
+    ierr = MatDestroy(&Bd);CHKERRQ(ierr);
+    ierr = MatChop(AinvBd, eps);CHKERRQ(ierr);
+    ierr = MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);CHKERRQ(ierr);
+    ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
+    ierr = MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
+    ierr = MatDestroy(&AinvB);CHKERRQ(ierr);
+  } else {
+    Mat Ainvd, Ainv;
+
+    ierr = PCComputeExplicitOperator(pc, &Ainvd);CHKERRQ(ierr);
+    ierr = MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);CHKERRQ(ierr);
+    ierr = MatDestroy(&Ainvd);CHKERRQ(ierr);
+#if 0
+    /* Symmetric version */
+    ierr = MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
+#else
+    /* Nonsymmetric version */
+    ierr = MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
+#endif
+    ierr = MatDestroy(&Ainv);CHKERRQ(ierr);
+  }
+
+  ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
+  ierr = MatView(*S, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+  ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  if (D) {
+    MatInfo info;
+
+    ierr = MatGetInfo(D, MAT_GLOBAL_SUM, &info);CHKERRQ(ierr);
+    if (info.nz_used) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented");
+  }
+  PetscFunctionReturn(0);
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "MatGetSchurComplement_Basic"
 /* Developer Notes: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */

File src/ksp/pc/impls/fieldsplit/fieldsplit.c

        http://chess.cs.umd.edu/~elman/papers/tax.pdf
 */
 
-const char *const PCFieldSplitSchurPreTypes[] = {"SELF","A11","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
+const char *const PCFieldSplitSchurPreTypes[] = {"SELF","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0};
 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0};
 
 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
   switch (jac->schurpre) {
   case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
   case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
-  case PC_FIELDSPLIT_SCHUR_PRE_USER:   /* Use a user-provided matrix if it is given, otherwise diagonal block */
+  case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
+  case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
   default:
     return jac->schur_user ? jac->schur_user : jac->pmat[1];
   }
       ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);break;
     case PC_FIELDSPLIT_SCHUR_PRE_A11:
       ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);break;
+    case PC_FIELDSPLIT_SCHUR_PRE_FULL:
+      ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr);break;
     case PC_FIELDSPLIT_SCHUR_PRE_USER:
       if (jac->schur_user) {
         ierr = PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr);
       ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr);
       ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr);
       ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr);
-      ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
         PC pcschur;
         ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr);
         ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr);
         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
+      } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
+        ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr);
       }
+      ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
       ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr);
       ierr = KSPSetOptionsPrefix(jac->kspschur,         Dprefix);CHKERRQ(ierr);
       /* propogate DM */
 -   userpre - matrix to use for preconditioning, or NULL
 
     Options Database:
-.     -pc_fieldsplit_schur_precondition <self,user,a11> default is a11
+.     -pc_fieldsplit_schur_precondition <self,user,a11,full> default is a11
 
     Notes:
 $    If ptype is
 $             to this function).
 $        a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original
 $             matrix associated with the Schur complement (i.e. A11)
+$        full then the preconditioner uses the exact Schur complement (this is expensive)
 $        self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
 $             The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
 $             preconditioner
                               been supplied explicitly by -pc_fieldsplit_%d_fields
 .   -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
 .   -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting
-.   -pc_fieldsplit_schur_precondition <self,user,a11> - default is a11
+.   -pc_fieldsplit_schur_precondition <self,user,a11,full> - default is a11
 .   -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver
 
 -    Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_

File src/mat/impls/dense/seq/dense.c

   Mat            B;
   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
   PetscErrorCode ierr;
-  PetscInt       i;
-  PetscInt       *rows;
-  MatScalar      *aa = a->v;
+  PetscInt       i, j;
+  PetscInt       *rows, *nnz;
+  MatScalar      *aa = a->v, *vals;
 
   PetscFunctionBegin;
   ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
   ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr);
-  ierr = MatSeqAIJSetPreallocation(B,A->cmap->n,NULL);CHKERRQ(ierr);
-
-  ierr = PetscMalloc1(A->rmap->n,&rows);CHKERRQ(ierr);
-  for (i=0; i<A->rmap->n; i++) rows[i] = i;
-
-  for (i=0; i<A->cmap->n; i++) {
-    ierr = MatSetValues(B,A->rmap->n,rows,1,&i,aa,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr);
+  for (j=0; j<A->cmap->n; j++) {
+    for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
+    aa += a->lda;
+  }
+  ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr);
+  aa = a->v;
+  for (j=0; j<A->cmap->n; j++) {
+    PetscInt numRows = 0;
+    for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
+    ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr);
     aa  += a->lda;
   }
-  ierr = PetscFree(rows);CHKERRQ(ierr);
+  ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr);
   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 

File src/vec/is/utils/vsectionis.c

   Level: developer
 
   Notes: Typical calling sequence
-       PetscSectionCreate(MPI_Comm,PetscSection *);
-       PetscSectionSetChart(PetscSection,low,high);
-       PetscSectionSetDof(PetscSection,point,numdof);
-       PetscSectionSetUp(PetscSection);
-       PetscSectionGetOffset(PetscSection,point,PetscInt *);
-       PetscSectionDestroy(PetscSection);
+$       PetscSectionCreate(MPI_Comm,PetscSection *);
+$       PetscSectionSetNumFields(PetscSection, numFields);
+$       PetscSectionSetChart(PetscSection,low,high);
+$       PetscSectionSetDof(PetscSection,point,numdof);
+$       PetscSectionSetUp(PetscSection);
+$       PetscSectionGetOffset(PetscSection,point,PetscInt *);
+$       PetscSectionDestroy(PetscSection);
 
        The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
        recommended they not be used in user codes unless you really gain something in their use.