dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

Anders Logg wrote:
> On Wed, Dec 02, 2009 at 08:53:12PM +0100, Marie Rognes wrote:
>   
> > A patch for setting slepc problem types. (Hopefully. Depending on
> > how well 'bzr send' worked.)
> >
> >     
>
> It didn't work. Try again. :-)
>   

Trying manually: see attachment.

--
Marie

# Bazaar merge directive format 2 (Bazaar 0.90)
# revision_id: meg@xxxxxxxxx-20091202194945-arvlq614lkei4419
# target_branch: bzr+ssh://bazaar.launchpad.net/~dolfin-\
#   core/dolfin/main/
# testament_sha1: 674ec601229653a601f274286db5e8439142f0fd
# timestamp: 2009-12-02 21:40:21 +0100
# base_revision_id: logg@xxxxxxxxx-20091202095634-mhvgsb6ooni710jo
# 
# Begin patch
=== modified file 'demo/la/eigensolver/python/demo.py'
--- demo/la/eigensolver/python/demo.py	2009-10-09 09:28:03 +0000
+++ demo/la/eigensolver/python/demo.py	2009-12-02 19:49:45 +0000
@@ -1,4 +1,6 @@
-"This simple program illustrates the use of the SLEPc eigenvalue solver."
+"""
+This program illustrates the use of the SLEPc eigenvalue solver for
+both standard and generalized eigenvalue problems."""
 
 __author__ = "Kristian B. Oelgaard (k.b.oelgaard@xxxxxxxxxx)"
 __date__ = "2007-11-28 -- 2009-10-09"
@@ -6,6 +8,7 @@
 __license__  = "GNU LGPL Version 2.1"
 
 # Modified by Anders Logg, 2008.
+# Modified by Marie Rognes, 2009.
 
 from dolfin import *
 import numpy
@@ -14,7 +17,7 @@
 if not has_la_backend("PETSc"):
     print "DOLFIN has not been configured with PETSc. Exiting."
     exit()
-    
+
 if not has_slepc():
     print "DOLFIN has not been configured with SLEPc. Exiting."
     exit()
@@ -23,23 +26,49 @@
 parameters["linear_algebra_backend"] = "PETSc"
 
 # Build stiftness matrix
-mesh = UnitSquare(64, 64)
+mesh = UnitSquare(4, 4)
 V = FunctionSpace(mesh, "CG", 1)
 v = TestFunction(V)
 u = TrialFunction(V)
 A = PETScMatrix()
 assemble(dot(grad(v), grad(u))*dx, tensor=A)
 
-# Compute the first n eigenvalues
-n = 10
+# Build mass matrix
+M = PETScMatrix()
+assemble(v*u*dx, tensor=M)
+
+# Compute the n largest eigenvalues of A x = \lambda x
+n = 3
 esolver = SLEPcEigenSolver()
-esolver.parameters["spectrum"] = "largest magnitude"
 esolver.solve(A, n)
 
-print "Eigenvalue solver converged in %d iterations" % (esolver.get_iteration_number(),)
-
 # Display eigenvalues
 for i in range(n):
     lr, lc, r_vec, c_vec = esolver.get_eigenpair(i)
     print "Eigenvalue " + str(i) + ": " + str(lr)
     print "Eigenvector " + str(i) + ": " + str(r_vec.array())
+
+# Compute the eigenvalues close to 10 of A x = \lambda M x
+print
+print "Computing eigenvalues of generalized problem"
+esolver = SLEPcEigenSolver()
+esolver.parameters["solver"] = "krylov-schur"
+esolver.parameters["spectral_transform"] = "shift-and-invert"
+esolver.parameters["spectral_shift"] = 10.0
+esolver.solve(A, M, n)
+
+m = esolver.get_number_converged()
+print "Number of converged eigenvalues: %d" % m
+for i in range(m):
+    r, c = esolver.get_eigenvalue(i)
+    print "Real(Eigenvalue) " + str(i) + ": " + str(r)
+
+
+# Test a different eigensolver type
+print
+print "Testing lanczos eigenvalue solver (not converging)"
+esolver = SLEPcEigenSolver()
+esolver.parameters["problem_type"] = "hermitian"
+esolver.parameters["solver"] = "lanczos"
+esolver.solve(M)
+

=== modified file 'dolfin/la/SLEPcEigenSolver.cpp'
--- dolfin/la/SLEPcEigenSolver.cpp	2009-10-12 18:07:27 +0000
+++ dolfin/la/SLEPcEigenSolver.cpp	2009-12-02 19:49:45 +0000
@@ -6,7 +6,7 @@
 // Modified by Marie Rognes, 2009.
 //
 // First added:  2005-08-31
-// Last changed: 2009-10-12
+// Last changed: 2009-12-02
 
 #ifdef HAS_SLEPC
 
@@ -142,9 +142,6 @@
   const uint nn = static_cast<int>(n);
   EPSSetDimensions(eps, nn, PETSC_DECIDE, PETSC_DECIDE);
 
