Tools for efficient sparse matrix computation

Abstract

To generate computationally efficient computer program code for carrying out computer computations on matrix organized input data, a program first is written in a relatively high-level language which includes programmer specifiable constructs for manipulating matrices and parts thereof; and which permits certain of the constructs to be annotated to specify programmer selected data structures and programmer selected operations on those data structures. This high-level program then is translated into a relatively high-level program into a relatively low-level language using low-level language routines that provide a compilable representation of the program, including all programmer selected data structures and all programmer selected operations on said data structures then the low-level representation of the program is compiled to generate computer executable code for implementing the program, including all programmer selected data structures and all programmer selected operations on said data structures.

Claims

What is claimed:

1. A method for generating computationally efficient computer program code for carrying out computer computations on a plurality of different matrix data types, said method comprising the steps of

writing a program in a relatively high-level language which includes programmer specifiable constructs for manipulating said matrix data types and parts thereof; at least certain of said constructs permitting annotations specifying programmer selected data structures for alternative implementations of matrix data types and programmer selected operations on said data structures for program implementation of the respective constructs;

translating the high-level program into a relatively low-level language using low-level language routines that provide a compilable representation of the program, including all programmer selected data structures and all programmer selected operations on said data structures;

compiling the low-level language representation of said program to generate computer executable code for implementing said program, including all programmer selected data structures and all programmer selected operations on said data structures.

2. The method of claim 2 wherein said low-level language is C++.

3. The method of claims 1 or 2 wherein said data structures include a sparse matrix data structures.

4. The method of any of claims 1-2 further including the step of executing said executable code on a computer for performing said computation.

Description

BACKGROUND OF THE INVENTION

When reality is modeled by computation, sparse matrices are often the link between the continuous physical world and the finite algorithmic one. Engineers and computational scientists who are not experts in sparse matrix algorithms must write programs that compute with sparse matrices for such applications as domain-specific preconditioning, numerical optimization, and prototyping modeling methods. There are three common approaches to writing such programs:

Use an existing numerical library. This is the easiest and safest course for the non-expert; but a library is unlikely to have the functionality or performance tuning required by a special-purpose application.

Write directly in a low-level computer language like C++ or Fortran. This can give flexibility and good performance; but developing low-level sparse matrix code is difficult, slow, and error-prone, and the non-expert is unlikely to be able to use the best available techniques effectively.

Write in a high-level interpreted matrix language like Matlab (as described, for example, in "Matlab Reference Guide", The MathWorks, Inc., Natick, Mass., 1992). This gives flexibility and ease of programming; but performance is unlikely to be acceptable because of the overhead of interpreting the code and because of inappropriate data structure choices, and inefficient data structure assesses especially if the application requires much vector- and scalar-level manipulation.

Thus, improved tools are needed for sparse matrix programming. The design goal is to provide a flexible, easy-to-use programming system that produces executable code that runs much faster than interpreted Matlab, and not much slower than benchmark Fortran or C library codes. The principal focus is on programs written at the vector (BLAS-1) and scalar (BLAS-0) level. Here efficiency requires choosing data representations based on more context than a single statement.

SUMMARY OF THE INVENTION

To achieve the design goal, this invention provides:

A flexible suite of data structures for representing sparse matrices and their constituents, including operations for manipulating and transforming the sparse matrices. This suite of data structures and transformations can be implemented in any low-level computer language like Fortran or C++. The implementation that is described herein is in the form of a C++ library.

A high-level language for describing matrix computations, coupled with an intermediate-level language for determining how the high-level operations should be implemented in terms of the data structure suite. This allows a programmer to describe sparse matrix computations at a high-level (comparable to Matlab) while, with a small added effort, getting performance comparable to tuned lower level code (comparable to Fortran or C++). In the illustrated embodiment of the invention, these intermediate-level descriptions are realized as annotations to Matlab; thus this embodiment employs a computer language called Annotated Matlab (or AML). A sparse matrix program written in AML can be translated into C++ either mechanically (by a compiler) or manually (by a programmer).

