[PYTHON MATRIX-SIG] Release 1.0alpha1 now available

[Janne Sinkkonen]

Konrad Hinsen      <hinsen@physik.rwth-aachen.de> writes:

> As soon as I am back in business, I will update my LAPACK interface
> for the new version. I don't pretend that it is the last word on
> linear algebra libraries (some problems, like application to different
> types of matrices (symmetric etc.) have been completely ignored
> until now), but it has worked quite well for me in the past, and I have
> a couple of application algorithms that use it.

Well, I need your LAPACK interface right now... see below. :)

By the way, dot(a,b) of the package Numeric now seems to reduce over
the same axis for both a and b. While this seems sensible and
efficient, especially if the axis is -1, the behaviour looks strange
for a newcomer who maybe expects to find a matrix multiplication
function. Should there be a function behaving like the ordinary matrix
multiplication for 2D arrays, say dot2(), so that dot(a,b) would be
equal to dot2(a,transpose(b)) (for 2D arrays)? On the other hand, this
kind of behaviour does not mean anything sensible for N-dimensional
arrays, or does it?

# This module contains high-level Python interfaces to the
# LAPACK library.
#
# Written by Konrad Hinsen <hinsenk@ere.umontreal.ca>
# last revision: 1996-3-13
#
# Modified 1996-3-19 by Doug Heisterkamp
#        Change type of pivots from "i" to "l" in solveLinearEquations
#
# Modified 1996-8-11 by Janne Sinkkonen <janne@iki.fi>
#        Made compatible with NumPy 1.0alpha1. A bug in eigenvalues()
#        (triggered by complex arguments) corrected. dot() is now
#        _dot2() and unitMatrix() is _unitMatrix().

# import Numeric,umath
import Numeric
import copy

transpose=Numeric.transpose
zeros=Numeric.zeros
cp=copy.copy

# Error object

LinAlgError = 'LinAlgError'

# Helper routines

_lapack_type = {'f': 0, 'd': 1, 'F': 2, 'D': 3}
_lapack_letter = ['s', 'd', 'c', 'z']
_array_kind = {'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1}
_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
_array_type = [['f', 'd'], ['F', 'D']]

def _commonType(*arrays):
    kind = 0
    precision = 0
    for a in arrays:
	t = a.typecode()
	kind = max(kind, _array_kind[t])
	precision = max(precision, _array_precision[t])
    return _array_type[kind][precision]

def _castCopyAndTranspose(type, *arrays):
    cast_arrays = ()
    for a in arrays:
	if a.typecode() == type:
	    cast_arrays = cast_arrays + (cp(transpose(a)),)
	else:
	    cast_arrays = cast_arrays + (transpose(a).asType(type),)
    if len(cast_arrays) == 1:
	return cast_arrays[0]
    else:
	return cast_arrays

def _findLapackRoutine(name, type):
    name = _lapack_letter[_lapack_type[type]] + name
    module = 'lapack_' + name[:3]
    exec 'import ' + module
    module = eval(module)
    return getattr(module, name)

def _assertRank2(*arrays):
    for a in arrays:
	if len(a.shape) != 2:
	    raise LinAlgError, 'Array must be two-dimensional'

def _assertSquareness(*arrays):
    for a in arrays:
	if max(a.shape) != min(a.shape):
	    raise LinAlgError, 'Array must be square'
    


# unit matrix

def _unitMatrix(n):
    i = Numeric.arange(n)
    return Numeric.equal(i[:,Numeric.NewAxis], i)

# A dot function behaving like matrix multiplication (for 2D arrays)

def _dot2(a1, a2):
    return(Numeric.dot(a1,transpose(a2)))




# Singular value decomposition

def singularValueDecomposition(a):
    _assertRank2(a)
    n = a.shape[1]
    m = a.shape[0]
    t =_commonType(a)
    real_t = _array_type[0][_array_precision[t]]
    lapack_routine = _findLapackRoutine('gesvd', t)
    a = _castCopyAndTranspose(t, a)
    s = zeros((min(n,m),), real_t)
    u = zeros((m,m,), t)
    vt =zeros((n,n,), t)
    lwork = 10*max(m,n) # minimum value: max(3*min(m,n)+max(m,n),5*min(m,n)-4)
    work = zeros((lwork,), t)
    if _array_kind[t] == 1: # Complex routines take different arguments
	rwork = zeros((max(3*min(m,n),5*min(m,n)-4),), real_t)
	results = lapack_routine('A', 'A', m, n, a, m, s, u, m, vt, n,
				 work, lwork, rwork, 0)
    else:
	results = lapack_routine('A', 'A', m, n, a, m, s, u, m, vt, n,
				 work, lwork, 0)
    if results['info'] > 0:
	raise LinAlgError, 'SVD did not converge'
    return cp(transpose(u)), s, cp(transpose(vt))


# Generalized inverse

def generalizedInverse(a):
    u, s, vt = singularValueDecomposition(a)
    m = u.shape[0]
    n = vt.shape[0]
    cutoff = 1.e-10*Numeric.maximum.reduce(s)
    sm = zeros((n,m,),s.typecode())
    for i in range(min(n,m)):
	if s[i] > cutoff:
	    sm[i,i] = 1./s[i]
    return _dot2(_dot2(transpose(vt), sm), transpose(u))


# Linear equations
	
def solveLinearEquations(a, b):
    one_eq = len(b.shape) == 1
    if one_eq:
	b = b[:, Numeric.NewAxis]
    _assertRank2(a, b)
    _assertSquareness(a)
    n_eq = a.shape[0]
    n_rhs = b.shape[1]
    if n_eq != b.shape[0]:
	raise LinAlgError, 'Incompatible dimensions'
    t =_commonType(a, b)
    lapack_routine = _findLapackRoutine('gesv', t)
    a, b = _castCopyAndTranspose(t, a, b)
    pivots = zeros((n_eq,), 'l')
    results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
    if results['info'] > 0:
	raise LinAlgError, 'Singular matrix'
    if one_eq:
	return cp(Numeric.ravel(b))
    else:
	return cp(transpose(b))


# Matrix inversion

def inverse(a):
    return solveLinearEquations(a, _unitMatrix(a.shape[0]))




# Eigenvalues

def eigenvalues(a):
    _assertRank2(a)
    _assertSquareness(a)
    t =_commonType(a)
    real_t = _array_type[0][_array_precision[t]]
    lapack_routine = _findLapackRoutine('geev', t)
    a = _castCopyAndTranspose(t, a)
    n = a.shape[0]
    dummy = Numeric.zeros((1,), t)
    lwork = max(1,6*n) # minimum value: max(1,3*n) real, max(1,2*n) complex
    work = Numeric.zeros((lwork,), t)
    if _array_kind[t] == 1: # Complex routines take different arguments
	w = Numeric.zeros((n,), t)
	rwork = Numeric.zeros((n,),real_t)
	results = lapack_routine('N', 'N', n, a, n, w,
				 dummy, 1, dummy, 1, work, lwork, rwork, 0)
    else:
	wr = Numeric.zeros((n,), t)
	wi = Numeric.zeros((n,), t)
	results = lapack_routine('N', 'N', n, a, n, wr, wi,
				 dummy, 1, dummy, 1, work, lwork, 0)
	if Numeric.logicalAnd.reduce(Numeric.equal(wi, 0.)):
	    w = wr
	else:
	    w = wr+1j*wi
    if results['info'] > 0:
	raise LinAlgError, 'Eigenvalues did not converge'
    return w

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================