diff -Nru krypy-0.1ubuntu1/.travis.yml krypy-0.1.1/.travis.yml --- krypy-0.1ubuntu1/.travis.yml 2013-02-05 11:23:26.000000000 +0000 +++ krypy-0.1.1/.travis.yml 2013-02-07 15:37:36.000000000 +0000 @@ -2,14 +2,9 @@ python: - "2.7" - "3.2" +virtualenv: + system_site_packages: true before_install: - # install numpy for python 3.2 (currently not installed in travis) - - "if [[ $TRAVIS_PYTHON_VERSION == '3.2' ]]; then sudo apt-get install python3-numpy 2>&1 | tail -n 2; fi" - - "if [[ $TRAVIS_PYTHON_VERSION == '3.2' ]]; then ln -s /usr/lib/python3/dist-packages/numpy ~/virtualenv/python3.2/lib/python3.2/site-packages/; fi" - # install scipy for python 3.2 - - "if [[ $TRAVIS_PYTHON_VERSION == '3.2' ]]; then sudo apt-get install python3-scipy 2>&1 | tail -n 2; fi" - - "if [[ $TRAVIS_PYTHON_VERSION == '3.2' ]]; then ln -s /usr/lib/python3/dist-packages/scipy ~/virtualenv/python3.2/lib/python3.2/site-packages/; fi" - # install scipy for python 2.7 - - "if [[ $TRAVIS_PYTHON_VERSION == '2.7' ]]; then sudo apt-get install python-scipy 2>&1 | tail -n 2; fi" - - "if [[ $TRAVIS_PYTHON_VERSION == '2.7' ]]; then ln -s /usr/lib/python2.7/dist-packages/scipy ~/virtualenv/python2.7/lib/python2.7/site-packages/; fi" -script: "python -m krypy.tests.tests" + - sudo apt-get update -qq + - sudo apt-get install -qq python-numpy python-scipy python3-numpy python3-scipy +script: python -m krypy.tests.tests diff -Nru krypy-0.1ubuntu1/README.md krypy-0.1.1/README.md --- krypy-0.1ubuntu1/README.md 2013-02-05 11:23:26.000000000 +0000 +++ krypy-0.1.1/README.md 2013-02-07 15:37:36.000000000 +0000 @@ -1,9 +1,13 @@ # KryPy -KryPy is a Python 3 module for Krylov subspace methods for the solution of linear algebraic systems. This includes enhanced versions of CG, MINRES and GMRES as well as methods for the efficient solution of sequences of linear systems. +KryPy is a Python (versions 2 and 3) module for Krylov subspace methods for the solution of linear algebraic systems. This includes enhanced versions of CG, MINRES and GMRES as well as methods for the efficient solution of sequences of linear systems. -# Dependencies +# Installing +### Ubuntu +There's an [Ubuntu PPA](https://launchpad.net/~andrenarchy/+archive/python) with a package of the Python 2 version. +### Installing from source +KryPy has the following dependencies: * NumPy * SciPy diff -Nru krypy-0.1ubuntu1/README.txt krypy-0.1.1/README.txt --- krypy-0.1ubuntu1/README.txt 2013-02-05 11:23:26.000000000 +0000 +++ krypy-0.1.1/README.txt 2013-02-07 15:37:36.000000000 +0000 @@ -1,9 +1,13 @@ # KryPy -KryPy is a Python 3 module for Krylov subspace methods for the solution of linear algebraic systems. This includes enhanced versions of CG, MINRES and GMRES as well as methods for the efficient solution of sequences of linear systems. +KryPy is a Python (versions 2 and 3) module for Krylov subspace methods for the solution of linear algebraic systems. This includes enhanced versions of CG, MINRES and GMRES as well as methods for the efficient solution of sequences of linear systems. -# Dependencies +# Installing +### Ubuntu +There's an [Ubuntu PPA](https://launchpad.net/~andrenarchy/+archive/python) with a package of the Python 2 version. +### Installing from source +KryPy has the following dependencies: * NumPy * SciPy diff -Nru krypy-0.1ubuntu1/debian/changelog krypy-0.1.1/debian/changelog --- krypy-0.1ubuntu1/debian/changelog 2013-02-05 17:23:11.000000000 +0000 +++ krypy-0.1.1/debian/changelog 2013-02-08 16:48:48.000000000 +0000 @@ -1,3 +1,9 @@ +krypy (0.1.1) quantal; urgency=low + + * new upstream version + + -- André Gaul Fri, 08 Feb 2013 17:48:23 +0100 + krypy (0.1ubuntu1) quantal; urgency=low * fix dependency diff -Nru krypy-0.1ubuntu1/krypy/linsys.py krypy-0.1.1/krypy/linsys.py --- krypy-0.1ubuntu1/krypy/linsys.py 2013-02-05 11:23:26.000000000 +0000 +++ krypy-0.1.1/krypy/linsys.py 2013-02-07 15:37:36.000000000 +0000 @@ -193,36 +193,38 @@ # initial relative residual norm relresvec = [norm_MMlr0 / norm_MMlb] + xk = x0.copy() + info = 0 # compute error? if exact_solution is not None: errvec = [utils.norm(exact_solution - x0, inner_product = inner_product)] - # Allocate and initialize the 'large' memory blocks. - if return_basis or full_reortho: - Vfull = numpy.c_[MMlr0 / norm_MMlr0, numpy.zeros([N,maxiter], dtype=cdtype)] - Pfull = numpy.c_[Mlr0 / norm_MMlr0, numpy.zeros([N,maxiter], dtype=cdtype)] - Hfull = numpy.zeros((maxiter+1,maxiter), dtype=numpy.float) - # Last and current Lanczos vector: - V = numpy.c_[numpy.zeros(N, dtype=cdtype), MMlr0 / norm_MMlr0] - # M*v[i] = P[1], M*v[i-1] = P[0] - P = numpy.c_[numpy.zeros(N, dtype=cdtype), Mlr0 / norm_MMlr0] - # Necessary for efficient update of yk: - W = numpy.c_[numpy.zeros(N, dtype=cdtype), numpy.zeros(N)] - # some small helpers - ts = 0.0 # (non-existing) first off-diagonal entry (corresponds to pi1) - y = [norm_MMlr0, 0] # first entry is (updated) residual - G2 = numpy.eye(2) # old givens rotation - G1 = numpy.eye(2) # even older givens rotation ;) - k = 0 + # initialize only if needed + if relresvec[-1] > tol: + # Allocate and initialize the 'large' memory blocks. + if return_basis or full_reortho: + Vfull = numpy.c_[MMlr0 / norm_MMlr0, numpy.zeros([N,maxiter], dtype=cdtype)] + Pfull = numpy.c_[Mlr0 / norm_MMlr0, numpy.zeros([N,maxiter], dtype=cdtype)] + Hfull = numpy.zeros((maxiter+1,maxiter), dtype=numpy.float) + # Last and current Lanczos vector: + V = numpy.c_[numpy.zeros(N, dtype=cdtype), MMlr0 / norm_MMlr0] + # M*v[i] = P[1], M*v[i-1] = P[0] + P = numpy.c_[numpy.zeros(N, dtype=cdtype), Mlr0 / norm_MMlr0] + # Necessary for efficient update of yk: + W = numpy.c_[numpy.zeros(N, dtype=cdtype), numpy.zeros(N)] + # some small helpers + ts = 0.0 # (non-existing) first off-diagonal entry (corresponds to pi1) + y = [norm_MMlr0, 0] # first entry is (updated) residual + G2 = numpy.eye(2) # old givens rotation + G1 = numpy.eye(2) # even older givens rotation ;) + k = 0 - # resulting approximation is xk = x0 + Mr*yk - yk = numpy.zeros((N,1), dtype=cdtype) - xk = x0.copy() - info = 0 + # resulting approximation is xk = x0 + Mr*yk + yk = numpy.zeros((N,1), dtype=cdtype) - if timer: - times['setup'][0] = time.time()-start + if timer: + times['setup'][0] = time.time()-start # -------------------------------------------------------------------------- # Lanczos + MINRES iteration @@ -433,10 +435,11 @@ # -------------------------------------------------------------------------- def _compute_explicit_xk(H, V, y): '''Compute approximation xk to the solution.''' - yy = numpy.linalg.solve(H, y) - u = utils.apply(Mr, numpy.dot(V, yy)) - xk = x0 + u - return xk + if (H.shape[0]>0): + yy = numpy.linalg.solve(H, y) + u = utils.apply(Mr, numpy.dot(V, yy)) + return x0+u + return x0 # -------------------------------------------------------------------------- def _compute_explicit_residual( xk ): '''Compute residual explicitly.''' @@ -458,6 +461,7 @@ # get memory for working variables V = numpy.zeros([N, maxiter+1], dtype=cdtype) # Arnoldi basis H = numpy.zeros([maxiter+1, maxiter], dtype=cdtype) # Hessenberg matrix + y = numpy.zeros( (maxiter+1,1), dtype=cdtype ) if M is not None: P = numpy.zeros([N,maxiter+1], dtype=cdtype) # V=M*P @@ -491,14 +495,15 @@ if exact_solution is not None: errvec = [utils.norm(exact_solution - x0, inner_product = inner_product)] - V[:, [0]] = MMlr0 / norm_MMlr0 - if M is not None: - P[:, [0]] = Mlr0 / norm_MMlr0 - # Right hand side of projected system: - y = numpy.zeros( (maxiter+1,1), dtype=cdtype ) - y[0] = norm_MMlr0 - # Givens rotations: - G = [] + # initialize only if needed + if relresvec[-1] > tol: + V[:, [0]] = MMlr0 / norm_MMlr0 + if M is not None: + P[:, [0]] = Mlr0 / norm_MMlr0 + # Right hand side of projected system: + y[0] = norm_MMlr0 + # Givens rotations: + G = [] k = 0 while relresvec[-1] > tol and k < maxiter: @@ -514,10 +519,8 @@ H[i, k] += inner_product(V[:, [i]], z)[0,0] z -= H[i, k] * V[:, [i]] Mz = utils.apply(M, z); - H[k+1, k] = utils.norm(z, Mz, inner_product=inner_product) - if M is not None: - P[:, [k+1]] = z / H[k+1, k] - V[:, [k+1]] = Mz / H[k+1, k] + norm_Mz = utils.norm(z, Mz, inner_product=inner_product) + H[k+1,k] = norm_Mz if return_basis: Horig[0:k+2, [k]] = H[0:k+2, [k]] @@ -559,9 +562,16 @@ % (k+1, relresvec[-1], tol, norm_ur)) info = 1 else: - warnigs.warn('Iter %d: Expl. res = %e >= tol = %e > upd. res = %e.' \ + warnings.warn('Iter %d: Expl. res = %e >= tol = %e > upd. res = %e.' \ % (k+1, relresvec[-1], tol, norm_ur)) + if relresvec[-1] > tol: + if norm_Mz < 1e-14: + warnings.warn('subdiagonal element is (close to) zero (%e) => breakdown in iteration %d' % (norm_Mz, k)) + if M is not None: + P[:, [k+1]] = z / norm_Mz + V[:, [k+1]] = Mz / norm_Mz + k += 1 xk = _compute_explicit_xk(H[:k,:k], V[:,:k], y[:k]) diff -Nru krypy-0.1ubuntu1/krypy/tests/tests.py krypy-0.1.1/krypy/tests/tests.py --- krypy-0.1ubuntu1/krypy/tests/tests.py 2013-02-05 11:23:26.000000000 +0000 +++ krypy-0.1.1/krypy/tests/tests.py 2013-02-07 15:37:36.000000000 +0000 @@ -16,52 +16,160 @@ for p in itertools.product(*d.values()): yield dict(zip(d.keys(), p)) -def get_spd_matrix(): +def get_operators(A): + return [ A, + LinearOperator(A.shape, lambda x: numpy.dot(A,x), + dtype=numpy.double), + csr_matrix(A) + ] +def get_vecs(v): + return [ v, numpy.reshape(v, (v.shape[0],)) ] + +def test_linsys_spd(): + # build spd diag matrix a = numpy.array(range(1,11)) - a[-1] = 1e2 - return numpy.diag(a), numpy.diag(1/a) + a[-1] = 1.e2 + A, Ainv = numpy.diag(a), numpy.diag(1./a) + + # solution + x = numpy.ones((10,1)) -def get_spd_precon(): + # preconditioner + m = numpy.array(range(1,11)) + m[-1] = 1. + M, Minv = numpy.diag(m), numpy.diag(1./m) + params_adds = [ + { 'M': [ None, Ainv ] + get_operators(Minv) }, + { 'Ml': [ None, Ainv ] + get_operators(Minv) }, + { 'Mr': [ None, Ainv ] + get_operators(Minv) } + ] + solvers = [ krypy.linsys.cg, krypy.linsys.minres, krypy.linsys.gmres ] + for case in produce_cases(A, x, params_adds, solvers): + yield case + +def test_linsys_symm_indef(): + # build symm indef diag matrix a = numpy.array(range(1,11)) - a[-1] = 1 - return numpy.diag(a), numpy.diag(1/a) + a[-1] = -1e2 + A, Ainv = numpy.diag(a), numpy.diag(1./a) -def test_linsys_spd(): - A, Ainv = get_spd_matrix() + # solution x = numpy.ones((10,1)) - b = numpy.dot(A,x) - M, Minv = get_spd_precon() - params = { - 'A': [ A, - LinearOperator(A.shape, lambda x: numpy.dot(A,x), - dtype=numpy.double), - csr_matrix(A) - ], - 'b': [ b, - numpy.reshape(b, (b.shape[0],))], - 'x0': [None, - numpy.zeros((10,1)), - numpy.ones((10,1)), - x], - 'tol': [1e-14, 1e-5, 1e-2], - 'maxiter': [15], - 'M': [ None, - Minv, - LinearOperator(Minv.shape, lambda x: numpy.dot(Minv,x), - dtype=numpy.double), - csr_matrix(Minv) - ], - 'exact_solution': [None, x] + + # preconditioner + m = numpy.array(range(1,11)) + m[-1] = 1. + M, Minv = numpy.diag(m), numpy.diag(1./m) + params_adds = [ + { 'M': [ None ] + get_operators(Minv) }, + { 'Ml': [ None, Ainv ] + get_operators(Minv) }, + { 'Mr': [ None, Ainv ] + get_operators(Minv) } + ] + solvers = [ krypy.linsys.minres, krypy.linsys.gmres ] + for case in produce_cases(A, x, params_adds, solvers): + yield case + +def test_linsys_nonsymm(): + # build nonsymm matrix + a = numpy.array(range(1,11)) + a[-1] = -1e2 + A = numpy.diag(a) + A[0,-1] = 1e2 + Ainv = numpy.linalg.inv(A) + + # solution + x = numpy.ones((10,1)) + + # preconditioner + m = numpy.array(range(1,11)) + m[-1] = 1. + M, Minv = numpy.diag(m), numpy.diag(1./m) + params_adds = [ + { 'M': [ None ] + get_operators(Minv) }, + { 'Ml': [ None, Ainv ] + get_operators(Minv) }, + { 'Mr': [ None, Ainv ] + get_operators(Minv) } + ] + solvers = [ krypy.linsys.gmres ] + for case in produce_cases(A, x, params_adds, solvers): + yield case + +def produce_cases(A, x, params_adds, solvers): + b = numpy.dot(A, x) + params_base = { + 'A': get_operators(A), + 'b': get_vecs(b), + 'x0': [ None, numpy.zeros(b.shape), x ] + get_vecs(numpy.ones(b.shape)), + 'tol': [ 1e-14, 1e-5, 1e-2 ], + 'maxiter': [ 15 ], + 'M': [ None ], + 'Ml': [ None ], + 'Mr': [ None ], + 'inner_product': [ krypy.utils.ip ], + 'exact_solution': [ None ] + get_vecs(x) } - for param in dictproduct(params): - yield linsys_cg, param -def linsys_cg(params): - ret = krypy.linsys.cg(**params) + for params_add in params_adds: + params = dict( list(params_base.items()) + list(params_add.items())) + for solver, param in itertools.product(solvers, dictproduct(params)): + yield run_case, solver, param + +def run_case(solver, params): + ret = solver(**params) + + # pick out the interesting data + A = params['A'] + b = krypy.utils.shape_vec(params['b']) + xk = krypy.utils.shape_vec(ret['xk']) + + # maxiter respected? assert( len(ret['relresvec'])-1 <= params['maxiter'] ) + + # tolerance reached (if not near machine eps)? if params['tol']>1e-15: assert( ret['relresvec'][-1] <= params['tol'] ) + # final residual norm correct? + # relresvec[-1] == ||M*Ml*(b-A*xk))||_{M^{-1}} / ||M*Ml*b||_{M^{-1}} + # compute residual norm + rk = b - krypy.utils.apply( A, xk) + Mlrk = krypy.utils.apply( params['Ml'], rk ) + MMlrk = krypy.utils.apply( params['M'], Mlrk ) + norm_MMlrk = krypy.utils.norm( Mlrk, MMlrk, inner_product=params['inner_product'] ) + # compute rhs norm + Mlb = krypy.utils.apply( params['Ml'], b ) + MMlb = krypy.utils.apply( params['M'], Mlb ) + norm_MMlb = krypy.utils.norm( Mlb, MMlb, inner_product=params['inner_product'] ) + # finally: the assertion + assert( abs(ret['relresvec'][-1] - norm_MMlrk/norm_MMlb) <= 1e-15 ) + + # final error norm correct? + # (if exact_solution was provided) + if params['exact_solution'] is not None: + assert( abs(ret['errvec'][-1] - + krypy.utils.norm( krypy.utils.shape_vec(params['exact_solution']) + - krypy.utils.shape_vec(ret['xk']))) <= 1e-15 ) + + # if the preconditioner is the inverse, then check if convergence + # occured after the first iteration + if isinstance(A, numpy.ndarray) and \ + isinstance(params['M'], numpy.ndarray) and \ + numpy.linalg.norm( numpy.eye(*A.shape)- numpy.dot(A, params['M']) ) < 1e-15: + assert( len(ret['relresvec'])<=2 ) + + # 0 iterations if initial guess was good enough? + if params['x0'] is not None: + r0 = b - krypy.utils.apply( A, krypy.utils.shape_vec(params['x0'])) + Mlr0 = krypy.utils.apply( params['Ml'], r0 ) + MMlr0 = krypy.utils.apply( params['M'], Mlr0 ) + norm_MMlr0 = krypy.utils.norm( Mlr0, MMlr0, inner_product=params['inner_product'] ) + if norm_MMlr0/norm_MMlb < params['tol']: + assert( len(ret['relresvec'])==1 ) + + # has gmres found the solution after max N iterations? + # (cg or minres may take longer because of roundoff errors) + if solver==krypy.linsys.gmres: + assert ( len(ret['relresvec'])-1 <= params['b'].shape[0] ) + if __name__ == '__main__': import nose nose.main() diff -Nru krypy-0.1ubuntu1/setup.py krypy-0.1.1/setup.py --- krypy-0.1ubuntu1/setup.py 2013-02-05 11:23:26.000000000 +0000 +++ krypy-0.1.1/setup.py 2013-02-08 16:48:07.000000000 +0000 @@ -8,7 +8,7 @@ setup( name = 'krypy', packages = ['krypy'], - version = '0.1', + version = '0.1.1', description = 'Krylov subspace methods for linear algebraic systems', long_description = read('README.md'), author = 'André Gaul',