BRIEF DESCRIPTION OF THE DRAWINGS

Further objects and advantages of this invention will become evident when the following detailed description is read in conjunction with the attached drawings, in which:

FIG. 1 outlines a process for creating an efficient sparse matrix program in accordance with this invention;

FIG. 2 illustrates a standard representation of sparse matrices and vectors in a compressed column format;

FIG. 3 illustrates an ordered SPA data structure;

FIG. 4 illustrates how this invention expresses a low-level operation of iterating through nonzeros, which cannot be naturally expressed in a standard high-level language like Matlab;

FIG. 5 the technique this invention provides for associating a matrix or vector with a permutation to achieve automatic referencing;

FIG. 6 shows a high-level description of an LU factorization;

FIG. 7 illustrates the update step of the LU factorization shown in FIG. 6;

FIG. 8 shows a version of the LU factorization written in Matlab;

FIG. 9 shows the LU factorization written, in accordance with this invention, in AML;

FIG. 10 shows the C++ library calls into which the AML code in FIG. 9 is transformed; and

FIG. 11 shows the results of an experiment to measure the efficiency of the executable code generated in accordance with this invention.

DETAILED DESCRIPTION OF THE INVENTION

While the invention is described in some detail hereinbelow with specific reference to a particular embodiment, it is to be understood that there is no intent to limit it to this embodiment. On the contrary, the intent is to cover all modifications, alternatives and equivalents falling within the spirit and scope of the invention as defined by the appended claims.

FIG. 1 outlines a process for creating an efficient sparse matrix program in accordance with this invention. The user (scientist, engineer, or programmer) writes a program in AML, as at 21. This program is then translated (either by a compiler 22 or a programmer) into calls 23 to a C++ library 24 that implements the suite of low-level data structures and operations that are described hereinbelow. This C++ program is in turn translated into executable machine code by a standard C++ compiler 25, such as g++, thereby enabling the program to be executed on a standard programmed computer.

The following Description is organized as follows: Section 1 describes a suitable suite of low-level data structures and transformations for implementing this invention. Section 2 then illustrates the use of the intermediate-level representations, using as an example a program to perform LU factorization with partial pivoting. Section 3, in turn, provides a detailed specification for the AML language that realizes the intermediate-level representation.

1.0 Low-level data structures and transformations

The low-level data structures and operations are outlined below.

The low-level data structures come with the expected operations on scalars (called BLAS-0 operations) and operations on vectors (called BLAS-1 operations). The library does not include operations between vectors and matrices (BLAS-2) and between matrices and matrices (BLAS-3), because they will be handled by specifying them at the AML level, and then translating them into C++ as above.

1.1 Full matrices and vectors

These are the standard dense-matrix representations of a matrix as a two-dimensional array, and of a vector as a one-dimensional array.

1.2 Sparse matrices and vectors

Sparse matrices and vectors are represented in a compressed column format, as shown in FIG. 2. This is a standard representation for sparse matrices, as described for example in George and Liu, Computer Solution of Large Sparse Positive Definite Systems, Prentice-Hall 1981. Within each column, nonzeros and their associated row index are stored consecutively in increasing order of row index. Three separate arrays give, for each column, a start pointer, an end pointer, and a headroom pointer. Thus it is not necessary for columns to be contiguous, though they will be by default. Explicit zeros are allowed in the data structure.

An option is for a sparse matrix to be row-traversable, that is, to have an auxiliary data structure that allows efficient sequencing through its rows. Another option is to require the columns to be contiguous, so that the headroom pointers can be eliminated. Further, if successive columns are required to start where the previous column ended, then the column end pointers can also be eliminated. Another option is not to store the nonzeros in increasing order of row index.

1.3 Sparse accumulator

