Changeset 16112

Timestamp:

Jan 17, 2008, 11:03:13 AM (18 years ago)

Author:

eugene

Message:

adding notes from kahan

File:

• trunk/Ohana/src/libohana/src/gaussj.c (modified) (6 diffs)

trunk/Ohana/src/libohana/src/gaussj.c

@@ -1 +1 @@
 # include <ohana.h>
 
-// Gauss-Jordan elimination using full pivots based on Press et al's description.  Major
-// modifications to conform to C indexing, use a boolean to track the completed pivot rows
-// and catch the singular matrix early on.  Also, much cleaner control loops than their
-// implementation.  XXX this really needs to check on round-off errors (see version by
+// Gauss-Jordan elimination using full pivots based on Press et al's description.  Substantially
+// reworked for Ohana: major modifications to conform to C indexing, use a boolean to track the
+// completed pivot rows and catch the singular matrix early on.  Also, much cleaner control loops
+// than their implementation.  XXX this really needs to check on round-off errors (see version by
 // William Kahan
 int dgaussjordan (double **A, double **B, int N, int M) {
@@ -104 +104 @@
   int *rowIndex;
   int *pivot;
-  int maxcol, maxrow;
-  float maxval, tmpval;
   
   int diag, col, row;
@@ -114 +112 @@
   memset (pivot, 0, N*sizeof(int));
 
-  // init these so gcc does not complain
-  maxval = 0.0;
-  maxrow = maxcol = 0;
-
   // we loop along the matrix diagonal, but we do not operate on the diagonal elements in
   // order instead, we are looking for the current max element and operating on that
@@ -125 +119 @@
   for (diag = 0; diag < N; diag++) {
 
-    maxval = 0.0;
+    float maxval = 0.0;
+    int maxrow = 0;
+    int maxcol = 0;
 
     // search for the next pivot
@@ -158 +154 @@
 
     /* rescale by pivot reciprocal */
-    tmpval = 1.0 / A[maxcol][maxcol];
+    float tmpval = 1.0 / A[maxcol][maxcol];
     A[maxcol][maxcol] = 1.0;
     for (col = 0; col < N; col++) A[maxcol][col] *= tmpval;
@@ -191 +187 @@
   return (FALSE);
 }
+
+
+/* Gauss-Jordan Inversion from William Kahan in Basic
+500 ' Gauss-Jordan Matrix Inversion     X = A^(-1) in IBM PC BASIC
+510 ' including checks for excessive    growth despite row-pivoting,
+520 '          and adjustments for zero pivots to avoid .../0 .
+530 ' DIM A(N,N), X(N,N), P(N) ...      are assumed.
+540     DEFINT I-N ' ... integer variables; the rest are REAL.
+550   '
+560   ' First determine levels of roundoff and over/underflow.
+570     UFL = 5.9E-39 ' ... = max{ under, 1/over}flow thresholds.
+580        G=4 : G=G/3 : G=G-1      ' ... = 1/3 + roundoff in 4/3
+590     EPS = ABS( ((G+G) - 1) + G ) ' ... = roundoff level.
+600     G = 1 ' ... will record pivot-growth factor
+610   '
+620   ' Copy A to X and record each column's biggest element.
+630     FOR J=1 TO N : P(J)=0
+640          FOR I=1 TO N : T = A(I,J) : X(I,J) = T : T = ABS(T)
+650                IF T > P(J) THEN P(J) = T
+660                NEXT I : NEXT J
+670   '
+680     FOR K=1 TO N :' ... perform elimination upon column K .
+690          Q=0 : J=K : ' ... search for Kth pivot ...
+700          FOR I=K TO N
+710                T=ABS(X(I,K)) : IF T>Q THEN Q=T : J=I
+720                NEXT I
+730          IF Q=0 THEN Q = EPS*P(K) + UFL : X(K,K)=Q
+740          IF P(K)>0 THEN Q=Q/P(K) : IF Q>G THEN G=Q
+750          IF G<=8*K THEN GOTO 790
+760      PRINT "Growth factor g = ";G;" exceeds ";8*K;" ; try"
+770      PRINT "moving A's column ";K;" to col. 1 to reduce g ."
+780            STOP ' ... or go back to re-order A's columns.
+790         P(K)=J ' ... record pivotal row exchange, if any.
+800         IF J=K THEN GOTO 830 ' ... Don't bother to swap.
+810             FOR L=1 TO N : Q=X(J,L) : X(J,L)=X(K,L)
+820                              X(K,L)=Q : NEXT L
+830         Q = X(K,K) : X(K,K) = 1
+840         FOR J=1 TO N : X(K,J) = X(K,J)/Q : NEXT J
+850         FOR I=1 TO N : IF I=K THEN GOTO 890
+860              Q = X(I,K) : X(I,K) = 0
+870              FOR J=1 TO N
+880                   X(I,J) = X(I,J) - X(K,J)*Q : NEXT J
+890              NEXT I : NEXT K
+900 '
+910 FOR K=N-1 TO 1 STEP -1 ' ... unswap columns of X
+920         J=P(K) : IF J=K THEN GOTO 950
+930         FOR I=1 TO N : Q=X(I,K) : X(I,K)=X(I,J)
+940                          X(I,J)=Q : NEXT I
+950         NEXT K
+960 RETURN
+*/