Commits

Jed Brown committed a02b86f

Mat: cache eigenvalue estimates for use in Chebyshev

Comments (0)

Files changed (4)

include/petsc-private/matimpl.h

   PetscViewerFormat      viewformatonassembly;
   PetscBool              checksymmetryonassembly,checknullspaceonassembly;
   PetscReal              checksymmetrytol;
-  };
+  PetscInt      neig_estimates_alloc;
+  PetscReal     *eig_estimates;
+};
 
 PETSC_INTERN PetscErrorCode MatAXPY_Basic(Mat,PetscScalar,Mat,MatStructure);
 PETSC_INTERN PetscErrorCode MatAXPY_BasicWithPreallocation(Mat,Mat,PetscScalar,Mat,MatStructure);
 #define MatPreallocateFinalize(dnz,onz) 0;_4_ierr = PetscFree2(dnz,onz);CHKERRQ(_4_ierr);}
 
 
+PETSC_EXTERN PetscErrorCode MatSetEigenvalueEstimates(Mat,PetscInt,const PetscReal*,const PetscReal*);
+PETSC_EXTERN PetscErrorCode MatGetEigenvalueEstimates(Mat,PetscInt*,const PetscReal**,const PetscReal**);
 
 /* Routines unique to particular data structures */
 PETSC_EXTERN PetscErrorCode MatShellGetContext(Mat,void *);

src/ksp/ksp/impls/cheby/cheby.c

   PetscErrorCode ierr;
   PetscInt       n,neig;
   PetscReal      *re,*im,min,max;
+  Mat            Amat;
 
   PetscFunctionBegin;
   ierr = KSPGetIterationNumber(kspest,&n);CHKERRQ(ierr);
   ierr = PetscMalloc2(n,PetscReal,&re,n,PetscReal,&im);CHKERRQ(ierr);
   ierr = KSPComputeEigenvalues(kspest,n,re,im,&neig);CHKERRQ(ierr);
+  ierr = KSPGetOperators(kspest,&Amat,NULL,NULL);CHKERRQ(ierr);
+  ierr = MatSetEigenvalueEstimates(Amat,neig,re,im);CHKERRQ(ierr);
   min  = PETSC_MAX_REAL;
   max  = PETSC_MIN_REAL;
   for (n=0; n<neig; n++) {
   PetscFunctionBegin;
   ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
+  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
 
   if (cheb->kspest && !cheb->estimate_current) {
+    PetscInt neig;
+    const PetscReal *re;
     PetscReal max,min;
-    Vec       X,B;
-    
-    if (hybrid && purification) {
-      X = ksp->vec_sol;
-    } else {
-      X = ksp->work[0];
-    }
 
-    if (cheb->random) {
-      B    = ksp->work[1];
-      ierr = VecSetRandom(B,cheb->random);CHKERRQ(ierr);
+    ierr = MatGetEigenvalueEstimates(Amat,&neig,&re,NULL);CHKERRQ(ierr);
+    if (neig) {
+      min  = PETSC_MAX_REAL;
+      max  = PETSC_MIN_REAL;
+      for (i=0; i<neig; i++) {
+        min = PetscMin(min,re[i]);
+        max = PetscMax(max,re[i]);
+      }
     } else {
-      B = ksp->vec_rhs;
-    }
-    ierr = KSPSolve(cheb->kspest,B,X);CHKERRQ(ierr);
-    if (hybrid) {
-      cheb->its = 0; /* initialize Chebyshev iteration associated to kspest */
-      ierr      = KSPSetInitialGuessNonzero(cheb->kspest,PETSC_TRUE);CHKERRQ(ierr); 
-    } else if (ksp->guess_zero) {
-      ierr = VecZeroEntries(X);CHKERRQ(ierr);
+      Vec       X,B;
+
+      if (hybrid && purification) {
+        X = ksp->vec_sol;
+      } else {
+        X = ksp->work[0];
+      }
+
+      if (cheb->random) {
+        B    = ksp->work[1];
+        ierr = VecSetRandom(B,cheb->random);CHKERRQ(ierr);
+      } else {
+        B = ksp->vec_rhs;
+      }
+      ierr = KSPSolve(cheb->kspest,B,X);CHKERRQ(ierr);
+      if (hybrid) {
+        cheb->its = 0; /* initialize Chebyshev iteration associated to kspest */
+        ierr      = KSPSetInitialGuessNonzero(cheb->kspest,PETSC_TRUE);CHKERRQ(ierr);
+      } else if (ksp->guess_zero) {
+        ierr = VecZeroEntries(X);CHKERRQ(ierr);
+      }
+      ierr = KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);CHKERRQ(ierr);
     }
-    ierr = KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);CHKERRQ(ierr);
 
     cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
     cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
   }
 
   ksp->its = 0;
