1 Introduction

This document describes a set of matrix abstractions that can be used for writing sparse matrix computations in a generic programming (GP) style. A progression of abstractions is used, with the first being the most abstract, and the subsequent exposing increasing amounts of information in order to allow progressively more specialized implementation to be developed.

This can be easily be done in C++ using abstract classes and templates. In Section 2, we will introduce a very basic matrix abstraction that can be used for generic matrix programming in C++. In the subsequent sections, we will introduce refinements to this model that allow a programmer to express more specialized and efficient implementations, without violating the basic generic programming spirit.

We envision these abstraction being used in one of two ways. First, these abstractions can be used by a programmer directly in order to implement a generic sparse matrix computation library. These abstractions can also be used as targets of a program restructuring tool, like our sparse compiler technology. In this note, we will sketch the techniques required to implement such a tool.

2 Base abstraction

We believe that, at its most basic level, a matrix can be describes in terms of two properties,

• a set of entries whose indices fall within a particular set of bounds, and

• a method for accessing individual indices in order to read or write their associated element.

Both dense and sparse matrices can be modeled using these two properties.

Abstract matrices are (almost?) always rectangular in shape. Their bounds can be conveniently expressed as the lower and upper bounding index on each dimension. Note that perhaps only a subset of the indices within these rectangular bounds are actually stored; this is a separate issue that is discussed in Section 6.

When accessing an individual matrix element, it may be permissible either to read or write the element, but not both. We leave this issue to particular implementations of this interface.

Figure 1 shows an abstract C++ class definition that provides this interface. Particular storage formats are subclasses of this abstract matrix class. Figure 2 shows an implementation of coordinate format using this scheme.

Given this interface, matrix computations can be expressed in a GP style in the usual manner. That is, the code can be written as a set of for or do loops, inside of which individual matrix elements are accessed. Modulo some array notation shorthand, this is essentially the type of code that can be found in a matrix computation textbook, such as [1]. We call this a control-centric

template <class T>
class matrix {
public:
    typedef T basetype;
    // The enveloping bounds of the matrix.
    virtual int rows() = 0;
    virtual int cols() = 0;
    // A read access to the matrix.
    virtual T rval(int i, int j) = 0;
    // An access suitable for writing.
    virtual T &lval(int i, int j) = 0;
    // destructors, constructors, etc.
};

Figure 1: The base sparse matrix abstraction

description of algorithms. Generic algorithms for matrix-matrix multiplication and forward solve are shown in Figure 3 and Figure 4, respectively.

3 Enumeration method

3.1 The problem

If we were to instantiate the algorithms in Figure 3 and Figure 4 for any common sparse matrix format, the performance of the resulting code would be abysmal. The fundamental reason is that the control-centric versions of these programs access entries regardless of whether or not they are zero or non-zero. In the case of sparse matrices, this can result in one or more factors of slowdown.

A fundamentally more efficient way of expressing sparse matrix programs is in a data-centric manner. In this approach, the non-zero entries of a sparse matrix are enumerated directly, and, for each entry, the iterations of the original control-centric code that accessed that entry are executed. That is, instead of enumerating the iterations and accessing the data, the data is enumerated directly and then the corresponding iterations are executed.

3.2 The data abstraction

In order to write GP’s in a data-centric manner, a method for enumerating the matrix entries is required. Figure 5 how the base matrix class can be refined to include such a method.

3.3 The program transformations

Program transformations to produce data-centric programs from control-centric programs are described in [5]. Using these techniques, the control-centric ma-

template <class T>
class coo_matrix : public matrix<T> {
protected:
int m, n, nnz;
vector<int> indx, jndx;
vector<T> vals;
public:
coo_matrix(...) { ...; }
virtual int rows() { return m; }
virtual int cols() { return n; }
virtual T rval(int i, int j) {
    for (int k = 0; k < nnz; k++)
        if (i == indx[k] && j == jndx[k])
            return vals[k];
    // (i,j) not stored -- implicitly 0.
    return 0;
}
virtual T &lval(int i, int j) {
    for (int k = 0; k < nnz; k++)
        if (i == indx[k] && j == jndx[k])
            return vals[k];
    // (i,j) not stored -- insert a new location.
    indx.push_back(i);
    jndx.push_back(j);
    vals.push_back(0);
    return vals.back();
}
}

