(DEFINE-FILE-INFO READTABLE "XCL" PACKAGE (DEFPACKAGE "LINPACK" ( USE "LISP" "PT" "BLAS") (
NICKNAMES "LP") ( PREFIX-NAME "LP")))
(il:filecreated " 6-Dec-88 18:30:14" il:{qv}<idl>next>linpack-ops.\;15 21087
il:|changes| il:|to:| (il:functions svd-factor qr-ols qr-q-transpose-y least-squares-regress lu-factor lu-solve lu-inverse cholesky-factor qr-factor)
(il:vars il:linpack-opscoms) (il:variables stack qr-example qr-response square cholesky-example) (xcl:file-environments "LINPACK-OPS")
il:|previous| il:|date:| " 5-Dec-88 19:05:02" il:{qv}<idl>next>linpack-ops.\;14)
; Copyright (c) 1988 by Xerox Corporation. All rights reserved.
(il:prettycomprint il:linpack-opscoms)
(il:rpaqq il:linpack-opscoms ((il:types square-matrix) (il:functions square-matrix-p) (il:* il:|;;| "LU decomposition ") (il:functions lu-factor lu-solve lu-inverse) (il:functions linear-solve matrix-inverse) (il:* il:|;;| "Cholesky decomposition") (il:functions cholesky-factor) (il:* il:|;;| "QR decomposition") (il:functions qr-factor qr-q-transpose-y qr-ols) (il:functions least-squares-regress) (il:* il:|;;| "SVD factorization") (il:functions svd-factor) (il:* il:|;;| "Debugging") (il:variables stack cholesky-example qr-example qr-response) (il:functions matrix-times) (eval-when (load) (il:p (export (quote (square-matrix square-matrix-p lu-factor lu-solve lu-inverse linear-solve matrix-inverse cholesky-factor qr-factor qr-q-transpose-y qr-ols least-squares-regress)) (find-package "LINPACK")))) (xcl:file-environments "LINPACK-OPS")))
(deftype square-matrix nil (quote (satisfies square-matrix-p)))
(defun square-matrix-p (x) (and (arrayp x) (eq 2 (array-rank x)) (let ((n (array-dimension x 0))) (and (eq n (array-dimension x 1)) n))))
(il:* il:|;;| "LU decomposition ")
(defun lu-factor (matrix &optional pivot-vector factors) (il:* il:|;;;| "Computes the LU decomposition of the square N x N matrix MATRIX by Gauss elimination with row pivoting. If FACTORS is provided, it will be overwritten with the packed factors. FACTORS may be eq to MATRIX. If PIVOT-VECTOR is provided it will be overwritten with the pivot permutation. ") (il:* il:|;;;| "Returns two values, the packed factors and a pivot permutation, which together define the result.") (il:* il:|;;;| "Based on LINPACK algorithm SGESL") (let ((n (square-matrix-p matrix))) (il:* il:|;;| "Arg Checks") (if (null n) (error "Not a square matrix: ~s" matrix)) (if (null pivot-vector) (setq pivot-vector (make-array n)) (if (not (and (vectorp pivot-vector) (eq (length pivot-vector) n))) (error "Bad PIVOT-VECTOR vector: ~s" pivot-vector))) (if (null factors) (setq factors (make-array (array-dimensions matrix) :element-type (array-element-type matrix))) (if (not (eq (square-matrix-p factors) n)) (error "Bad FACTORS matrix: ~s" factors))) (il:* il:|;;| "Copy MATRIX to FACTORMATRIX") (if (not (eq matrix factors)) (copy-array matrix factors)) (il:* il:|;;| "Compute the LU decomposition of MATRIX") (let ((factors-vector (make-array (* n n) :element-type (array-element-type factors) :displaced-to factors)) pivot-index) (dotimes (i (1- n)) (il:* il:|;;| "Find pivot index") (setq pivot-index (+ i (max-abs-index factors-vector :offset (+ i (* n i)) :skip n :count (- n i)))) (setf (aref pivot-vector i) pivot-index) (when (not (zerop (aref factors pivot-index i))) (if (not (eq pivot-index i)) (il:* il:|;;| "Interchange ") (rotatef (aref factors pivot-index i) (aref factors i i))) (il:* il:|;;| "Compute Multpliers") (ax (- (/ (aref factors i i))) factors-vector :offset (+ i (* (1+ i) n)) :skip n :count (- n (1+ i))) (il:* il:|;;| "Row eliminate with column indexing") (do* ((i+1 (1+ i)) (n*i+1 (* n i+1)) (x-offset (+ i n*i+1)) (count (- n i+1)) (j i+1 (1+ j))) ((eq j n)) (if (not (eq pivot-index i)) (il:* il:|;;| "Interchange") (rotatef (aref factors pivot-index j) (aref factors i j))) (ax+y (aref factors i j) factors-vector factors-vector :x-offset x-offset :x-skip n :y-offset (+ j n*i+1) :y-skip n :count count))))) (il:* il:|;;| "No row elimination on last column") (setf (aref pivot-vector (1- n)) (1- n)) (values factors pivot-vector)))
(defun lu-solve (lu-matrix pivot-vector c-vector &optional solution) (il:* il:|;;;| "Solves to system A x = b. LU-MATRIX and PIVOT-VECTOR should be the outputs of LUFACTOR. CVECTOR should to the RHS of the system. Returns SOLUTION") (il:* il:|;;;| "lifted from LINPACK SGESL") (let ((n (square-matrix-p lu-matrix))) (if (null n) (error "Not a square matrix: ~s" lu-matrix)) (if (not (eq (length pivot-vector) n)) (error "PIVOT-VECTOR not of correct length: ~s" pivot-vector)) (if (not (eq (length c-vector) n)) (error "C-VECTOR not of correct length: ~s" c-vector)) (if (null solution) (setq solution (make-array n :element-type (array-element-type c-vector))) (if (not (eq (length solution) n)) (error "SOLUTION not of correct length: ~s" solution))) (il:* il:|;;| "Copy C-VECTOR to SOLUTION") (if (not (eq c-vector solution)) (copy-array c-vector solution)) (let ((lu-vector (make-array (* n n) :element-type (array-element-type lu-matrix) :displaced-to lu-matrix))) (il:* il:|;;| "First solve L*y = b") (dotimes (k (1- n)) (let ((pivot-index (aref pivot-vector k))) (if (not (eq pivot-index k)) (il:* il:|;;| "Interchange") (rotatef (aref solution k) (aref solution pivot-index))) (ax+y (aref solution k) lu-vector solution :x-offset (+ k (* n (1+ k))) :x-skip n :y-offset (1+ k) :y-skip 1 :count (- n (1+ k))))) (il:* il:|;;| "Then solve U*x = y") (do ((k (1- n) (1- k))) ((< k 0)) (setf (aref solution k) (/ (aref solution k) (aref lu-matrix k k))) (ax+y (- (aref solution k)) lu-vector solution :x-offset k :x-skip n :y-offset 0 :y-skip 1 :count k)) solution)))
(defun lu-inverse (lu-matrix pivot-vector &optional solution) (il:* il:|;;;| "Solves to system A x = b. LU-MATRIX and PIVOT-VECTOR should be the outputs of LUFACTOR. CVECTOR should to the RHS of the system. Returns SOLUTION") (il:* il:|;;;| "lifted from LINPACK SGESL") (let ((n (square-matrix-p lu-matrix))) (if (null n) (error "Not a square matrix: ~s" lu-matrix)) (if (not (eq (length pivot-vector) n)) (error "PIVOT-VECTOR not of correct length: ~s" pivot-vector)) (if (null solution) (setq solution (make-array (list n n) :element-type (array-element-type lu-matrix))) (if (not (eq (square-matrix-p solution) n)) (error "SOLUTION matrix not of correct size: ~s" solution))) (il:* il:|;;| "Copy LU-MATRIX to SOLUTION") (if (not (eq lu-matrix solution)) (copy-array lu-matrix solution)) (let ((solution-vector (make-array (* n n) :element-type (array-element-type solution) :displaced-to solution))) (il:* il:|;;| "first compute INVERSE (U)") (dotimes (i n) (setf (aref solution i i) (/ (aref solution i i))) (ax (- (aref solution i i)) solution-vector :offset i :skip n :count i) (do ((j (1+ i) (1+ j)) temp) ((eq j n)) (setq temp (aref solution i j)) (il:* il:|;;| "This baroqueness because we don't know the appropriate zero for this array") (setf (aref solution i j) (* 0 temp)) (ax+y temp solution-vector solution-vector :x-offset i :x-skip n :y-offset j :y-skip n :count (1+ i)))) (il:* il:|;;| "Form INVERSE (U) *INVERSE (L)") (let ((temp-vector (make-array n :element-type (array-element-type solution)))) (do ((k (- n 2) (1- k))) ((< k 0)) (do ((i (1+ k) (1+ i))) ((eq i n)) (shiftf (aref temp-vector i) (aref solution i k) (* 0 (aref solution i k)))) (do ((j (1+ k) (1+ j))) ((eq j n)) (ax+y (aref temp-vector j) solution-vector solution-vector :x-offset j :x-skip n :y-offset k :y-skip n :count n)) (if (not (eq k (aref pivot-vector k))) (swap-x-and-y solution-vector solution-vector :x-offset k :x-skip n :y-offset (aref pivot-vector k) :y-skip n :count n)))) solution)))
(defun linear-solve (matrix constants) (multiple-value-bind (lu-factors pivot-vector) (lu-factor matrix) (lu-solve lu-factors pivot-vector constants)))
(defun matrix-inverse (matrix) (multiple-value-bind (lu-factors pivot-vector) (lu-factor matrix) (lu-inverse lu-factors pivot-vector)))
(il:* il:|;;| "Cholesky decomposition")
(defun cholesky-factor (matrix &optional factor) (il:* il:|;;;| "Based on LINPACK algorithm SCHDC") (il:* il:|;;;| "Returns the upper-triangular CHOLESKY factor for the symmetric positive-definite matrix MATRIX") (let ((n (square-matrix-p matrix))) (il:* il:|;;;| "Arg Checks") (if (null n) (error "Not a square matrix: ~s" matrix)) (if (null factor) (setq factor (make-array (array-dimensions matrix) :element-type (array-element-type matrix))) (if (not (eq (square-matrix-p factor) n)) (error "Bad FACTOR matrix: ~s" factor))) (il:* il:|;;;| "Copy MATRIX to result") (if (not (eq matrix factor)) (copy-array matrix factor)) (il:* il:|;;;| "Compute the cholesky decomposition of FACTOR") (let ((factor-vector (make-array (* n n) :element-type (array-element-type factor) :displaced-to factor)) (work (make-array n :element-type (array-element-type factor)))) (dotimes (i n) (if (<= (aref factor i i) 0) (error "Zero pivot element at index: ~s" i)) (setf (aref work i) (sqrt (aref factor i i))) (setf (aref factor i i) (aref work i)) (do ((j (1+ i) (1+ j))) ((eq j n)) (setf (aref factor i j) (/ (aref factor i j) (aref work i))) (setf (aref work j) (aref factor i j)) (ax+y (- (aref factor i j)) work factor-vector :x-offset (1+ i) :x-skip 1 :y-offset (+ j (* n (1+ i))) :y-skip n :count (- j i))))) (il:* il:|;;| "Zero out lower diagonal components") (do ((i 1 (1+ i))) ((eq i n)) (do ((j 0 (1+ j))) ((eq j i)) (setf (aref factor i j) 0.0))) factor))
(il:* il:|;;| "QR decomposition")
(defun qr-factor (matrix &optional qr-aux factors) (il:* il:|;;;| "Computes the QR decomposition of the N x P matrix MATRIX by Householder reduction. If FACTORS is provided, it will be overwritten with the packed factors. FACTORS may be eq to MATRIX. If QR-AUX is provided it will be overwritten with the auxillary information required to recover the orthogonal factors.") (il:* il:|;;;| "Returns two values, the packed factors and a QR-AUX, which together define the result.") (il:* il:|;;;| "Lifted from LINPACK algorithmSQRDC") (let ((n (array-dimension matrix 0)) (p (array-dimension matrix 1))) (il:* il:|;;| "Arg Checks") (if (null qr-aux) (setq qr-aux (make-array p :element-type (array-element-type matrix))) (unless (and (eq 1 (array-rank qr-aux)) (eq p (array-total-size qr-aux))) (error "QR-AUX not of size: ~s" p))) (if (null factors) (setq factors (make-array (array-dimensions matrix) :element-type (array-element-type matrix))) (unless (equal (array-dimensions matrix) (array-dimensions factors)) (error "Bad FACTORS matrix: ~s" factors))) (il:* il:|;;| "Copy MATRIX to FACTORMATRIX") (unless (eq matrix factors) (copy-array matrix factors)) (il:* il:|;;| "Compute the QR decomposition of FACTORMATRIX") (fill-array qr-aux 0.0) (let ((factors-vector (make-array (* n p) :element-type (array-element-type factors) :displaced-to factors))) (dotimes (l (min (1- n) p)) (il:* il:|;;| "Compute the Householder transformation for column L") (let ((nrmxl (l2-norm factors-vector :offset (+ l (* p l)) :skip p :count (- n l)))) (when (> nrmxl 0.0) (if (< (aref factors l l) 0.0) (setq nrmxl (- nrmxl))) (ax (/ 1.0 nrmxl) factors-vector :offset (+ l (* p l)) :skip p :count (- n l)) (setf (aref factors l l) (+ 1.0 (aref factors l l))) (il:* il:|;;| "apply the transform to the remaining columns") (do ((j (1+ l) (1+ j))) ((eq j p)) (let ((temp (- (/ (<xy> factors-vector factors-vector :x-offset (+ l (* p l)) :x-skip p :y-offset (+ j (* p l)) :y-skip p :count (- n l)) (aref factors l l))))) (ax+y temp factors-vector factors-vector :x-offset (+ l (* p l)) :x-skip p :y-offset (+ j (* p l)) :y-skip p :count (- n l)))) (setf (aref qr-aux l) (aref factors l l)) (setf (aref factors l l) (- nrmxl)))))) (values factors qr-aux)))
(defun qr-q-transpose-y (qr-matrix qr-aux y &optional product) (il:* il:|;;;| "COMPUTE (TRANS Q) * Y given a QR factorization described by QR-MATRIX and QR-AUX where Y is an N vector") (il:* il:|;;;| "Product may be eq to Y.") (il:* il:|;;;| "Lifted from LINPACK algorithm SQRSL") (let ((n (array-dimension qr-matrix 0)) (p (array-dimension qr-matrix 1))) (il:* il:|;;| "Arg Checks") (if (not (and (eq 1 (array-rank y)) (eq n (array-total-size y)))) (error "Size mismatch: ~s" y)) (if (null product) (setq product (make-array n :element-type (array-element-type y))) (unless (and (eq 1 (array-rank product)) (eq n (array-total-size product))) (error "Size mismatch: ~s" product))) (il:* il:|;;| "Copy Y to product") (if (not (eq y product)) (copy-array y product)) (let ((qr-vector (make-array (* n p) :element-type (array-element-type qr-matrix) :displaced-to qr-matrix))) (dotimes (j (min p n)) (unless (zerop (aref qr-aux j)) (let ((temp (aref qr-matrix j j))) (setf (aref qr-matrix j j) (aref qr-aux j)) (ax+y (- (/ (<xy> qr-vector product :x-offset (+ j (* p j)) :x-skip p :y-offset j :y-skip 1 :count (- n j)) (aref qr-matrix j j))) qr-vector product :x-offset (+ j (* p j)) :x-skip p :y-offset j :y-skip 1 :count (- n j)) (setf (aref qr-matrix j j) temp))))) product))
(defun qr-ols (qr-matrix qr-aux y &key q-transpose-y coefficients residuals y-hat) (il:* il:|;;| "Lifted from LINPACK algorithm SQRSL") (let ((n (array-dimension qr-matrix 0)) (p (array-dimension qr-matrix 1))) (macrolet ((check-or-make-vector (v size) (il:bquote (if (null (il:\\\, v)) (setq (il:\\\, v) (make-array (il:\\\, size) :element-type (array-element-type y))) (check-vector (il:\\\, v) (il:\\\, size))))) (check-vector (v size) (il:bquote (unless (and (eq 1 (array-rank (il:\\\, v))) (eq (il:\\\, size) (array-total-size (il:\\\, v)))) (error "Size mismatch: ~s" (il:\\\, v)))))) (il:* il:|;;| "Arg Checks") (check-vector y n) (check-or-make-vector q-transpose-y n) (check-or-make-vector coefficients p) (let ((qr-vector (make-array (* n p) :element-type (array-element-type qr-matrix) :displaced-to qr-matrix))) (il:* il:|;;| "Compute TRANS (Q) * Y") (qr-q-transpose-y qr-matrix qr-aux y q-transpose-y) (il:* il:|;;| "Compute Coefficients") (dotimes (i p) (il:* il:|;;| "Set up computation of Coefficients") (setf (aref coefficients i) (aref q-transpose-y i))) (do ((j (1- p) (1- j))) ((< j 0)) (when (zerop (aref qr-matrix j j)) (error "Singular matrix at index: ~s" j)) (setf (aref coefficients j) (/ (aref coefficients j) (aref qr-matrix j j))) (unless (eq j 0) (ax+y (- (aref coefficients j)) qr-vector coefficients :x-offset j :x-skip p :y-offset 0 :y-skip 1 :count j))) (macrolet ((compute-y-component (component) (il:bquote (do ((j (1- (min p (1- n))) (1- j))) ((< j 0)) (unless (zerop (aref qr-matrix j j)) (let ((temp (aref qr-matrix j j))) (setf (aref qr-matrix j j) (aref qr-aux j)) (ax+y (- (/ (<xy> qr-vector (il:\\\, component) :x-offset (+ j (* p j)) :x-skip p :y-offset j :y-skip 1 :count (- n j)) (aref qr-matrix j j))) qr-vector (il:\\\, component) :x-offset (+ j (* p j)) :x-skip p :y-offset j :y-skip 1 :count (- n j)) (setf (aref qr-matrix j j) temp))))))) (il:* il:|;;| "Compute RSD") (when residuals (check-vector residuals n) (il:* il:|;;| "Set up computation of RSD") (dotimes (i n) (setf (aref residuals i) (if (< i p) 0.0 (aref q-transpose-y i)))) (compute-y-component residuals)) (when y-hat (check-vector y-hat n) (il:* il:|;;| "Set up computation of Y-HAT") (dotimes (i n) (setf (aref y-hat i) (if (< i p) (aref q-transpose-y i) 0.0))) (il:* il:|;;| "Compute YHAT") (compute-y-component y-hat)))) (values coefficients residuals y-hat))))
(defun least-squares-regress (y x &key coefficients residuals y-hat) (il:* il:|;;;| "LEAST-SQUARES-REGRESS calculates the least squares (multiple) regression of Y on X. ") (multiple-value-bind (qrmatrix qraux) (qr-factor x) (qr-ols qrmatrix qraux y :coefficients coefficients :residuals residuals :y-hat y-hat)))
(il:* il:|;;| "SVD factorization")
(defun svd-factor (matrix &optional s-vector u-matrix v-matrix) (il:* il:|;;;| "Singular-value decomposition by means of orthogonalization by plane rotations. ") (il:* il:|;;;| "Taken from Nash and Shlien: \"Partial svd algorithms.\"") (il:* il:|;;;| " On entry MATRIX contains the M by N matrix to be decomposed.") (il:* il:|;;;| " SVECTOR must be a vector of length N and VMATRIX must be a square N by N matrix. ") (il:* il:|;;;| "On return UMATRIX has been overwritten by the left singular vectors, SVECTOR contains the singular values, and VMATRIX contains the right singular vectors.") (let ((m (array-dimension matrix 0)) (n (array-dimension matrix 1))) (il:* il:|;;| "Args checks") (macrolet ((check-array (array rank dimensions) (il:bquote (let ((dims (il:\\\, dimensions))) (if (null (il:\\\, array)) (setq (il:\\\, array) (make-array dims :element-type (array-element-type matrix))) (unless (and (eq (il:\\\, rank) (array-rank (il:\\\, array))) (equal dims (array-dimensions (il:\\\, array)))) (error "Size mismatch: ~s" (il:\\\, array)))))))) (check-array s-vector 1 (list n)) (check-array u-matrix 2 (list m n)) (check-array v-matrix 2 (list n n))) (il:* il:|;;| "Copy MATRIX to UMATRIX") (copy-array matrix u-matrix) (il:* il:|;;| "Initialize VMATRIX to identity matrix.") (fill-array v-matrix 0.0) (dotimes (i n) (setf (aref v-matrix i i) 1.0)) (il:* il:|;;| "Start the computation") (let ((u-vector (make-array (* m n) :element-type (array-element-type u-matrix) :displaced-to u-matrix)) (v-vector (make-array (* n n) :element-type (array-element-type v-matrix) :displaced-to v-matrix)) (count n) (eps 1.0E-6) (sweep-limit (max (truncate n 4) 6)) (sweep-count 0)) (il:* il:|;;| "The main loop: repeatedly sweep over all pairs of columns in U, rotating as needed, until no rotations in a complete sweep are effective. Check the opportunity for rank reduction at the conclusion of each sweep.") (loop (let ((rcount (truncate (* count (1- count)) 2))) (if (> (incf sweep-count) sweep-limit) (error "Number of sweeps exceeds sweep limit: ~s" sweep-count)) (dotimes (j (1- count)) (do ((k (1+ j) (1+ k))) ((eq k count)) (let ((p (<xy> u-vector u-vector :x-offset j :x-skip n :y-offset k :y-skip n :count m)) (q (<xy> u-vector u-vector :x-offset j :x-skip n :y-offset j :y-skip n :count m)) (r (<xy> u-vector u-vector :x-offset k :x-skip n :y-offset k :y-skip n :count m)) c s v) (setf (aref s-vector j) q) (setf (aref s-vector k) r) (cond ((< q r) (setq q (- (/ q r) 1.0)) (setq p (/ p r)) (setq v (sqrt (+ (* 4.0 p p) (* q q)))) (setq s (sqrt (* 0.5 (- 1.0 (/ q v))))) (if (< p 0.0) (setq s (- s))) (setq c (/ p (* v s))) (givens-rotate c s u-vector u-vector :x-offset j :x-skip n :y-offset k :y-skip n :count m) (givens-rotate c s v-vector v-vector :x-offset j :x-skip n :y-offset k :y-skip n :count n)) ((or (<= (* q r) (* eps eps)) (<= (* (/ p q) (/ p r)) eps)) (decf rcount)) (t (setq r (- 1.0 (/ r q))) (setq p (/ p q)) (setq v (sqrt (+ (* 4.0 p p) (* r r)))) (setq c (sqrt (* 0.5 (+ 1.0 (/ r v))))) (setq s (/ p (* v c))) (givens-rotate c s u-vector u-vector :x-offset j :x-skip n :y-offset k :y-skip n :count m) (givens-rotate c s v-vector v-vector :x-offset j :x-skip n :y-offset k :y-skip n :count n)))))) (loop (unless (and (>= count 3) (<= (/ (aref s-vector (1- count)) (+ (aref s-vector 0) eps)) eps)) (return nil)) (decf count)) (unless (> rcount 0) (return nil)))) (il:* il:|;;| "Finish the decomposition by returning all N singular values, and by normalizing those columns of UMATRIX judged to be non-zero.") (dotimes (j n) (let ((q (sqrt (aref s-vector j)))) (setf (aref s-vector j) q) (if (<= j count) (ax (/ 1.0 q) u-vector :offset j :skip n :count m)))) (values s-vector u-matrix v-matrix))))
(il:* il:|;;| "Debugging")
(defparameter stack (make-array (quote (21 3)) :element-type (quote single-float) :initial-contents (quote ((80.0 27.0 89.0) (80.0 27.0 88.0) (75.0 25.0 90.0) (62.0 24.0 87.0) (62.0 22.0 87.0) (62.0 23.0 87.0) (62.0 24.0 93.0) (62.0 24.0 93.0) (58.0 23.0 87.0) (58.0 18.0 80.0) (58.0 18.0 89.0) (58.0 17.0 88.0) (58.0 18.0 82.0) (58.0 19.0 93.0) (50.0 18.0 89.0) (50.0 18.0 86.0) (50.0 19.0 72.0) (50.0 19.0 79.0) (50.0 20.0 80.0) (56.0 20.0 82.0) (70.0 20.0 91.0)))))
(defparameter cholesky-example (make-array (quote (2 2)) :initial-contents (quote ((2 -2) (-2 5)))))
(defparameter qr-example (make-array (quote (5 2)) :element-type (quote single-float) :initial-contents (quote ((12.0 5.0) (15.0 2.0) (10.0 7.0) (23.0 5.0) (16.0 6.0)))))
(defparameter qr-response (make-array 5 :element-type (quote single-float) :initial-contents (with-collection (dotimes (i 5) (collect (+ (* 1.0 (aref qr-example i 0)) (* 10.0 (aref qr-example i 1))))))))
(defun matrix-times (a b) (ea:inner-product (quote +) (quote *) a b))
(eval-when (load)
(export (quote (square-matrix square-matrix-p lu-factor lu-solve lu-inverse linear-solve matrix-inverse cholesky-factor qr-factor qr-q-transpose-y qr-ols least-squares-regress)) (find-package "LINPACK"))
)
(xcl:define-file-environment "LINPACK-OPS" :package (xcl:defpackage "LINPACK" (:use "LISP" "PT" "BLAS") (:nicknames "LP") (:prefix-name "LP")) :readtable "XCL" :compiler :compile-file)
(il:putprops il:linpack-ops il:copyright ("Xerox Corporation" 1988))
(il:declare\: il:dontcopy
(il:filemap (nil)))
il:stop