Simple Sparse Linear Solver

LUdecomposition.cc

// 
// Sparse LU Decomposition (Gilbert-Peierls)
//
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <functional>
#include <algorithm>
#include <cmath>
#include <map>

using namespace std;

#define ALL(c) c.begin(), c.end()
#define FOR(i,c) for(typeof(c.begin())i=c.begin();i!=c.end();++i)
#define FORR(i,c) for(typeof(c.rbegin())i=c.rbegin();i!=c.rend();++i)
#define REP(i,n) for(int i=0;i<n;++i)
#define fst first
#define snd second

// Simplified Compressed Row Storage
typedef pair<int, double> Entry;
typedef vector<Entry> Row;
typedef vector<Row> Matrix;
const double EPS = 1e-8;

// LU decomposition without reordering (Gilbert-Peierls)
//
// Note: A is LU-decomposable <=> all principal minors are nonsingular
map<int,double> Lsolve(const Matrix &L, const Row &b) {
  map<int, double> x;
  FOR(e, b) x[e->fst] = e->snd;
  FOR(e, x) FOR(l, L[e->fst])
    x[l->fst] -= l->snd * e->snd;
  return x;
}
Row Usolve(Matrix U, map<int, double> b) {
  FORR(e, b) FORR(u, U[e->fst]) 
    if (u->fst == e->fst) e->snd /= u->snd;
    else b[u->fst] -= u->snd * e->snd;
  Row x;
  FOR(e, b) x.push_back(*e);
  return x;
}
pair<Matrix, Matrix> LUdecomposition(Matrix A) {
  int n = A.size();
  Matrix L(n), U(n);
  REP(k, n) {
    map<int, double> s = Lsolve(L, A[k]);
    map<int, double>::iterator j = s.find(k);
    double D = (j++)->snd;
    U[k].assign(s.begin(), j);
    L[k].assign(j, s.end());
    FOR(l, L[k]) l->snd /= D;
  }
  return make_pair(L, U);
}


int main() {
  int n = 100;
  Matrix A(n);
  REP(i, n) {
    if (i == n-1) A[i].push_back( Entry(0,   1) );
    if (i-2 >= 0) A[i].push_back( Entry(i-2, 3) );
    if (i-1 >= 0) A[i].push_back( Entry(i-1, 4) );
                  A[i].push_back( Entry(i  , 5) );
    if (i+1 <  n) A[i].push_back( Entry(i+1, 2) );
    if (i   == 0) A[i].push_back( Entry(n-1, 1) );
  }
  // M = |5 4 3     1|
  //     |2 5 4 3    |
  //     |  2 5 4 3  |
  //     |    2 5 4 3|
  //     |      2 5 4|
  //     |1       2 5|

  Row x;
  REP(i, n) x.push_back( Entry(i, i+1) );
  // x = [1,2,...,n]'

  pair<Matrix,Matrix> LU = LUdecomposition(A);
  x = Usolve( LU.snd, Lsolve( LU.fst, x ) );

  FOR(e, x) printf("%d %lf\n", e->fst, e->snd);
}