Add GMRES example

examples/README.txt

@@ -73,6 +73,14 @@ This directory includes three implementations of the algorithm:
   light culling and shading.
 
 
+GMRES
+=====
+
+An implementation of the generalized minimal residual method for solving
+sparse matrix equations.
+(http://en.wikipedia.org/wiki/Generalized_minimal_residual_method)
+
+
 Mandelbrot
 ==========
 

examples/common.mk

@@ -1,16 +1,22 @@
 
 TASK_CXX=../tasksys.cpp
 TASK_LIB=-lpthread
-TASK_OBJ=tasksys.o
+TASK_OBJ=objs/tasksys.o
 
 CXX=g++
 CXXFLAGS=-Iobjs/ -O2 -m64
+CC=gcc
+CCFLAGS=-Iobjs/ -O2 -m64
+
 LIBS=-lm $(TASK_LIB) -lstdc++
 ISPC=ispc -O2 --arch=x86-64 $(ISPC_FLAGS)
 ISPC_OBJS=$(addprefix objs/, $(ISPC_SRC:.ispc=)_ispc.o $(ISPC_SRC:.ispc=)_ispc_sse2.o \
 	$(ISPC_SRC:.ispc=)_ispc_sse4.o $(ISPC_SRC:.ispc=)_ispc_avx.o)
 ISPC_HEADER=objs/$(ISPC_SRC:.ispc=_ispc.h)
-CPP_OBJS=$(addprefix objs/, $(CPP_SRC:.cpp=.o) $(TASK_OBJ))
+
+CPP_OBJS=$(addprefix objs/, $(CPP_SRC:.cpp=.o))
+CC_OBJS=$(addprefix objs/, $(CC_SRC:.c=.o))
+OBJS=$(CPP_OBJS) $(CC_OBJS) $(TASK_OBJ) $(ISPC_OBJS)
 
 default: $(EXAMPLE)
 
@@ -26,12 +32,15 @@ objs/%.cpp objs/%.o objs/%.h: dirs
 clean:
 	/bin/rm -rf objs *~ $(EXAMPLE) $(EXAMPLE)-sse4 $(EXAMPLE)-generic16
 
-$(EXAMPLE): $(CPP_OBJS) $(ISPC_OBJS)
+$(EXAMPLE): $(OBJS)
 	$(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS)
 
 objs/%.o: %.cpp dirs $(ISPC_HEADER)
 	$(CXX) $< $(CXXFLAGS) -c -o $@
 
+objs/%.o: %.c dirs $(ISPC_HEADER)
+	$(CC) $< $(CCFLAGS) -c -o $@
+
 objs/%.o: ../%.cpp dirs
 	$(CXX) $< $(CXXFLAGS) -c -o $@
 

examples/gmres/Makefile

@@ -0,0 +1,8 @@
+
+EXAMPLE=gmres
+CPP_SRC=algorithm.cpp main.cpp matrix.cpp
+CC_SRC=mmio.c
+ISPC_SRC=matrix.ispc
+ISPC_TARGETS=sse2,sse4-x2,avx-x2
+
+include ../common.mk

examples/gmres/algorithm.cpp

@@ -0,0 +1,231 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+/*===========================================================================*\
+|* Includes
+\*===========================================================================*/
+#include "algorithm.h"
+#include "stdio.h"
+#include "debug.h"
+
+
+/*===========================================================================*\
+|* GMRES
+\*===========================================================================*/
+/* upper_triangular_right_solve:
+ * ----------------------------
+ * Given upper triangular matrix R and rhs vector b, solve for
+ * x.  This "solve" ignores the rows, columns of R that are greater than the
+ * dimensions of x.
+ */
+void upper_triangular_right_solve (const DenseMatrix &R, const Vector &b, Vector &x) 
+{
+    // Dimensionality check
+    ASSERT(R.rows() >= b.size());
+    ASSERT(R.cols() >= x.size());
+    ASSERT(b.size() >= x.size());
+
+    int max_row = x.size() - 1;
+
+    // first solve step:
+    x[max_row] = b[max_row] / R(max_row, max_row);
+
+    for (int row = max_row - 1; row >= 0; row--) {
+        double xi = b[row];
+        for (int col = max_row; col > row; col--)
+            xi -= x[col] * R(row, col);
+        x[row] = xi / R(row, row);
+    }
+}
+
+/* create_rotation (used in gmres):
+ * -------------------------------
+ * Construct a Givens rotation to zero out the lowest non-zero entry in a partially
+ * factored Hessenburg matrix.  Note that the previous Givens rotations should be
+ * applied to this column before creating a new rotation.
+ */
+void create_rotation (const DenseMatrix &H, size_t col, Vector &Cn, Vector &Sn) 
+{
+    double a = H(col,     col);
+    double b = H(col + 1, col);
+    double r;
+
+    if (b == 0) {
+        Cn[col] = copysign(1, a);
+        Sn[col] = 0;
+    } 
+    else if (a == 0) {
+        Cn[col] = 0;
+        Sn[col] = copysign(1, b);
+    }
+    else {
+        r       = sqrt(a*a + b*b);
+        Sn[col] = -b / r;
+        Cn[col] =  a / r;
+    }
+}
+
+/* Applies the 'col'th Givens rotation stored in vectors Sn and Cn to the 'col'th 
+ * column of the DenseMatrix M.  (Previous columns don't need the rotation applied b/c
+ * presumeably, the first col-1 columns are already upper triangular, and so their
+ * entries in the col and col+1 rows are 0.)
+ */
+void apply_rotation (DenseMatrix &H, size_t col, Vector &Cn, Vector &Sn) 
+{
+    double c = Cn[col];
+    double s = Sn[col];
+    double tmp    = c * H(col, col) - s * H(col+1, col);
+    H(col+1, col) = s * H(col, col) + c * H(col+1, col);
+    H(col,   col) = tmp;
+}
+
+/* Applies the 'col'th Givens rotation to the vector.
+ */
+void apply_rotation (Vector &v, size_t col, Vector &Cn, Vector &Sn) 
+{
+    double a = v[col];
+    double b = v[col + 1];
+
+    double c = Cn[col];
+    double s = Sn[col];
+
+    v[col]     = c * a - s * b;
+    v[col + 1] = s * a + c * b;
+}
+
+/* Applies the first 'col' Givens rotations to the newly-created column
+ * of H.  (Leaves other columns alone.)
+ */
+void update_column (DenseMatrix &H, size_t col, Vector &Cn, Vector &Sn) 
+{
+    for (int i = 0; i < col; i++) {
+        double c    = Cn[i];
+        double s    = Sn[i];
+        double t    = c * H(i,col) - s * H(i+1,col);
+        H(i+1, col) = s * H(i,col) + c * H(i+1,col);
+        H(i,   col) = t;
+    }
+}
+
+/* After a new column has been added to the hessenburg matrix, factor it back into
+ * an upper-triangular matrix by:
+ * - applying the previous Givens rotations to the new column
+ * - computing the new Givens rotation to make the column upper triangluar
+ * - applying the new Givens rotation to the column, and
+ * - applying the new Givens rotation to the solution vector
+ */
+void update_qr_decomp (DenseMatrix &H, Vector &s, size_t col, Vector &Cn, Vector &Sn)
+{
+    update_column(  H, col, Cn, Sn);
+    create_rotation(H, col, Cn, Sn);
+    apply_rotation( H, col, Cn, Sn);
+    apply_rotation( s, col, Cn, Sn);
+}
+
+void gmres (const Matrix &A, const Vector &b, Vector &x, int num_iters, double max_err)  
+{
+    DEBUG_PRINT("gmres starting!\n");
+    x.zero();
+
+    ASSERT(A.rows() == A.cols());
+    DenseMatrix Qstar(num_iters + 1, A.rows());
+    DenseMatrix H(num_iters + 1, num_iters);
+
+    // arrays for storing parameters of givens rotations
+    Vector Sn(num_iters);
+    Vector Cn(num_iters);
+
+    // array for storing the rhs projected onto the hessenburg's column space
+    Vector G(num_iters+1);
+    G.zero();
+
+    double beta = b.norm();
+    G[0] = beta;
+
+    // temp vector, stores Aqi
+    Vector w(A.rows());
+
+    w.copy(b);
+    w.normalize();
+    Qstar.set_row(0, w);
+
+    int iter = 0;
+    Vector temp(A.rows(), false);
+    double rel_err;
+
+    while (iter < num_iters) 
+    {
+        // w = Aqi
+        Qstar.row(iter, temp);
+        A.multiply(temp, w);
+
+        // construct ith column of H, i+1th row of Qstar:        
+        for (int row = 0; row <= iter; row++) {
+            Qstar.row(row, temp);
+            H(row, iter) = temp.dot(w);
+            w.add_ax(-H(row, iter), temp);
+        }
+
+        H(iter+1, iter) = w.norm();
+        w.divide(H(iter+1, iter));
+        Qstar.set_row(iter+1, w);
+
+        update_qr_decomp (H, G, iter, Cn, Sn);
+
+        rel_err = fabs(G[iter+1] / beta);
+
+        if (rel_err < max_err)
+            break;
+
+        if (iter % 100 == 0)
+            DEBUG_PRINT("Iter %d: %f err\n", iter, rel_err);
+
+        iter++;
+    }
+
+    if (iter == num_iters) {
+        fprintf(stderr, "Error: gmres failed to converge in %d iterations (relative err: %f)\n", num_iters, rel_err);
+        exit(-1);
+    }
+
+    // We've reached an acceptable solution (?):
+
+    DEBUG_PRINT("gmres completed in %d iterations (rel. resid. %f, max %f)\n", num_iters, rel_err, max_err);
+    Vector y(iter+1);
+    upper_triangular_right_solve(H, G, y);
+    for (int i = 0; i < iter + 1; i++) {
+        Qstar.row(i, temp);
+        x.add_ax(y[i], temp);
+    }
+}

examples/gmres/algorithm.h

@@ -0,0 +1,50 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+#ifndef __ALGORITHM_H__
+#define __ALGORITHM_H__
+
+#include "matrix.h"
+
+
+/* Generalized Minimal Residual Method:
+ * -----------------------------------
+ * Takes a square matrix and an rhs and uses GMRES to find an estimate for x.
+ * The specified error is relative.
+ */
+void gmres (const Matrix &A, const Vector &b, Vector &x, int num_iters, double err);
+
+
+
+#endif

examples/gmres/debug.h

@@ -0,0 +1,55 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+#ifndef __DEBUG_H__
+#define __DEBUG_H__
+
+#include <cassert>
+
+
+/**************************************************************\
+| Macros
+\**************************************************************/
+#define DEBUG
+
+#ifdef DEBUG
+#define ASSERT(expr) assert(expr)
+#define DEBUG_PRINT(...) printf(__VA_ARGS__)
+#else
+#define ASSERT(expr)
+#define DEBUG_PRINT(...)
+#endif
+
+
+#endif

examples/gmres/main.cpp

@@ -0,0 +1,79 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+#include "matrix.h"
+#include "algorithm.h"
+#include "util.h"
+#include <cmath>
+#include "../timing.h"
+
+
+int main (int argc, char **argv) 
+{
+    if (argc < 4) {
+        printf("usage: %s <input-matrix> <input-rhs> <output-file>\n", argv[0]);
+        return -1;
+    }
+
+    double gmres_cycles;
+
+    DEBUG_PRINT("Loading A...\n");
+    Matrix *A = CRSMatrix::matrix_from_mtf(argv[1]);
+    if (A == NULL) 
+        return -1;
+    DEBUG_PRINT("... size: %lu\n", A->cols());
+
+    DEBUG_PRINT("Loading b...\n");
+    Vector *b = Vector::vector_from_mtf(argv[2]);
+    if (b == NULL)
+        return -1;
+
+    Vector x(A->cols());
+    DEBUG_PRINT("Beginning gmres...\n");
+    gmres(*A, *b, x, A->cols() / 2, .01);
+
+    // Write result out to file
+    x.to_mtf(argv[argc-1]);
+
+    // Compute residual (double-check)
+#ifdef DEBUG
+    Vector bprime(b->size());
+    A->multiply(x, bprime);
+    Vector resid(bprime.size(), &(bprime[0]));
+    resid.subtract(*b);
+    DEBUG_PRINT("residual error check: %lg\n", resid.norm() / b->norm());
+#endif
+    // Print profiling results
+    DEBUG_PRINT("-- Total mcycles to solve : %.03f --\n", gmres_cycles);
+}

examples/gmres/matrix.cpp

@@ -0,0 +1,246 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+/**************************************************************\
+| Includes
+\**************************************************************/
+#include "matrix.h"
+#include "matrix_ispc.h"
+
+extern "C" {
+#include "mmio.h"
+}
+
+/**************************************************************\
+| DenseMatrix methods
+\**************************************************************/
+void DenseMatrix::multiply (const Vector &v, Vector &r) const 
+{
+    // Dimensionality check
+    ASSERT(v.size() == cols());
+    ASSERT(r.size() == rows());
+
+    for (int i = 0; i < rows(); i++)
+        r[i] = v.dot(entries + i * num_cols);
+}
+
+const Vector *DenseMatrix::row (size_t row) const {
+    return new Vector(num_cols, entries + row * num_cols, true);
+}
+
+void DenseMatrix::row (size_t row, Vector &r) {
+    r.entries = entries + row * cols();
+    r._size   = cols();
+}
+
+void DenseMatrix::set_row(size_t row, const Vector &v) 
+{
+    ASSERT(v.size() == num_cols);
+    memcpy(entries + row * num_cols, v.entries, num_cols * sizeof(double));
+}
+
+
+/**************************************************************\
+| CRSMatrix Methods
+\**************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <vector>
+#include <algorithm>
+
+
+struct entry {
+    int row;
+    int col;
+    double val;
+};
+
+bool compare_entries(struct entry i, struct entry j) {
+    if (i.row < j.row)
+        return true;
+    if (i.row > j.row)
+        return false;
+
+    return i.col < j.col;
+}
+
+#define ERR_OUT(...) { fprintf(stderr, __VA_ARGS__); return NULL; }
+
+CRSMatrix *CRSMatrix::matrix_from_mtf (char *path) {
+    FILE *f;
+    MM_typecode matcode;
+
+    int m, n, nz;
+
+    if ((f = fopen(path, "r")) == NULL) 
+        ERR_OUT("Error: %s does not name a valid/readable file.\n", path);
+
+    if (mm_read_banner(f, &matcode) != 0)
+        ERR_OUT("Error: Could not process Matrix Market banner.\n");
+
+    if (mm_is_complex(matcode)) 
+        ERR_OUT("Error: Application does not support complex numbers.\n")
+
+    if (mm_is_dense(matcode))
+        ERR_OUT("Error: supplied matrix is dense (should be sparse.)\n");
+
+    if (!mm_is_matrix(matcode))
+        ERR_OUT("Error: %s does not encode a matrix.\n", path)
+
+    if (mm_read_mtx_crd_size(f, &m, &n, &nz) != 0)
+        ERR_OUT("Error: could not read matrix size from file.\n");
+
+    if (m != n)
+        ERR_OUT("Error: Application does not support non-square matrices.");
+
+    std::vector<struct entry> entries;
+    entries.resize(nz);
+
+    for (int i = 0; i < nz; i++) {
+        fscanf(f, "%d %d %lg\n", &entries[i].row, &entries[i].col, &entries[i].val);
+        // Adjust from 1-based to 0-based
+        entries[i].row--;
+        entries[i].col--;
+    }
+
+    sort(entries.begin(), entries.end(), compare_entries);
+
+    CRSMatrix *M = new CRSMatrix(m, n, nz);
+    int cur_row = -1;
+    for (int i = 0; i < nz; i++) {
+        while (entries[i].row > cur_row)
+            M->row_offsets[++cur_row] = i;
+        M->entries[i] = entries[i].val;
+        M->columns[i] = entries[i].col;
+    }
+
+    return M;
+}
+
+Vector *Vector::vector_from_mtf (char *path) {
+    FILE *f;
+    MM_typecode matcode;
+
+    int m, n, nz;
+
+    if ((f = fopen(path, "r")) == NULL) 
+        ERR_OUT("Error: %s does not name a valid/readable file.\n", path);
+
+    if (mm_read_banner(f, &matcode) != 0)
+        ERR_OUT("Error: Could not process Matrix Market banner.\n");
+
+    if (mm_is_complex(matcode)) 
+        ERR_OUT("Error: Application does not support complex numbers.\n")
+
+    if (mm_is_dense(matcode)) {
+        if (mm_read_mtx_array_size(f, &m, &n) != 0)
+            ERR_OUT("Error: could not read matrix size from file.\n");
+    } else {
+        if (mm_read_mtx_crd_size(f, &m, &n, &nz) != 0)
+            ERR_OUT("Error: could not read matrix size from file.\n");
+    }
+    if (n != 1)
+        ERR_OUT("Error: %s does not describe a vector.\n", path);
+
+    Vector *x = new Vector(m);
+
+    if (mm_is_dense(matcode)) {
+        double val;
+        for (int i = 0; i < m; i++) {
+            fscanf(f, "%lg\n", &val);
+            (*x)[i] = val;
+        }
+    }
+    else {
+        x->zero();
+        double val;
+        int row;
+        int col;
+        for (int i = 0; i < nz; i++) {
+            fscanf(f, "%d %d %lg\n", &row, &col, &val);
+            (*x)[row-1] = val;
+        }
+    }
+    return x;
+}
+
+#define ERR(...) { fprintf(stderr, __VA_ARGS__); exit(-1); }
+
+void Vector::to_mtf (char *path) {
+    FILE *f;
+    MM_typecode matcode;
+
+    mm_initialize_typecode(&matcode);
+    mm_set_matrix(&matcode);
+    mm_set_real(&matcode);
+    mm_set_dense(&matcode);
+    mm_set_general(&matcode);
+
+    if ((f = fopen(path, "w")) == NULL)
+        ERR("Error: cannot open/write to %s\n", path);
+
+    mm_write_banner(f, matcode);
+    mm_write_mtx_array_size(f, size(), 1);
+    for (int i = 0; i < size(); i++)
+        fprintf(f, "%lg\n", entries[i]);
+
+    fclose(f);
+}
+
+void CRSMatrix::multiply (const Vector &v, Vector &r) const
+{
+    ASSERT(v.size() == cols());
+    ASSERT(r.size() == rows());
+
+    for (int row = 0; row < rows(); row++) 
+    {
+        int row_offset = row_offsets[row];
+        int next_offset = ((row + 1 == rows()) ? _nonzeroes : row_offsets[row + 1]);
+
+        double sum = 0;
+        for (int i = row_offset; i < next_offset; i++)
+        {
+            sum += v[columns[i]] * entries[i];
+        }
+        r[row] = sum;
+    }
+}
+
+void CRSMatrix::zero ( ) 
+{
+    entries.clear();
+    row_offsets.clear();
+    columns.clear();
+    _nonzeroes = 0;
+}

examples/gmres/matrix.h

@@ -0,0 +1,279 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+#ifndef __MATRIX_H__
+#define __MATRIX_H__
+
+/**************************************************************\
+| Includes
+\**************************************************************/
+#include <cstring> // size_t
+#include <cstdlib> // malloc, memcpy, etc.
+#include <cmath>   // sqrt
+#include <vector>
+
+#include "debug.h"
+#include "matrix_ispc.h"
+
+
+class DenseMatrix;
+/**************************************************************\
+| Vector class
+\**************************************************************/
+class Vector {
+ public:
+    static Vector *vector_from_mtf(char *path);
+    void to_mtf (char *path);
+
+    Vector(size_t size, bool alloc_mem=true) 
+        {
+            shared_ptr = false;
+            _size      = size;
+			
+            if (alloc_mem)
+                entries = (double *) malloc(sizeof(double) * _size);
+            else {
+                shared_ptr = true;
+                entries    = NULL;
+            }
+        }
+
+    Vector(size_t size, double *content, bool share_ptr=false) 
+        {
+            _size = size;
+            if (share_ptr) {
+                entries = content;
+                shared_ptr = true;
+            }
+            else {
+                shared_ptr = false;
+                entries = (double *) malloc(sizeof(double) * _size);
+                memcpy(entries, content, sizeof(double) * _size);
+            }
+        }
+
+    ~Vector() { if (!shared_ptr) free(entries); }
+
+    const double & operator [] (size_t index) const 
+    { 
+        ASSERT(index < _size); 
+        return *(entries + index); 
+    }
+
+    double &operator [] (size_t index) 
+    {
+        ASSERT(index < _size);
+        return *(entries + index);
+    }
+
+    bool operator == (const Vector &v) const 
+    {
+        if (v.size() != _size)
+            return false;
+
+        for (int i = 0; i < _size; i++)
+            if (entries[i] != v[i])
+                return false;
+
+        return true;
+    }
+
+    size_t size() const {return _size; }
+
+    double dot (const Vector &b) const 
+    {
+        ASSERT(b.size() == this->size());
+        return ispc::vector_dot(entries, b.entries, size());
+    }
+
+    double dot (const double * const b) const 
+    {
+        return ispc::vector_dot(entries, b, size());
+    }
+
+    void zero () 
+    {
+        ispc::zero(entries, size()); 
+    }
+
+    double norm () const { return sqrtf(dot(entries)); }
+
+    void normalize () { this->divide(this->norm()); }
+
+    void add (const Vector &a) 
+    {
+        ASSERT(size() == a.size());
+        ispc::vector_add(entries, a.entries, size());
+    }
+
+    void subtract (const Vector &s)
+    {
+        ASSERT(size() == s.size());
+        ispc::vector_sub(entries, s.entries, size());
+    }
+
+    void multiply (double scalar) 
+    {
+        ispc::vector_mult(entries, scalar, size());
+    }
+
+    void divide (double scalar) 
+    {
+        ispc::vector_div(entries, scalar, size());
+    }
+
+    // Note: x may be longer than *(this)
+    void add_ax (double a, const Vector &x) {
+        ASSERT(x.size() >= size());
+        ispc::vector_add_ax(entries, a, x.entries, size());
+    }
+
+    // Note that copy only copies the first size() elements of the
+    // supplied vector, i.e. the supplied vector can be longer than
+    // this one.  This is useful in least squares calculations.
+    void copy (const Vector &other) {
+        ASSERT(other.size() >= size());
+        memcpy(entries, other.entries, size() * sizeof(double));
+    }
+
+    friend class DenseMatrix;
+
+ private:
+    size_t  _size;
+    bool     shared_ptr;
+    double  *entries;
+};
+
+
+/**************************************************************\
+| Matrix base class
+\**************************************************************/
+class Matrix {
+    friend class Vector;
+	
+ public:
+    Matrix(size_t size_r, size_t size_c) 
+        { 
+            num_rows = size_r; 
+            num_cols = size_c; 
+        }
+    ~Matrix(){}
+
+    size_t rows() const { return num_rows; }
+    size_t cols() const { return num_cols; }
+
+    virtual void multiply (const Vector &v, Vector &r) const = 0;
+    virtual void zero () = 0;
+
+ protected:
+    size_t num_rows;
+    size_t num_cols;
+};
+
+/**************************************************************\
+| DenseMatrix class
+\**************************************************************/
+class DenseMatrix : public Matrix { 
+    friend class Vector;
+
+ public:
+ DenseMatrix(size_t size_r, size_t size_c) : Matrix(size_r, size_c) 
+        {
+            entries = (double *) malloc(size_r * size_c * sizeof(double));
+        }
+
+ DenseMatrix(size_t size_r, size_t size_c, const double *content) : Matrix (size_r, size_c)
+        {
+            entries = (double *) malloc(size_r * size_c * sizeof(double));
+            memcpy(entries, content, size_r * size_c * sizeof(double));
+        }
+
+    virtual void multiply (const Vector &v, Vector &r) const;
+
+    double &operator () (unsigned int r, unsigned int c)
+    {
+        return *(entries + r * num_cols + c);
+    }
+
+    const double &operator () (unsigned int r, unsigned int c) const
+    {
+        return *(entries + r * num_cols + c);			
+    }
+
+    const Vector *row(size_t row) const;
+    void          row(size_t row, Vector &r);
+    void      set_row(size_t row, const Vector &v);
+
+    virtual void zero() { ispc::zero(entries, rows() * cols()); }
+
+    void copy (const DenseMatrix &other) 
+    {
+        ASSERT(rows() == other.rows());
+        ASSERT(cols() == other.cols());
+        memcpy(entries, other.entries, rows() * cols() * sizeof(double));
+    }
+
+ private:
+    double *entries;
+    bool shared_ptr;
+};
+
+/**************************************************************\
+| CSRMatrix (compressed row storage, a sparse matrix format)
+\**************************************************************/
+class CRSMatrix : public Matrix { 
+ public:
+    CRSMatrix (size_t size_r, size_t size_c, size_t nonzeroes) :
+    Matrix(size_r, size_c) 
+        {
+            _nonzeroes = nonzeroes;
+            entries.resize(nonzeroes);
+            columns.resize(nonzeroes);
+            row_offsets.resize(size_r);
+        }
+
+    virtual void multiply(const Vector &v, Vector &r) const;
+
+    virtual void zero();
+
+    static CRSMatrix *matrix_from_mtf (char *path);
+
+ private:
+    unsigned int        _nonzeroes;
+    std::vector<double>  entries;
+    std::vector<int>     row_offsets;
+    std::vector<int>     columns;
+};
+
+#endif

examples/gmres/matrix.ispc

@@ -0,0 +1,122 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+/**************************************************************\
+| General
+\**************************************************************/
+export void zero (uniform double data[],
+                  uniform int size)
+{
+    foreach (i = 0 ... size)
+        data[i] = 0.0;
+}
+
+
+/**************************************************************\
+| Vector helpers
+\**************************************************************/
+export void vector_add (uniform double a[], 
+                        const uniform double b[], 
+                        const uniform int size) 
+{
+    foreach (i = 0 ... size)
+        a[i] += b[i];
+}
+
+export void vector_sub (uniform double a[], 
+                        const uniform double b[], 
+                        const uniform int size) 
+{
+    foreach (i = 0 ... size)
+        a[i] -= b[i];
+}
+
+export void vector_mult (uniform double a[],
+                         const uniform double b,
+                         const uniform int size)
+{
+    foreach (i = 0 ... size)
+        a[i] *= b;
+}
+
+export void vector_div (uniform double a[],
+                        const uniform double b,
+                        const uniform int size)
+{
+    foreach (i = 0 ... size)
+        a[i] /= b;
+}
+
+export void vector_add_ax (uniform double r[],
+                           const uniform double a,
+                           const uniform double x[],
+                           const uniform int    size)
+{
+    foreach (i = 0 ... size)
+        r[i] += a * x[i];
+}
+
+export uniform double vector_dot (const uniform double a[],
+                                  const uniform double b[],
+                                  const uniform int size)
+{
+    varying double sum = 0.0;
+    foreach (i = 0 ... size)
+        sum += a[i] * b[i];
+    return reduce_add(sum);
+}
+
+/**************************************************************\
+| Matrix helpers
+\**************************************************************/
+export void sparse_multiply (const uniform double entries[],
+                             const uniform double columns[],
+                             const uniform double row_offsets[],
+                             const uniform int rows,
+                             const uniform int cols,
+                             const uniform int nonzeroes,
+                             const uniform double v[],
+                             uniform double r[]) 
+{
+    foreach (row = 0 ... rows) {
+        int row_offset = row_offsets[row];
+        int next_offset = ((row + 1 == rows) ? nonzeroes : row_offsets[row+1]);
+
+        double sum = 0;
+        for (int j = row_offset; j < next_offset; j++)
+            sum += v[columns[j]] * entries[j];
+        r[row] = sum;
+    }
+}
+

examples/gmres/mmio.c

@@ -0,0 +1,511 @@
+/* 
+*   Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <ctype.h>
+
+#include "mmio.h"
+
+int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
+                double **val_, int **I_, int **J_)
+{
+    FILE *f;
+    MM_typecode matcode;
+    int M, N, nz;
+    int i;
+    double *val;
+    int *I, *J;
+ 
+    if ((f = fopen(fname, "r")) == NULL)
+            return -1;
+ 
+ 
+    if (mm_read_banner(f, &matcode) != 0)
+    {
+        printf("mm_read_unsymetric: Could not process Matrix Market banner ");
+        printf(" in file [%s]\n", fname);
+        return -1;
+    }
+ 
+ 
+ 
+    if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
+            mm_is_sparse(matcode)))
+    {
+        fprintf(stderr, "Sorry, this application does not support ");
+        fprintf(stderr, "Market Market type: [%s]\n",
+                mm_typecode_to_str(matcode));
+        return -1;
+    }
+ 
+    /* find out size of sparse matrix: M, N, nz .... */
+ 
+    if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
+    {
+        fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
+        return -1;
+    }
+ 
+    *M_ = M;
+    *N_ = N;
+    *nz_ = nz;
+ 
+    /* reseve memory for matrices */
+ 
+    I = (int *) malloc(nz * sizeof(int));
+    J = (int *) malloc(nz * sizeof(int));
+    val = (double *) malloc(nz * sizeof(double));
+ 
+    *val_ = val;
+    *I_ = I;
+    *J_ = J;
+ 
+    /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
+    /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
+    /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
+ 
+    for (i=0; i<nz; i++)
+    {
+        fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
+        I[i]--;  /* adjust from 1-based to 0-based */
+        J[i]--;
+    }
+    fclose(f);
+ 
+    return 0;
+}
+
+int mm_is_valid(MM_typecode matcode)
+{
+    if (!mm_is_matrix(matcode)) return 0;
+    if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
+    if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
+    if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || 
+                mm_is_skew(matcode))) return 0;
+    return 1;
+}
+
+int mm_read_banner(FILE *f, MM_typecode *matcode)
+{
+    char line[MM_MAX_LINE_LENGTH];
+    char banner[MM_MAX_TOKEN_LENGTH];
+    char mtx[MM_MAX_TOKEN_LENGTH]; 
+    char crd[MM_MAX_TOKEN_LENGTH];
+    char data_type[MM_MAX_TOKEN_LENGTH];
+    char storage_scheme[MM_MAX_TOKEN_LENGTH];
+    char *p;
+
+
+    mm_clear_typecode(matcode);  
+
+    if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) 
+        return MM_PREMATURE_EOF;
+
+    if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, 
+        storage_scheme) != 5)
+        return MM_PREMATURE_EOF;
+
+    for (p=mtx; *p!='\0'; *p=tolower(*p),p++);  /* convert to lower case */
+    for (p=crd; *p!='\0'; *p=tolower(*p),p++);  
+    for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
+    for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
+
+    /* check for banner */
+    if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
+        return MM_NO_HEADER;
+
+    /* first field should be "mtx" */
+    if (strcmp(mtx, MM_MTX_STR) != 0)
+        return  MM_UNSUPPORTED_TYPE;
+    mm_set_matrix(matcode);
+
+
+    /* second field describes whether this is a sparse matrix (in coordinate
+            storgae) or a dense array */
+
+
+    if (strcmp(crd, MM_SPARSE_STR) == 0)
+        mm_set_sparse(matcode);
+    else
+    if (strcmp(crd, MM_DENSE_STR) == 0)
+            mm_set_dense(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+    
+
+    /* third field */
+
+    if (strcmp(data_type, MM_REAL_STR) == 0)
+        mm_set_real(matcode);
+    else
+    if (strcmp(data_type, MM_COMPLEX_STR) == 0)
+        mm_set_complex(matcode);
+    else
+    if (strcmp(data_type, MM_PATTERN_STR) == 0)
+        mm_set_pattern(matcode);
+    else
+    if (strcmp(data_type, MM_INT_STR) == 0)
+        mm_set_integer(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+    
+
+    /* fourth field */
+
+    if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
+        mm_set_general(matcode);
+    else
+    if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
+        mm_set_symmetric(matcode);
+    else
+    if (strcmp(storage_scheme, MM_HERM_STR) == 0)
+        mm_set_hermitian(matcode);
+    else
+    if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
+        mm_set_skew(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+        
+
+    return 0;
+}
+
+int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
+{
+    if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
+        return MM_COULD_NOT_WRITE_FILE;
+    else 
+        return 0;
+}
+
+int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
+{
+    char line[MM_MAX_LINE_LENGTH];
+    int num_items_read;
+
+    /* set return null parameter values, in case we exit with errors */
+    *M = *N = *nz = 0;
+
+    /* now continue scanning until you reach the end-of-comments */
+    do 
+    {
+        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
+            return MM_PREMATURE_EOF;
+    }while (line[0] == '%');
+
+    /* line[] is either blank or has M,N, nz */
+    if (sscanf(line, "%d %d %d", M, N, nz) == 3)
+        return 0;
+        
+    else
+    do
+    { 
+        num_items_read = fscanf(f, "%d %d %d", M, N, nz); 
+        if (num_items_read == EOF) return MM_PREMATURE_EOF;
+    }
+    while (num_items_read != 3);
+
+    return 0;
+}
+
+
+int mm_read_mtx_array_size(FILE *f, int *M, int *N)
+{
+    char line[MM_MAX_LINE_LENGTH];
+    int num_items_read;
+    /* set return null parameter values, in case we exit with errors */
+    *M = *N = 0;
+	
+    /* now continue scanning until you reach the end-of-comments */
+    do 
+    {
+        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
+            return MM_PREMATURE_EOF;
+    }while (line[0] == '%');
+
+    /* line[] is either blank or has M,N, nz */
+    if (sscanf(line, "%d %d", M, N) == 2)
+        return 0;
+        
+    else /* we have a blank line */
+    do
+    { 
+        num_items_read = fscanf(f, "%d %d", M, N); 
+        if (num_items_read == EOF) return MM_PREMATURE_EOF;
+    }
+    while (num_items_read != 2);
+
+    return 0;
+}
+
+int mm_write_mtx_array_size(FILE *f, int M, int N)
+{
+    if (fprintf(f, "%d %d\n", M, N) != 2)
+        return MM_COULD_NOT_WRITE_FILE;
+    else 
+        return 0;
+}
+
+
+
+/*-------------------------------------------------------------------------*/
+
+/******************************************************************/
+/* use when I[], J[], and val[]J, and val[] are already allocated */
+/******************************************************************/
+
+int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
+        double val[], MM_typecode matcode)
+{
+    int i;
+    if (mm_is_complex(matcode))
+    {
+        for (i=0; i<nz; i++)
+            if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
+                != 4) return MM_PREMATURE_EOF;
+    }
+    else if (mm_is_real(matcode))
+    {
+        for (i=0; i<nz; i++)
+        {
+            if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
+                != 3) return MM_PREMATURE_EOF;
+
+        }
+    }
+
+    else if (mm_is_pattern(matcode))
+    {
+        for (i=0; i<nz; i++)
+            if (fscanf(f, "%d %d", &I[i], &J[i])
+                != 2) return MM_PREMATURE_EOF;
+    }
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+    return 0;
+        
+}
+
+int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
+        double *real, double *imag, MM_typecode matcode)
+{
+    if (mm_is_complex(matcode))
+    {
+            if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
+                != 4) return MM_PREMATURE_EOF;
+    }
+    else if (mm_is_real(matcode))
+    {
+            if (fscanf(f, "%d %d %lg\n", I, J, real)
+                != 3) return MM_PREMATURE_EOF;
+
+    }
+
+    else if (mm_is_pattern(matcode))
+    {
+            if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
+    }
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+    return 0;
+        
+}
+
+
+/************************************************************************
+    mm_read_mtx_crd()  fills M, N, nz, array of values, and return
+                        type code, e.g. 'MCRS'
+
+                        if matrix is complex, values[] is of size 2*nz,
+                            (nz pairs of real/imaginary values)
+************************************************************************/
+
+int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, 
+        double **val, MM_typecode *matcode)
+{
+    int ret_code;
+    FILE *f;
+
+    if (strcmp(fname, "stdin") == 0) f=stdin;
+    else
+    if ((f = fopen(fname, "r")) == NULL)
+        return MM_COULD_NOT_READ_FILE;
+
+
+    if ((ret_code = mm_read_banner(f, matcode)) != 0)
+        return ret_code;
+
+    if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && 
+            mm_is_matrix(*matcode)))
+        return MM_UNSUPPORTED_TYPE;
+
+    if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
+        return ret_code;
+
+
+    *I = (int *)  malloc(*nz * sizeof(int));
+    *J = (int *)  malloc(*nz * sizeof(int));
+    *val = NULL;
+
+    if (mm_is_complex(*matcode))
+    {
+        *val = (double *) malloc(*nz * 2 * sizeof(double));
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+    else if (mm_is_real(*matcode))
+    {
+        *val = (double *) malloc(*nz * sizeof(double));
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+
+    else if (mm_is_pattern(*matcode))
+    {
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+
+    if (f != stdin) fclose(f);
+    return 0;
+}
+
+int mm_write_banner(FILE *f, MM_typecode matcode)
+{
+    char *str = mm_typecode_to_str(matcode);
+    int ret_code;
+
+    ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
+    free(str);
+    if (ret_code !=2 )
+        return MM_COULD_NOT_WRITE_FILE;
+    else
+        return 0;
+}
+
+int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
+        double val[], MM_typecode matcode)
+{
+    FILE *f;
+    int i;
+
+    if (strcmp(fname, "stdout") == 0) 
+        f = stdout;
+    else
+    if ((f = fopen(fname, "w")) == NULL)
+        return MM_COULD_NOT_WRITE_FILE;
+    
+    /* print banner followed by typecode */
+    fprintf(f, "%s ", MatrixMarketBanner);
+    fprintf(f, "%s\n", mm_typecode_to_str(matcode));
+
+    /* print matrix sizes and nonzeros */
+    fprintf(f, "%d %d %d\n", M, N, nz);
+
+    /* print values */
+    if (mm_is_pattern(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d\n", I[i], J[i]);
+    else
+    if (mm_is_real(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
+    else
+    if (mm_is_complex(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i], 
+                        val[2*i+1]);
+    else
+    {
+        if (f != stdout) fclose(f);
+        return MM_UNSUPPORTED_TYPE;
+    }
+
+    if (f !=stdout) fclose(f);
+
+    return 0;
+}
+  
+
+/**
+*  Create a new copy of a string s.  mm_strdup() is a common routine, but
+*  not part of ANSI C, so it is included here.  Used by mm_typecode_to_str().
+*
+*/
+char *mm_strdup(const char *s)
+{
+	int len = strlen(s);
+	char *s2 = (char *) malloc((len+1)*sizeof(char));
+	return strcpy(s2, s);
+}
+
+char  *mm_typecode_to_str(MM_typecode matcode)
+{
+    char buffer[MM_MAX_LINE_LENGTH];
+    char *types[4];
+	char *mm_strdup(const char *);
+    int error =0;
+
+    /* check for MTX type */
+    if (mm_is_matrix(matcode)) 
+        types[0] = MM_MTX_STR;
+    else
+        error=1;
+
+    /* check for CRD or ARR matrix */
+    if (mm_is_sparse(matcode))
+        types[1] = MM_SPARSE_STR;
+    else
+    if (mm_is_dense(matcode))
+        types[1] = MM_DENSE_STR;
+    else
+        return NULL;
+
+    /* check for element data type */
+    if (mm_is_real(matcode))
+        types[2] = MM_REAL_STR;
+    else
+    if (mm_is_complex(matcode))
+        types[2] = MM_COMPLEX_STR;
+    else
+    if (mm_is_pattern(matcode))
+        types[2] = MM_PATTERN_STR;
+    else
+    if (mm_is_integer(matcode))
+        types[2] = MM_INT_STR;
+    else
+        return NULL;
+
+
+    /* check for symmetry type */
+    if (mm_is_general(matcode))
+        types[3] = MM_GENERAL_STR;
+    else
+    if (mm_is_symmetric(matcode))
+        types[3] = MM_SYMM_STR;
+    else 
+    if (mm_is_hermitian(matcode))
+        types[3] = MM_HERM_STR;
+    else 
+    if (mm_is_skew(matcode))
+        types[3] = MM_SKEW_STR;
+    else
+        return NULL;
+
+    sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
+    return mm_strdup(buffer);
+
+}

examples/gmres/mmio.h

@@ -0,0 +1,135 @@
+/* 
+*   Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+
+#ifndef MM_IO_H
+#define MM_IO_H
+
+#define MM_MAX_LINE_LENGTH 1025
+#define MatrixMarketBanner "%%MatrixMarket"
+#define MM_MAX_TOKEN_LENGTH 64
+
+typedef char MM_typecode[4];
+
+#include <stdio.h>
+
+char *mm_typecode_to_str(MM_typecode matcode);
+
+int mm_read_banner(FILE *f, MM_typecode *matcode);
+int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz);
+int mm_read_mtx_array_size(FILE *f, int *M, int *N);
+
+int mm_write_banner(FILE *f, MM_typecode matcode);
+int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz);
+int mm_write_mtx_array_size(FILE *f, int M, int N);
+
+
+/********************* MM_typecode query fucntions ***************************/
+
+#define mm_is_matrix(typecode)	((typecode)[0]=='M')
+
+#define mm_is_sparse(typecode)	((typecode)[1]=='C')
+#define mm_is_coordinate(typecode)((typecode)[1]=='C')
+#define mm_is_dense(typecode)	((typecode)[1]=='A')
+#define mm_is_array(typecode)	((typecode)[1]=='A')
+
+#define mm_is_complex(typecode)	((typecode)[2]=='C')
+#define mm_is_real(typecode)		((typecode)[2]=='R')
+#define mm_is_pattern(typecode)	((typecode)[2]=='P')
+#define mm_is_integer(typecode) ((typecode)[2]=='I')
+
+#define mm_is_symmetric(typecode)((typecode)[3]=='S')
+#define mm_is_general(typecode)	((typecode)[3]=='G')
+#define mm_is_skew(typecode)	((typecode)[3]=='K')
+#define mm_is_hermitian(typecode)((typecode)[3]=='H')
+
+int mm_is_valid(MM_typecode matcode);		/* too complex for a macro */
+
+
+/********************* MM_typecode modify fucntions ***************************/
+
+#define mm_set_matrix(typecode)	((*typecode)[0]='M')
+#define mm_set_coordinate(typecode)	((*typecode)[1]='C')
+#define mm_set_array(typecode)	((*typecode)[1]='A')
+#define mm_set_dense(typecode)	mm_set_array(typecode)
+#define mm_set_sparse(typecode)	mm_set_coordinate(typecode)
+
+#define mm_set_complex(typecode)((*typecode)[2]='C')
+#define mm_set_real(typecode)	((*typecode)[2]='R')
+#define mm_set_pattern(typecode)((*typecode)[2]='P')
+#define mm_set_integer(typecode)((*typecode)[2]='I')
+
+
+#define mm_set_symmetric(typecode)((*typecode)[3]='S')
+#define mm_set_general(typecode)((*typecode)[3]='G')
+#define mm_set_skew(typecode)	((*typecode)[3]='K')
+#define mm_set_hermitian(typecode)((*typecode)[3]='H')
+
+#define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
+									(*typecode)[2]=' ',(*typecode)[3]='G')
+
+#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
+
+
+/********************* Matrix Market error codes ***************************/
+
+
+#define MM_COULD_NOT_READ_FILE	11
+#define MM_PREMATURE_EOF		12
+#define MM_NOT_MTX				13
+#define MM_NO_HEADER			14
+#define MM_UNSUPPORTED_TYPE		15
+#define MM_LINE_TOO_LONG		16
+#define MM_COULD_NOT_WRITE_FILE	17
+
+
+/******************** Matrix Market internal definitions ********************
+
+   MM_matrix_typecode: 4-character sequence
+
+				    ojbect 		sparse/   	data        storage 
+						  		dense     	type        scheme
+
+   string position:	 [0]        [1]			[2]         [3]
+
+   Matrix typecode:  M(atrix)  C(oord)		R(eal)   	G(eneral)
+						        A(array)	C(omplex)   H(ermitian)
+											P(attern)   S(ymmetric)
+								    		I(nteger)	K(kew)
+
+ ***********************************************************************/
+
+#define MM_MTX_STR		"matrix"
+#define MM_ARRAY_STR	"array"
+#define MM_DENSE_STR	"array"
+#define MM_COORDINATE_STR "coordinate" 
+#define MM_SPARSE_STR	"coordinate"
+#define MM_COMPLEX_STR	"complex"
+#define MM_REAL_STR		"real"
+#define MM_INT_STR		"integer"
+#define MM_GENERAL_STR  "general"
+#define MM_SYMM_STR		"symmetric"
+#define MM_HERM_STR		"hermitian"
+#define MM_SKEW_STR		"skew-symmetric"
+#define MM_PATTERN_STR  "pattern"
+
+
+/*  high level routines */
+
+int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
+		 double val[], MM_typecode matcode);
+int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
+		double val[], MM_typecode matcode);
+int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *img,
+			MM_typecode matcode);
+
+int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
+                double **val_, int **I_, int **J_);
+
+
+
+#endif

examples/gmres/util.h

@@ -0,0 +1,53 @@
+/*
+  Copyright (c) 2012, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
+*/
+
+
+#ifndef __UTIL_H__
+#define __UTIL_H__
+
+#include <stdio.h>
+#include "matrix.h"
+
+
+inline void printMatrix (DenseMatrix &M, const char *name) {
+    printf("Matrix %s:\n", name);
+    for (int row = 0; row < M.rows(); row++) {
+        printf("row %2d: ", row + 1);
+        for (int col = 0; col < M.cols(); col++)
+            printf("%6f ", M(row, col));
+        printf("\n");
+    }
+    printf("\n");
+}
+
+#endif