-  ierr     = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
   maxit    = ksp->max_it;
 
   /* These three point to the three active solutions, we

src/mat/interface/matrix.c

   if ((*A)->ops->destroy) {
     ierr = (*(*A)->ops->destroy)(*A);CHKERRQ(ierr);
   }
+  ierr = PetscFree((*A)->eig_estimates);CHKERRQ(ierr);
   ierr = MatNullSpaceDestroy(&(*A)->nullsp);CHKERRQ(ierr);
   ierr = MatNullSpaceDestroy(&(*A)->nearnullsp);CHKERRQ(ierr);
   ierr = PetscLayoutDestroy(&(*A)->rmap);CHKERRQ(ierr);
   ierr   = PetscLogEventEnd(MAT_TransposeColoringCreate,mat,0,0,0);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
+
+#undef __FUNCT__
+#define __FUNCT__ "MatGetEigenvalueEstimateKey"
+static PetscErrorCode MatGetEigenvalueEstimateKey(Mat A,PetscInt *key)
+{
+  static PetscInt _key = -1;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if (_key == -1) {
+    ierr = PetscObjectComposedDataRegister(&_key);CHKERRQ(ierr);
+  }
+  *key = _key;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSetEigenvalueEstimates"
+/*@
+   MatSetEigenvalueEstimates - set estimated eigenvalues to be used by various analysis components
+
+   Logically Collective
+
+   Input Arguments:
++  A - the matrix
+.  n - number of eigenvalue estimates
+.  re - real parts of n eigenvalue estimates
+-  im - imaginary parts of n eigenvalue estimates (or NULL)
+
+   Level: developer
+
+.seealso: MatGetEigenvalueEstimate(), KSPChebyshevSetEigenvalues()
+@*/
+PetscErrorCode MatSetEigenvalueEstimates(Mat A,PetscInt n,const PetscReal *re,const PetscReal *im)
+{
+  PetscErrorCode ierr;
+  PetscInt key;
+
+  PetscFunctionBegin;
+  if (n > A->neig_estimates_alloc) {
+    ierr = PetscFree(A->eig_estimates);CHKERRQ(ierr);
+    ierr = PetscMalloc(n*2*sizeof(PetscReal),&A->eig_estimates);CHKERRQ(ierr);
+    A->neig_estimates_alloc = n;
+  }
+  ierr = PetscMemcpy(A->eig_estimates,re,n*sizeof(PetscReal));CHKERRQ(ierr);
+  if (im) {
+    ierr = PetscMemcpy(A->eig_estimates+n,im,n*sizeof(PetscReal));CHKERRQ(ierr);
+  } else {
+    ierr = PetscMemzero(A->eig_estimates+n,n*sizeof(PetscReal));CHKERRQ(ierr);
+  }
+  ierr = MatGetEigenvalueEstimateKey(A,&key);CHKERRQ(ierr);
+  ierr = PetscObjectComposedDataSetInt((PetscObject)A,key,n);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatGetEigenvalueEstimates"
+/*@
+   MatGetEigenvalueEstimates - get estimated eigenvalues to be used by various analysis components
+
+   Logically Collective
+
+   Input Arguments:
+.  A - the matrix
+
+   Output Arguments:
++  n - number of eigenvalue estimates (zero if none)
+.  re - real parts of eigenvalue estimates
+-  im - imaginary parts of eigenvalue estimates
+
+   Level: developer
+
+.seealso: MatSetEigenvalueEstimates(), KSPChebyshevGetEigenvalues()
+@*/
+PetscErrorCode MatGetEigenvalueEstimates(Mat A,PetscInt *n,const PetscReal **re,const PetscReal **im)
+{
+  PetscErrorCode ierr;
+  PetscInt key,neig;
+  PetscBool flg;
+
+  PetscFunctionBegin;
+  ierr = MatGetEigenvalueEstimateKey(A,&key);CHKERRQ(ierr);
+  ierr = PetscObjectComposedDataGetInt((PetscObject)A,key,neig,flg);CHKERRQ(ierr);
+  if (flg) {
+    *n = neig;;
+    *re = A->eig_estimates;
+    if (im) *im = A->eig_estimates + neig;
+  } else {
+    *n = 0;
+    *re = NULL;
+    if (im) *im = NULL;
+  }
+  PetscFunctionReturn(0);
+}