The sparse accumulator (or "SPA") is a representation of a sparse vector that permits random access to individual elements in constant time, and sequencing through the nonzero elements, in constant time per element accessed. The SPA may be implemented as an array of values, a parallel array of "occupied" flags, and a list of occupied row indices. Such a representation was described for example in Gilbert, Moler, and Schreiber, "Sparse matrices in Matlab: Design and implementation," SIAM Journal on Matrix Analysis 13: 333-356, 1992.

An option is for a SPA to be ordered, which means that it permits sequencing through the nonzero elements in increasing order of row index even if fill is being added at the same time. Such sequencing takes logarithmic time per element. This requires adding an ordered data structure of occupied row indices to the data structure. The ordered SPA representation and the operations are described in a commonly assigned, concurrently filed U.S. patent application Ser. No. 08/573,078 of Gilbert et al. (which is hereby incorporated by reference). An ordered SPA data structure is illustrated in FIG. 3.

1.4 Iterating through nonzeros

A commonly required operation is to iterate through the nonzero elements of a sparse vector. The compressed column format for a sparse vector (section 1.2) permits this operation, simply by sequencing through the segments of the arrays of elements and row indices that represent the given vector. The SPA permits this operation, by sequencing through the list of occupied row indices. (This sequence may not correspond to increasing row number.) The ordered SPA permits this operation, by repeatedly finding the minimum element of the heap of row indices.

As FIG. 4 shows, the low-level operation of iterating through nonzeros cannot be written in a natural fashion in a high-level language like Matlab. (The naive loop iterates through both zeros and nonzeros; the Matlab loop is inefficient, complicated, and difficult to program correctly.) FIG. 4 also shows our intermediate-level construction "for nzs . . . in" that specifies this iteration in AML.

1.5 Permutations

Another commonly required operation is to permute the rows or columns of a matrix or vector. Rather than doing this by explicitly reordering the elements, we can simply record the permutation as an array of integers. Then, for example, two rows of a matrix can be exchanged implicitly just by exchanging two entries in the permutation vector.

If a matrix has an associated permutation vector, then every access to an element of the matrix must be made indirectly through the permutation. As shown in FIG. 5, our intermediate-level representation includes a way to associate a matrix or vector with a permutation so that the indirect reference is made automatically. The AML construction for this is "view . . . through".

2.0 Example: Sparse LU factorization

This section describes a typical process for writing an efficient sparse matrix program using the tools provided by this invention. The program is one to compute the LU factorization, with partial pivoting, of a sparse matrix. The algorithm used is a "left-looking" column-oriented algorithm; see for example Gilbert and Peierls, "Sparse partial pivoting in time proportional to arithmetic operations," SIAM Journal on Scientific and Statistical Computing 9: 862-874, 1988 for a description of left-looking LU factorization algorithms.

FIG. 6 shows a high-level description of the LU factorization. The main loop is over the columns of the matrix; each iteration computes one column of the factors L and U. The most computationally expensive part of the computation is the Update step, which is shown in more detail in FIG. 7. The update step is itself a loop over columns to the left of the current column. Although the loop runs over all values of k from 1 to j-1, the update is empty if the element A(k,j) is zero. Thus, for efficiency, the loop should only run over the nonzero elements of column j of A. Since this column is changing during the loop, it will be necessary to use an ordered SPA representation for it.

FIG. 8 shows a version of the LU factorization written in Matlab. This implementation is inefficient both because it is executed by an interpreter, and because there is no way to specify that the active column should be represented and manipulated as a SPA. It is also difficult to read and error-prone to write, because of the complicated expression for iteration over nonzeros and the need to mention the permutation P explicitly whenever one of the matrix or vector variables is used.

FIG. 9 shows the LU factorization written in AML. Here the program is simpler to understand and easier to write correctly, because the permutation is specified only in the view statements, and the iteration through nonzeros is made explicit. Also, the program is much more efficient, because AML allows the specification of a SPA. (In fact, this will be an ordered SPA, because of the "in order" clause in the iteration.)

