From 0a968002501102b211f1dcbb97c276743ba24844 Mon Sep 17 00:00:00 2001
From: Adam Szalkowski <adam@ipd.info.uni-karlsruhe.de>
Date: Wed, 31 May 2006 18:39:11 +0000
Subject: [PATCH] added a routine to solve a system of linear equations using
 Gauss-Jordan elimination

[r7843]
---
 ir/adt/Makefile.in   |   5 +-
 ir/adt/gaussjordan.c | 130 +++++++++++++++++++++++++++++++++++++++++++
 ir/adt/gaussjordan.h |  12 ++++
 3 files changed, 145 insertions(+), 2 deletions(-)
 create mode 100644 ir/adt/gaussjordan.c
 create mode 100644 ir/adt/gaussjordan.h

diff --git a/ir/adt/Makefile.in b/ir/adt/Makefile.in
index b9f9178b2..e8efe058d 100644
--- a/ir/adt/Makefile.in
+++ b/ir/adt/Makefile.in
@@ -16,12 +16,13 @@ topdir = ../..
 subdir := ir/adt
 disable_libiberty := @disable_libiberty@
 
-INSTALL_HEADERS_ADT = pset.h set.h pmap.h eset.h hashptr.h array.h pdeq.h iterator.h align.h fourcc.h util.h plist.h bipartite.h xmalloc.h
+INSTALL_HEADERS_ADT = pset.h set.h pmap.h eset.h hashptr.h array.h pdeq.h iterator.h align.h fourcc.h util.h plist.h \
+		      bipartite.h xmalloc.h gaussjordan.h
 
 SOURCES = Makefile.in 			        \
 	array.c obst.h pdeq.c   \
 	set.c pmap.c eset.c iterator.c xmalloc.h \
-	plist.c bipartite.c
+	plist.c bipartite.c gaussjordan.c
 
 ifeq ($(disable_libiberty),no)
 SOURCES += xmalloc.c
diff --git a/ir/adt/gaussjordan.c b/ir/adt/gaussjordan.c
new file mode 100644
index 000000000..52d804032
--- /dev/null
+++ b/ir/adt/gaussjordan.c
@@ -0,0 +1,130 @@
+/**
+ * taken from: http://coastal.er.usgs.gov/rvm/toolkit/rvmlibv2.c
+ * unknown license
+ */
+
+/*------------------------------------------------------*/
+/* gauss.c                                              */
+/*                                                      */
+/* 20 jan 89                                            */
+/*                                                      */
+/* Now does full pivoting--row and column swapping.     */
+/* Requires keeping track of which variable corresponds */
+/* to each vector position.                             */
+/*                                                      */
+/* 18 jan 89                                            */
+/* paul o'neill                                         */
+/*                                                      */
+/* from rob holman's pascal program--geo.pas            */
+/*                                                      */
+/* Gauss-Jordan procedure to solve even-determined      */
+/* square matrix equation Ax = vec,where A is N x N     */
+/* and vec is N x 1.  This is taken roughly from        */
+/* Menke's book (Geophysical Data Analysis: Discrete    */
+/* Inverse Theory--1984 p210 ), but performs actual     */
+/* row switching to simplify the programming.           */
+/* Partial pivoting is used.                            */
+/*                                                      */
+/* A[][] is a square matrix, N x N                      */
+/* vec[] is N x 1 of the matrix                         */
+/* nsize is the size of the equation system             */
+/*                                                      */
+/* returns 0 if successful                              */
+/* returns -1 if ill-conditioned matrix                 */
+/*------------------------------------------------------*/
+
+#include <math.h>
+#include <stdlib.h>
+
+#define SMALL 0.00001
+
+int firm_gaussjordansolve(double *A, double *vec, int nsize)
+{
+  int i, j, row, col, col2, biggest_r, biggest_c;
+  double big, temp, sum;
+  double *x = malloc(nsize * sizeof(double));
+  double *scramvec = malloc(nsize * sizeof(double));
+
+#define _A(row,col) A[(row)*nsize + (col)]
+  /* init x[] */
+  for(i=0;i<nsize;i++)
+    x[i] = i;
+
+  /* triangularize A */
+  /* ie A has zeros below it's diagonal */
+  for(col=0;col<nsize-1;col++) {
+    big = 0;
+    /* find the largest left in LRH box */
+    for(row=col;row<nsize;row++) {
+      for(col2=col;col2<nsize;col2++) {
+	temp = fabs(_A(row,col2));
+	if (temp > big) {
+	  biggest_r = row;
+	  biggest_c = col2;
+	  big = temp; /* largest element left */
+	}
+      }
+    }
+    if (big < SMALL) {
+      return (-1);
+    }
+    /* swap rows */
+    for(i=0;i<nsize;i++) {
+      temp = _A(col,i);
+      _A(col,i) = _A(biggest_r,i);
+      _A(biggest_r,i) = temp;
+    }
+    /* swap vec elements */
+    temp = vec[col];
+    vec[col] = vec[biggest_r];
+    vec[biggest_r] = temp;
+
+    /* swap columns */
+    for(i=0;i<nsize;i++) {
+      temp = _A(i,col);
+      _A(i,col) = _A(i,biggest_c);
+      _A(i,biggest_c) = temp;
+    }
+    /* swap unknowns */
+    temp = x[col];
+    x[col] = x[biggest_c];
+    x[biggest_c] = temp;
+
+    /* partially annihilate this col */
+    /* zero columns below diag */
+    for(row=col+1;row<nsize;row++) {
+
+      /* changes during calc */
+      temp = _A(row,col) / _A(col,col);
+
+      /* annihilates A[][] */
+      for(i=col;i<nsize;i++)
+	_A(row,i) = _A(row,i) - temp * _A(col,i);
+
+      /* same op on vec */
+      vec[row] = vec[row] - temp * vec[col];
+    }
+  }
+  /* back-solve, store solution in srcamvec */
+  scramvec[nsize - 1] = vec[nsize - 1] / _A(nsize - 1,nsize - 1);
+
+  /* answer needs sorting */
+  for(i=nsize-2;i>=0;i--) {
+    sum = 0;
+    for(j=i+1;j<nsize;j++)
+      sum = sum + _A(i,j) * scramvec[j];
+    scramvec[i] = (vec[i] - sum) / _A(i,i);
+  }
+  /* reorder unknowns--return in vec */
+  for(i=0;i<nsize;i++) {
+    j = x[i];
+    vec[j] = scramvec[i];
+  }
+
+  free(x);
+  free(scramvec);
+
+  return (0);
+}
+
+#undef _A
diff --git a/ir/adt/gaussjordan.h b/ir/adt/gaussjordan.h
new file mode 100644
index 000000000..98b0a01a0
--- /dev/null
+++ b/ir/adt/gaussjordan.h
@@ -0,0 +1,12 @@
+#ifndef FIRM_GAUSSJORDAN_H_
+#define FIRM_GAUSSJORDAN_H_
+
+/**
+ * solves a system of linear equations and returns 0 if successful
+ *
+ * @param A    the linear equations as matrix
+ * @param b    the result vector, will contain the result if successful
+ */
+int firm_gaussjordansolve(double *A, double *b, int nsize);
+
+#endif
-- 
2.20.1