Figure 2: An implementation of coordinate storage

// C = alpha * A * B + beta * C
template<class T>
void matrix_multip(T alpha, matrix<T> &A, matrix<T> &B,
                    T beta, matrix<T> &C)
{
    assert(C.rows() == A.rows());
    assert(A.cols() == B.rows());
    assert(C.cols() == B.cols());
    for (int i = 0; i < C.rows(); i++)
        for (int j = 0; j < C.cols(); j++) {
            T sum = C.rval(i,j);
            sum *= beta;
            for (int k = 0; k < A.cols(); k++)
                sum += alpha * A.rval(i,k) * B.rval(k,j);
            C.lval(i,j) = sum;
        }
}

Figure 3: Generic algorithm for matrix multiplication

// X = L^-1 B
template<class T>
void forward_solve(matrix<T> &L, matrix<T> &B, matrix<T> &X)
{
    assert(L.rows() == L.cols() && L.cols() == B.rows()
           && B.rows() == X.rows());
    assert(B.cols() == X.cols());

    for (int i = 0; i < B.rows(); i++)
        for (int k = 0; k < B.cols(); k++)
            X.lval(i,k) = B.rval(i,k);
    for (int i=0; i<L.rows(); i++)
        for (int k = 0; k < B.cols(); k++) {
            for (int j=0; j<i-1; j++)
                X.lval(i,k) -= L.rval(i,j)*X.rval(j,k);
            X.lval(i,k) /= L.rval(i,i);
        }
}

Figure 4: Generic algorithm for forward solve

template <class T>
class enumerable_matrix
: public matrix<T> {
protected:
class iterator_state {
friend class iterator;
protected:
virtual int row_index() = 0;
virtual int col_index() = 0;
virtual T& value() = 0;
virtual bool more_p() = 0;
virtual void advance() = 0;
// constructors, destructors, etc.
};
public:
class iterator {
friend class enumerable_matrix;
private:
iterator_state *p;
protected:
iterator(iterator_state *_p) p(_p) { };
public:
int row_index() { return p->row_index(); }
int col_index() { return p->col_index(); }
T& operator *() { return p->value(); }
operator bool() { return p->more_p(); }
void operator (++)() { p->advance(); }
// destructors, etc.
};
virtual iterator enum() = 0;
// constructors, destructors, etc.
};

Figure 5: The enumerable matrix abstraction trix multiplication template code shown in Figure 3 can be used to generate a specialized data-centric template code that can be used with subclasses of the enumerable matrix class. This derived code is shown in Figure 6.

// C = alpha * A * B + beta * C
template<class T>
void matrix_multip(T alpha, enumerable_matrix<T> &A,
              matrix<T> &B, T beta, matrix<T> &C)
{
assert(C.rows() == A.rows());
assert(A.cols() == B.rows());
assert(C.cols() == B.cols());
for (int i = 0; i < C.rows(); i++)
   for (int j = 0; j < C.cols(); j++)
       C.lval(i,j) *= beta;
for (enumerable_matrix<T>::iterator it = A.enum();
    it; it++) {
   int i = it.row_index(), k = it.col_index(); T A_ik = *it;
   for (int j = 0; j < C.cols(); j++)
       C.rval(i,j) += alpha * A_ik * B.rval(k,j);
}
}

Figure 6: Data-centric matrix multiplication

In addition to the actual data-centric transformations, there are two changes that need to be made to a conventional C++ compiler to accommodate this abstraction,

• The compiler must be “aware” of the enumerable matrix class, as it is only subclasses of this abstract class that are guaranteed to have an enum method with the “known” semantics. If the compiler were not “aware” of the enumerable matrix class, then it would have to analyze arbitrary user code to determine the existence of a suitable enumeration method.

• The data-centric specialization must be inserted somewhere in the compilation process. One possibility would be to introduce it into the instantiation logic. That is, whenever a template function takes a matrix as an argument, the compiler would attempt to make the instantiated function data-centric. Another possibility would be to do the data-centric transformations on the template code directly. However, there is not enough information to type check template code, so performing program transformations may not be possible at this stage.

This question of when to apply program transformations to generic algorithms arises throughout this document.

4 Distance/direction

4.1 The problem

