--- /dev/null
+
+Copyright (c) 2008 NANRI <southly@gmail.com>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(in-package :cl-jama)
+
+(defclass Cholesky ()
+ ((L_)
+ (isspd)))
+
+(defmethod is-spd ((A Cholesky))
+ (with-slots (isspd) A
+ isspd))
+
+(defmethod get-L ((A Cholesky))
+ (with-slots (L_) A
+ (copy-array L_)))
+
+(defmethod initialize-instance :after ((obj Cholesky) &key matrix)
+ (with-slots (L_ isspd) obj
+ (let ((m (array-dimension matrix 0))
+ (n (array-dimension matrix 1)))
+ (setf isspd (eql m n))
+ (unless (eql m n)
+ (setf L_ nil)
+ (return-from initialize-instance nil))
+ (setf L_ (make-array (list n n)))
+ ;;
+ (loop :for j :from 0 :below n
+ :do (let ((d 0.0))
+ (loop :for k :from 0 :below j
+ :do (let ((s 0.0))
+ (loop :for i :from 0 :below k
+ :do (incf s (* (aref L_ k i) (aref L_ j i))))
+ (setf s (/ (- (aref matrix j k) s) (aref L_ k k)))
+ (setf (aref L_ j k) s)
+ (setf d (+ d (* s s)))
+ (setf isspd (and isspd (eql (aref matrix k j) (aref matrix j k))))))
+ (setf d (- (aref matrix j j) d))
+ (setf isspd (and isspd (> d 0.0)))
+ (setf (aref L_ j j) (sqrt (if (> d 0.0) d 0.0)))
+ (loop :for k :from (1+ j) :below n
+ :do (setf (aref L_ j k) 0.0)))))))
+
+(defmethod solve ((A Cholesky) (b vector))
+ (with-slots (L_) A
+ (unless (eql (array-dimension L_ 0)
+ (array-dimension b 0))
+ (return-from solve nil))
+ (let ((n (array-dimension L_ 0))
+ (x (copy-array b)))
+ (loop :for k :from 0 :below n
+ :do (loop :for i :from 0 :below k
+ :do (decf (aref x k) (* (aref x i) (aref L_ k i))))
+ (setf (aref x k) (/ (aref x k) (aref L_ k k))))
+ (loop :for k :from (1- n) :downto 0
+ :do (loop :for i :from (1+ k) :below n
+ :do (decf (aref x k) (* (aref x i) (aref L_ i k))))
+ (setf (aref x k) (/ (aref x k) (aref L_ k k))))
+ x)))
+
+(defmethod solve ((A Cholesky) (B array))
+ (with-slots (L_) A
+ (unless (eql (array-dimension L_ 0)
+ (array-dimension B 0))
+ (return-from solve nil))
+ (let ((n (array-dimension L_ 0))
+ (X (copy-array B))
+ (nx (array-dimension B 1)))
+ (loop :for j :from 0 :below nx
+ :do (loop :for k :from 0 :below n
+ :do (loop :for i :from 0 :below k
+ :do (decf (aref X k j) (* (aref X i j) (aref L_ k i))))
+ (setf (aref X k j) (/ (aref X k j) (aref L_ k k)))))
+ (loop :for j :from 0 :below nx
+ :do (loop :for k :from (1- n) :downto 0
+ :do (loop :for i :from (1+ k) :below n
+ :do (decf (aref X k j) (* (aref X i j) (aref L_ i k))))
+ (setf (aref X k j) (/ (aref X k j) (aref L_ k k)))))
+ X)))
--- /dev/null
+(asdf:defsystem :cl-jama
+ :description "The CL port of JAMA/C++ Linear Algebra Package"
+ :version "1.2.5-1"
+ :author "NANRI <southly@gmail.com>"
+ :licence "The MIT License"
+ :depends-on ()
+ :components ((:file "package")
+ (:file "utils" :depends-on ("package"))
+ (:file "lu" :depends-on ("utils"))
+ (:file "cholesky" :depends-on ("utils"))
+ (:file "qr" :depends-on ("utils"))
+ (:file "svd" :depends-on ("utils"))
+ (:file "eigenvalue" :depends-on ("utils"))
+ ))
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(in-package :cl-jama)
+
+(declaim (optimize (debug 3) (safety 3)))
+
+(defclass Eigenvalue ()
+ ((n)
+ (issymmetric)
+ (d)
+ (e)
+ (V)
+ (H)
+ (ort)
+ (cdivi)
+ (cdivr)))
+
+(defmethod tred2 ((A Eigenvalue))
+ (with-slots (n d e V) A
+ (loop :for j :from 0 :below n
+ :do (setf (aref d j) (aref V (1- n) j)))
+ ;;
+ (loop :for i :from (1- n) :above 0
+ :for scale = 0.0 :then 0.0
+ :for h = 0.0 :then 0.0
+ :do (loop :for k :from 0 :below i
+ :do (setf scale (+ scale (abs (aref d k)))))
+ :do (cond ((zerop scale)
+ (setf (aref e i) (aref d (1- i)))
+ (loop :for j :from 0 :below i
+ :do (setf (aref d j) (aref V (1- i) j)
+ (aref V i j) 0.0
+ (aref V j i) 0.0)))
+ (t
+ (loop :for k :from 0 :below i
+ :do (setf (aref d k) (/ (aref d k) scale))
+ (incf h (* (aref d k) (aref d k))))
+ (let ((f (aref d (1- i)))
+ (g (sqrt h))
+ hh)
+ (when (plusp f)
+ (setf g (- g)))
+ (setf (aref e i) (* scale g)
+ h (- h (* f g))
+ (aref d (1- i)) (- f g))
+ (loop :for j :from 0 :below i
+ :do (setf (aref e j) 0.0))
+ ;;
+ (loop :for j :from 0 :below i
+ :do (setf f (aref d j)
+ (aref V j i) f
+ g (+ (aref e j) (* (aref V j j) f)))
+ :do (loop :for k :from (1+ j) :to (1- i)
+ :do (incf g (* (aref V k j) (aref d k)))
+ (incf (aref e k) (* (aref V k j) f)))
+ :do (setf (aref e j) g))
+ (setf f 0.0)
+ (loop :for j :from 0 :below i
+ :do (setf (aref e j) (/ (aref e j) h))
+ (incf f (* (aref e j) (aref d j))))
+ (setf hh (/ f (+ h h)))
+ (loop :for j :from 0 :below i
+ :do (decf (aref e j) (* hh (aref d j))))
+ (loop :for j :from 0 :below i
+ :do (setf f (aref d j)
+ g (aref e j))
+ :do (loop :for k :from j :to (1- i)
+ :do (decf (aref V k j) (+ (* f (aref e k)) (* g (aref d k)))))
+ :do (setf (aref d j) (aref V (1- i) j)
+ (aref V i j) 0.0)))))
+ :do (setf (aref d i) h))
+ ;;
+ (loop :for i :from 0 :below (1- n)
+ :do (setf (aref V (1- n) i) (aref V i i)
+ (aref V i i) 1.0)
+ (let ((h (aref d (1+ i))))
+ (when (not (zerop h))
+ (loop :for k :from 0 :to i
+ :do (setf (aref d k) (/ (aref V k (1+ i)) h)))
+ (loop :for j :from 0 :to i
+ :do (let ((g 0.0))
+ (loop :for k :from 0 :to i
+ :do (incf g (* (aref V k (1+ i)) (aref V k j))))
+ (loop :for k :from 0 :to i
+ :do (decf (aref V k j) (* g (aref d k))))))))
+ (loop :for k :from 0 :to i
+ :do (setf (aref V k (1+ i)) 0.0)))
+ (loop :for j :from 0 :below n
+ :do (setf (aref d j) (aref V (1- n) j)
+ (aref V (1- n) j) 0.0))
+ (setf (aref V (1- n) (1- n)) 1.0
+ (aref e 0) 0.0)))
+
+(defmethod tql2 ((A Eigenvalue))
+ (with-slots (n d e V) A
+ (loop :for i :from 1 :below n
+ :do (setf (aref e (1- i)) (aref e i)))
+ (setf (aref e (1- n)) 0.0)
+ (let ((f 0.0)
+ (tst1 0.0)
+ (eps (expt 2.0 -52.0)))
+ (loop :for l :from 0 :below n
+ :do (setf tst1 (max tst1 (+ (abs (aref d l)) (abs (aref e l)))))
+ :do (let ((m l))
+ (loop :while (< m n)
+ :thereis (<= (abs (aref e m)) (* eps tst1))
+ :do (incf m))
+ (when (> m l)
+ (let ((iter 0))
+ (loop :do (setf iter (1+ iter))
+ :do (let* ((g (aref d l))
+ (p (/ (- (aref d (1+ l)) g) (* 2.0 (aref e l))))
+ (r (hypot p 1.0))
+ dl1 h c c2 c3 el1 s s2)
+ (when (minusp p)
+ (setf r (- r)))
+ (setf (aref d l) (/ (aref e l) (+ p r))
+ (aref d (1+ l)) (* (aref e l) (+ p r))
+ dl1 (aref d (1+ l))
+ h (- g (aref d l)))
+ (loop :for i :from (+ l 2) :below n
+ :do (decf (aref d i) h))
+ (setf f (+ f h))
+ ;;
+ (setf p (aref d m)
+ c 1.0
+ c2 c
+ c3 c
+ el1 (aref e (1+ l))
+ s 0.0
+ s2 0.0)
+ (loop :for i :from (1- m) :downto l
+ :do (setf c3 c2
+ c2 c
+ s2 s
+ g (* c (aref e i))
+ h (* c p)
+ r (hypot p (aref e i))
+ (aref e (1+ i)) (* s r)
+ s (/ (aref e i) r)
+ c (/ p r)
+ p (- (* c (aref d i)) (* s g))
+ (aref d (1+ i)) (+ h (* s (+ (* c g) (* s (aref d i))))))
+ :do (loop :for k :from 0 :below n
+ :do (setf h (aref V k (1+ i))
+ (aref V k (1+ i)) (+ (* s (aref V k i)) (* c h))
+ (aref V k i) (- (* c (aref V k i)) (* s h)))))
+ (setf p (/ (* (- s) s2 c3 el1 (aref e l)) dl1)
+ (aref e l) (* s p)
+ (aref d l) (* c p)))
+ :while (> (abs (aref e l)) (* eps tst1)))))
+ (setf (aref d l) (+ (aref d l) f)
+ (aref e l) 0.0)))
+ ;;
+ (loop :for i :from 0 :below (1- n)
+ :for k = i :then i
+ :for p = (aref d i) :then (aref d i)
+ :do (loop :for j :from (1+ i) :below n
+ :do (when (< (aref d j) p)
+ (setf k j
+ p (aref d j))))
+ (unless (eql k i)
+ (setf (aref d k) (aref d i)
+ (aref d i) p)
+ (loop :for j :from 0 :below n
+ :do (setf p (aref V j i)
+ (aref V j i) (aref V j k)
+ (aref V j k) p)))))))
+
+(defmethod orthes ((A Eigenvalue))
+ (with-slots (n ort H V) A
+ (let ((low 0)
+ (high (1- n)))
+ (loop :for m :from (1+ low) :to (1- high)
+ :do (let ((scale 0.0))
+ (loop :for i :from m :to high
+ :do (setf scale (+ scale (abs (aref H i (1- m))))))
+ (unless (zerop scale)
+ (let ((rh 0.0) g)
+ (loop :for i :from high :downto m
+ :do (setf (aref ort i) (/ (aref H i (1- m)) scale))
+ (incf rh (* (aref ort i) (aref ort i))))
+ (setf g (sqrt rh))
+ (when (plusp (aref ort m))
+ (setf g (- g)))
+ (setf rh (- rh (* (aref ort m) g))
+ (aref ort m) (- (aref ort m) g))
+ (loop :for j :from m :below n
+ :do (let ((f 0.0))
+ (loop :for i :from high :downto m
+ :do (incf f (* (aref ort i) (aref H i j))))
+ (setf f (/ f rh))
+ (loop :for i :from m :to high
+ :do (decf (aref H i j) (* f (aref ort i))))))
+ (loop :for i :from 0 :to high
+ :do (let ((f 0.0))
+ (loop :for j :from high :downto m
+ :do (incf f (* (aref ort j) (aref H i j))))
+ (setf f (/ f rh))
+ (loop :for j :from m :to high
+ :do (decf (aref H i j) (* f (aref ort j))))))
+ (setf (aref ort m) (* scale (aref ort m))
+ (aref H m (1- m)) (* scale g))))))
+ ;;
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below n
+ :do (setf (aref V i j) (if (eql i j) 1.0 0.0))))
+
+ (loop :for m :from (1- high) :downto (1+ low)
+ :do (unless (zerop (aref H m (1- m)))
+ (loop :for i :from (1+ m) :to high
+ :do (setf (aref ort i) (aref H i (1- m))))
+ (loop :for j :from m :to high
+ :do (let ((g 0.0))
+ (loop :for i :from m :to high
+ :do (incf g (* (aref ort i) (aref V i j))))
+ ;;
+ (setf g (/ (/ g (aref ort m)) (aref H m (1- m))))
+ (loop :for i :from m :to high
+ :do (incf (aref V i j) (* g (aref ort i)))))))))))
+
+
+(defmethod cdiv ((A Eigenvalue) xr xi yr yi)
+ (with-slots (cdivr cdivi) A
+ (let (r d)
+ (if (> (abs yr) (abs yi))
+ (setf r (/ yi yr)
+ d (+ yr (* r yi))
+ cdivr (/ (+ xr (* r xi)) d)
+ cdivi (/ (- xi (* r xr)) d))
+ (setf r (/ yr yi)
+ d (+ yi (* r yr))
+ cdivr (/ (+ (* r xr) xi) d)
+ cdivi (/ (- (* r xi) xr) d))))))
+
+(defmethod hqr2 ((A Eigenvalue))
+ (with-slots (d e H V cdivr cdivi) A
+ (let* ((nn (slot-value A 'n))
+ (n (1- nn))
+ (low 0)
+ (high (1- nn))
+ (eps (expt 2.0 -52.0))
+ (exshift 0.0)
+ (p 0.0) (q 0.0) (r 0.0) (s 0.0) (z 0.0)
+ rt w x y
+ ;;
+ (norm 0.0)
+ (iter 0))
+ (loop :for i :from 0 :below nn
+ :do (when (or (< i low) (> i high))
+ (setf (aref d i) (aref H i i)
+ (aref e i) 0.0))
+ :do (loop :for j :from (max (1- i) 0) :below nn
+ :do (setf norm (+ norm (abs (aref H i j))))))
+ ;;
+ (loop :while (> n low)
+ :do (break "n = ~S~%~S~%~S" n d e)
+ :do (let ((l n))
+ (loop :while (> l low)
+ :do (setf s (+ (abs (aref H (1- l) (1- l))) (abs (aref H l l))))
+ (when (zerop s)
+ (setf s norm))
+ :thereis (< (abs (aref H l (1- l))) (* eps s))
+ :do (decf l))
+ ;;
+ (cond ((eql l n)
+ (setf (aref H n n) (+ (aref H n n) exshift)
+ (aref d n) (aref H n n)
+ (aref e n) 0.0)
+ (decf n)
+ (setf iter 0))
+ ((eql l (1- n))
+ (setf w (* (aref H n (1- n)) (aref H (1- n) n))
+ p (/ (- (aref H (1- n) (1- n)) (aref H n n)) 2.0)
+ q (+ (* p p) w)
+ z (sqrt (abs q))
+ (aref H n n) (+ (aref H n n) exshift)
+ (aref H (1- n) (1- n)) (+ (aref H (1- n) (1- n)) exshift)
+ x (aref H n n))
+ (cond ((>= q 0)
+ (if (>= p 0)
+ (setf z (+ p z))
+ (setf z (- p z)))
+ (setf (aref d (1- n)) (+ x z)
+ (aref d n) (aref d (1- n)))
+ (unless (zerop z)
+ (setf (aref d n) (- x (/ w z))))
+ (setf (aref e (1- n)) 0.0
+ (aref e n) 0.0
+ x (aref H n (1- n))
+ s (+ (abs x) (abs z))
+ p (/ x s)
+ q (/ z s)
+ r (sqrt (+ (* p p) (* q q)))
+ p (/ p r)
+ q (/ q r))
+ ;;
+ (loop :for j :from (1- n) :below nn
+ :do (setf z (aref H (1- n) j)
+ (aref H (1- n) j) (+ (* q z) (* p (aref H n j)))
+ (aref H n j) (- (* q (aref H n j)) (* p z))))
+ ;;
+ (loop :for i :from 0 :to n
+ :do (setf z (aref H i (1- n))
+ (aref H i (1- n)) (+ (* q z) (* p (aref H i n)))
+ (aref H i n) (- (* q (aref H i n)) (* p z))))
+ ;;
+ (loop :for i :from low :to high
+ :do (setf z (aref V i (1- n))
+ (aref V i (1- n)) (+ (* q z) (* p (aref V i n)))
+ (aref V i n) (- (* q (aref V i n)) (* p z)))))
+ (t
+ (setf (aref d (1- n)) (+ x p)
+ (aref d n) (+ x p)
+ (aref e (1- n)) z
+ (aref e n) (- z))))
+ (setf n (- n 2)
+ iter 0))
+ (t
+ (setf x (aref H n n)
+ y 0.0
+ w 0.0)
+ (when (< l n)
+ (setf y (aref H (1- n) (1- n))
+ w (* (aref H n (1- n)) (aref H (1- n) n))))
+ ;;
+ (when (eql iter 10)
+ (incf exshift x)
+ (loop :for i :from low :to n
+ :do (decf (aref H i i) x))
+ (setf s (+ (abs (aref H n (1- n))) (abs (aref H (1- n) (- n 2))))
+ x (* 0.75 s)
+ y (* 0.75 s)
+ w (* -0.4375 s s)))
+ ;;
+ (when (eql iter 30)
+ (setf s (/ (- y x) 2.0)
+ s (+ (* s s) w))
+ (when (plusp s)
+ (setf s (sqrt s))
+ (when (< y x)
+ (setf s (- s)))
+ (setf s (- x (/ w (+ (/ (- y x) 2.0) s))))
+ (loop :for i :from low :to n
+ :do (decf (aref H i i) s))
+ (incf exshift s)
+ (setf x 0.964
+ y 0.964
+ w 0.964)))
+ (setf iter (+ iter 1))
+ ;;
+ (let ((m (- n 2)))
+ (loop :while (>= m l)
+ :do (setf z (aref H m m)
+ r (- x z)
+ s (- y z)
+ p (+ (/ (- (* r s) w) (aref H (1+ m) m)) (aref H m (1+ m)))
+ q (- (aref H (1+ m) (1+ m)) z r s)
+ r (aref H (+ m 2) (1+ m))
+ s (+ (abs p) (abs q) (abs r))
+ p (/ p s)
+ q (/ q s)
+ r (/ r s))
+ :thereis (eql m l)
+ :thereis (< (* (abs (aref H m (1- m))) (+ (abs q) (abs r)))
+ (* eps (* (abs p) (+ (abs (aref H (1- m) (1- m))) (abs z) (abs (aref H (1+ m) (1+ m)))))))
+ :do (decf m))
+ (loop :for i :from (+ m 2) :to n
+ :do (setf (aref H i (- i 2)) 0.0)
+ (when (> i (+ m 2))
+ (setf (aref H i (- i 3)) 0.0)))
+ ;;
+ (loop :for k :from m :to (1- n)
+ :do (let ((notlast (not (eql k (1- n)))))
+ (unless (eql k m)
+ (setf p (aref H k (1- k))
+ q (aref H (1+ k) (1- k))
+ r (if notlast (aref H (+ k 2) (1- k)) 0.0)
+ x (+ (abs p) (abs q) (abs r)))
+ (unless (zerop x)
+ (setf p (/ p x)
+ q (/ q x)
+ r (/ r x))))
+ (when (zerop x)
+ (return))
+ (setf s (sqrt (+ (* p p) (* q q) (* r r))))
+ (when (minusp p)
+ (setf s (- s)))
+ (unless (zerop s)
+ (cond ((not (eql k m))
+ (setf (aref H k (1- k)) (* (- s) x)))
+ ((not (eql l m))
+ (setf (aref H k (1- k)) (- (aref H k (1- k))))))
+ (setf p (+ p s)
+ x (/ p s)
+ y (/ q s)
+ z (/ r s)
+ q (/ q p)
+ r (/ r p))
+ ;;
+ (loop :for j :from k :below nn
+ :do (setf p (+ (aref H k j) (* q (aref H (1+ k) j))))
+ :when notlast
+ :do (setf p (+ p (* r (aref H (+ k 2) j)))
+ (aref H (+ k 2) j) (- (aref H (+ k 2) j) (* p z)))
+ :end
+ :do (setf (aref H k j) (- (aref H k j) (* p x))
+ (aref H (1+ k) j) (- (aref H (1+ k) j) (* p y))))
+ ;;
+ (loop :for i :from 0 :to (min n (+ k 3))
+ :do (setf p (+ (* x (aref H i k)) (* y (aref H i (1+ k)))))
+ :when notlast
+ :do (setf p (+ p (* z (aref H i (+ k 2))))
+ (aref H i (+ k 2)) (- (aref H i (+ k 2)) (* p r)))
+ :end
+ :do (setf (aref H i k) (- (aref H i k) p)
+ (aref H i (1+ k)) (- (aref H i (1+ k)) (* p q))))
+ ;;
+ (loop :for i :from low :to high
+ :do (setf p (+ (* x (aref V i k)) (* y (aref V i (1+ k)))))
+ :when notlast
+ :do (setf p (+ p (* z (aref V i (+ k 2))))
+ (aref V i (+ k 2)) (- (aref V i (+ k 2)) (* p r)))
+ :end
+ :do (setf (aref V i k) (- (aref V i k) p)
+ (aref V i (1+ k)) (- (aref V i (1+ k)) (* p q))))))))))))
+ (when (zerop norm)
+ (return-from hqr2 nil))
+ (loop :for n :from (1- nn) :downto 0
+ :do (setf p (aref d n)
+ q (aref e n))
+ (cond ((zerop q)
+ (let ((l n))
+ (setf (aref H n n) 1.0)
+ (loop :for i :from (1- n) :downto 0
+ :do (setf w (- (aref H i i) p)
+ r 0.0)
+ (loop :for j :from l :to n
+ :do (setf r (+ r (* (aref H i j) (aref H j n)))))
+ (cond ((minusp (aref e i))
+ (setf z w
+ s r))
+ (t
+ (setf l i)
+ (cond ((zerop (aref e i))
+ (if (not (zerop w))
+ (setf (aref H i n) (/ (- r) w))
+ (setf (aref H i n) (/ (- r) (* eps norm)))))
+ (t
+ (setf x (aref H i (1+ i))
+ y (aref H (1+ i) i)
+ q (+ (* (- (aref d i) p) (- (aref d i) p)) (* (aref e i) (aref e i)))
+ rt (/ (- (* x s) (* z r)) q)
+ (aref H i n) rt)
+ (if (> (abs x) (abs z))
+ (setf (aref H (1+ i) n) (/ (- (- r) (* w rt)) x))
+ (setf (aref H (1+ i) n) (/ (- (- s) (* y rt)) z)))
+ (setf rt (abs (aref H i n)))
+ (when (> (* (* eps rt) rt) 1)
+ (loop :for j :from i :to n
+ :do (setf (aref H j n) (/ (aref H j n) rt)))))))))))
+ ((< q 0)
+ (let ((l (1- n)))
+ (cond ((> (abs (aref H n (1- n))) (abs (aref H (1- n) n)))
+ (setf (aref H (1- n) (1- n)) (/ q (aref H n (1- n)))
+ (aref H (1- n) n) (/ (- (- (aref H n n) p)) (aref H n (1- n)))))
+ (t
+ (cdiv A 0.0 (- (aref H (1- n) n)) (- (aref H (1- n) (1- n)) p) q)
+ (setf (aref H (1- n) (1- n)) cdivr
+ (aref H (1- n) n) cdivi)))
+ (setf (aref H n (1- n)) 0.0
+ (aref H n n) 1.0)
+ (loop :for i :from (- n 2) :downto 0
+ :do (let (ra sa vr vi)
+ (setf ra 0.0
+ sa 0.0)
+ (loop :for j :from l :to n
+ :do (setf ra (+ ra (* (aref H i j) (aref H j (1- n))))
+ sa (+ sa (* (aref H i j) (aref H j n)))))
+ (setf w (- (aref H i i) p))
+ (cond ((minusp (aref e i))
+ (setf z w
+ r ra
+ s sa))
+ (t
+ (setf l i)
+ (cond ((zerop (aref e i))
+ (cdiv A (- ra) (- sa) w q)
+ (setf (aref H i (1- n)) cdivr
+ (aref H i n) cdivi))
+ (t
+ ;;
+ (setf x (aref H i (1+ i))
+ y (aref H (1+ i) i)
+ vr (- (+ (* (- (aref d i) p) (- (aref d i) p)) (* (aref e i) (aref e i))) (* q q))
+ vi (* (- (aref d i) p) 2.0 q))
+ (cond ((> (abs x) (+ (abs z) (abs q)))
+ (setf (aref H (1+ i) (1- n)) (/ (+ (- (- ra) (* w (aref H i (1- n)))) (* q (aref H i n))) x)
+ (aref H (1+ i) n) (/ (- (- sa) (* w (aref H i n)) (* q (aref H i (1- n)))) x)))
+ (t
+ (cdiv A (- (- r) (* y (aref H i (1- n)))) (- (- s) (* y (aref H i n))) z q)
+ (setf (aref H (1+ i) (1- n)) cdivr
+ (aref H (1+ i) n) cdivi)))))
+ ;;
+ (setf rt (max (abs (aref H i (1- n))) (abs (aref H i n))))
+ (when (> (* (* eps rt) rt) 1)
+ (loop :for j :from i :to n
+ :do (setf (aref H j (1- n)) (/ (aref H j (1- n)) rt)
+ (aref H j n) (/ (aref H j n) rt))))))))
+ ;;
+ (loop :for i :from 0 :below nn
+ :when (or (< i low) (> i high))
+ :do (loop :for j :from i :below nn
+ :do (setf (aref V i j) (aref H i j))))
+ ;;
+ (loop :for j :from (1- nn) :downto low
+ :do (loop :for i :from low :to high
+ :do (setf z 0.0)
+ (loop :for k :from low :to (min j high)
+ :do (setf z (+ z (* (aref V i k) (aref H k j)))))
+ (setf (aref V i j) z))))))))))
+
+(defmethod initialize-instance :after ((obj Eigenvalue) &key matrix)
+ (with-slots (n issymmetric d e V H ort) obj
+ (setf n (array-dimension matrix 1)
+ V (make-array (list n n))
+ d (make-array (list n))
+ e (make-array (list n))
+ issymmetric t)
+ (loop :for j :from 0 :below n
+ :while issymmetric
+ :do (loop :for i :from 0 :below n
+ :while issymmetric
+ :do (setf issymmetric (= (aref matrix i j) (aref matrix j i)))))
+ (cond (issymmetric
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below n
+ :do (setf (aref V i j) (aref matrix i j))))
+ ;;
+ (tred2 obj)
+ ;;
+ (tql2 obj))
+ (t
+ (setf H (make-array (list n n))
+ ort (make-array (list n)))
+ (loop :for j :from 0 :below n
+ :do (loop :for i :from 0 :below n
+ :do (setf (aref H i j) (aref matrix i j))))
+ ;;
+ (orthes obj)
+ ;;
+ (hqr2 obj)))))
+
+(defmethod get-V ((A Eigenvalue))
+ (with-slots (V) A
+ (copy-array V)))
+
+(defmethod get-RealEigenvalues ((A Eigenvalue))
+ (with-slots (d) A
+ (copy-array d)))
+
+(defmethod get-ImagEigenvalues ((A Eigenvalue))
+ (with-slots (e) A
+ (copy-array e)))
+
+(defmethod get-D ((A Eigenvalue))
+ (with-slots (n d e) A
+ (let ((D_ (make-array (list n n))))
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below n
+ :do (setf (aref D_ i j) 0.0))
+ (setf (aref D_ i i) (aref d i))
+ (cond ((plusp (aref e i))
+ (setf (aref D_ i (1+ i)) (aref e i)))
+ ((minusp (aref e i))
+ (setf (aref D_ i (1- i)) (aref e i)))))
+ D_)))
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(in-package :cl-jama)
+
+(defclass LU ()
+ ((LU_)
+ (m)
+ (n)
+ (pivsign)
+ (piv)))
+
+(defun permute-copy-matrix (A piv j0 j1)
+ (let ((X (make-array (list (array-dimension piv 0) (1+ (- j1 j0))))))
+ (loop :for i :from 0 :below (array-dimension piv 0)
+ :do (loop :for j :from j0 :to j1
+ :do (setf (aref X i (- j j0)) (aref A (aref piv i) j))))
+ X))
+
+(defun permute-copy-vector (A piv)
+ (unless (eql (array-dimension piv 0) (array-dimension A 0))
+ (return-from permute-copy-vector nil))
+ (let ((x (make-array (list (array-dimension piv 0)))))
+ (loop :for i :from 0 :below (array-dimension piv 0)
+ :do (setf (aref x i) (aref A (aref piv i))))
+ x))
+
+(defmethod initialize-instance :after ((obj LU) &key matrix)
+ (with-slots (LU_ m n pivsign piv) obj
+ (setf LU_ (copy-array matrix))
+ (setf m (array-dimension matrix 0))
+ (setf n (array-dimension matrix 1))
+ (setf piv (make-array (list (array-dimension matrix 0))))
+ ;; Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+ (loop :for i :from 0 :below m
+ :do (setf (aref piv i) i))
+ (setf pivsign 1)
+ (let ((LUcolj (make-array (list m)))
+ p)
+ ;; Outer loop.
+ (loop :for j :from 0 :below n
+ ;; Make a copy of the j-th column to localize references.
+ :do (loop :for i :from 0 :below m
+ :do (setf (aref LUcolj i) (aref LU_ i j)))
+ ;; Apply previous transformations.
+ :do (loop :for i :from 0 :below m
+ :do (decf (aref LUcolj i) (loop :for k :from 0 :below (min i j)
+ :sum (* (aref LU_ i k) (aref LUcolj k))))
+ (setf (aref LU_ i j) (aref LUcolj i)))
+ ;; Find pivot and exchange if necessary.
+ :do (setf p j)
+ :do (loop :for i :from (1+ j) :below m
+ :when (> (abs (aref LUcolj i)) (abs (aref LUcolj p)))
+ :do (setf p i))
+ :do (unless (eql p j)
+ (loop :for k :from 0 :below n
+ :do (psetf (aref LU_ p k) (aref LU_ j k)
+ (aref LU_ j k) (aref LU_ p k)))
+ (psetf (aref piv p) (aref piv j)
+ (aref piv j) (aref piv p))
+ (setf pivsign (- pivsign)))
+ ;; Compute multipliers.
+ :do (when (and (< j m) (not (zerop (aref LU_ j j))))
+ (loop :for i :from (1+ j) :below m
+ :do (setf (aref LU_ i j) (/ (aref LU_ i j) (aref LU_ j j)))))))))
+
+(defmethod is-nonsingular ((A LU))
+ (with-slots (LU_ n) A
+ (loop :for j :from 0 :below n
+ :when (eql (aref LU_ j j) 0)
+ :do (return-from is-nonsingular nil)))
+ t)
+
+(defmethod get-L ((A LU))
+ (with-slots (LU_ m n) A
+ (let ((L_ (make-array (list m n))))
+ (loop :for i :from 0 :below m
+ :do (loop :for j :from 0 :below n
+ :if (> i j)
+ :do (setf (aref L_ i j) (aref LU_ i j))
+ :else
+ :if (= i j)
+ :do (setf (aref L_ i j) 1.0)
+ :else
+ :do (setf (aref L_ i j) 0.0)))
+ L_)))
+
+(defmethod get-U ((A LU))
+ (with-slots (LU_ n) A
+ (let ((U_ (make-array (list n n))))
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below n
+ :if (<= i j)
+ :do (setf (aref U_ i j) (aref LU_ i j))
+ :else
+ :do (setf (aref U_ i j) 0.0)))
+ U_)))
+
+(defmethod get-pivot ((A LU))
+ (with-slots (piv) A
+ (copy-array piv)))
+
+(defmethod det ((A LU))
+ (with-slots (m n pivsign LU_) A
+ (when (/= m n)
+ (return-from det 0))
+ (apply #'* pivsign (loop :for j :from 0 :below n
+ :collect (aref LU_ j j)))))
+
+(defmethod solve ((A LU) (B array))
+ (with-slots (LU_ m n piv) A
+ (unless (eql m (array-dimension B 0))
+ (return-from solve nil))
+ (unless (is-nonsingular A)
+ (return-from solve nil))
+ (let* ((nx (array-dimension B 1))
+ (X (permute-copy-matrix B piv 0 (1- nx))))
+ (loop :for k :from 0 :below n
+ :do (loop :for i :from (1+ k) :below n
+ :do (loop :for j :from 0 :below nx
+ :do (decf (aref X i j) (* (aref X k j) (aref LU_ i k))))))
+ (loop :for k :from (1- n) :downto 0
+ :do (loop :for j :from 0 :below nx
+ :do (setf (aref X k j) (/ (aref X k j) (aref LU_ k k))))
+ :do (loop :for i :from 0 :below k
+ :do (loop :for j :from 0 :below nx
+ :do (decf (aref X i j) (* (aref X k j) (aref LU_ i k))))))
+ X)))
+
+(defmethod solve ((A LU) (b vector))
+ (with-slots (LU_ m n piv) A
+ (unless (eql m (array-dimension b 0))
+ (return-from solve nil))
+ (unless (is-nonsingular A)
+ (return-from solve nil))
+ (let ((x (permute-copy-vector b piv)))
+ (loop :for k :from 0 :below n
+ :do (loop :for i :from (1+ k) :below n
+ :do (decf (aref x i) (* (aref x k) (aref LU_ i k)))))
+ (loop :for k :from (1- n) :downto 0
+ :do (setf (aref x k) (/ (aref x k) (aref LU_ k k)))
+ (loop :for i :from 0 :below k
+ :do (decf (aref x i) (* (aref x k) (aref LU_ i k)))))
+ x)))
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(defpackage :cl-jama
+ (:use :common-lisp)
+ (:nicknames :jama)
+ (:export :copy-array :transpose :multiply
+ :LU :is-nonsingular :get-L :get-U :get-pivot :det :solve
+ :Cholesky :get-L :is-spd
+ :QR :is-FullRank :get-Householder :get-R :get-Q
+ :SVD :get-V :get-singular-values :get-S :norm2 :norm2-of-condition-number :rank
+ :Eigenvalue :get-RealEigenvalues :get-ImagEigenvalues :get-D)
+ (:documentation ""))
+
+(in-package :cl-jama)
+
+(defgeneric cdiv (A xr xi yr yi)
+ )
+(defgeneric det (A)
+ )
+(defgeneric get-D (A)
+ )
+(defgeneric get-Householder (A)
+ )
+(defgeneric get-L (A)
+ )
+(defgeneric get-Q (A)
+ )
+(defgeneric get-R (A)
+ )
+(defgeneric get-RealEigenvalues (A)
+ )
+(defgeneric get-ImagEigenvalues (A)
+ )
+(defgeneric get-S (A)
+ )
+(defgeneric get-U (A)
+ )
+(defgeneric get-V (A)
+ )
+(defgeneric get-pivot (A)
+ )
+(defgeneric get-singular-values (A)
+ )
+(defgeneric hqr2 (A)
+ )
+(defgeneric is-FullRank (A)
+ )
+(defgeneric is-nonsingular (A)
+ )
+(defgeneric is-spd (A)
+ )
+(defgeneric norm2 (A)
+ )
+(defgeneric norm2-of-condition-number (A)
+ )
+(defgeneric orthes (A)
+ )
+(defgeneric rank (A)
+ )
+(defgeneric solve (A B)
+ )
+(defgeneric tql2 (A)
+ )
+(defgeneric tred2 (A)
+ )
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(in-package :cl-jama)
+
+(defclass QR ()
+ ((QR_)
+ (m)
+ (n)
+ (Rdiag)))
+
+(defmethod initialize-instance :after ((obj QR) &key matrix)
+ (with-slots (QR_ m n Rdiag) obj
+ (setf QR_ (copy-array matrix))
+ (setf m (array-dimension matrix 0))
+ (setf n (array-dimension matrix 1))
+ (setf Rdiag (make-array (list n)))
+ ;;
+ (loop :for k :from 0 :below n
+ :do (let ((nrm 0.0))
+ (loop :for i :from k :below m
+ :do (setf nrm (hypot nrm (aref QR_ i k))))
+ (unless (zerop nrm)
+ (when (minusp (aref QR_ k k))
+ (setf nrm (- nrm)))
+ (loop :for i :from k :below m
+ :do (setf (aref QR_ i k) (/ (aref QR_ i k) nrm)))
+ (incf (aref QR_ k k) 1.0)
+ ;;
+ (loop :for j :from (1+ k) :below n
+ :do (let ((s 0.0))
+ (loop :for i :from k :below m
+ :do (incf s (* (aref QR_ i k) (aref QR_ i j))))
+ (setf s (/ (- s) (aref QR_ k k)))
+ (loop :for i :from k :below m
+ :do (incf (aref QR_ i j) (* s (aref QR_ i k)))))))
+ (setf (aref Rdiag k) (- nrm))))))
+
+(defmethod is-FullRank ((A QR))
+ (with-slots (n Rdiag) A
+ (loop :for j :from 0 :below n
+ :when (zerop (aref Rdiag j))
+ :do (return-from is-FullRank nil))
+ t))
+
+(defmethod get-Householder ((A QR))
+ (with-slots (QR_ m n) A
+ (let ((H (make-array (list m n))))
+ (loop :for i :from 0 :below m
+ :do (loop :for j :from 0 :below n
+ :if (>= i j)
+ :do (setf (aref H i j) (aref QR_ i j))
+ :else
+ :do (setf (aref H i j) 0.0)))
+ H)))
+
+(defmethod get-R ((A QR))
+ (with-slots (QR_ n Rdiag) A
+ (let ((R (make-array (list n n))))
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below n
+ :if (< i j)
+ :do (setf (aref R i j) (aref QR_ i j))
+ :else
+ :if (= i j)
+ :do (setf (aref R i j) (aref Rdiag i))
+ :else
+ :do (setf (aref R i j) 0.0)))
+ R)))
+
+(defmethod get-Q ((A QR))
+ (with-slots (QR_ m n) A
+ (let ((Q (make-array (list m n))))
+ (loop :for k :from (1- n) :downto 0
+ :do (loop :for i :from 0 :below m
+ :do (setf (aref Q i k) 0.0))
+ (setf (aref Q k k) 1.0)
+ (loop :for j :from k :below n
+ :do (unless (zerop (aref QR_ k k))
+ (let ((s 0.0))
+ (loop :for i :from k :below m
+ :do (incf s (* (aref QR_ i k) (aref Q i j))))
+ (setf s (/ (- s) (aref QR_ k k)))
+ (loop :for i :from k :below m
+ :do (incf (aref Q i j) (* s (aref QR_ i k))))
+ ))))
+ Q)))
+
+(defmethod solve ((A QR) (b vector))
+ (with-slots (QR_ m n Rdiag) A
+ (unless (eql m (array-dimension b 0))
+ (return-from solve nil))
+ (unless (is-FullRank A)
+ (return-from solve nil))
+ (let ((x (copy-array b))
+ (x_ (make-array (list n))))
+ (loop :for k :from 0 :below n
+ :do (let ((s 0.0))
+ (loop :for i :from k :below m
+ :do (incf s (* (aref QR_ i k) (aref x i))))
+ (setf s (/ (- s) (aref QR_ k k)))
+ (loop :for i :from k :below m
+ :do (incf (aref x i) (* s (aref QR_ i k))))))
+ (loop :for k :from (1- n) :downto 0
+ :do (setf (aref x k) (/ (aref x k) (aref Rdiag k)))
+ (loop :for i :from 0 :below k
+ :do (decf (aref x i) (* (aref x k) (aref QR_ i k)))))
+ (loop :for i :from 0 :below n
+ :do (setf (aref x_ i) (aref x i)))
+ x_)))
+
+(defmethod solve ((A QR) (B array))
+ (with-slots (QR_ m n Rdiag) A
+ (unless (eql m (array-dimension B 0))
+ (return-from solve nil))
+ (unless (is-FullRank A)
+ (return-from solve nil))
+ (let* ((nx (array-dimension B 1))
+ (X (copy-array B))
+ (X_ (make-array (list n nx))))
+ (loop :for k :from 0 :below n
+ :do (loop :for j :from 0 :below nx
+ :do (let ((s 0.0))
+ (loop :for i :from k :below m
+ :do (incf s (* (aref QR_ i k) (aref X i j))))
+ (setf s (/ (- s) (aref QR_ k k)))
+ (loop :for i :from k :below m
+ :do (incf (aref X i j) (* s (aref QR_ i k)))))))
+ (loop :for k :from (1- n) :downto 0
+ :do (loop :for j :from 0 :below nx
+ :do (setf (aref X k j) (/ (aref X k j) (aref Rdiag k))))
+ (loop :for i :from 0 :below k
+ :do (loop :for j :from 0 :below nx
+ :do (decf (aref X i j) (* (aref X k j) (aref QR_ i k))))))
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below nx
+ :do (setf (aref X_ i j) (aref X i j))))
+ X_)))
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(in-package :cl-jama)
+
+(defclass SVD ()
+ ((U)
+ (V)
+ (s)
+ (m)
+ (n)))
+
+(defmethod initialize-instance :after ((obj SVD) &key matrix)
+ (with-slots (U V s m n) obj
+ (let (nu e work A (wantu t) (wantv t) nct nrt p pp eps iter)
+ (setf m (array-dimension matrix 0))
+ (setf n (array-dimension matrix 1))
+ (setf nu (min m n))
+ (setf s (make-array (list (min (1+ m) n))))
+ (setf U (make-array (list m nu) :initial-element 0.0))
+ (setf V (make-array (list n n)))
+ (setf e (make-array (list n)))
+ (setf work (make-array (list m)))
+ (setf A (copy-array matrix))
+ ;;
+ (setf nct (min (1- m) n))
+ (setf nrt (max 0 (min (- n 2) m)))
+ (loop :for k :from 0 :below (max nct nrt)
+ :do (when (< k nct)
+ ;;
+ (setf (aref s k) 0)
+ (loop :for i :from k :below m
+ :do (setf (aref s k) (hypot (aref s k) (aref A i k))))
+ (unless (zerop (aref s k))
+ (when (minusp (aref A k k))
+ (setf (aref s k) (- (aref s k))))
+ (loop :for i :from k :below m
+ :do (setf (aref A i k) (/ (aref A i k) (aref s k))))
+ (incf (aref A k k) 1.0))
+ (setf (aref s k) (- (aref s k))))
+ (loop :for j :from (1+ k) :below n
+ :do (when (and (< k nct) (not (zerop (aref s k))))
+ ;;
+ (let ((rt 0.0))
+ (loop :for i :from k :below m
+ :do (incf rt (* (aref A i k) (aref A i j))))
+ (setf rt (/ (- rt) (aref A k k)))
+ (loop :for i :from k :below m
+ :do (incf (aref A i j) (* rt (aref A i k))))))
+ (setf (aref e j) (aref A k j)))
+ (when (and wantu (< k nct))
+ (loop :for i :from k :below m
+ :do (setf (aref U i k) (aref A i k))))
+ (when (< k nrt)
+ ;;
+ (setf (aref e k) 0)
+ (loop :for i :from (1+ k) :below n
+ :do (setf (aref e k) (hypot (aref e k) (aref e i))))
+ (unless (zerop (aref e k))
+ (when (minusp (aref e (1+ k)))
+ (setf (aref e k) (- (aref e k))))
+ (loop :for i :from (1+ k) :below n
+ :do (setf (aref e i) (/ (aref e i) (aref e k))))
+ (incf (aref e (1+ k)) 1.0))
+ (setf (aref e k) (- (aref e k)))
+ (when (and (< (1+ k) m) (not (zerop (aref e k))))
+ ;;
+ (loop :for i :from (1+ k) :below m
+ :do (setf (aref work i) 0.0))
+ (loop :for j :from (1+ k) :below n
+ :do (loop :for i :from (1+ k) :below m
+ :do (incf (aref work i) (* (aref e j) (aref A i j)))))
+ (loop :for j :from (1+ k) :below n
+ :do (let ((rt (/ (- (aref e j)) (aref e (1+ k)))))
+ (loop :for i :from (1+ k) :below m
+ :do (incf (aref A i j) (* rt (aref work i)))))))
+ (when wantv
+ ;;
+ (loop :for i :from (1+ k) :below n
+ :do (setf (aref V i k) (aref e i))))))
+ ;;
+ (setf p (min n (1+ m)))
+ (when (< nct n)
+ (setf (aref s nct) (aref A nct nct)))
+ (when (< m p)
+ (setf (aref s (1- p)) 0.0))
+ (when (< (1+ nrt) p)
+ (setf (aref e nrt) (aref A nrt (1- p))))
+ (setf (aref e (1- p)) 0.0)
+ ;;
+ (when wantu
+ (loop :for j :from nct :below nu
+ :do (loop :for i :from 0 :below m
+ :do (setf (aref U i j) 0.0))
+ (setf (aref U j j) 1.0))
+ (loop :for k :from (1- nct) :downto 0
+ :do (cond ((not (zerop (aref s k)))
+ (loop :for j :from (1+ k) :below nu
+ :do (let ((rt 0.0))
+ (loop :for i :from k :below m
+ :do (incf rt (* (aref U i k) (aref U i j))))
+ (setf rt (/ (- rt) (aref U k k)))
+ (loop :for i :from k :below m
+ :do (incf (aref U i j) (* rt (aref U i k))))))
+ (loop :for i :from k :below m
+ :do (setf (aref U i k) (- (aref U i k))))
+ (setf (aref U k k) (+ 1.0 (aref U k k)))
+ (loop :for i :from 0 :below (1- k)
+ :do (setf (aref U i k) 0.0)))
+ (t
+ (loop :for i :from 0 :below m
+ :do (setf (aref U i k) 0.0))
+ (setf (aref U k k) 1.0)))))
+ ;;
+ (when wantv
+ (loop :for k :from (1- n) :downto 0
+ :do (when (and (< k nrt) (not (zerop (aref e k))))
+ (loop :for j :from (1+ k) :below nu
+ :do (let ((rt 0.0))
+ (loop :for i :from (1+ k) :below n
+ :do (incf rt (* (aref V i k) (aref V i j))))
+ (setf rt (/ (- rt) (aref V (1+ k) k)))
+ (loop :for i :from (1+ k) :below n
+ :do (incf (aref V i j) (* rt (aref V i k)))))))
+ (loop :for i :from 0 :below n
+ :do (setf (aref V i k) 0.0))
+ (setf (aref V k k) 1.0)))
+ ;;
+ (setf pp (1- p))
+ (setf iter 0)
+ (setf eps (expt 2.0 -52.0))
+ (loop :while (plusp p)
+ :do (let ((k 0) (kase 0))
+ ;;
+ (loop :for lk :from (- p 2) :downto -1
+ :when (= lk -1)
+ :return (setf k lk)
+ :end
+ :when (<= (abs (aref e lk)) (* eps (+ (abs (aref s lk)) (abs (aref s (1+ lk))))))
+ :return (setf (aref e lk) 0.0
+ k lk)
+ :end
+ :finally (setf k lk)) ; XXX
+ (cond ((= k (- p 2))
+ (setf kase 4))
+ (t
+ (let ((ks 0))
+ (loop :for lks :from (1- p) :downto k
+ :when (= lks k)
+ :return (setf ks lks)
+ :end
+ :when (<= (abs (aref s lks)) (* eps (+ (if (/= lks p) (abs (aref e lks)) 0.0)
+ (if (/= lks (1+ k)) (abs (aref e (1- lks))) 0.0))))
+ :return (setf (aref s lks) 0.0
+ ks lks)
+ :end
+ :finally (setf ks lks)) ; XXX
+ (cond ((= ks k)
+ (setf kase 3))
+ ((= ks (1- p))
+ (setf kase 1))
+ (t
+ (setf kase 2)
+ (setf k ks))))))
+ (incf k)
+ ;;
+ ;; switch
+ (case kase
+ (1 (let ((f (aref e (- p 2))))
+ (setf (aref e (- p 2)) 0.0)
+ (loop :for j :from (- p 2) :downto k
+ :for rt = (hypot (aref s j) f) :then (hypot (aref s j) f)
+ :for cs = (/ (aref s j) rt) :then (/ (aref s j) rt)
+ :for sn = (/ f rt) :then (/ f rt)
+ :do (setf (aref s j) rt)
+ (unless (eql j k)
+ (setf f (* (- sn) (aref e (1- j))))
+ (setf (aref e (1- j)) (* cs (aref e (1- j)))))
+ (when wantv
+ (loop :for i :from 0 :below n
+ :do (setf rt (+ (* cs (aref V i j)) (* sn (aref V i (1- p)))))
+ (setf (aref V i (1- p)) (+ (* (- sn) (aref V i j)) (* cs (aref V i (1- p)))))
+ (setf (aref V i j) rt))))))
+ (2 (let ((f (aref e (1- k))))
+ (setf (aref e (1- k)) 0.0)
+ (loop :for j :from k :below p
+ :for rt = (hypot (aref s j) f) :then (hypot (aref s j) f)
+ :for cs = (/ (aref s j) rt) :then (/ (aref s j) rt)
+ :for sn = (/ f rt) :then (/ f rt)
+ :do (setf (aref s j) rt)
+ (setf f (* (- sn) (aref e j)))
+ (setf (aref e j) (* cs (aref e j)))
+ (when wantu
+ (loop :for i :from 0 :below m
+ :do (setf rt (+ (* cs (aref U i j)) (* sn (aref U i (1- k)))))
+ (setf (aref U i (1- k)) (+ (* (- sn) (aref U i j)) (* cs (aref U i (1- k)))))
+ (setf (aref U i j) rt))))))
+ (3 (let* ((scale (max (abs (aref s (- p 1))) (abs (aref s (- p 2))) (abs (aref e (- p 2)))
+ (abs (aref s k)) (abs (aref e k))))
+ (sp (/ (aref s (1- p)) scale))
+ (spm1 (/ (aref s (- p 2)) scale))
+ (epm1 (/ (aref e (- p 2)) scale))
+ (sk (/ (aref s k) scale))
+ (ek (/ (aref e k) scale))
+ (b (/ (+ (* (+ spm1 sp) (- spm1 sp)) (* epm1 epm1)) 2.0))
+ (c (* sp epm1 sp epm1))
+ (shift 0.0)
+ f g)
+ (when (or (not (zerop b)) (not (zerop c)))
+ (setf shift (sqrt (+ (* b b) c)))
+ (when (minusp b)
+ (setf shift (- shift)))
+ (setf shift (/ c (+ b shift))))
+ (setf f (+ (* (+ sk sp) (- sk sp)) shift))
+ (setf g (* sk ek))
+ ;;
+ (loop :for j :from k :below (1- p)
+ :do (let* ((rt (hypot f g))
+ (cs (/ f rt))
+ (sn (/ g rt)))
+ (when (/= j k)
+ (setf (aref e (1- j)) rt))
+ (setf f (+ (* cs (aref s j)) (* sn (aref e j)))
+ (aref e j) (- (* cs (aref e j)) (* sn (aref s j)))
+ g (* sn (aref s (1+ j)))
+ (aref s (1+ j)) (* cs (aref s (1+ j))))
+ (when wantv
+ (loop :for i :from 0 :below n
+ :do (setf rt (+ (* cs (aref V i j)) (* sn (aref V i (1+ j))))
+ (aref V i (1+ j)) (+ (* (- sn) (aref V i j)) (* cs (aref V i (1+ j))))
+ (aref V i j) rt)))
+ (setf rt (hypot f g))
+ (setf cs (/ f rt))
+ (setf sn (/ g rt))
+ (setf (aref s j) rt)
+ (setf f (+ (* cs (aref e j)) (* sn (aref s (1+ j)))))
+ (setf (aref s (1+ j)) (+ (* (- sn) (aref e j)) (* cs (aref s (1+ j)))))
+ (setf g (* sn (aref e (1+ j))))
+ (setf (aref e (1+ j)) (* cs (aref e (1+ j))))
+ (when (and wantu (< j (1- m)))
+ (loop :for i :from 0 :below m
+ :do (setf rt (+ (* cs (aref U i j)) (* sn (aref U i (1+ j))))
+ (aref U i (1+ j)) (+ (* (- sn) (aref U i j)) (* cs (aref U i (1+ j))))
+ (aref U i j) rt)))))
+ (setf (aref e (- p 2)) f)
+ (setf iter (1+ iter))))
+ (4 (when (<= (aref s k) 0.0)
+ (setf (aref s k) (if (minusp (aref s k)) (- (aref s k)) 0.0))
+ (when wantv
+ (loop :for i :from 0 :to pp
+ :do (setf (aref V i k) (- (aref V i k))))))
+ (loop :while (< k pp)
+ :thereis (>= (aref s k) (aref s (1+ k)))
+ :do (psetf (aref s k) (aref s (1+ k))
+ (aref s (1+ k)) (aref s k))
+ :when (and wantv (< k (1- n)))
+ :do (loop :for i :from 0 :below n
+ :do (psetf (aref V i (1+ k)) (aref V i k)
+ (aref V i k) (aref V i (1+ k))))
+ :end
+ :when (and wantu (< k (1- m)))
+ :do (loop :for i :from 0 :below m
+ :do (psetf (aref U i (1+ k)) (aref U i k)
+ (aref U i k) (aref U i (1+ k))))
+ :do (incf k))
+ (setf iter 0)
+ (decf p))))))))
+
+(defmethod get-U ((A SVD))
+ (with-slots (U m n) A
+ (let* ((minm (min (1+ m) n))
+ (X (make-array (list m minm))))
+ (loop :for i :from 0 :below m
+ :do (loop :for j :from 0 :below minm
+ :do (setf (aref X i j) (aref U i j))))
+ X)))
+
+(defmethod get-V ((A SVD))
+ (with-slots (V) A
+ (copy-array V)))
+
+(defmethod get-singular-values ((A SVD))
+ (with-slots (s) A
+ (copy-array s)))
+
+(defmethod get-S ((A SVD))
+ (with-slots (n s) A
+ (let ((X (make-array (list n n))))
+ (loop :for i :from 0 :below n
+ :do (loop :for j :from 0 :below n
+ :do (setf (aref X i j) 0.0))
+ (setf (aref X i i) (aref s i)))
+ X)))
+(defmethod norm2 ((A SVD))
+ (with-slots (s) A
+ (aref s 0)))
+
+(defmethod norm2-of-condition-number ((A SVD))
+ (with-slots (m n s) A
+ (/ (aref s 0) (aref s (1- (min m n))))))
+
+(defmethod rank ((A SVD))
+ (with-slots (m n s) A
+ (let* ((eps (expt 2.0 -52.0))
+ (tol (* (max m n) (aref s 0) eps)))
+ (loop :for i :from 0 :below (array-dimension s 0)
+ :count (> (aref s 0) tol)))))
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(in-package :cl-jama-test)
+(rem-all-tests)
+
+(defparameter *chol.1* (make-instance 'Cholesky
+ :matrix (make-matrix '((5.49 2.68 0.0 0.0)
+ (2.68 5.63 -2.39 0.0)
+ (0.0 -2.39 2.60 -2.22)
+ (0.0 0.0 -2.22 5.17)))))
+(deftest cholesky.1
+ (with-output-to-string (os)
+ (dump-vector (solve *chol.1* (make-vector '(22.09 9.31 -5.24 11.83))) os))
+" 5.0000 -2.0000 -3.0000 1.0000
+")
+
+(deftest cholesky.2
+ (with-output-to-string (os)
+ (dump-matrix (get-L *chol.1*) os))
+" 2.3431 0.0000 0.0000 0.0000
+ 1.1438 2.0789 0.0000 0.0000
+ 0.0000 -1.1497 1.1306 0.0000
+ 0.0000 0.0000 -1.9635 1.1465
+")
+
+(deftest cholesky.3
+ (with-output-to-string (os)
+ (dump-matrix (solve *chol.1* (make-matrix '((22.09 5.10)
+ (9.31 30.81)
+ (-5.24 -25.82)
+ (11.83 22.90))))
+ os))
+" 5.0000 -2.0000
+ -2.0000 6.0000
+ -3.0000 -1.0000
+ 1.0000 4.0000
+")
+
+(defparameter *chol.2* (make-instance 'Cholesky
+ :matrix (make-matrix '((4.16 -3.12 0.56 -0.10)
+ (-3.12 5.03 -0.83 1.18)
+ (0.56 -0.83 0.76 0.34)
+ (-0.10 1.18 0.34 1.18)))))
+(deftest cholesky.4
+ (with-output-to-string (os)
+ (dump-vector (solve *chol.2* (make-vector '(8.70 -13.35 1.89 -4.14))) os))
+" 1.0000 -1.0000 2.0000 -3.0000
+")
+
+(deftest cholesky.5
+ (with-output-to-string (os)
+ (dump-matrix (get-L *chol.2*) os))
+" 2.0396 0.0000 0.0000 0.0000
+ -1.5297 1.6401 0.0000 0.0000
+ 0.2746 -0.2500 0.7887 0.0000
+ -0.0490 0.6737 0.6617 0.5347
+")
+
+(deftest cholesky.6
+ (with-output-to-string (os)
+ (dump-matrix (solve *chol.2* (make-matrix '((8.70 8.30)
+ (-13.35 2.13)
+ (1.89 1.61)
+ (-4.14 5.00))))
+ os))
+" 1.0000 4.0000
+ -1.0000 3.0000
+ 2.0000 2.0000
+ -3.0000 1.0000
+")
+
+(do-tests)
--- /dev/null
+
+(asdf:defsystem :cl-jama-test
+ :depends-on (:cl-jama :rt)
+ :components ((:file "package")
+ (:file "utils" :depends-on ("package"))
+ (:file "lu" :depends-on ("utils"))
+ (:file "cholesky" :depends-on ("utils"))
+ (:file "qr" :depends-on ("utils"))
+ (:file "svd" :depends-on ("utils"))
+ (:file "eigenvalue" :depends-on ("utils"))
+ ))
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(in-package :cl-jama-test)
+(rem-all-tests)
+
+(defparameter *A.1* (make-matrix '((4.0 1.0 1.0)
+ (1.0 2.0 3.0)
+ (1.0 3.0 6.0))))
+(defparameter *eig.1* (make-instance 'Eigenvalue
+ :matrix *A.1*))
+(defparameter *D.1* (get-D *eig.1*))
+(defparameter *V.1* (get-V *eig.1*))
+
+(deftest eig.1
+ (with-output-to-string (os)
+ (dump-matrix *D.1* os))
+" 0.3451 0.0000 0.0000
+ 0.0000 3.5956 0.0000
+ 0.0000 0.0000 8.0593
+")
+
+(deftest eig.2
+ (with-output-to-string (os)
+ (dump-matrix *V.1* os))
+" 0.1195 -0.9406 0.3178
+ -0.8856 0.0436 0.4623
+ 0.4487 0.3367 0.8278
+")
+
+(deftest eig.3
+ (equal (with-output-to-string (os)
+ (dump-matrix (multiply *A.1* *V.1*) os))
+ (with-output-to-string (os)
+ (dump-matrix (multiply *V.1* *D.1*) os)))
+ T)
+
+(defparameter *A.2* (make-matrix '((0.0 1.0 0.0 0.0)
+ (1.0 0.0 2.0d-7 0.0)
+ (0.0 -2.0d-7 0.0 1.0)
+ (0.0 0.0 1.0 0.0))))
+(defparameter *eig.2* (make-instance 'Eigenvalue
+ :matrix *A.2*))
+(defparameter *D.2* (get-D *eig.2*))
+(defparameter *V.2* (get-V *eig.2*))
+
+(deftest eig.4
+ (equal (with-output-to-string (os)
+ (dump-matrix (multiply *A.2* *V.2*) os))
+ (with-output-to-string (os)
+ (dump-matrix (multiply *V.2* *D.2*) os)))
+ T)
+
+(do-tests)
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(in-package :cl-jama-test)
+
+(rem-all-tests)
+
+(defparameter *lu.1* (make-instance 'LU :matrix (make-matrix '((-0.23 2.54 -3.66 0.0)
+ (-6.98 2.46 -2.73 -2.13)
+ (0.0 2.56 2.46 4.07)
+ (0.0 0.0 -4.78 -3.82)))))
+(deftest lu.1
+ (with-output-to-string (os)
+ (dump-vector (solve *lu.1* (make-vector '(4.42 27.13 -6.14 10.50))) os))
+" -2.0000 3.0000 1.0000 -4.0000
+")
+
+(deftest lu.2
+ (with-output-to-string (os)
+ (dump-matrix (solve *lu.1* (make-matrix '((4.42 -36.01)
+ (27.13 -31.67)
+ (-6.14 -1.16)
+ (10.50 -25.82))))
+ os))
+" -2.0000 1.0000
+ 3.0000 -4.0000
+ 1.0000 7.0000
+ -4.0000 -2.0000
+")
+
+(defparameter *lu.2* (make-instance 'LU :matrix (make-matrix '((1.80 2.88 2.05 -0.89)
+ (5.25 -2.95 -0.95 -3.80)
+ (1.58 -2.69 -2.90 -1.04)
+ (-1.11 -0.66 -0.59 0.80)))))
+(deftest lu.3
+ (with-output-to-string (os)
+ (dump-vector (solve *lu.2* (make-vector '(9.52 24.35 0.77 -6.22))) os))
+" 1.0000 -1.0000 3.0000 -5.0000
+")
+
+(defparameter *lu.3* (make-instance 'lu :matrix (make-matrix '(( 1.80 2.88 2.05 -0.89)
+ (525.00 -295.00 -95.00 -380.00)
+ ( 1.58 -2.69 -2.90 -1.04)
+ ( -1.11 -0.66 -0.59 0.80)))))
+(deftest lu.4
+ (with-output-to-string (os)
+ (dump-matrix (solve *lu.3* (make-matrix '((9.52 18.47)
+ (2435.00 225.00)
+ (0.77 -13.28)
+ (-6.22 -6.21))))
+ os))
+" 1.0000 3.0000
+ -1.0000 2.0000
+ 3.0000 4.0000
+ -5.0000 1.0000
+")
+
+(do-tests)
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(defpackage :cl-jama-test
+ (:use :cl :cl-jama :rt))
+
+(in-package :cl-jama-test)
+
+(rem-all-tests)
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(in-package :cl-jama-test)
+(rem-all-tests)
+
+(defparameter *a.1* (make-matrix '((1.0d0 2.0d0) (3.0d0 4.0d0))))
+(defparameter *qr.1* (make-instance 'QR :matrix *a.1*))
+(defparameter *q.1* (get-Q *qr.1*))
+(defparameter *r.1* (get-R *qr.1*))
+
+(deftest qr.1
+ (equal (with-output-to-string (os)
+ (dump-matrix (multiply *q.1* *r.1*) os))
+ (with-output-to-string (os)
+ (dump-matrix *a.1* os)))
+ t)
+
+(deftest qr.2
+ (let ((X (multiply *q.1* (transpose *q.1*)))
+ (Y #(1 0 0 1))
+ (eps 1.0d10))
+ (reduce #'eql
+ (loop :for i :from 0 :below (array-total-size X)
+ :collect (< (abs (- (row-major-aref X i) (aref Y i))) eps))))
+ T)
+
+(deftest qr.3
+ (let ((X (multiply (transpose *q.1*) *q.1*))
+ (Y #(1 0 0 1))
+ (eps 1.0d10))
+ (reduce #'eql
+ (loop :for i :from 0 :below (array-total-size X)
+ :collect (< (abs (- (row-major-aref X i) (aref Y i))) eps))))
+ T)
+
+(defparameter *qr.2* (make-instance 'QR
+ :matrix (make-matrix '((-0.57d0 -1.28d0 -0.39d0 0.25d0)
+ (-1.93d0 1.08d0 -0.31d0 -2.14d0)
+ ( 2.30d0 0.24d0 0.40d0 -0.35d0)
+ (-1.93d0 0.64d0 -0.66d0 0.08d0)
+ ( 0.15d0 0.30d0 0.15d0 -2.13d0)
+ (-0.02d0 1.03d0 -1.43d0 0.50d0)))))
+
+(deftest qr.4
+ (with-output-to-string (os)
+ (dump-vector (solve *qr.2* (make-vector '(-2.67d0 -0.55d0 3.34d0 -0.77d0 0.48d0 4.10d0)))
+ os))
+" 1.5339 1.8707 -1.5241 0.0392
+")
+
+(do-tests)
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(in-package :cl-jama-test)
+(rem-all-tests)
+
+(defparameter *hilb* (let ((A (make-array '(3 3))))
+ (loop :for i :from 0 :below 3
+ :do (loop :for j :from 0 :below 3
+ :do (setf (aref A i j) (/ 1.0d0 (+ i j 1)))))
+ A))
+
+(defparameter *svd.1* (make-instance 'SVD :matrix *hilb*))
+(defparameter *u.1* (get-U *svd.1*))
+(defparameter *s.1* (get-S *svd.1*))
+(defparameter *v.1* (get-V *svd.1*))
+
+(deftest svd.1
+ (with-output-to-string (os)
+ (dump-matrix *u.1* os 11 5))
+" 0.82704 -0.54745 0.12766
+ 0.45986 0.52829 -0.71375
+ 0.32330 0.64901 0.68867
+")
+
+(deftest svd.2
+ (with-output-to-string (os)
+ (dump-matrix *s.1* os 11 5))
+" 1.40832 0.00000 0.00000
+ 0.00000 0.12233 0.00000
+ 0.00000 0.00000 0.00269
+")
+
+(deftest svd.3
+ (with-output-to-string (os)
+ (dump-matrix *v.1* os 11 5))
+" 0.82704 -0.54745 0.12766
+ 0.45986 0.52829 -0.71375
+ 0.32330 0.64901 0.68867
+")
+
+(deftest svd.4
+ (equal (with-output-to-string (os)
+ (dump-matrix (multiply *u.1* (multiply *s.1* (transpose *v.1*))) os 11 5))
+ (with-output-to-string (os)
+ (dump-matrix *hilb* os 11 5)))
+ T)
+
+(deftest svd.5
+ (let ((X (multiply *u.1* (transpose *u.1*)))
+ (Y #(1 0 0 0 1 0 0 0 1))
+ (eps 1.0d10))
+ (reduce #'eql
+ (loop :for i :from 0 :below (array-total-size X)
+ :collect (< (abs (- (row-major-aref X i) (aref Y i))) eps))))
+ T)
+
+(deftest svd.6
+ (let ((X (multiply *v.1* (transpose *v.1*)))
+ (Y #(1 0 0 0 1 0 0 0 1))
+ (eps 1.0d10))
+ (reduce #'eql
+ (loop :for i :from 0 :below (array-total-size X)
+ :collect (< (abs (- (row-major-aref X i) (aref Y i))) eps))))
+ T)
+
+(defparameter *svd.2* (make-instance 'SVD
+ :matrix (make-matrix '(( 2.27d0 -1.54d0 1.15d0 -1.94d0)
+ ( 0.28d0 -1.67d0 0.94d0 -0.78d0)
+ (-0.48d0 -3.09d0 0.99d0 -0.21d0)
+ ( 1.07d0 1.22d0 0.79d0 0.63d0)
+ (-2.35d0 2.93d0 -1.45d0 2.30d0)
+ ( 0.62d0 -7.39d0 1.03d0 -2.57d0)))))
+
+(deftest svd.7
+ (with-output-to-string (os)
+ (dump-vector (get-singular-values *svd.2*) os))
+" 9.9966 3.6831 1.3569 0.5000
+")
+
+(do-tests)
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+
+(in-package :cl-jama-test)
+
+(defun make-matrix (x)
+ (let ((m (length x))
+ (n (length (first x))))
+ (make-array (list m n) :initial-contents x)))
+
+(defun make-vector (x)
+ (make-array (list (length x)) :initial-contents x))
+
+(defun dump-matrix (x &optional (s *standard-output*) (w 11) (d 4))
+ (loop :for i :from 0 :below (array-dimension x 0)
+ :do (loop :for j :from 0 :below (array-dimension x 1)
+ :do (format s "~V,VF" w d (aref x i j)))
+ (terpri s))
+ (values))
+
+(defun dump-vector (x &optional (s *standard-output*))
+ (loop :for i :from 0 :below (array-dimension x 0)
+ :do (format s "~11,4F" (aref x i)))
+ (terpri s)
+ (values))
+
+(rem-all-tests)
+
+(deftest make-matrix.1
+ (make-matrix '((1.23 -2.34 3.45) (-4.56 5.67 -6.78)))
+ #2A((1.23 -2.34 3.45) (-4.56 5.67 -6.78)))
+
+(deftest make-vector.1
+ (make-vector '(9.87 -8.76 7.65))
+ #(9.87 -8.76 7.65))
+
+(deftest dump-matrix.1
+ (with-output-to-string (os) (dump-matrix #2A((1.23 -2.34 3.45) (-4.56 5.67 -6.78)) os))
+" 1.2300 -2.3400 3.4500
+ -4.5600 5.6700 -6.7800
+")
+
+(deftest dump-vector.1
+ (with-output-to-string (os) (dump-vector #(9.87 -8.76 7.65) os))
+" 9.8700 -8.7600 7.6500
+")
+
+(deftest copy-array.1
+ (copy-array #(1.23 -2.34 3.45))
+ #(1.23 -2.34 3.45))
+
+(deftest copy-array.2
+ (copy-array #2A((-4.56 5.67 -6.78) (7.89 -8.90 9.01)))
+ #2A((-4.56 5.67 -6.78) (7.89 -8.90 9.01)))
+
+(deftest transpose.1
+ (transpose #2A((-4.56 5.67 -6.78) (7.89 -8.90 9.01)))
+ #2A((-4.56 7.89) (5.67 -8.90) (-6.78 9.01)))
+
+(deftest multiply.1
+ (multiply #2A((1 1 1) (2 2 2) (3 3 3)) #2A((1 2 3) (1 2 3) (1 2 3)))
+ #2A((3 6 9) (6 12 18) (9 18 27)))
+
+(do-tests)
--- /dev/null
+;;; -*- mode: lisp; coding: utf-8 -*-
+;;;
+;;; Copyright (c) 2008 NANRI <southly@gmail.com>
+;;;
+;;; Permission is hereby granted, free of charge, to any person obtaining a copy
+;;; of this software and associated documentation files (the "Software"), to deal
+;;; in the Software without restriction, including without limitation the rights
+;;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+;;; copies of the Software, and to permit persons to whom the Software is
+;;; furnished to do so, subject to the following conditions:
+;;;
+;;; The above copyright notice and this permission notice shall be included in
+;;; all copies or substantial portions of the Software.
+;;;
+;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+;;; THE SOFTWARE.
+
+(in-package :cl-jama)
+
+(defun copy-array (A)
+ (let ((X (make-array (array-dimensions A))))
+ (loop :for i :from 0 :below (array-total-size A)
+ :do (setf (row-major-aref X i) (row-major-aref A i)))
+ X))
+
+(defun hypot (a b)
+ (cond ((> (abs a) (abs b))
+ (* (abs a) (sqrt (1+ (expt (/ b a) 2)))))
+ ((not (zerop b))
+ (* (abs b) (sqrt (1+ (expt (/ a b) 2)))))
+ (t
+ 0.0)))
+
+(defun transpose (X)
+ (let* ((m (array-dimension X 0))
+ (n (array-dimension X 1))
+ (Y (make-array (list n m))))
+ (loop :for i :from 0 :below m
+ :do (loop :for j :from 0 :below n
+ :do (setf (aref Y j i) (aref X i j))))
+ Y))
+
+(defun multiply (A B)
+ (unless (eql (array-dimension A 1) (array-dimension B 0))
+ (return-from multiply nil))
+ (let* ((C (make-array (list (array-dimension A 0) (array-dimension B 1)))))
+ (loop :for i :from 0 :below (array-dimension A 0)
+ :do (loop :for j :from 0 :below (array-dimension B 1)
+ :do (let ((sum 0))
+ (loop :for k :from 0 :below (array-dimension A 1)
+ :do (incf sum (* (aref A i k) (aref B k j))))
+ (setf (aref C i j) sum))))
+ C))
OSDN Git Service