-  // Set algorithm type (Hermitian matrix)
-  //EPSSetProblemType(eps, EPS_NHEP);
-
   // Set parameters from PETSc parameter database
   EPSSetFromOptions(eps);
 
@@ -175,6 +172,7 @@
 //-----------------------------------------------------------------------------
 void SLEPcEigenSolver::read_parameters()
 {
+  set_problem_type(parameters["problem_type"]);
   set_spectrum(parameters["spectrum"]);
   set_solver(parameters["solver"]);
   set_tolerance(parameters["tolerance"], parameters["maximum_iterations"]);
@@ -184,6 +182,23 @@
 }
 //-----------------------------------------------------------------------------
 
+void SLEPcEigenSolver::set_problem_type(std::string type)
+{
+  // Do nothing if default type is specified
+  if (type == "default")
+    return;
+
+  if (type == "hermitian")
+    EPSSetProblemType(eps, EPS_HEP);
+  else if (type == "non_hermitian")
+    EPSSetProblemType(eps, EPS_NHEP);
+  else if (type == "gen_hermitian")
+    EPSSetProblemType(eps, EPS_GHEP);
+  else if (type == "gen_non_hermitian")
+    EPSSetProblemType(eps, EPS_GNHEP);
+  else
+    error("Unknown problem type: \"%s\".", type.c_str());
+}
 
 void SLEPcEigenSolver::set_spectral_transform(std::string transform,
                                               double shift)
@@ -202,7 +217,6 @@
     error("Unknown transform: \"%s\".", transform.c_str());
 }
 