The data-centric approach cannot always be used. In particular, if there are loop carried dependencies on the control-centric code, and if there are no guarantees about the order in which matrix entries are visited by an enumeration method, then the data-centric approach cannot be used because it cannot be guaranteed that the original loop’s dependencies are satisfied. If, however, certain guarantees could be made about the order of enumeration, then it might be possible to show that the dependencies were satisfied.

For instance, the control-centric version of triangular solve code shown in shown in Figure 4 has a dependence carried by the outermost loop. If a matrix can be enumerated such that its entries are visited in (row, column) lexicographical order, then the row-orientied control-centric loops can be used to generate a legal data-centric program. Similarly, if a matrix can be enumerated in (column, row) lexicographical order, then a column-oriented version can be used.

4.2 The data abstraction

An annotation on an enumeration method could be used to convey any guarantees about the order in which matrix entries are visited. Figure 7 shows such an annotation added to the coo matrix class.

template <class T>
class coo_matrix : public enumerable_matrix<T> {
    // ...
public:
    // ...
    class iterator_state :
        public coo_matrix<T>::iterator_state
    { ... };
#pragma ordered (0,+) (+,*)
    iterator enum() { ... };
    // ...
};

Figure 7: Coordinate storage with an ordered enumerator

4.3 The program transformations

Given the ordering information obtained from the annotation on the enumeration method, the compiler can easily determine whether the resulting datacentric code is legal. This is handled within the existing data-centric program transformation framework.

Figure 8 shows the forward solve template code specialized using the datacentric transformations. Notice the annotation on the template indicating what ordering must be satisfied by a matrix class in order for the specialized template code to be applicable. This condition must be checked by the template instantator.

// X = L^-1 B
#pragma L.enum() is ordered (0,+) (+,*)
template<class T>
void forward_solve(enumerable_matrix<T> &L, matrix<T> &B,
                   matrix<T> &X)
{
    assert(L.rows() == L.cols() && L.cols() == B.rows()
           && B.rows() == X.rows());
    assert(B.cols() == X.cols());

    for (int i = 0; i < B.rows(); i++)
        for (int k = 0; k < B.cols(); k++)
            X.lval(i,k) = B.rval(i,k);
    for (enumerable_matrix<T>::iterator it = L.enum();
         it; it++) {
        int i = it.row_index(), j = it.col_index(); T L_ij = *it;
        if (i > j)
            for (int k = 0; k < B.cols(); k++)
                X.lval(i,k) -= L.rval(i,j)*X.rval(j,k);
        else if (i == j)
            for (int k = 0; k < B.cols(); k++)
                X.lval(i,k) /= L.rval(i,i);
        else
            ; // nothing
    }
}

Figure 8: Data-centric forward solve

5 Combining operator

5.1 The problem

We have ignored one minor detail to this point: what happens when a sparse matrix storage format contains multiple values for the same (i, j) entry? This case arises, for instance, when the format stores unassembled portions of a sparse matrix. In this case, the multiple values are to be added to produce the final (i, j) entry. We call the operator used to combine duplicate entries the combining operator.

If a matrix format contains duplicate entries, the data-centric approach cannot be applied in general. It can be applied, however, if it can be shown that visiting the duplicate entries will not affect the final result. For instance, matrixmultiplication can safely by made data-centric if it is known that duplicate entries either do not occur or can be combining using the addition operator to produce the final result.

5.2 The data abstraction

A combining operator can be associated with each enumeration method. The absence of such an annotation could be interpreted to mean that a matrix is known not to have duplicate elements. Figure 9 shows the coordinate storage format augmented with an assertion indicating that the combining operator for duplicate elements is +.

template <class T>
class coo_matrix : public enumerable_matrix<T> {
    // ...
public:
#pragma order (0,+) (+,*)
#pragma combining (+)
    iterator enum() { ... };
};

Figure 9: Coordinate storage with a combining operator

5.3 The program transformations

As part of the data-centric transformations, the compiler must ascertain whether the transformed code is safe in the presence of duplicates. This situation is discussed in a little more detail in [8]. Figure 10 shows the result of applying the data-centric transformations to the matrix-multiplication template. The resulting template is annotated to only be instantiated for formats that contain no duplicate entries or whose combining operator is +.

