;;; -*- Mode: Lisp -*-
;;;
;;; $Header: /home/gene/library/website/docsrc/linear/RCS/numerical-analysis.lisp,v 395.1 2008/04/20 17:25:55 gene Exp $
;;;
;;; Copyright (c) 2007 Gene Michael Stover.  All rights reserved.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU Lesser General Public License as
;;; published by the Free Software Foundation; either version 2 of the
;;; License, or (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public
;;; License along with this program; if not, write to the Free Software
;;; Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
;;; USA
;;;

(defpackage "COM.CYBERTIGGYR.GENE.NUMERICAL-ANALYSIS"
  (:use "COMMON-LISP")
  (:import-from "CYBERTIGGYR-TEST" "DEFTEST")
  (:export "ABSERR"
	   "ABSERRP"
	   "GAUSS0"
	   "GAUSS1"
	   "GAUSS2"
	   "POLY"
	   "RELERR"
	   "RELERRP"))
(in-package "COM.CYBERTIGGYR.GENE.NUMERICAL-ANALYSIS")

;;;
;;; Most (many? all?) of these algorithms are from
;;; /Numerical Analysis/, fourth edition, by Burden & Faires.
;;;
;;; Documentation is at <http://cybertiggyr.com/gene/linear/>.
;;;

(deftest test0000 ()
  "Null test.  Always succeeds."
  'test0000)

(defun dump-matrix (a &optional (strm *standard-output*))
  (declare (type (array real (* *)) a))
  (format strm "~&Dump matrix:")
  (let ((rcount (array-dimension a 0))
	(ccount (array-dimension a 1)))
    (declare (type (integer 0) rcount ccount))
    (do ((row 0 (1+ row)))
	((>= row rcount))
	(declare (type (integer 0) row))
	(format strm "~&[~A]" row)
	(do ((col 0 (1+ col)))
	    ((>= col ccount))
	    (declare (type (integer 0) col))
	    (format strm " ~A" (aref a row col)))))
  strm)

(defun copy-array2 (a)
  "Return a copy of a 2-D array.  The copy will be simple: not
adjustable & no fill pointer."
  (declare (type (array real (* *)) a))
  (do ((b (make-array (array-dimensions a) :element-type 'real
		      :initial-element 0))
       (rcount (array-dimension a 0))
       (r 0 (1+ r))
       (ccount (array-dimension a 1)))
      ((>= r rcount) b)
      (declare (type (array real (* *)) b))
      (do ((c 0 (1+ c)))
	  ((>= c ccount))
	  (setf (aref b r c) (aref a r c)))))

(defun diff-matrix (predicate a b)
  "A & B must have the same dimensions.  Both must be 2-D.  If any
element differs between A & B (corresponding elements, same
coordinates), return their coordinates in a list.  Otherwise, return
NIL.  PREDICATE returns true when elements differ."
  (declare (type function predicate) (type (array real (* *)) a b))
  (assert (equal (array-dimensions a) (array-dimensions b)))
  (let ((coords nil)
	(rcount (array-dimension a 0))
	(ccount (array-dimension a 1)))
    (do ((row 0 (1+ row)))
	((or coords (>= row rcount)) coords)
	(declare (type (integer 0) row))
	(do ((col 0 (1+ col)))
	    ((or coords (>= col ccount)) coords)
	    (declare (type (integer 0) col))
	    (when (funcall predicate (aref a row col) (aref b row col))
	      (setq coords (list row col)))))))

(defun equal-array2 (a b)
  "Return true if & only if A & B are both augmented matrices of the
same size & whose elements are equal."
  (and (typep a '(array real (* *)))
       (typep b '(array real (* *)))
       (equal (array-dimensions a) (array-dimensions b))
       (let ((a0 a) (b0 b))
	 (declare (type (array real (* *)) a0 b0))
	 (do ((nrow (array-dimension a0 0))
	      (ncol (array-dimension a0 1))
	      (i 0 (1+ i))
	      (found-difference nil))
	     ((or found-difference (>= i nrow)) (not found-difference))
	     (do ((j 0 (1+ j)))
		 ((or found-difference (>= j ncol)))
		 (unless (eql (aref a0 i j) (aref b0 i j))
		   (setf found-difference t)))))))

(defun abserr (y y0)
  "Return the absolute error in Y as an approximation of Y0."
  (declare (type real y y0))
  (abs (- y y0)))

(deftest test0010 () (zerop (abserr 0 0.0)))
(deftest test0011 () (eql (abserr  1  2) 1))
(deftest test0012 () (eql (abserr -1  2) 3))
(deftest test0013 () (eql (abserr -1 -2) 1))
(deftest test0014 () (eql (abserr  1 -2) 3))

(defun abserrp (y y0 range)
  "True if & only if the absolute error of Y as an approximation
of Y0 is less than or equal to RANGE."
  (declare (type real y y0 range))
  (<= (abserr y y0) range))

(deftest test0021 () (abserrp 0 0 0))
(deftest test0022 () (abserrp 0 0 10))
(deftest test0023 () (abserrp 2 3 1))
(deftest test0024 () (abserrp 2 3 2))
(deftest test0025 () (not (abserrp 2 3 0.5)))

(defun relerr (y y0)
  "Return the relative error in Y as an approximation of Y0.
The relative error on Y0 = 0 is 0 if Y = 0 or the absolute value
of Y otherwise."
  (declare (type real y y0))
  (abs (/ (- y y0) (if (zerop y0) 1 y0))))

(deftest test0031 () (zerop (relerr 0 0)))
(deftest test0032 () (eql (relerr 1 0) 1))
(deftest test0033 () (zerop (relerr 1 1)))
(deftest test0034 () (zerop (relerr 1 1.0)))
(deftest test0035 () (eql (relerr 1 2) 1/2))
(deftest test0036 () (eql (relerr 1.0 2.0) 0.5))

(defun relerrp (y y0 range)
  "True if & only if the relative error of Y as an approximation
of Y0 is less than or equal to RANGE."
  (declare (type real y y0 range))
  (<= (relerr y y0) range))

(deftest test0041 () (relerrp 0 0 0))
(deftest test0042 () (not (relerrp 1 0 0)))
(deftest test0043 () (relerrp 2 3 1))
(deftest test0044 () (relerrp 2 3 1/3))
(deftest test0045 () (not (relerrp 2 3 1/4)))
(deftest test0046 () (not (relerrp 2 3 0.25)))
(deftest test0047 () (relerrp 2 3 0.34))

(defun is-augmented-matrix (x)
  "It's an augmented matrix if it's a 2-d array whose number of
rows N is one less than the number of columns N+1."
  (and (typep x '(array real (* *)))
       (= (1+ (array-dimension x 0)) (array-dimension x 1))))

(deftest test0101 ()
  "Verify that NIL is not an augmented matrix."
  (not (is-augmented-matrix nil)))

(deftest test0102 ()
  "Test that a 2-by-3 matrix is an augmented matrix.  We use defaults
for adjustable & fill pointer."
  (is-augmented-matrix (make-array '(2 3) :element-type 'real)))

(deftest test0103 ()
  "Test that a 2-by-3 matrix, which we explicitly say is not
adjustable & doesn't have a fill pointer, is an augmented
matrix."
  (is-augmented-matrix
   (make-array '(2 3) :element-type 'real :adjustable nil
	       :fill-pointer nil)))

(deftest test0104 ()
  "Test that a 2-by-3 matrix, which we explicitly say is adjustable
but doesn't have a fill pointer, is an augmented matrix."
  (is-augmented-matrix (make-array '(2 3) :element-type 'real :adjustable t
				   :fill-pointer nil)))

(deftest test0105 ()
  "Test that a 2-by-2 matrix is not an augmented matrix."
  (not (is-augmented-matrix (make-array '(2 2) :element-type 'real))))

(deftest test0106 ()
  "Test that a 2-by-4 matrix (one too many columns) is not an
augmented matrix."
  (not (is-augmented-matrix (make-array '(2 4) :element-type 'real))))

(deftest test0107 ()
  "Test that a 2-by-10 matrix (way too many columns) is not an
augmented matrix."
  (not (is-augmented-matrix (make-array '(2 10) :element-type 'real))))

(deftest test0108 ()
  "Test that a 10-by-3 matrix (way too many rows) is not an
augmented matrix."
  (not (is-augmented-matrix (make-array '(10 3) :element-type 'real))))

(defun swap-rows! (i j a)
  "Swap rows I & J in matrix A.  Modifies A.  Returns A.  It is not
an error if I = J."
  (declare (type (integer 0) i j) (type (array t (* *)) a))
  (assert (< i (array-dimension a 0)))
  (assert (< j (array-dimension a 0)))
  (unless (= i j)
    (do ((columns (array-dimension a 1))
	 (c 0 (1+ c)))
	((>= c columns))
	(declare (type (integer 0) c))
	(rotatef (aref a i c) (aref a j c))))
  a)

(deftest test0121 ()
  "Test SWAP-ROWS! on a hard-coded 2-by-3 matrix.  The rows are not
the same row."
  (let ((a (make-array '(2 3) :element-type 'real)))
    (declare (type (simple-array real (2 3)) a))
    (setf (aref a 0 0) 0
	  (aref a 0 1) 1
	  (aref a 0 2) 2
	  (aref a 1 0) 3
	  (aref a 1 1) 4
	  (aref a 1 2) 5)
    (swap-rows! 0 1 a)
    (and (= 3 (aref a 0 0))
	 (= 4 (aref a 0 1))
	 (= 5 (aref a 0 2))
	 (= 0 (aref a 1 0))
	 (= 1 (aref a 1 1))
	 (= 2 (aref a 1 2)))))

(deftest test0122 ()
  "Test SWAP-ROWS! on a hard-coded 2-by-3 matrix.  We swap the same
row.  Ideally, SWAP-ROWS! will do nothing, but the main thing is
that the resulting matrix is correct."
  (let ((a (make-array '(2 3) :element-type 'real)))
    (declare (type (simple-array real (2 3)) a))
    (setf (aref a 0 0) 0
	  (aref a 0 1) 1
	  (aref a 0 2) 2
	  (aref a 1 0) 3
	  (aref a 1 1) 4
	  (aref a 1 2) 5)
    (swap-rows! 0 0 a)
    (and (= 0 (aref a 0 0))
	 (= 1 (aref a 0 1))
	 (= 2 (aref a 0 2))
	 (= 3 (aref a 1 0))
	 (= 4 (aref a 1 1))
	 (= 5 (aref a 1 2)))))

(deftest test0123 ()
  "Test SWAP-ROWS! on a hard-coded 3-by-2 matrix.  We swap the
first & last rows."
  (let ((a (make-array '(3 2) :element-type 'real)))
    (declare (type (simple-array real (3 2)) a))
    (setf (aref a 0 0) 0
	  (aref a 0 1) 1
	  (aref a 1 0) 2
	  (aref a 1 1) 3
	  (aref a 2 0) 4
	  (aref a 2 1) 5)
    (swap-rows! 0 2 a)
    (and (= (aref a 0 0) 4)
	 (= (aref a 0 1) 5)
	 (= (aref a 1 0) 2)
	 (= (aref a 1 1) 3)
	 (= (aref a 2 0) 0)
	 (= (aref a 2 1) 1))))

(defun gauss0-select-q (i a)
  "Corresponds to step 2 in algorithm 6.1.  Instead of calling the
value P as in algorithm 6.1, we call it Q.  Return NIL if we can
not find an adequate Q."
  (declare (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  (do ((n (array-dimension a 0))        ; number of rows
       (q i (1+ q)))
      ((or (>= q n) (not (zerop (aref a q i))))
       (if (< q n) q -1))))

(deftest test0151 ()
  "Test that GAUSS0-SELECT-Q works for this hard-coded, 2-by-3
input set:
  ( 0 1 2
    3 4 5 )
where I = 0.  The correct output value is 1, the subscript of
the second row, because the element at (0,I) is zero, & the
element at (1,I) is the first non-zero element if we scan from
the first row to the last."
  (let ((a (make-array '(2 3) :element-type 'real)))
    (declare (type (simple-array real (* *)) a))
    (setf (aref a 0 0) 0
	  (aref a 0 1) 1
	  (aref a 0 2) 2
	  (aref a 1 0) 3
	  (aref a 1 1) 4
	  (aref a 1 2) 5)
    (= 1 (gauss0-select-q 0 a))))

(deftest test0152 ()
  "Like TEST0151 but I = 1.  We expect the same result."
  (let ((a (make-array '(2 3) :element-type 'real)))
    (declare (type (simple-array real (* *)) a))
    (setf (aref a 0 0) 0
	  (aref a 0 1) 1
	  (aref a 0 2) 2
	  (aref a 1 0) 3
	  (aref a 1 1) 4
	  (aref a 1 2) 5)
    (= 1 (gauss0-select-q 1 a))))

(deftest test0153 ()
  "Test that GAUSS0-SELECT-Q returns -1 for this hard-coded, 2-by-3
input set:
  ( 0 1 2
    0 4 5 )
where I = 0.  The correct output value is -1 because A[J][I] = 0
for J in (0, 1)."
  (let ((a (make-array '(2 3) :element-type 'real)))
    (declare (type (simple-array real (* *)) a))
    (setf (aref a 0 0) 0
	  (aref a 0 1) 1
	  (aref a 0 2) 2
	  (aref a 1 0) 0
	  (aref a 1 1) 4
	  (aref a 1 2) 5)
    (let ((x (gauss0-select-q 0 a))
	  (expect -1))
      (unless (= expect x)
	(format t "~&~A: ~A returned ~A.  Expected ~A." 'test0153
		'gauss0-select-q x expect))
      (= expect x))))

(defun gauss0-update-rows! (i a)
  "Steps 4, 5, & 6 from Algorithm 6.1  Multiplies all the
rows after I by M, & subtracts.  I dunno.  Not descriptive.  Just
see algorithm 6.1, stpes 4, 5, & 6 & deal with it.  Fuck it."
  (declare (type (integer 0) i) (type (array real (* *))))
  (assert (is-augmented-matrix a))
  (let ((n (array-dimension a 0))
	(columns (array-dimension a 1)))
    (declare (type (integer 0) n))
    (when (< i (1- n))
      (do ((j (1+ i) (1+ j)))
	  ((>= j n))
	  (declare (type (integer 0) j))
	  (let ((m (/ (aref a j i) (aref a i i))))
	    (declare (type real m))
	    (do ((k 0 (1+ k)))
		((>= k columns))
		(declare (type (integer 0) k))
		(decf (aref a j k) (* m (aref a i k))))))))
  a)

(defun gauss0-eliminate! (a)
  "Perform the \"elimination process\" for Gaussian Elimination with
Backward Substituion, algorithm 6.1.  Modifies the matrix A."
  (do* ((n (array-dimension a 0))
	(i 0 (1+ i))
	(q (gauss0-select-q 0 a) (gauss0-select-q i a)))
       ((or (>= i n) (= -1 q))
	(if (= -1 q) nil a))
       (when (/= q i) (swap-rows! q i a))
       (setq a (gauss0-update-rows! i a)))
  a)

(deftest test0171 ()
  "Test GAUSS0-ELIMINATE! on a hard-coded augmented matrix from
page 319 of the Burden & Faires book."
  (let ((a (make-array '(4 5) :element-type 'real))
	(expect (make-array '(4 5) :element-type 'real)))
    (declare (type (simple-array real (4 5)) a expect))
    (setf (aref a 0 0)   1	      ; Input matrix.  Hard-coded from
	  (aref a 0 1)  -1              ; page 319 of the book.
	  (aref a 0 2)   2
	  (aref a 0 3)  -1
	  (aref a 0 4)  -8
	  (aref a 1 0)   2
	  (aref a 1 1)  -2
	  (aref a 1 2)   3
	  (aref a 1 3)  -3
	  (aref a 1 4) -20
	  (aref a 2 0)   1
	  (aref a 2 1)   1
	  (aref a 2 2)   1
	  (aref a 2 3)   0
	  (aref a 2 4)  -2
	  (aref a 3 0)   1
	  (aref a 3 1)  -1
	  (aref a 3 2)   4
	  (aref a 3 3)   3
	  (aref a 3 4)   4)
    (setf (aref expect 0 0)    1
	  (aref expect 0 1)   -1
	  (aref expect 0 2)    2
	  (aref expect 0 3)   -1
	  (aref expect 0 4)   -8
	  (aref expect 1 0)    0
	  (aref expect 1 1)    2
	  (aref expect 1 2)   -1
	  (aref expect 1 3)    1
	  (aref expect 1 4)    6
	  (aref expect 2 0)    0
	  (aref expect 2 1)    0
	  (aref expect 2 2)   -1
	  (aref expect 2 3)   -1
	  (aref expect 2 4)   -4
	  (aref expect 3 0)    0
	  (aref expect 3 1)    0
	  (aref expect 3 2)    0
	  (aref expect 3 3)    2
	  (aref expect 3 4)    4)
    (setq a (gauss0-eliminate! a))
    (let ((diff (diff-matrix #'/= a expect)))
      (when diff
	(format t "~&~A:" 'test0171)
	(format t " Matrices differ at ~A. " diff)
	(format t " A's element is ~A. " (apply #'aref a diff))
	(format t " EXPECT's element is ~A." (apply #'aref expect diff))
	(format t "~&Matrix A is")
	(dump-matrix a)
	(format t "~&Matrix EXPECT is ")
	(dump-matrix expect))
      (not diff))))

(defun gauss0-sum (i a x)
  (declare (type (integer 0) i) (type (array real (* *)) a)
	   (type (array real (*)) x))
  (assert (is-augmented-matrix a))
  (do ((n (array-dimension a 0))
       (sum 0 (+ sum (* (aref a i j) (aref x j))))
       (j (1+ i) (1+ j)))
      ((>= j n) sum)
      (declare (type (integer 0) n) (type real sum) (type (integer 0) j))))

(defun gauss0-substitute-backward (a)
  (declare (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  (let* ((n (array-dimension a 0))
	 (n1 (1- n))
	 (x (make-array n :element-type 'real :initial-element 0)))
    (declare (type (integer 0) n) (type (array real (*)) x))
    (setf (aref x n1) (/ (aref a n1 n) (aref a n1 n1)))
    (do ((i (1- n1) (1- i)))
	((< i 0))
	(declare (type integer i))
	(setf (aref x i) (/ (- (aref a i n) (gauss0-sum i a x))
			    (aref a i i))))
    x))

(defun gauss0 (a)
  "Return a vector which approximates the solutions to the linear
system in A.  A is an ``augmented matrix'' according to the book,
page 321.  Uses Gaussian Elimination with Backward Substitution.
That's algorithm 6.1 (page 321) in the book.  Return NIL if there
is no unique solution."
  (declare (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  ;; Our N is the subscript of the last row.
  (let ((n (1- (array-dimension a 0))))
    ;; Easiest way to ensure we don't modify A is to copy it.  We'll
    ;; bind the copy to A again.
    (setq a (copy-array2 a))
    (setq a (gauss0-eliminate! a))
    (cond ((null a) nil)                ; No unique solution
	  ((zerop (aref a n n))         ; No unique solution
	   ;; In algorithm 6.1, this is step 7: a[n][n] = 0.
	   ;; There is no unique solution.
	   nil)
	  (t (gauss0-substitute-backward a)))))

(deftest test0201 ()
  "Test GAUSS0 on a hard-coded input.  The input is from page
319 of the book."
  (let ((a (make-array '(4 5) :element-type 'real)))
    (declare (type (simple-array real (4 5)) a))
    (setf (aref a 0 0)   1	      ; Input matrix.  Hard-coded from
	  (aref a 0 1)  -1              ; page 319 of the book.
	  (aref a 0 2)   2
	  (aref a 0 3)  -1
	  (aref a 0 4)  -8
	  (aref a 1 0)   2
	  (aref a 1 1)  -2
	  (aref a 1 2)   3
	  (aref a 1 3)  -3
	  (aref a 1 4) -20
	  (aref a 2 0)   1
	  (aref a 2 1)   1
	  (aref a 2 2)   1
	  (aref a 2 3)   0
	  (aref a 2 4)  -2
	  (aref a 3 0)   1
	  (aref a 3 1)  -1
	  (aref a 3 2)   4
	  (aref a 3 3)   3
	  (aref a 3 4)   4)
    (let ((x (gauss0 a)))
      (declare (type (simple-array real (4)) x))
      (let ((is-good (and (= -7 (aref x 0))
			  (=  3 (aref x 1))
			  (=  2 (aref x 2))
			  (=  2 (aref x 3)))))
	(unless is-good
	  (format t "~&~A:" 'test0201)
	  (format t " X is ~A. " x)
	  (format t " Expected 2, 2, 3, -7."))
	is-good))))

(deftest test0210 (&optional (n 100))
  "Construct a large linear system, use GAUSS0 to solve it, check
the result."
  (declare (type (integer 0) n))
  (let* ((x (make-array n :element-type 'real :initial-element 0))
	 (a (make-array (list n (1+ n)) :element-type 'real
			:initial-element 0)))
    (declare (type (simple-array real (*)) x)
	     (type (simple-array real (* *)) a))
    (do ((i 0 (1+ i)))
	((>= i n))
	(setf (aref x i) (1+ (random 100.0)))
	(do ((j 0 (1+ j)))
	    ((>= j n))
	    (setf (aref a i j) (1+ (random 100.0)))))
    (do ((i 0 (1+ i)))
	((>= i n))
	(setf (aref a i n)
	      (loop for j from 0 to (1- n) sum (* (aref a i j) (aref x j)))))
    (let ((y (gauss0 a)))
      (declare (type (array real (*)) y))
      (let ((is-good t))
	(do ((i 0 (1+ i)))
	    ((>= i n))
	    (unless (relerrp (aref y i) (aref x i) 5/100)
	      (format t "~&X[~A] is ~A.  Y[~A] is ~A. "
		      i (aref x i) i (aref y i))
	      (format t " Relative error is ~A." (relerr (aref y i)
							 (aref x i)))
	      (setf is-good nil)))
	is-good))))

(defun make-nrow (a)
  "Return an NROW vector as used by the Gaussian Elimination functions which
simulate row swaps by using indirection through NROW.  NROW has as
many elements as A has rows.  NROW's values are the subscripts of A's
rows.  In other words, if A has N rows, then NROW will have 0, 1, ...
N-1."
  (declare (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  (let ((n (array-dimension a 0)))
    (do ((nrow (make-array n :element-type '(integer 0)))
	 (i 0 (1+ i)))
	((>= i n) nrow)
	(declare (type (array (integer 0) (*)) nrow) (type (integer 0) i)))))
fixme

(deftest test0221 ()
  "Test MAKE-NROW on an matrix of 1 row.  The matrix has two columns
because we are using \"augmented matrices\"; they have N rows & N+1
columns."
  (let ((input (make-array '(1 2) :element-type 'real :initial-element 0))
	(expect (make-array 1 :element-type '(integer 0)
			    :initial-contents '(0))))
    (let ((x (make-nrow input)))
      (equal-array2 x expect))))

(defun gaussn-eliminate! (a select-q)
  "Perform the \"elimination process\" for Gaussian Elimination with
Backward Substituion, algorithm 6.1.  Modifies the matrix A."
  (declare (type (array real (* *)) a) (type function select-q))
  (assert (is-augmented-matrix a))
  (do* ((n (array-dimension a 0))
	(i 0 (1+ i))
	(q (gauss0-select-q 0 a) (gauss0-select-q i a)))
       ((or (>= i n) (= -1 q))
	(if (= -1 q) nil a))
       (when (/= q i) (swap-rows! q i a))
       (setq a (gauss0-update-rows! i a)))
  a)

(defun gaussn (a select-q)
  "Perform Gaussian Elimination in which the SELECT-Q function is a
parameter.  This means that whether you get no column pivoting,
maximal column pivoting, or scaled column pivoting depends on the
SELECT-Q function you provide."
  (declare (type (array real (* *)) a) (type function select-q))
  (assert (is-augmented-matrix a))
  ;; Our N is the subscript of the last row.
  (let ((n (1- (array-dimension a 0)))
	(nrow (make-nrow a)))
    ;; Easiest way to ensure we don't modify A is to copy it.  We'll
    ;; bind the copy to A again.
    (setq a (copy-array2 a))
    (setq a (gaussn-eliminate! a select-q))
    (cond ((null a) nil)                ; No unique solution
	  ((zerop (aref a n n))         ; No unique solution
	   ;; In algorithm 6.1, this is step 7: a[n][n] = 0.
	   ;; There is no unique solution.
	   nil)
	  (t (gaussn-substitute-backward a)))))

(defun gauss1 (a)
  "Like GAUSS0 in that it performas Gaussian Elimination with Backward
Substitution, but implemented with the more generic functions which we
use for the later to Gaussian techniques."
  (declare (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  (gaussn a #'gauss1-select-q))

(defun gauss2-absmax-row (i nrow a)
  "Return an index into NROW.  This is the predicate from step 3
of Algorithm 6.2 on page 330.  You'll have to see the book
understand this predicate.  Sorry."
  (declare (type (integer 0) i) (type (array (integer 0) (*)) nrow)
	   (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  (reduce #'max
	  (do ((lst nil (cons (abs (aref a (aref nrow j) i)) lst))
	       (n (array-dimension a 0)) ; number of rows
	       (j i (1+ j)))
	      ((>= j n) lst)
	      (declare (type list lst) (type (integer 0) n j)))))

(deftest test0215 ()
  "Test GAUSS2-ABSMAX-ROW with a hard-coded input.  The test
array has only one row, so the result should be 0."
  (let ((a (make-array '(1 2) :element-type 'real))
	(nrow (make-array 1 :element-type '(integer 0)))
	(expect 0))
    (setf (aref a 0 0) -1
	  (aref a 0 1) -1)
    (setf (aref nrow 0) 0)
    (let ((x (gauss2-absmax-row 0 nrow a)))
      (unless (eql x expect)
	(format t "~&~A: ~A returned ~A.  Expected ~A." 'test0215
		'gauss2-absmax-row x expect))
      (eql x expect))))

(defun gauss2-select-q (i a)
  "Corresponds to step 2 in algorithm 6.1.  Instead of calling the
value P as in algorithm 6.1, we call it Q.  Return NIL if we can
not find an adequate Q."
  (declare (type (integer 0) i) (type (array real (* *)) a))
  (assert (is-augmented-matrix a))
  (do ((n (array-dimension a 0))        ; number of rows
       (q i (1+ q)))
      ((or (>= q n) (not (zerop (aref a q i))))
       (if (< q n) q -1))
      (declare (type (integer 0) n q))))

(defun gauss2-eliminate! (a)
  "Perform the \"elimination process\" for Gaussian Elimination with
Backward Substituion, algorithm 6.1.  Modifies the matrix A."
;; fixme
  (gauss0-eliminate! a))

(defun gauss2 (a)
  "Return a vector which approximates the solutions to the linear
system in A.  A is an ``augmented matrix'' according to the book,
page 330.  Uses Gaussian Elimination with Maximal Column Pivoting.
That's algorithm 6.2 (page 330) in the book.  Return NIL if there
is no unique solution.
    Very similar to GAUSS0 except that it uses a vector to
simulate row exchanges, & of course it has a different way of
selecting the row Q to exchange."
  ;; Our N is the subscript of the last row.
  (let ((n (1- (array-dimension a 0))))
    ;; Easiest way to ensure we don't modify A is to copy it.  We'll
    ;; bind the copy to A again.
    (setq a (copy-array2 a))
    (setq a (gauss2-eliminate! a))
    (cond ((null a) nil)                ; No unique solution
	  ((zerop (aref a n n))         ; No unique solution
	   ;; In algorithm 6.1, this is step 7: a[n][n] = 0.
	   ;; There is no unique solution.
	   nil)
	  (t (gauss0-substitute-backward a)))))

(deftest test0251 ()
  "Test GAUSS2 on a hard-coded input.  The input is from page
319 of the book."
  (let ((a (make-array '(4 5) :element-type 'real)))
    (declare (type (simple-array real (4 5)) a))
    (setf (aref a 0 0)   1	      ; Input matrix.  Hard-coded from
	  (aref a 0 1)  -1              ; page 319 of the book.
	  (aref a 0 2)   2
	  (aref a 0 3)  -1
	  (aref a 0 4)  -8
	  (aref a 1 0)   2
	  (aref a 1 1)  -2
	  (aref a 1 2)   3
	  (aref a 1 3)  -3
	  (aref a 1 4) -20
	  (aref a 2 0)   1
	  (aref a 2 1)   1
	  (aref a 2 2)   1
	  (aref a 2 3)   0
	  (aref a 2 4)  -2
	  (aref a 3 0)   1
	  (aref a 3 1)  -1
	  (aref a 3 2)   4
	  (aref a 3 3)   3
	  (aref a 3 4)   4)
    (let ((x (gauss2 a)))
      (declare (type (simple-array real (4)) x))
      (let ((is-good (and (= -7 (aref x 0))
			  (=  3 (aref x 1))
			  (=  2 (aref x 2))
			  (=  2 (aref x 3)))))
	(unless is-good
	  (format t "~&~A:" 'test0201)
	  (format t " X is ~A. " x)
	  (format t " Expected 2, 2, 3, -7."))
	is-good))))

(deftest test0260 (&optional (n 100))
  "Construct a large linear system, use GAUSS2 to solve it, check
the result."
  (declare (type (integer 0) n))
  (let* ((x (make-array n :element-type 'real :initial-element 0))
	 (a (make-array (list n (1+ n)) :element-type 'real
			:initial-element 0)))
    (declare (type (simple-array real (*)) x)
	     (type (simple-array real (* *)) a))
    (do ((i 0 (1+ i)))
	((>= i n))
	(setf (aref x i) (1+ (random 100.0)))
	(do ((j 0 (1+ j)))
	    ((>= j n))
	    (setf (aref a i j) (1+ (random 100.0)))))
    (do ((i 0 (1+ i)))
	((>= i n))
	(setf (aref a i n)
	      (loop for j from 0 to (1- n) sum (* (aref a i j) (aref x j)))))
    (let ((y (gauss2 a)))
      (declare (type (array real (*)) y))
      (let ((is-good t))
	(do ((i 0 (1+ i)))
	    ((>= i n))
	    (unless (relerrp (aref y i) (aref x i) 2/100)
	      (format t "~&X[~A] is ~A.  Y[~A] is ~A. "
		      i (aref x i) i (aref y i))
	      (format t " Relative error is ~A." (relerr (aref y i)
							 (aref x i)))
	      (setf is-good nil)))
	is-good))))

;;; --- end of file ---