[Scipy-svn] r4058 - trunk/scipy/sparse/linalg/eigen/lobpcg
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Mar 31 13:35:37 EDT 2008
Author: wnbell
Date: 2008-03-31 12:35:30 -0500 (Mon, 31 Mar 2008)
New Revision: 4058
Removed:
trunk/scipy/sparse/linalg/eigen/lobpcg/README
Modified:
trunk/scipy/sparse/linalg/eigen/lobpcg/__init__.py
trunk/scipy/sparse/linalg/eigen/lobpcg/info.py
trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
trunk/scipy/sparse/linalg/eigen/lobpcg/setup.py
Log:
updated docstrings
Deleted: trunk/scipy/sparse/linalg/eigen/lobpcg/README
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/README 2008-03-31 14:43:29 UTC (rev 4057)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/README 2008-03-31 17:35:30 UTC (rev 4058)
@@ -1,12 +0,0 @@
-Pure SciPy implementation of Locally Optimal Block Preconditioned Conjugate
-Gradient Method (LOBPCG), see
-http://www-math.cudenver.edu/~aknyazev/software/BLOPEX/
-
-lobpcg.py code was written by Robert Cimrman. Many thanks belong to Andrew
-Knyazev, the author of the algorithm, for lots of advice and support.
-
-The algorithm of LOBPCG is described in detail in:
-
-A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
-
-A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov, Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc (2007). http://arxiv.org/abs/0705.2626
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/__init__.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/__init__.py 2008-03-31 14:43:29 UTC (rev 4057)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/__init__.py 2008-03-31 17:35:30 UTC (rev 4058)
@@ -1,4 +1,4 @@
-"LOBPCG"
+"""LOBPCG eigensolver"""
from info import __doc__
import lobpcg
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/info.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/info.py 2008-03-31 14:43:29 UTC (rev 4057)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/info.py 2008-03-31 17:35:30 UTC (rev 4058)
@@ -1,88 +1,108 @@
"""
-The algorithm of LOBPCG is described in detail in:
+Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG)
-A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal
-Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific
-Computing 23 (2001), no. 2,
-pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
+LOBPCG is a preconditioned eigensolver for large symmetric positive definite
+(SPD) generalized eigenproblems.
-A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov, Block Locally
-Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc
-(2007). http://arxiv.org/abs/0705.2626
-
Call the function lobpcg - see help for lobpcg.lobpcg. See also lobpcg.as2d,
which can be used in the preconditioner (example below)
-Example:
+Acknowledgements
+----------------
- # Solve A x = lambda B x with constraints and preconditioning.
+lobpcg.py code was written by Robert Cimrman. Many thanks belong to Andrew
+Knyazev, the author of the algorithm, for lots of advice and support.
- n = 100
- vals = [nm.arange( n, dtype = nm.float64 ) + 1]
+Examples
+--------
- # Matrix A.
- operatorA = spdiags( vals, 0, n, n )
-
- # Matrix B
- operatorB = nm.eye( n, n )
-
- # Constraints.
- Y = nm.eye( n, 3 )
-
- # Initial guess for eigenvectors, should have linearly independent
- # columns. Column dimension = number of requested eigenvalues.
- X = sc.rand( n, 3 )
-
- # Preconditioner - inverse of A.
- ivals = [1./vals[0]]
- def precond( x ):
+>>> # Solve A x = lambda B x with constraints and preconditioning.
+>>> n = 100
+>>> vals = [nm.arange( n, dtype = nm.float64 ) + 1]
+>>> # Matrix A.
+>>> operatorA = spdiags( vals, 0, n, n )
+>>> # Matrix B
+>>> operatorB = nm.eye( n, n )
+>>> # Constraints.
+>>> Y = nm.eye( n, 3 )
+>>> # Initial guess for eigenvectors, should have linearly independent
+>>> # columns. Column dimension = number of requested eigenvalues.
+>>> X = sc.rand( n, 3 )
+>>> # Preconditioner - inverse of A.
+>>> ivals = [1./vals[0]]
+>>> def precond( x ):
invA = spdiags( ivals, 0, n, n )
y = invA * x
if sp.issparse( y ):
y = y.toarray()
-
+
return as2d( y )
- # Alternative way of providing the same preconditioner.
- #precond = spdiags( ivals, 0, n, n )
- tt = time.clock()
- eigs, vecs = lobpcg( X, operatorA, operatorB, blockVectorY = Y,
- operatorT = precond,
- residualTolerance = 1e-4, maxIterations = 40,
- largest = False, verbosityLevel = 1 )
- print 'solution time:', time.clock() - tt
- print eigs
+>>>
+>>> # Alternative way of providing the same preconditioner.
+>>> #precond = spdiags( ivals, 0, n, n )
+>>>
+>>> tt = time.clock()
+>>> eigs, vecs = lobpcg( X, operatorA, operatorB, blockVectorY = Y,
+>>> operatorT = precond,
+>>> residualTolerance = 1e-4, maxIterations = 40,
+>>> largest = False, verbosityLevel = 1 )
+>>> print 'solution time:', time.clock() - tt
+>>> print eigs
-Usage notes:
-Notation: n - matrix size, m - number of required eigenvalues (smallest or
-largest)
-1) The LOBPCG code internally solves eigenproblems of the size 3m on every
- iteration by calling the"standard" eigensolver, so if m is not small enough
- compared to n, it does not make sense to call the LOBPCG code, but rather
- one should use the "standard" eigensolver, e.g. symeig function in this
- case. If one calls the LOBPCG algorithm for 5m>n, it will most likely break
- internally, so the code tries to call symeig instead.
+Notes
+-----
- It is not that n should be large for the LOBPCG to work, but rather the
- ratio n/m should be large. It you call the LOBPCG code with m=1 and n=10, it
- should work, though n is small. The method is intended for extremely large
- n/m, see e.g., reference [28] in http://arxiv.org/abs/0705.2626
+In the following ``n`` denotes the matrix size and ``m`` the number
+of required eigenvalues (smallest or largest).
-2) The convergence speed depends basically on two factors:
+The LOBPCG code internally solves eigenproblems of the size 3``m`` on every
+iteration by calling the "standard" dense eigensolver, so if ``m`` is not
+small enough compared to ``n``, it does not make sense to call the LOBPCG
+code, but rather one should use the "standard" eigensolver, e.g. symeig
+function in this case. If one calls the LOBPCG algorithm for 5``m``>``n``,
+it will most likely break internally, so the code tries to call symeig
+instead.
- a) how well relatively separated the seeking eigenvalues are from the rest of
- the eigenvalues. One can try to vary m to make this better.
+It is not that n should be large for the LOBPCG to work, but rather the
+ratio ``n``/``m`` should be large. It you call the LOBPCG code with ``m``=1
+and ``n``=10, it should work, though ``n`` is small. The method is intended
+for extremely large ``n``/``m``, see e.g., reference [28] in
+http://arxiv.org/abs/0705.2626
- b) how "well conditioned" the problem is. This can be changed by using proper
- preconditioning. For example, a rod vibration test problem (under tests
- directory) is ill-conditioned
- for large n, so convergence will be slow, unless efficient preconditioning is
- used. For this specific problem, a good simple preconditioner function would
- be a linear solve for A, which is easy to code since A is tridiagonal.
+The convergence speed depends basically on two factors:
+1. How well relatively separated the seeking eigenvalues are from the rest of
+ the eigenvalues. One can try to vary ``m`` to make this better.
+
+2. How well conditioned the problem is. This can be changed by using proper
+ preconditioning. For example, a rod vibration test problem (under tests
+ directory) is ill-conditioned for large ``n``, so convergence will be
+ slow, unless efficient preconditioning is used. For this specific problem,
+ a good simple preconditioner function would be a linear solve for A, which
+ is easy to code since A is tridiagonal.
+
+
+References
+----------
+A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal
+Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific
+Computing 23 (2001), no. 2,
+pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
+
+A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov, Block Locally
+Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc
+(2007). http://arxiv.org/abs/0705.2626
+
+A. V. Knyazev's C and MATLAB implementations:
+http://www-math.cudenver.edu/~aknyazev/software/BLOPEX/
+
"""
+
+__docformat__ = "restructuredtext en"
+
postpone_import = 1
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-03-31 14:43:29 UTC (rev 4057)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-03-31 17:35:30 UTC (rev 4058)
@@ -138,59 +138,61 @@
residualTolerance = None, maxIterations = 20,
largest = True, verbosityLevel = 0,
retLambdaHistory = False, retResidualNormsHistory = False ):
- """
- LOBPCG solves symmetric partial eigenproblems using preconditioning.
+ """LOBPCG solves symmetric partial eigenproblems using preconditioning.
- Required input:
+ TODO write in terms of Ax=lambda B x
- blockVectorX - initial approximation to eigenvectors, full or sparse matrix
- n-by-blockSize
+ Parameters
+ ----------
+ blockVectorX : array_like
+ initial approximation to eigenvectors shape=(n,blockSize)
+ operatorA : {dense matrix, sparse matrix, LinearOperator}
+ the linear operator of the problem, usually a sparse matrix
+ often called the "stiffness matrix"
- operatorA - the operator of the problem, can be given as a matrix or as an
- M-file
+ Returns
+ -------
+ (lambda,blockVectorV) : tuple of arrays
+ blockVectorX and lambda are computed blockSize eigenpairs, where
+ blockSize=size(blockVectorX,2) for the initial guess blockVectorX
+ if it is full rank.
+ Other Parameters
+ ----------------
+ operatorB : {dense matrix, sparse matrix, LinearOperator}
+ the right hand side operator in a generalized eigenproblem.
+ by default, operatorB = Identity
+ often called the "mass matrix"
+ operatorT : {dense matrix, sparse matrix, LinearOperator}
+ preconditioner to operatorA; by default operatorT = Identity
+ operatorT should approximate the inverse of operatorA
+ blockVectorY : array_like
+ n-by-sizeY matrix of constraints, sizeY < n
+ The iterations will be performed in the (operatorB-) orthogonal
+ complement of the column-space of blockVectorY.
+ blockVectorY must be full rank.
- Optional input:
+ residualTolerance : scalar
+ solver tolerance. default: residualTolerance=n*sqrt(eps)
+ maxIterations: integer
+ maximum number of iterations
+ by default: maxIterations=min(n,20)
+ largest : boolean
+ when True, solve for the largest eigenvalues, otherwise the smallest
+ verbosityLevel : integer
+ controls solver output. default: verbosityLevel = 0.
+ retLambdaHistory : boolean
+ whether to return eigenvalue history
+ retResidualNormsHistory : boolean
+ whether to return history of residual norms
- operatorB - the second operator, if solving a generalized eigenproblem; by
- default, or if empty, operatorB = I.
- operatorT - preconditioner; by default, operatorT = I.
-
-
- Optional constraints input:
-
- blockVectorY - n-by-sizeY matrix of constraints, sizeY < n. The iterations
- will be performed in the (operatorB-) orthogonal complement of the
- column-space of blockVectorY. blockVectorY must be full rank.
-
-
- Optional scalar input parameters:
-
- residualTolerance - tolerance, by default, residualTolerance=n*sqrt(eps)
-
- maxIterations - max number of iterations, by default, maxIterations =
- min(n,20)
-
- largest - when true, solve for the largest eigenvalues, otherwise for the
- smallest
-
- verbosityLevel - by default, verbosityLevel = 0.
-
- retLambdaHistory - return eigenvalue history
-
- retResidualNormsHistory - return history of residual norms
-
- Output:
-
- blockVectorX and lambda are computed blockSize eigenpairs, where
- blockSize=size(blockVectorX,2) for the initial guess blockVectorX if it is
- full rank.
-
+ Notes
+ -----
If both retLambdaHistory and retResidualNormsHistory are True, the
- return tuple has the flollowing order:
+ return tuple has the following format:
+ (lambda, blockVectorV, lambda history, residual norms history)
- lambda, blockVectorX, lambda history, residual norms history
"""
failureFlag = True
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/setup.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/setup.py 2008-03-31 14:43:29 UTC (rev 4057)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/setup.py 2008-03-31 17:35:30 UTC (rev 4058)
@@ -1,6 +1,4 @@
#!/usr/bin/env python
-from os.path import join
-import sys
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
More information about the Scipy-svn
mailing list