// C = alpha * A * B + beta * C
#pragma A.enum() is combining (none,+)
template<class T>
void matrix_multip(T alpha, enumerable_matrix<T> &A,
                   matrix<T> &B, T beta, matrix<T> &C)
{
    ...
}

Figure 10: Data-centric unassembled matrix multiplication

6 Bounds on stored indices

6.1 The problem

The data-centric transformation requires that, for each matrix entry generated, a bounds check be done to determine that the value is accessed by an iteration of the original control-centric code. In the case of dense matrices and linear loop bounds, these bounds checks can often be eliminated or hoisted using Fourier-Motzkin elimination (FME). However, in the of GP programs written using the abstractions discussed above, there is no way of “pushing” these bounds into a matrix’s enumeration method, because that code may not be present, or may be arbitrarily complex.

Suppose that the forward solve algorithm is specialized (by the user, presumably) for the case when the diagonal of the matrix is assumed to be identity (all 1’s). One such implementation of this algorithm is given in Figure 11. This code can be made data-centric (as in Figure 8), but bounds checks will have to be introduced to ensure that only the matrix entries in the strictly lower triangle of the matrix are used. If a particular storage format is able to convey to the compiler that it only stored the strictly lower triangle of the matrix, then the compiler can safely eliminate these checks.

6.2 The data abstraction

In order to eliminate the bounds checks when they are not needed, a sparse matrix format must provide information about what set of indices it actually stores. That is, while any instance of the abstract matrix class has bounds that are specified as [1 … m, 1 … n], a particular matrix implemenation may store only a predetermined subset of these indices. If this subset can be exposed to the program, then the same FME techniques can be used to optimize or eliminate the bounds tests.

Each matrix implementation can provide information about what subset of the declared bounds are actually stored within the matrix. Figure 12 shows an implementation of the coordinate format with the assertion that it only stores the strictly lower triangle of a matrix.

// X = L^-1 B
// L is strictly lower triangular. A unit diagonal is assumed.
template<class T>
void forward_solve_unitd(matrix<T> &L, matrix<T> &B,
                         matrix<T> &X)
{
    assert(L.rows() == L.cols() && L.cols() == B.rows()
           && B.rows() == X.rows());
    assert(B.cols() == X.cols());
    for (int i = 0; i < B.rows(); i++)
        for (int k = 0; k < B.cols(); k++)
            X.lval(i,k) = B.rval(i,k);
    for (int i=0; i<L.rows(); i++)
        for (int k = 0; k < B.cols(); k++)
            for (int j=0; j<i-1; j++)
                X.lval(i,k) -= L.rval(i,j)*X.rval(j,k);
}

Figure 11: Forward solve with implicit unit diagonal

template <class T>
class coo_matrix : public enumerable_matrix<T> {
    // ...
public:
#pragma order (0,+) (+,*)
#pragma constrained row_index() > col_index()
    iterator enum() { ... };
};

Figure 12: Strictly lower triangular coordinate storage

The difference between a matrix’s declared bounds and stored bounds is subtle. Why, for instance, can’t the declared bounds be eliminated and only the stored bounds used? Consider a storage format in which only the m diagonal elements are stored. In order to implement a generic program for matrix-vector multiplication, it is necessary to know whether the declared bounds of this matrix are (m, m), (m, n), or even (n, m).

6.3 The program transformations

FME and “gisting” ([10]) can be used during the data-centric program transformation to eliminate redundant bounds checks. Figure 13 shows a data-centric version of the code in Figure 11, with the assertion that it only be specialized for formats that store the strictly lower-triangle of a matrix.

// X = L^-1 B
// L is strictly lower triangular. A unit diagonal is assumed.
#pragma L.enum() is ordered (0,+) (+,*)
template<class T>
void forward_solve_unitd(enumerable_matrix<T> &L, matrix<T> &B,
                         matrix<T> &X)
{
    assert(L.rows() == L.cols() && L.cols() == B.rows()
           && B.rows() == X.rows());
    assert(B.cols() == X.cols());
    for (int i = 0; i < B.rows(); i++)
        for (int k = 0; k < B.cols(); k++)
            X.lval(i,k) = B.rval(i,k);
    for (enumerable_matrix<T>::iterator it = L.enum();
         it; it++) {
        int i = it.row_index(), j = it.col_index(); T L_ij = *it;
        for (int k = 0; k < B.cols(); k++)
            X.lval(i,k) -= L.rval(i,j)*X.rval(j,k);
    }
}

