Add the code that generates lapack_lite from LAPACK sources.

This is from Numeric. I think we're still using the same generated sources.

5 files changed, 711 insertions, 0 deletions

diff --git a/numpy/linalg/lapack_lite/README b/numpy/linalg/lapack_lite/README
new file mode 100644
index 000000000..1c97d4d1e
--- /dev/null
+++ b/numpy/linalg/lapack_lite/README
@@ -0,0 +1,40 @@
+Regenerating lapack_lite source
+===============================
+
+:Author: David M. Cooke <cookedm@physics.mcmaster.ca>
+
+The ``numpy/linalg/blas_lite.c``, ``numpy/linalg/dlapack_lite.c``, and
+``numpy/linalg/zlapack_lite.c`` are ``f2c``'d versions of the LAPACK routines
+required by the ``LinearAlgebra`` module, and wrapped by the ``lapack_lite``
+module. The scripts in this directory can be used to create these files
+automatically from a directory of LAPACK source files.
+
+You'll need `Plex 1.1.4`_ installed to do the appropriate scrubbing.
+
+.. _Plex 1.1.4: http://www.cosc.canterbury.ac.nz/~greg/python/Plex/
+
+The routines that ``lapack_litemodule.c`` wraps are listed in
+``wrapped_routines``, along with a few exceptions that aren't picked up
+properly. Assuming that you have an unpacked LAPACK source tree in
+``~/LAPACK``, you generate the new routines in a directory ``new-lite/`` with::
+
+$ python ./make_lite.py wrapped_routines ~/LAPACK new-lite/
+
+This will grab the right routines, with dependencies, put them into the
+appropiate ``blas_lite.f``, ``dlapack_lite.f``, or ``zlapack_lite.f`` files,
+run ``f2c`` over them, then do some scrubbing similiar to that done to
+generate the CLAPACK_ distribution.
+
+.. _CLAPACK: http://netlib.org/clapack/index.html
+
+The versions in Numeric CVS as of 2005-04-12 use the LAPACK source from the
+`Debian package lapack3`_, version 3.0.20000531a-6. It was found that these
+(being regularly maintained) worked better than the patches to the last
+released version of LAPACK available at the LAPACK_ page.
+
+.. _Debian package lapack3: http://packages.debian.org/unstable/libs/lapack3
+.. _LAPACK: http://netlib.org/lapack/index.html
+
+A slightly-patched ``f2c`` was used to add parentheses around ``||`` expressions
+and the arguments to ``<<`` to silence gcc warnings. Edit
+the ``src/output.c`` in the ``f2c`` source to do this.
diff --git a/numpy/linalg/lapack_lite/clapack_scrub.py b/numpy/linalg/lapack_lite/clapack_scrub.py
new file mode 100644
index 000000000..df00851c3
--- /dev/null
+++ b/numpy/linalg/lapack_lite/clapack_scrub.py
@@ -0,0 +1,276 @@
+#!/usr/bin/env python2.4
+
+import sys, os
+from cStringIO import StringIO
+import re
+
+from Plex import *
+from Plex.Traditional import re as Re
+
+class MyScanner(Scanner):
+    def __init__(self, info, name='<default>'):
+        Scanner.__init__(self, self.lexicon, info, name)
+
+    def begin(self, state_name):
+#        if self.state_name == '':
+#            print '<default>'
+#        else:
+#            print self.state_name
+        Scanner.begin(self, state_name)
+
+def sep_seq(sequence, sep):
+    pat = Str(sequence[0])
+    for s in sequence[1:]:
+        pat += sep + Str(s)
+    return pat
+
+def runScanner(data, scanner_class, lexicon=None):
+    info = StringIO(data)
+    outfo = StringIO()
+    if lexicon is not None:
+        scanner = scanner_class(lexicon, info)
+    else:
+        scanner = scanner_class(info)
+    while 1:
+        value, text = scanner.read()
+        if value is None:
+            break
+        elif value is IGNORE:
+            pass
+        else:
+            outfo.write(value)
+    return outfo.getvalue(), scanner
+
+class LenSubsScanner(MyScanner):
+    """Following clapack, we remove ftnlen arguments, which f2c puts after
+    a char * argument to hold the length of the passed string. This is just
+    a nuisance in C.
+    """
+    def __init__(self, info, name='<ftnlen>'):
+        MyScanner.__init__(self, info, name)
+        self.paren_count = 0
+
+    def beginArgs(self, text):
+        if self.paren_count == 0:
+            self.begin('args')
+        self.paren_count += 1
+        return text
+
+    def endArgs(self, text):
+        self.paren_count -= 1
+        if self.paren_count == 0:
+            self.begin('')
+        return text
+
+    digits = Re('[0-9]+')
+    iofun = Re(r'\([^;]*;')
+    decl = Re(r'\([^)]*\)[,;'+'\n]')
+    any = Re('[.]*')
+    S = Re('[ \t\n]*')
+    cS = Str(',') + S
+    len_ = Re('[a-z][a-z0-9]*_len')
+
+    iofunctions = Str("s_cat", "s_copy", "s_stop", "s_cmp",
+                      "i_len", "do_fio", "do_lio") + iofun
+
+    # Routines to not scrub the ftnlen argument from
+    keep_ftnlen = (Str('ilaenv_') | Str('s_rnge')) + Str('(')
+
+    lexicon = Lexicon([
+        (iofunctions,                           TEXT),
+        (keep_ftnlen,                           beginArgs),
+        State('args', [
+            (Str(')'),   endArgs),
+            (Str('('),   beginArgs),
+            (AnyChar,    TEXT),
+        ]),
+        (cS+Re(r'[1-9][0-9]*L'),                IGNORE),
+        (cS+Str('ftnlen')+Opt(S+len_),          IGNORE),
+        (cS+sep_seq(['(', 'ftnlen', ')'], S)+S+digits,      IGNORE),
+        (Bol+Str('ftnlen ')+len_+Str(';\n'),    IGNORE),
+        (cS+len_,                               TEXT),
+        (AnyChar,                               TEXT),
+    ])
+
+def scrubFtnlen(source):
+    return runScanner(source, LenSubsScanner)[0]
+
+def cleanSource(source):
+    # remove whitespace at end of lines
+    source = re.sub(r'[\t ]+\n', '\n', source)
+    # remove comments like .. Scalar Arguments ..
+    source = re.sub(r'(?m)^[\t ]*/\* *\.\. .*?\n', '', source)
+    # collapse blanks of more than two in-a-row to two
+    source = re.sub(r'\n\n\n\n+', r'\n\n\n', source)
+    return source
+
+class LineQueue(object):
+    def __init__(self):
+        object.__init__(self)
+        self._queue = []
+
+    def add(self, line):
+        self._queue.append(line)
+
+    def clear(self):
+        self._queue = []
+
+    def flushTo(self, other_queue):
+        for line in self._queue:
+            other_queue.add(line)
+        self.clear()
+
+    def getValue(self):
+        q = LineQueue()
+        self.flushTo(q)
+        s = ''.join(q._queue)
+        self.clear()
+        return s
+
+class CommentQueue(LineQueue):
+    def __init__(self):
+        LineQueue.__init__(self)
+
+    def add(self, line):
+        if line.strip() == '':
+            LineQueue.add(self, '\n')
+        else:
+            line = '  ' + line[2:-3].rstrip() + '\n'
+            LineQueue.add(self, line)
+
+    def flushTo(self, other_queue):
+        if len(self._queue) == 0:
+            pass
+        elif len(self._queue) == 1:
+            other_queue.add('/*' + self._queue[0][2:].rstrip() + ' */\n')
+        else:
+            other_queue.add('/*\n')
+            LineQueue.flushTo(self, other_queue)
+            other_queue.add('*/\n')
+        self.clear()
+
+# This really seems to be about 4x longer than it needs to be
+def cleanComments(source):
+    lines = LineQueue()
+    comments = CommentQueue()
+    def isCommentLine(line):
+        return line.startswith('/*') and line.endswith('*/\n')
+
+    blanks = LineQueue()
+    def isBlank(line):
+        return line.strip() == ''
+
+    def SourceLines(line):
+        if isCommentLine(line):
+            comments.add(line)
+            return HaveCommentLines
+        else:
+            lines.add(line)
+            return SourceLines
+    def HaveCommentLines(line):
+        if isBlank(line):
+            blanks.add('\n')
+            return HaveBlankLines
+        elif isCommentLine(line):
+            comments.add(line)
+            return HaveCommentLines
+        else:
+            comments.flushTo(lines)
+            lines.add(line)
+            return SourceLines
+    def HaveBlankLines(line):
+        if isBlank(line):
+            blanks.add('\n')
+            return HaveBlankLines
+        elif isCommentLine(line):
+            blanks.flushTo(comments)
+            comments.add(line)
+            return HaveCommentLines
+        else:
+            comments.flushTo(lines)
+            blanks.flushTo(lines)
+            lines.add(line)
+            return SourceLines
+
+    state = SourceLines
+    for line in StringIO(source):
+        state = state(line)
+    comments.flushTo(lines)
+    return lines.getValue()
+
+def removeHeader(source):
+    lines = LineQueue()
+
+    def LookingForHeader(line):
+        m = re.match(r'/\*[^\n]*-- translated', line)
+        if m:
+            return InHeader
+        else:
+            lines.add(line)
+            return LookingForHeader
+    def InHeader(line):
+        if line.startswith('*/'):
+            return OutOfHeader
+        else:
+            return InHeader
+    def OutOfHeader(line):
+        if line.startswith('#include "f2c.h"'):
+            pass
+        else:
+            lines.add(line)
+        return OutOfHeader
+
+    state = LookingForHeader
+    for line in StringIO(source):
+        state = state(line)
+    return lines.getValue()
+
+def replaceDlamch(source):
+    """Replace dlamch_ calls with appropiate macros"""
+    def repl(m):
+        s = m.group(1)
+        return dict(E='EPSILON', P='PRECISION', S='SAFEMINIMUM',
+                    B='BASE')[s[0]]
+    source = re.sub(r'dlamch_\("(.*?)"\)', repl, source)
+    source = re.sub(r'^\s+extern.*? dlamch_.*?;$(?m)', '', source)
+    return source
+
+# do it
+
+def scrubSource(source, nsteps=None, verbose=False):
+    steps = [
+             ('scrubbing ftnlen', scrubFtnlen),
+             ('remove header', removeHeader),
+             ('clean source', cleanSource),
+             ('clean comments', cleanComments),
+             ('replace dlamch_() calls', replaceDlamch),
+            ]
+
+    if nsteps is not None:
+        steps = steps[:nsteps]
+
+    for msg, step in steps:
+        if verbose:
+            print msg
+        source = step(source)
+
+    return source
+
+if __name__ == '__main__':
+    filename = sys.argv[1]
+    outfilename = os.path.join(sys.argv[2], os.path.basename(filename))
+    fo = open(filename, 'r')
+    source = fo.read()
+    fo.close()
+
+    if len(sys.argv) > 3:
+        nsteps = int(sys.argv[3])
+    else:
+        nsteps = None
+
+    source = scrub_source(source, nsteps, verbose=True)
+
+    writefo = open(outfilename, 'w')
+    writefo.write(source)
+    writefo.close()
+
diff --git a/numpy/linalg/lapack_lite/fortran.py b/numpy/linalg/lapack_lite/fortran.py
new file mode 100644
index 000000000..7be986a8e
--- /dev/null
+++ b/numpy/linalg/lapack_lite/fortran.py
@@ -0,0 +1,114 @@
+import re
+import itertools
+
+def isBlank(line):
+    return not line
+def isLabel(line):
+    return line[0].isdigit()
+def isComment(line):
+    return line[0] != ' '
+def isContinuation(line):
+    return line[5] != ' '
+
+COMMENT, STATEMENT, CONTINUATION = 0, 1, 2
+def lineType(line):
+    """Return the type of a line of Fortan code."""
+    if isBlank(line):
+        return COMMENT
+    elif isLabel(line):
+        return STATEMENT
+    elif isComment(line):
+        return COMMENT
+    elif isContinuation(line):
+        return CONTINUATION
+    else:
+        return STATEMENT
+
+class LineIterator(object):
+    """LineIterator(iterable)
+
+    Return rstrip()'d lines from iterable, while keeping a count of the
+    line number in the .lineno attribute.
+    """
+    def __init__(self, iterable):
+        object.__init__(self)
+        self.iterable = iter(iterable)
+        self.lineno = 0
+    def __iter__(self):
+        return self
+    def next(self):
+        self.lineno += 1
+        line = self.iterable.next()
+        line = line.rstrip()
+        return line
+
+class PushbackIterator(object):
+    """PushbackIterator(iterable)
+
+    Return an iterator for which items can be pushed back into.
+    Call the .pushback(item) method to have item returned as the next
+    value of .next().
+    """
+    def __init__(self, iterable):
+        object.__init__(self)
+        self.iterable = iter(iterable)
+        self.buffer = []
+
+    def __iter__(self):
+        return self
+
+    def next(self):
+        if self.buffer:
+            return self.buffer.pop()
+        else:
+            return self.iterable.next()
+
+    def pushback(self, item):
+        self.buffer.append(item)
+
+def fortranSourceLines(fo):
+    """Return an iterator over statement lines of a Fortran source file.
+
+    Comment and blank lines are stripped out, and continuation lines are
+    merged.
+    """
+    numberingiter = LineIterator(fo)
+    # add an extra '' at the end
+    with_extra = itertools.chain(numberingiter, [''])
+    pushbackiter = PushbackIterator(with_extra)
+    for line in pushbackiter:
+        t = lineType(line)
+        if t == COMMENT:
+            continue
+        elif t == STATEMENT:
+            lines = [line]
+            # this is where we need the extra '', so we don't finish reading
+            # the iterator when we don't want to handle that
+            for next_line in pushbackiter:
+                t = lineType(next_line)
+                if t == CONTINUATION:
+                    lines.append(next_line[6:])
+                else:
+                    pushbackiter.pushback(next_line)
+                    break
+            yield numberingiter.lineno, ''.join(lines)
+        else:
+            raise ValueError("jammed: continuation line not expected: %s:%d" %
+                             (fo.name, numberingiter.lineno))
+
+def getDependencies(filename):
+    """For a Fortran source file, return a list of routines declared as EXTERNAL
+    in it.
+    """
+    fo = open(filename)
+    external_pat = re.compile(r'^\s*EXTERNAL\s', re.I)
+    routines = []
+    for lineno, line in fortranSourceLines(fo):
+        m = external_pat.match(line)
+        if m:
+            names = line = line[m.end():].strip().split(',')
+            names = [n.strip().lower() for n in names]
+            names = [n for n in names if n]
+            routines.extend(names)
+    fo.close()
+    return routines
diff --git a/numpy/linalg/lapack_lite/make_lite.py b/numpy/linalg/lapack_lite/make_lite.py
new file mode 100755
index 000000000..dec0be017
--- /dev/null
+++ b/numpy/linalg/lapack_lite/make_lite.py
@@ -0,0 +1,264 @@
+#!/usr/bin/env python
+
+import sys, os
+import fortran
+import clapack_scrub
+
+try: set
+except NameError:
+    from sets import Set as set
+
+# Arguments to pass to f2c. You'll always want -A for ANSI C prototypes
+# Others of interest: -a to not make variables static by default
+#                     -C to check array subscripts
+F2C_ARGS = '-A'
+
+# The header to add to the top of the *_lite.c file. Note that dlamch_() calls
+# will be replaced by the macros below by clapack_scrub.scrub_source()
+HEADER = '''\
+/*
+NOTE: This is generated code. Look in Misc/lapack_lite for information on
+      remaking this file.
+*/
+#include "Numeric/f2c.h"
+
+#ifdef HAVE_CONFIG
+#include "config.h"
+#else
+extern doublereal dlamch_(char *);
+#define EPSILON dlamch_("Epsilon")
+#define SAFEMINIMUM dlamch_("Safe minimum")
+#define PRECISION dlamch_("Precision")
+#define BASE dlamch_("Base")
+#endif
+
+extern doublereal dlapy2_(doublereal *x, doublereal *y);
+
+'''
+
+class FortranRoutine:
+    """Wrapper for a Fortran routine in a file.
+    """
+    type = 'generic'
+    def __init__(self, name=None, filename=None):
+        self.filename = filename
+        if name is None:
+            root, ext = os.path.splitext(filename)
+            name = root
+        self.name = name
+        self._dependencies = None
+
+    def dependencies(self):
+        if self._dependencies is None:
+            deps = fortran.getDependencies(self.filename)
+            self._dependencies = [d.lower() for d in deps]
+        return self._dependencies
+
+class UnknownFortranRoutine(FortranRoutine):
+    """Wrapper for a Fortran routine for which the corresponding file
+    is not known.
+    """
+    type = 'unknown'
+    def __init__(self, name):
+        FortranRoutine.__init__(self, name=name, filename='<unknown>')
+
+    def dependencies(self):
+        return []
+
+class FortranLibrary:
+    """Container for a bunch of Fortran routines.
+    """
+    def __init__(self, src_dirs):
+        self._src_dirs = src_dirs
+        self.names_to_routines = {}
+
+    def _findRoutine(self, rname):
+        rname = rname.lower()
+        for s in self._src_dirs:
+            ffilename = os.path.join(s, rname + '.f')
+            if os.path.exists(ffilename):
+                return self._newFortranRoutine(rname, ffilename)
+        return UnknownFortranRoutine(rname)
+
+    def _newFortranRoutine(self, rname, filename):
+        return FortranRoutine(rname, filename)
+
+    def addIgnorableRoutine(self, rname):
+        """Add a routine that we don't want to consider when looking at
+        dependencies.
+        """
+        rname = rname.lower()
+        routine = UnknownFortranRoutine(rname)
+        self.names_to_routines[rname] = routine
+
+    def addRoutine(self, rname):
+        """Add a routine to the library.
+        """
+        self.getRoutine(rname)
+
+    def getRoutine(self, rname):
+        """Get a routine from the library. Will add if it's not found.
+        """
+        unique = []
+        rname = rname.lower()
+        routine = self.names_to_routines.get(rname, unique)
+        if routine is unique:
+            routine = self._findRoutine(rname)
+            self.names_to_routines[rname] = routine
+        return routine
+
+    def allRoutineNames(self):
+        """Return the names of all the routines.
+        """
+        return self.names_to_routines.keys()
+
+    def allRoutines(self):
+        """Return all the routines.
+        """
+        return self.names_to_routines.values()
+
+    def resolveAllDependencies(self):
+        """Try to add routines to the library to satisfy all the dependencies
+        for each routine in the library.
+
+        Returns a set of routine names that have the dependencies unresolved.
+        """
+        done_this = set()
+        last_todo = set()
+        while 1:
+            todo = set(self.allRoutineNames()) - done_this
+            if todo == last_todo:
+                break
+            for rn in todo:
+                r = self.getRoutine(rn)
+                deps = r.dependencies()
+                for d in deps:
+                    self.addRoutine(d)
+                done_this.add(rn)
+            last_todo = todo
+        return todo
+
+class LapackLibrary(FortranLibrary):
+    def _newFortranRoutine(self, rname, filename):
+        routine = FortranLibrary._newFortranRoutine(self, rname, filename)
+        if filename.find('BLAS') != -1:
+            routine.type = 'blas'
+        elif rname.startswith('z'):
+            routine.type = 'zlapack'
+        else:
+            routine.type = 'dlapack'
+        return routine
+
+    def allRoutinesByType(self, typename):
+        routines = [(r.name,r) for r in self.allRoutines() if r.type == typename]
+        routines.sort()
+        return [a[1] for a in routines]
+
+def printRoutineNames(desc, routines):
+    print desc
+    for r in routines:
+        print '\t%s' % r.name
+
+def getLapackRoutines(wrapped_routines, ignores, lapack_dir):
+    blas_src_dir = os.path.join(lapack_dir, 'BLAS', 'SRC')
+    if not os.path.exists(blas_src_dir):
+        blas_src_dir = os.path.join(lapack_dir, 'blas', 'src')
+    lapack_src_dir = os.path.join(lapack_dir, 'SRC')
+    if not os.path.exists(lapack_src_dir):
+        lapack_src_dir = os.path.join(lapack_dir, 'src')
+    library = LapackLibrary([blas_src_dir, lapack_src_dir])
+
+    for r in ignores:
+        library.addIgnorableRoutine(r)
+
+    for w in wrapped_routines:
+        library.addRoutine(w)
+
+    library.resolveAllDependencies()
+
+    return library
+
+def getWrappedRoutineNames(wrapped_routines_file):
+    fo = open(wrapped_routines_file)
+    routines = []
+    ignores = []
+    for line in fo:
+        line = line.strip()
+        if not line or line.startswith('#'):
+            continue
+        if line.startswith('IGNORE:'):
+            line = line[7:].strip()
+            ig = line.split()
+            ignores.extend(ig)
+        else:
+            routines.append(line)
+    return routines, ignores
+
+def dumpRoutineNames(library, output_dir):
+    for typename in ['unknown', 'blas', 'dlapack', 'zlapack']:
+        routines = library.allRoutinesByType(typename)
+        filename = os.path.join(output_dir, typename + '_routines.lst')
+        fo = open(filename, 'w')
+        for r in routines:
+            deps = r.dependencies()
+            fo.write('%s: %s\n' % (r.name, ' '.join(deps)))
+        fo.close()
+
+def concatenateRoutines(routines, output_file):
+    output_fo = open(output_file, 'w')
+    for r in routines:
+        fo = open(r.filename, 'r')
+        source = fo.read()
+        fo.close()
+        output_fo.write(source)
+    output_fo.close()
+
+class F2CError(Exception):
+    pass
+
+def runF2C(fortran_filename, output_dir):
+    # we're assuming no funny business that needs to be quoted for the shell
+    cmd = "f2c %s -d %s %s" % (F2C_ARGS, output_dir, fortran_filename)
+    rc = os.system(cmd)
+    if rc != 0:
+        raise F2CError
+
+def scrubF2CSource(c_file):
+    fo = open(c_file, 'r')
+    source = fo.read()
+    fo.close()
+    source = clapack_scrub.scrubSource(source, verbose=True)
+    fo = open(c_file, 'w')
+    fo.write(HEADER)
+    fo.write(source)
+    fo.close()
+
+def main():
+    if len(sys.argv) != 4:
+        print 'Usage: %s wrapped_routines_file lapack_dir output_dir' % \
+              (sys.argv[0],)
+        return
+    wrapped_routines_file = sys.argv[1]
+    lapack_src_dir = sys.argv[2]
+    output_dir = sys.argv[3]
+
+    wrapped_routines, ignores = getWrappedRoutineNames(wrapped_routines_file)
+    library = getLapackRoutines(wrapped_routines, ignores, lapack_src_dir)
+
+    dumpRoutineNames(library, output_dir)
+
+    for typename in ['blas', 'dlapack', 'zlapack']:
+        print 'creating %s_lite.c ...'  % typename
+        routines = library.allRoutinesByType(typename)
+        fortran_file = os.path.join(output_dir, typename+'_lite.f')
+        c_file = fortran_file[:-2] + '.c'
+        concatenateRoutines(routines, fortran_file)
+        try:
+            runF2C(fortran_file, output_dir)
+        except F2CError:
+            print 'f2c failed on %s' % fortran_file
+            break
+        scrubF2CSource(c_file)
+
+if __name__ == '__main__':
+    main()
diff --git a/numpy/linalg/lapack_lite/wrapped_routines b/numpy/linalg/lapack_lite/wrapped_routines
new file mode 100644
index 000000000..8b7153ba9
--- /dev/null
+++ b/numpy/linalg/lapack_lite/wrapped_routines
@@ -0,0 +1,17 @@
+dgeev
+zgeev
+dsyevd
+zheevd
+dgelsd
+zgelsd
+dgesv
+zgesv
+dgetrf
+zgetrf
+dpotrf
+zpotrf
+dgesdd
+zgesdd
+# need this b/c it's not properly declared as external in the BLAS source
+dcabs1
+IGNORE: dlamch