-
 void SLEPcEigenSolver::set_spectrum(std::string spectrum)
 {
   // Do nothing if default type is specified
@@ -235,7 +249,10 @@
   if (solver == "default")
     return;
 
-  // Choose solver
+  // Choose solver.
+
+  // (Note that lanczos will give PETSc error unless problem_type is
+  // set to 'hermitian' or 'gen_hermitian')
   if (solver == "power")
     EPSSetType(eps, EPSPOWER);
   else if (solver == "subspace")

=== modified file 'dolfin/la/SLEPcEigenSolver.h'
--- dolfin/la/SLEPcEigenSolver.h	2009-10-12 18:07:27 +0000
+++ dolfin/la/SLEPcEigenSolver.h	2009-12-02 19:49:45 +0000
@@ -6,7 +6,7 @@
 // Modified by Marie Rognes, 2009.
 //
 // First added:  2005-08-31
-// Last changed: 2009-10-12
+// Last changed: 2009-12-02
 
 #ifndef __SLEPC_EIGEN_SOLVER_H
 #define __SLEPC_EIGEN_SOLVER_H
@@ -112,6 +112,7 @@
     {
       Parameters p("slepc_eigenvalue_solver");
 
+      p.add("problem_type",       "default");
       p.add("spectrum",           "largest magnitude");
       p.add("solver",             "krylov-schur");
       p.add("tolerance",          1e-15);
@@ -130,6 +131,9 @@
     /// Callback for changes in parameter values
     void read_parameters();
 
+    // Set problem type (used for SLEPc internals)
+    void set_problem_type(std::string type);
+
     // Set spectral transform
     void set_spectral_transform(std::string transform, double shift);
 

# Begin bundle
IyBCYXphYXIgcmV2aXNpb24gYnVuZGxlIHY0CiMKQlpoOTFBWSZTWYyT+pcABDJfgFgQWv///37H
Xg6////6YAiJ9692ue7VAAQWVrZGhmbVUMJJSZqMmVPyMjVP0mTymap+iTag8I1GjQAAPUGimj0T
ENTBQDQ00yAAAAAPJADVTw1IaaeoBkaNAAANGgAAAAEiIgmITRqeVPZTap+qeU3pTyag9J6QAAPU
NAcZMmjEMTTAQMCaYIwTE000AGEEkQTQEyYQIGpT2mmKTamT1HpDDRHqAA+CiM8+J6l+EzFcbnxu
boGqhqH/nm6ZGUEkMLB3tYxNRh/rqQO6RG4lBRmU1MSjg9/5aq9YJDKkwGNiBC1gYDKw3Tj4a5Mi
hB0WdGio5kDVCeHyefmwG8RW6YBASkR7mSI9PdQPRBwjy7vPKH6IxtoYAy/Q5jJJhg0+gGZ6Ky8t
G7pptc9O2czh37opWus27rLXbOF9Ic41K41XVEq3lR3jKu3tsqCyY5JvmjL7uR9O828OiOw9O5ht
gKZxyR9J+2auZXVITBWQ4YUSC5xZUc0VSFiTDSwVQAlwtZTod0FPvkasGbsdZs150CUaHEuCVZCd
qSUGnGpiydNdGlGhMOWU66lV2eZurFLVid70NIZj0aMifLu+RaejKXtlhgFmIbEWOzNh2htosr96
WoZd+6sCsaMkeUFU8yiB9Fmw23zZZ7YIQxsHb5kyyAyRJnzkGHveEEPA7nVomIbhkCccDltvzB3P
hu5gzFeyddLfj/UOv3f5qTpJWC6efK6w0HmWW02FxrIEK5BquEYHrMav8OKDa/0NMduxtPQxzlUB
ogGRomPrS5r3H8HC/hjYoL19FiQehnSw8tdsE1QVUtAuRthyvtJQsEFZkovsgeIzGkEkrIqc0Mri
T9oyvLwJMptNlXtd68oy7zq7sVW1A2gf4ZkFltTUNCxWGRGjHt+3GGGA4w52bhGiIHdjwESY0kDm
TnG0nCzIiyUUDiHctWprWRLYmZxtOJZ8aEQ4oEzoTcq61OXEYnERI1hmMKqbyMC9pFEEpnEzcVCh
VaRx5oK1kNa1wxNVsbubkdyvIELbtTGY5ftyK4EhjyQWboEtnbKX1RTyQbkyHwKhDwbcI2sEczlo
4VRLSvs2IiN5MYNDjGIiG9q33RfbV+l4i57h7JjvatMIEiQzDlFxtJSOpDLJBjBrZ5WtigY/4IyB
yi7Vu20IG+ovzb5ZV8GTBMxN5Abkrq5D1au5uxzM+85qY5TYkYF5uOfLC4oYTjrGvWJA0ii8yMKj
fzK4kr98dAZo7ghvJueQ3Km/DeWaPyDEr3mltLzLEnVpmWfXOiCEy3RPFyMAsTKY6E6DQdC5aDqh
XARWVCKTW55ELzcZhd1QXFnH6tvLU1qwwCyRzha+ZqxUy3rLKAUdDXPTaGIoL8z09jS9QeK2NjaR
6efJtpvUrkwqjn19qs+UAqfDpoZJjbf6VC5z+wd4SH7oPE4owp7jzPOP6QUM7s2QXO22ufQFRoDY
Eykx4gf2RiGxFlpS3LkEz2ChsgVwbbrw1wQ4AiAt+eFya1RgluxcpBfaBId5xKHE5TQ6z/p8H7Lv
MyzFfFw+Gn7kEMiCdLWSTxWfwlGsbLclp8XQV8pIhFPRaPRtGMX8W9b6e1BnHoDxjs5Ut5wOrbw3
kxuGJGlAMl/peQhthcxOTzmOuJBwNeKFEXkxLTp0hrWePjmdM8ox2b1bBTKJifg4QbxYj88YS533
1RR4BK5FQ/gNZfSsKiEdHh0et9+fqGbi1GCJioHV1Zpr332Wpd1e03KqqTHUnsTmKLRcQFnlCuct
QvQgpjzFp0Xmu0mcwEW1a1yrsadh6ZTCkfKP495TNIqrGlizWwRitPb6cW9tS2Qn2EdbQZe9g5ce
1vhc5IOaS33g/MBuaT4o4polBQPEKBXhbKETGyac5vlYdSV2XtgVd/Wpv0z6eneZhpap9nREeGng
QFSVVJ4JZJQtRxTMZ9vDbs6kui8ZnUoSDvNnfvvS/wuKmuPj0gRbZj30GklicUyQ/mlyHJo9Pxas
gJS+BzTxc9iv7ZrYgBucih7PIZrmmS6t/qngtADD3k8m+/VIoVD6TPUHdMkz+gdhFWiVus5IjQMw
SbIUHOEsNkA6/WKWki1sDuYVWvXspncKqqxu4XDkBdMqXVe1cTrARzFX4QjyCgrowCJI4oD+ZuIN
3B5sJgCwLe+teIie1rnBZflkIpiFSyggYtBLdbGYac+yCROKaRWRdXnxNAZOWNjJJuisS54WBZaE
qeYuSC/LhNyU2rM1cJjbiVQVROkVg3SBRUWP+mnkgN+rox45shYiaG02i1oIWKeTqWT/ZwB5WG1B
kQL2NzTpkaF3QVnXqQbhmYGZA2rsWyhEY3blTluToymZk8R0itfDjeWFYZ/F8Ox6OMoUX/enxMDE
9v6mD1gmR1CoPUBVLC3ivLYPk1qO1epQZZBks6s5xHDu9/SXJsb7YgeiCGb1r5pSCiBpMwpDnrRL
plJQm02QWV5Z0CZNjN6zUUQZgqlmdQ+4GB9aCyLnxLBF14hr18KThukaI8Tg9YNzawX4X16radYX
WjT0zGY8KoKjbw6f+QGoUrOjf81jfUY4DxoRBIxpPoJobKjbGQ7R8BXRMiCrFVkt8HhQz2Y8ScST
Bm0EwfQVArtKol09KOKURW0dDeyL2pX4Fl19aYVyZewxcLSNbIswi1DJFHYTI/H+SD7uWaD58UBw
uXnuty9dThHK46p40QS0mQTCEwUWguWwjx91dljg1GBLv2Qypj3oGtMYL2gwEHopjFhrWI/4u5Ip
woSEZJ/UuA==