Figure 13: Data-centric forward solve with implicit unit diagonal

7 Linear views

7.1 The problem

Some sparse matrix storage formats are most conveniently expressed, not in terms of rows and columns, but in terms of some other set of axes. Diagonallyoriented formats (sparse diagonal, compressed diagonal) are examples of this. For instance, if the entries of a matrix can be enumerated lexicographically, w.r.t. a diagonal orientation of the matrix, there is no consise way of expressing this as distances or directions w.r.t the row and column axes of the matrix. The best that can be expressed is (*,*), which conveys no information.

7.2 The data abstraction

Each matrix implementation can provide a non-singular linear transformation that maps between the row and column indices of the base abstraction and the index space used within the implementation.

A mechanism for mapping between different axes of different matrix representations would solve two problems. First, it would allow storage format implementations to be easily reused with different orientations (e.g., a compressed row (CRS) matrix class could be used to implement a compressed diagonal (CDS) matrix class). Second, it would allow ordering and bounds annotations to be expressed w.r.t. one set of axes and to map to a different set of axes without losing information.

Novak calls these sorts of mappings “views” ([6], [7]).

Views for particular types of mappings can be implemented directly as C++ templates. Figure 14 shows a template for defining non-singular mappings between two sets of axes. Figure 15 shows how this template can define a compressed diagonal storage format in terms of a compressed row storage format.

7.3 The program transformations

The information from an instantiation of a non singular view can easily be incorporated by the compiler into the linear systems that it uses to reason about program transformations. The catch is that the compiler have to know about the non singular view template class and its intended use.

8 Permutation views

8.1 The problem

Non-singular mappings are not the only sort of mappings between index spaces found in sparse matrix computations. Permutation mappings are another.

// some_matrix_type B;
// non_singular_view<some_matrix_type,a,b,c,d> A(B);
//
//  -> (Ai,Aj)' = | a b ; c d | (Bi,Bj)'
template <class M, int a, int b, int c, int d>
class non_singular_view : public enumerable_matrix {
private:
static void apply_xform(int Bi, int Bj, int &Ai, int &Aj) {
    Ai = a*Bi + b*Bj;
    Aj = c*Bi + d*Bj;
}
static void apply_inv_xform(int Ai, int Aj, int &Bi, int &Bj)
    { ... }
M *B;
public:
non_singular_view(M *_B) : B(_B) { }
virtual M::basetype rval(int Ai, int Aj) {
    int Bi, Bj;
    apply_inv_xform(Ai, Aj, Bi, Bj);
    return B->rval(Bi,Bj);
}
virtual M::basetype &lval(int Ai, int Aj) { ... }
protected:
class iterator_state :
    public enumerable_matrix::iterator_state {
private:
    M::iterator it;
    iterator_state() { it = B->enum; }
    virtual int row_index() {
        int Bi = it->row_index(), Bj = it->col_index();
        int Ai, Aj;
        apply_xform(Bi,Bj,Ai,Aj);
        return Ai;
    }
    ...
};
virtual iterator enum() {
    return iterator(new iterator_state());
}
};

Figure 14: Non-singular view class

template <class T>
class cds_matrix :
    public non_singular_view<crs_matrix<T>, 1, -1, 0, 1>
{
public:
    cds_matrix(...) :
        non_singular_view<crs_matrix<T>, 1, -1, 0, 1>(
            new crs_matrix<T>(...))
        { }
};

Figure 15: A compressed diagonal matrix

8.2 The data abstraction

Each matrix implementation can provide either one or two permutation vectors that specify a permutation on either the rows or columns of the matrix.

See Figure 16.

8.3 The program transformations

Program transformations that optimize the use of permutation mappings are discusses in [2] and [9].

9 Nested matrix structures

9.1 The problem

There are many instances where it is advantageous to be able to access a matrix by either its rows, columns, or some other substructure. In order to motivate this, here are three examples,

