From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001
From: Eric Andersen <andersen@codepoet.org>
Date: Thu, 22 Nov 2001 14:04:29 +0000
Subject: Totally rework the math library, this time based on the MacOs X math
 library (which is itself based on the math lib from FreeBSD).  -Erik

---
 libm/double/simq.c | 180 -----------------------------------------------------
 1 file changed, 180 deletions(-)
 delete mode 100644 libm/double/simq.c

(limited to 'libm/double/simq.c')

diff --git a/libm/double/simq.c b/libm/double/simq.c
deleted file mode 100644
index 96d63e521..000000000
--- a/libm/double/simq.c
+++ /dev/null
@@ -1,180 +0,0 @@
-/*							simq.c
- *
- *	Solution of simultaneous linear equations AX = B
- *	by Gaussian elimination with partial pivoting
- *
- *
- *
- * SYNOPSIS:
- *
- * double A[n*n], B[n], X[n];
- * int n, flag;
- * int IPS[];
- * int simq();
- *
- * ercode = simq( A, B, X, n, flag, IPS );
- *
- *
- *
- * DESCRIPTION:
- *
- * B, X, IPS are vectors of length n.
- * A is an n x n matrix (i.e., a vector of length n*n),
- * stored row-wise: that is, A(i,j) = A[ij],
- * where ij = i*n + j, which is the transpose of the normal
- * column-wise storage.
- *
- * The contents of matrix A are destroyed.
- *
- * Set flag=0 to solve.
- * Set flag=-1 to do a new back substitution for different B vector
- * using the same A matrix previously reduced when flag=0.
- *
- * The routine returns nonzero on error; messages are printed.
- *
- *
- * ACCURACY:
- *
- * Depends on the conditioning (range of eigenvalues) of matrix A.
- *
- *
- * REFERENCE:
- *
- * Computer Solution of Linear Algebraic Systems,
- * by George E. Forsythe and Cleve B. Moler; Prentice-Hall, 1967.
- *
- */
-
-/*							simq	2 */
-
-#include <stdio.h>
-#define fabs(x) ((x) < 0 ? -(x) : (x))
-
-int simq( A, B, X, n, flag, IPS )
-double A[], B[], X[];
-int n, flag;
-int IPS[];
-{
-int i, j, ij, ip, ipj, ipk, ipn;
-int idxpiv, iback;
-int k, kp, kp1, kpk, kpn;
-int nip, nkp, nm1;
-double em, q, rownrm, big, size, pivot, sum;
-
-nm1 = n-1;
-if( flag < 0 )
-	goto solve;
-
-/*	Initialize IPS and X	*/
-
-ij=0;
-for( i=0; i<n; i++ )
-	{
-	IPS[i] = i;
-	rownrm = 0.0;
-	for( j=0; j<n; j++ )
-		{
-		q = fabs( A[ij] );
-		if( rownrm < q )
-			rownrm = q;
-		++ij;
-		}
-	if( rownrm == 0.0 )
-		{
-		printf("SIMQ ROWNRM=0");
-		return(1);
-		}
-	X[i] = 1.0/rownrm;
-	}
-
-/*							simq	3 */
-/*	Gaussian elimination with partial pivoting 	*/
-
-for( k=0; k<nm1; k++ )
-	{
-	big= 0.0;
-	idxpiv = 0;
-	for( i=k; i<n; i++ )
-		{
-		ip = IPS[i];
-		ipk = n*ip + k;
-		size = fabs( A[ipk] ) * X[ip];
-		if( size > big )
-			{
-			big = size;
-			idxpiv = i;
-			}
-		}
-
-	if( big == 0.0 )
-		{
-		printf( "SIMQ BIG=0" );
-		return(2);
-		}
-	if( idxpiv != k )
-		{
-		j = IPS[k];
-		IPS[k] = IPS[idxpiv];
-		IPS[idxpiv] = j;
-		}
-	kp = IPS[k];
-	kpk = n*kp + k;
-	pivot = A[kpk];
-	kp1 = k+1;
-	for( i=kp1; i<n; i++ )
-		{
-		ip = IPS[i];
-		ipk = n*ip + k;
-		em = -A[ipk]/pivot;
-		A[ipk] = -em;
-		nip = n*ip;
-		nkp = n*kp;
-		for( j=kp1; j<n; j++ )
-			{
-			ipj = nip + j;
-			A[ipj] = A[ipj] + em * A[nkp + j];
-			}
-		}
-	}
-kpn = n * IPS[n-1] + n - 1;	/* last element of IPS[n] th row */
-if( A[kpn] == 0.0 )
-	{
-	printf( "SIMQ A[kpn]=0");
-	return(3);
-	}
-
-/*							simq 4 */
-/*	back substitution	*/
-
-solve:
-ip = IPS[0];
-X[0] = B[ip];
-for( i=1; i<n; i++ )
-	{
-	ip = IPS[i];
-	ipj = n * ip;
-	sum = 0.0;
-	for( j=0; j<i; j++ )
-		{
-		sum += A[ipj] * X[j];
-		++ipj;
-		}
-	X[i] = B[ip] - sum;
-	}
-
-ipn = n * IPS[n-1] + n - 1;
-X[n-1] = X[n-1]/A[ipn];
-
-for( iback=1; iback<n; iback++ )
-	{
-/* i goes (n-1),...,1	*/
-	i = nm1 - iback;
-	ip = IPS[i];
-	nip = n*ip;
-	sum = 0.0;
-	for( j=i+1; j<n; j++ )
-		sum += A[nip+j] * X[j];
-	X[i] = (X[i] - sum)/A[nip+i];
-	}
-return(0);
-}
-- 
cgit v1.2.3