FIG. 10 shows the C++ library calls into which the AML code in FIG. 9 is transformed. This code is efficient, but is very difficult to read or to write. As described above, it is produced mechanically (by a programmer or a compiler) from the AML code.

FIG. 11 shows the results of an experiment to measure the efficiency of the executable code generated by our method, compared to several other codes for LU factorization. The four codes compared to ours are (1) the Fortran code from Gilbert and Peierls (cited above), which uses a complicated special-purpose depth-first search technique to replace the iteration through the ordered SPA; (2) the sparse LU factorization built into the core of Matlab, which is written in C and uses the same technique as Gilbert/Peierls; (3) the Matlab code from FIG. 8; and (4) an experimental code called SuperLU that uses very sophisticated cache optimization techniques, and is at present the fastest available LU factorization code. For each code, a bar shows the ratio of its running time to the running time of our executable code (on a Sun SparcStation 20 computer). Thus, a bar extending upward indicates that our code was faster; a bar extending downward indicates that our code was slower. The table shows timing results for each of seven sample sparse matrices.

The executable code derived from AML is faster than the interpreted Matlab code by a factor of about 100. The AML code ranges from about 3 times as fast to about 1/3 as fast as the benchmark Gilbert/Peierls and Matlab-core codes. In this experiment it is never worse than 1/4 times as fast as the best research code, SuperLU. This confirms that that the present invention attains the goal of being much faster than interpreted Matlab, and not much slower than benchmark Fortran or C codes.

3.0 Annotated Matlab

This section includes a detailed specification of the annotations that describe how to implement the matrix and vector data structures and operations. Adding these annotations to Matlab yields the AML language. As will be seen, a notion of lexical block is introduced into the Matlab programming environment to supplement its existing notion of function-wide scoping. Thus, it is to be understood that all statements that have statement lists inside (if, for, and while) define a block consisting of the statements inside. Further, a new begin/end statement also defines such a scope block.

    ______________________________________
               Syntax:
               begin <statement*> end
    ______________________________________

    ______________________________________
    Syntax:
    view <ArrayVariable> through <PermutationVariable*>
    <PermutationVariable*> is a list of permutation vectors, one
    corresponding to each dimension in the <ArrayVariable>.
    ______________________________________

    ______________________________________
    view <ArrayVariable> through <Permutation Variable*> as
    <OtherArrayVariable>
    ______________________________________

    ______________________________________
    Syntax:
    <Type> <Variable*>
    <Type>::= {int .vertline. real .vertline. complex} {scalar .vertline.
    vector .vertline. row vector .vertline. matrix} .vertline.
      permutation .vertline. range
    ______________________________________

    ______________________________________
             Syntax:
             <Representation> <Variable*>
             <Representation>::= {spa .vertline. sparse .vertline. full}
               {rows direct}
               {columns direct}
               {transposed}
               {ordered columns}
    ______________________________________

    ______________________________________
    Syntax:
    for nzs <variable> {in order} in <Vector> <Statement*> end
    ______________________________________

    ______________________________________
    Syntax:
    for (rows .vertline. columns) <variable> in <Array> <Statement*> end
    ______________________________________

    ______________________________________
             incremental <Array variable>›<index*>!
    ______________________________________

• Inventors
  Gilbert, John R.; Lamping, John O.; Mendhekar, Anurag; Shpeisman, Tatiana;
• Assignee
  Xerox Corporation (Stamford, CT)
• Published
  Jul-14-1998
• Current US Classes:
  708/520
  717/116
  717/147
  717/163
• Application #
  573707
• International Classes
  G06F 009/45
• Field of Search
  395/705 395/707 395/708 395/706 395/601 395/614 364/280.4 364/736.03
• Examiner
  Voeltz; Emanuel Todd
• US Patent References:
  5519868
  5526475