• Suppose that a matrix can be accessed, first, by enumerating its rows, and, then, by enumerating the non-zero entries within each row (e.g., CRS). An efficient method for implementing the matrix-vector multiplication, y + = A ∗ x, is, for each row, i, load y(i) into a register, update this register with the result of A(i,:) ∗ x(:), and write the register back to y(i). This plan is efficient only when a matrix can provide access to individual rows, A(i,:).

• An efficient method for performing the matrix-matrix multiplication, C += A ∗ B, involved three nested loops, which essentially, find matching nonempty rows, C(i,:) and A(i,:), find the matching row B(k,:) for each nonzero A(i, k), and then find matching non-zero entries, C(i, j) and B(k, j) (see [8] for details). This can only be done efficiently if the matrix formats for A, B, and C provide efficient access to non-empty rows.

// some_matrix_type B;
// row_perm_view<some_matrix_type> A(B,perm);
//
//  -> (Ai,Aj)' = (perm[Bi],Bj)'
template <class M>
class row_perm_view : enumerable_matrix {
private:
vector<int> perm;
vector<int> iperm;
M *B;
public:
row_perm_view(M *_B, vector<int> _perm)
    : B(_B), perm(_perm) {
    for (int i; i<perm.size(); i++)
        iperm[perm[i]] = i;
}
virtual M::basetype rval(int Ai, int Aj) {
    int Bi = iperm[Ai], Bj = Aj;
    return B->rval(Bi,Bj);
}
virtual M::basetype &lval(int Ai, int Aj) { ... }
protected:
class iterator_state :
    public enumerable_matrix::iterator_state  {
private:
    M::iterator it;
    iterator_state() { it = B->enum; }
    virtual int row_index() { return Bi[it->row_index()]; }
    ...
};
virtual iterator enum() {
    return iterator(new iterator_state());
}
};

Figure 16: Row permutation view class • In the BlockSolve library ([3], [4]), a matrix is reordered in such a way that

small, dense (i.e., few or no zero entries) subregions are present. Then the matrix is stored in carefully designed storage format, in which, these dense subregions are stored in small dense matrices. By carefully structuring the algorithms that access this matrix, level 2 and 3 BLAS operations can be exposed, which can greatly improve the performance of these codes.

In all of these cases, access to particular substructures of a matrix (rows, blocks, etc.) is required for the efficient implementation of matrix computations.

9.2 The data abstraction

There are several ways to view such matrix formats. Perhaps the most intuitive is as nested matrix structures. In particular, a 2-d matrix of doubles that provides access to rows can be thought of as a 1-d vector of 1-d vectors of doubles, where the “outer” vector stores each of the rows, and the “inner” vectors stores the non-zeros within each row. Similarly, a sparse, blocked matrix can be thought of as a 2-d sparse array of 2-d dense arrays, where the “outer” array stored each of the non-zero blocks, and the “inner” arrays stored the non-zeros within each block.

Each matrix implementation can provide a description of its nested substructures. Each nesting level should provide its own lookup and enumeration method.

9.3 The program transformations

The query optimization techniques developed for the existing sparse compiler ([11], [8]) can be used to convert a control-centric program into a set of nested “joins” or loops that exploit these structures. Lookup and enumeration methods for the entire matrix can be derived automatically from each of the individual lookup and enumeration methods can also be distilled from this work.

10 Perspective

10.1 The problem

As was mentioned above, it is possible to automatically derive a lookup or enumeration method for the entire matrix, given lookup and enumeration methods for each nesting level. However, there are cases where the derived methods are not as efficient as are required.

One example of this occurs in the Jagged Diagonal (JAD) storage format. In this format, the non-zeros of a matrix are compacted into “diagonals”, which are laid out as long vectors in memory. Enumerating the entries of a single diagonal corresponds to enumerating a set of non-zeros whose row indices are consecutive increasing, but whose column indices are essentially random. In addition, with some slightly complicated address arithmetic, it is possible to access the individual non-zeros within each row of the matrix.

The first approach is the most efficient way to enumerate the non-zeros of the matrix, as it provides unit stride access to memory. For an algorithm like matrix-vector multiplication, in which, there are no restrictions on the order in which entries are enumerated, this is the preferred enumeration method. However, for a code like triangular solve, there are restrictions on the order in which matrix entries are visited. In this case, the second method of enumeration is preferred.

