From: Barton W. <wil...@us...> - 2006-01-08 15:41:13
|
Update of /cvsroot/maxima/maxima/share/linearalgebra In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv5988/share/linearalgebra Modified Files: linalg.mac linalgcholesky.lisp lu.lisp mring.lisp Log Message: (1) change psqrt in *rationalfield* so that it returns the square root when the input is a perfect square (2) added ratdisrep to the crering psqrt function (3) fixed bugs in $zerofor and $identfor (lu.lisp) -- these functions needed to convert to maxima representation (4) added function full-matrix-map to lu.lisp (5) linalg now loads linalgcholesky (6) fixed bug in *bigfloatfield* psqrt Index: linalg.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/linalg.mac,v retrieving revision 1.10 retrieving revision 1.11 diff -u -d -r1.10 -r1.11 --- linalg.mac 7 Jan 2006 19:50:38 -0000 1.10 +++ linalg.mac 8 Jan 2006 15:41:01 -0000 1.11 @@ -33,6 +33,7 @@ eval_when([batch,load,translate,compile], load("polynomialp"), load("lu"), + load("linalgcholesky"), load("matrixexp"), load("linalg-utilities")); Index: linalgcholesky.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/linalgcholesky.lisp,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- linalgcholesky.lisp 8 Jan 2006 04:09:01 -0000 1.3 +++ linalgcholesky.lisp 8 Jan 2006 15:41:01 -0000 1.4 @@ -1,4 +1,4 @@ -;; Copyright 2005 by Barton Willis +;; Copyright 2005, 2006 by Barton Willis ;; This is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License, @@ -34,6 +34,7 @@ (freciprocal (mring-reciprocal fld))) (setq l ($zerofor m)) + (full-matrix-map l fconvert) (loop for k from 1 to n do (push k perm)) @@ -51,7 +52,6 @@ (setmatelem l lii i i) (if (funcall fzerop lii) (merror "Unable to find the Cholesky factorization")) (setq lii-inv (funcall fadjoint (funcall freciprocal lii))) - (loop for j from (+ i 1) to n do (setq acc (fn-melem m perm j i fconvert)) (loop for k from 1 to (- i 1) do @@ -59,7 +59,7 @@ (m-elem l perm j k) (funcall fadjoint (m-elem l perm i k)))))) (setmatelem l (funcall fmult acc lii-inv) j i))) - (matrix-map l n n (mring-mring-to-maxima fld)) + (full-matrix-map l (mring-mring-to-maxima fld)) l)) Index: lu.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/lu.lisp,v retrieving revision 1.5 retrieving revision 1.6 diff -u -d -r1.5 -r1.6 --- lu.lisp 1 Jan 2006 13:28:32 -0000 1.5 +++ lu.lisp 8 Jan 2006 15:41:01 -0000 1.6 @@ -31,6 +31,20 @@ (loop for j from 1 to c do (setmatelem m (funcall fn (nth j (nth i m))) i j)))) +;; Map the lisp function fn over the r by c Maxima matrix m. This function is +;; block matrix friendly. + +(defun full-matrix-map (mat fn) + (cond (($matrixp mat) + (let ((r) (c)) + (setq r ($matrix_size mat)) + (setq c ($second r)) + (setq r ($first r)) + (loop for i from 1 to r do + (loop for j from 1 to c do + (setmatelem mat (full-matrix-map (nth j (nth i mat)) fn) i j))))) + (t (funcall fn mat)))) + ;; Return the i,j entry of the Maxima matrix m. The rows of m have been permuted according ;; to the Maxima list p. @@ -275,8 +289,9 @@ (defun $identfor (mat &optional (fld-name '$generalring) (p 1) (q 1)) ($require_ring fld-name "$second" "$identfor") - (let ((add-id (funcall (mring-add-id (get fld-name 'ring)))) - (mul-id (funcall (mring-mult-id (get fld-name 'ring))))) + (let* ((fld (get fld-name 'ring)) + (add-id (funcall (mring-mring-to-maxima fld) (funcall (mring-add-id fld)))) + (mul-id (funcall (mring-mring-to-maxima fld) (funcall (mring-mult-id fld))))) (cond (($matrixp mat) ($require_square_matrix mat "$first" "$identfor") @@ -291,13 +306,17 @@ (defun $zerofor (mat &optional (fld-name '$generalring)) ($require_ring fld-name "$second" "$zerofor") - (let ((add-id (funcall (mring-add-id (get fld-name 'ring))))) + (let* ((fld (get fld-name 'ring)) + (add-id (funcall (mring-mring-to-maxima fld) (funcall (mring-add-id fld))))) (cond (($matrixp mat) - (let ((n ($first ($matrix_size mat))) (mc)) + (let ((r) (c) (mc)) + (setq r ($matrix_size mat)) + (setq c ($second r)) + (setq r ($first r)) (setq mc (copy-tree mat)) - (loop for i from 1 to n do - (loop for j from 1 to n do - (setf (nth j (nth i mc)) ($zerofor (nth j (nth i mat)))))) + (loop for i from 1 to r do + (loop for j from 1 to c do + (setf (nth j (nth i mc)) ($zerofor (nth j (nth i mat)) fld-name)))) mc)) (t add-id)))) Index: mring.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/mring.lisp,v retrieving revision 1.6 retrieving revision 1.7 diff -u -d -r1.6 -r1.7 --- mring.lisp 8 Jan 2006 04:09:01 -0000 1.6 +++ mring.lisp 8 Jan 2006 15:41:01 -0000 1.7 @@ -132,7 +132,12 @@ :mult #'* :sub #'- :negate #'- - :psqrt #'(lambda (s) nil) + :psqrt #'(lambda (s) (let ((x)) + (cond ((>= s 0) + (setq x (isqrt (numerator s))) + (setq x (/ x (isqrt (denominator s)))) + (if (= s (* x x)) x nil)) + (t nil)))) :add-id #'(lambda () 0) :mult-id #'(lambda () 1) :fzerop #'(lambda (s) (= s 0)) @@ -159,7 +164,7 @@ :mult #'mult :sub #'sub :negate #'(lambda (s) (mult -1 s)) - :psqrt #'(lambda (s) (if (member (csign s) `($pos $pz $zero)) (take '(%sqrt) s) nil)) + :psqrt #'(lambda (s) (if (member (csign ($ratdisrep s)) `($pos $pz $zero)) (take '(%sqrt) s) nil)) :add-id #'(lambda () 0) :mult-id #'(lambda () 1) :fzerop #'(lambda (s) (like s 0)) @@ -208,7 +213,7 @@ :mult #'(lambda (a b) ($rectform (mult a b))) :sub #'(lambda (a b) ($rectform (sub a b))) :negate #'(lambda (a) (mult -1 a)) - :psqrt #'(lambda (s) (mlsp s 0) nil (take '(%sqrt) s)) + :psqrt #'(lambda (s) (if (mlsp s 0) nil (take '(%sqrt) s))) :add-id #'(lambda () 0) :mult-id #'(lambda () 1) :fzerop #'(lambda (s) (like s bigfloatzero)) @@ -264,11 +269,11 @@ :add #'fp+ :div #'fp/ :rdiv #'fp/ - :reciprocal #'(lambda (s) (div (list 1 0) s)) + :reciprocal #'(lambda (s) (fp/ (list 1 0) s)) :mult #'fp* :sub #'fp- :negate #'(lambda (s) (list (- (first s)) (second s))) - :psqrt #'(lambda (s) (if (>= 0 (first s)) (list (cl:sqrt (first s)) (+ 1 (second s))) nil)) + :psqrt #'(lambda (s) (if (> (first s) 0) (list (cl:sqrt (first s)) (+ 1 (second s))) nil)) :add-id #'(lambda () (list 0 0)) :mult-id #'(lambda () (list 1 0)) :fzerop #'(lambda (s) (like (first s) 0)) |