For JAD, it makes sense to provide both a fast enumeration method for the entire matrix, as well as the nested matrix structure with its individual enumeration methods. Then a programmer or sparse compiler could chose the most efficient method on a per-program basis. There are other structures where having multiple ways of enumerating or looking up entries in a matrix make sense. Some storage formats can be accessed equally efficiently by rows or by columns, so providing two nested matrix structures makes sense. We refer to each of these different ways of accessing the same matrix structure as different perspectives on the matrix.

10.2 The data abstraction

Each matrix implementation may provide multiple perspectives, each of which can be an instance of the abstractions enumerated above.

10.3 The program transformations

11 Conclusions

In this note, we present a base abstraction of abstract matrices, with which matrix computations can be written in a GP style. Building upon this base abstraction, we propose the following refinements,

• Enumeration methods

• Distance/direction annotations on enumeration methods.

• Combining operator associated with enumeration methods.

• Bounds on stored indices

• Linear views on matrices

• Permutation views on matrices

• Nested matrix structures

• Perspectives on matrix structures

Here are some issues left unresolved by this note, • What language constructs are appropriate for expressing each of the dif-

ferent abstract levels? In some cases, it seems obvious what languages mechanisms can be used (e.g., matrix implementations as classes, lookup and enumeration methods as methods of each class). In other cases, it is not obvious at all (e.g., linear views, perspective).

• How can one build a compiler in which these and other abstraction mechanisms can be easily incorporated? The refinements presented in this note are simply those that we require in order to express efficient implementations of the algorithm that we are interested in. This list is by no means complete. It would be foolish to define and build a GP system that hard-wired these and did not provide a mechanism with which additional refinements could be defined.

References

[1] Gene H. Golub and Charles F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, MD, 2nd edition, 1989.

[2] John Irwin, Jean-Marc Loingtier, John Gilbert, Gregor Kiczales, John Lamping, Anurag Mendhekar, and Tatiana Shpeisman. Aspect-oriented programming of sparse matrix code. Technical Report SPL97-007 P9710045, Xerox Palo Alto Research Center, February 1997.

http://www.parc.xerox.com/spl/projects/aop/tr-aml.htm.

[3] Mark T. Jones and Paul E. Plassmann. The efficient parallel iterative solution of large sparse linear systems. In A. George, J. Gilbert, and J. W. H. Liu, editors, Graph Theory and Sparse Matrix Computation, volume 56 of IMA Volumes in Mathematics and Its Applications, pages 229–245. Springer-Verlag, 1993. Also MSC Preprint P314-0692.

[4] Mark T. Jones and Paul E. Plassmann. Blocksolve95 users manual: Scalable library software for the parallel solution of sparse linear systems. Technical Report ANL-95/48, Argonne National Laboratory, December 1995.

[5] Induprakas Kodukula, Nawaaz Ahmed, and Keshav Pingali. Datacentric multi-level blocking. In Proceedings of the ACM SIGPLAN ’97 Conference on Programming Language Design and Implementation, Las Vegas, NV, June 16–18, 1997.

[6] Jr. Novak, Gordon S. Creation of views for reuse of software with different data representations. IEEE Transactions on Software Engineering, 21(5):993–1005, December 1995.

http://www.cs.utexas.edu/users/novak/tose95.html.

[7] Jr. Novak, Gordon S. Software reuse by specialization of generic procedures through views. IEEE Transactions on Software Engineering, 23(7), July 1997.

http://www.cs.utexas.edu/users/novak/tose97.html.

[8] Paul Stodghill. A Relational Approach to the Automatic Generation of Sequential Sparse Matrix Codes. PhD thesis, Cornell University, July 1997. also as Technical Report CORNELLCS:TR97-1635.

http://cs-tr.cs.cornell.edu/Dienst/UI/1.0/Display/ncstrl.cornell/TR97-1635.

[9] William Pugh and Tatiana Shpeisman. Generation of the efficient code for sparse matrix computations. submitted to LCPC98, February 1, 1998.

[10] William Pugh. The omega test: a fast and practical integer programming algorithm for dependence analysis. Communications of the ACM, 8:102–114, August 1992.

[11] Vladimir Kotlyar, Keshav Pingali, and Paul Stodghill. A relational approach to the compilation of sparse matrix programs. In EU-ROPAR, 1997.