## maxima-commits

 ```Update of /cvsroot/maxima/maxima/share/linearalgebra In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv28340/share/linearalgebra Modified Files: eigens-by-jacobi.lisp linalg-utilities.lisp linalg.mac linalgcholesky.lisp lu.lisp matrixexp.lisp mring.lisp Added Files: test-eigens-by-jacobi.mac Log Message: (1) changing matrixexp to use \$require_square_matrix (did have its own require-square-matrix) (2) change order of arguments in full-matrix-map. (3) renamed my-matrix-map to matrix-map and moved to linalg-utilities. (4) changed infolevel values to \$debug, \$verbose, and \$silent. (5) moved the put for infolevel to linalg.mac. (6) moved function 'inform' to linalg-utilities. (7) changed argument order in require_symmetric_matrix. (8) array based eigens-by-jacobi. (9) changed add-id and mult-id for the bigfloatring. --- NEW FILE: test-eigens-by-jacobi.mac --- (listequalp(p,q) := block([listarith : true], listp(sort(p) - sort(q), lambda([x], is(abs(x) < 1.0d-10)))),0); 0\$ (matrixequalp(p,q) := matrixp(p-q, lambda([x], is(abs(x) < 1.0d-10))),0); 0\$ (m : matrix([1]),0); 0\$ eigens_by_jacobi(m,floatfield); [[1.0], matrix([1.0])]; eigens_by_jacobi(m,bigfloatfield); [[1.0b0], matrix([1.0b0])]; (m : matrix([0]),0); 0\$ eigens_by_jacobi(m,floatfield); [[0.0], matrix([1.0])]; eigens_by_jacobi(m,bigfloatfield); [[0.0b0], matrix([1.0b0])]; (m : matrix([%i]),0); 0\$ errcatch(eigens_by_jacobi(m,floatfield)); []\$ errcatch(eigens_by_jacobi(m,bigfloatfield)); []\$ (m : matrix([0,0],[0,0]),0); 0\$ eigens_by_jacobi(m,floatfield); [[0.0,0.0], matrix([1.0,0.0],[0.0,1.0])]\$ eigens_by_jacobi(m,bigfloatfield); [[0.0b0,0.0b0], matrix([1.0b0,0.0b0],[0.0b0,1.0b0])]\$ /*-------*/ (m : matrix([1,2],[2,1]),0); 0\$ (e : eigens_by_jacobi(m),0); 0\$ listequalp(first(e),[-1,3]); true\$ matrixequalp(transpose(second(e)) . m . second(e), apply('diag_matrix, first(e))); true\$ /*(e : eigens_by_jacobi(m,bigfloatfield),0); 0\$ listequalp(first(e),[-1,3]); true\$ matrixequalp(transpose(second(e)) . m . second(e), apply('diag_matrix, first(e))); true\$ */ /*-------*/ (m : matrix([14,32,50],[32,77,122],[50,122,194]),0); 0\$ (e : eigens_by_jacobi(m),0); 0\$ listequalp(first(e), map('rhs, allroots(charpoly(m,z)))); true\$ matrixequalp(transpose(second(e)) . m . second(e), apply('diag_matrix, first(e))); true\$ (remvalue(m,e),0); 0\$ (remfunction(matrixequalp, listequalp),0); 0\$ Index: eigens-by-jacobi.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/eigens-by-jacobi.lisp,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- eigens-by-jacobi.lisp 22 Jan 2006 04:29:10 -0000 1.1 +++ eigens-by-jacobi.lisp 28 Jan 2006 03:10:14 -0000 1.2 @@ -8,12 +8,7 @@ ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. (eval-when (:compile-toplevel :load-toplevel :execute) - (\$put '\$infolevel '\$debug '|\$linalg|) - (\$put '\$eigensbyjacobi 1 '\$version)) ;; Let's have version numbers 1,2,3,... - -(defun inform (level pck msg &rest arg) - (if (member level (member (\$get '\$infolevel pck) `(\$debug \$obnoxious \$chatty \$silent))) - (apply 'mtell `(,msg ,@arg)))) + (\$put '\$eigensbyjacobi 1 '\$version)) ;; Let's have version numbers 1,2,3,... ;; One sweep zeros each member of the matrix; for a n x n matrix, this requires n(n-1)/2 ;; Jacobi rotations. @@ -34,11 +29,11 @@ (merror "The field must either be 'floatfield' or 'bigfloatfield'")) (setq mm (mfuncall '\$mat_fullunblocker mm)) - (\$require_symmetric_matrix mm "\$first" "\$eigens_by_jacobi") + (\$require_symmetric_matrix mm "\$first" "\$eigens_by_jacobi") - (let* ((mat (copy-tree mm)) (g) (h) (sweeps 0) (rotations 0) (eps) (change) - (theta) (mpq) (c) (s) (tee) (tau) (d) (v (\$identfor mat)) (x) - (n (\$first (\$matrix_size mat))) (continue (> n 1)) + (let* ((mat) (g) (h) (sweeps 0) (rotations 0) (eps) (change) + (theta) (mpq) (c) (s) (tee) (tau) (d) (v) (x) (row) + (n (\$first (\$matrix_size mm))) (continue (> n 1)) (fld (\$require_ring fld-name "\$second" "\$eigens_by_jacobi")) (one (funcall (mring-mult-id fld))) (zero (funcall (mring-add-id fld)))) @@ -55,80 +50,92 @@ (fgreat (a b) (funcall (mring-great fld) a b)) (fmax (a b) (if (funcall (mring-great fld) a b) a b)) (fconvert (a) (funcall (mring-maxima-to-mring fld) a))) - - (matrix-map mat n n #'fconvert) + + (setq mat (make-array (list n n) :initial-contents (mapcar #'rest + (margs (matrix-map #'fconvert mm))))) + (setq v (make-array (list n n) :initial-element zero)) + (setq d (make-array n)) + (setq eps (if (eq fld-name '\$floatfield) double-float-epsilon (\$bfloat (div 1 (power 2 fpprec))))) - (loop for p from 1 to n do (push (array-elem mat p p) d)) - (setq d (reverse d)) - (push '(mlist) d) - + (decf n) + (loop for i from 0 to n do + (setf (aref v i i) one) + (setf (aref d i) (aref mat i i))) + (while continue (if (> sweeps 50) (merror "Exceeded maximum allowable number of Jacobi sweeps")) (incf sweeps) - - (setq change zero) - (loop for p from 1 to n do + (loop for p from 0 to n do (loop for q from (+ p 1) to n do - (setq mpq (array-elem mat p q)) + (setq mpq (aref mat p q)) (cond ((not (fzerop mpq)) (incf rotations) - (setq theta (fdiv (fsub (array-elem mat q q) (array-elem mat p p))(fmult 2 mpq))) + (setq theta (fdiv (fsub (aref mat q q) (aref mat p p))(fmult 2 mpq))) (setq tee (fdiv one (fadd (fabs theta) (fpsqrt (fadd one (fmult theta theta)))))) (if (fgreat 0 theta) (setq tee (fnegate tee))) (setq c (fdiv one (fpsqrt (fadd one (fmult tee tee))))) (setq s (fmult tee c)) (setq tau (fdiv s (fadd one c))) - (setmatelem mat zero p q) + (setf (aref mat p q) zero) - (loop for k from 1 to (- p 1) do - (setq g (array-elem mat k p)) - (setq h (array-elem mat k q)) - (setmatelem mat (fsub g (fmult s (fadd h (fmult g tau)))) k p) - (setmatelem mat (fadd h (fmult s (fsub g (fmult h tau)))) k q)) + (loop for k from 0 to (- p 1) do + (setq g (aref mat k p)) + (setq h (aref mat k q)) + (setf (aref mat k p) (fsub g (fmult s (fadd h (fmult g tau))))) + (setf (aref mat k q) (fadd h (fmult s (fsub g (fmult h tau)))))) (loop for k from (+ p 1) to (- q 1) do - (setq g (array-elem mat p k)) - (setq h (array-elem mat k q)) - (setmatelem mat (fsub g (fmult s (fadd h (fmult g tau)))) p k) - (setmatelem mat (fadd h (fmult s (fsub g (fmult h tau)))) k q)) + (setq g (aref mat p k)) + (setq h (aref mat k q)) + (setf (aref mat p k) (fsub g (fmult s (fadd h (fmult g tau))))) + (setf (aref mat k q) (fadd h (fmult s (fsub g (fmult h tau)))))) (loop for k from (+ q 1) to n do - (setq g (array-elem mat p k)) - (setq h (array-elem mat q k)) - (setmatelem mat (fsub g (fmult s (fadd h (fmult g tau)))) p k) - (setmatelem mat (fadd h (fmult s (fsub g (fmult h tau)))) q k)) + (setq g (aref mat p k)) + (setq h (aref mat q k)) + (setf (aref mat p k) (fsub g (fmult s (fadd h (fmult g tau))))) + (setf (aref mat q k) (fadd h (fmult s (fsub g (fmult h tau)))))) - (setmatelem mat (fsub (array-elem mat p p) (fmult tee mpq)) p p) - (setmatelem mat (fadd (array-elem mat q q) (fmult tee mpq)) q q) - - (loop for k from 1 to n do - (setq g (array-elem v k p)) - (setq h (array-elem v k q)) - (setmatelem v (fsub g (fmult s (fadd h (fmult g tau)))) k p) - (setmatelem v (fadd h (fmult s (fsub g (fmult h tau)))) k q)))))) + (setf (aref mat p p) (fsub (aref mat p p) (fmult tee mpq))) + (setf (aref mat q q) (fadd (aref mat q q) (fmult tee mpq))) + (loop for k from 0 to n do + (setq g (aref v k p)) + (setq h (aref v k q)) + (setf (aref v k p) (fsub g (fmult s (fadd h (fmult g tau))))) + (setf (aref v k q)(fadd h (fmult s (fsub g (fmult h tau)))))))))) (setq change zero) - (loop for i from 1 to n do - (setq x (array-elem mat i i)) - (setq change (fmax change (if (fgreat x eps) (fabs (fdiv (fsub (nth i d) x) x)) zero))) - (setf (nth i d) x)) + (loop for i from 0 to n do + (setq x (aref mat i i)) + (setq change (fmax change (if (fgreat (fabs x) eps) (fabs (fdiv (fsub (aref d i) x) x)) zero))) + (setf (aref d i) x)) (inform '\$debug '|\$linalg| "The largest percent change was ~:M~%" change) (setq continue (fgreat change eps))) - (inform '\$chatty '|\$linalg| "number of sweeps: ~:M~%" sweeps) - (inform '\$chatty '|\$linalg| "number of rotations: ~:M~%" rotations) - `((mlist) ,d ,v)))) + (inform '\$verbose '|\$linalg| "number of sweeps: ~:M~%" sweeps) + (inform '\$verbose '|\$linalg| "number of rotations: ~:M~%" rotations) + + (setq mm nil) + (loop for i from 0 to n do + (setq row nil) + (loop for j from 0 to n do + (push (aref v i j) row)) + (setq row (reverse row)) + (push '(mlist) row) + (push row mm)) + (setq mm (reverse mm)) + (push '(\$matrix) mm) + (setq d `((mlist) ,@(coerce d 'list))) + `((mlist) ,d ,mm)))) - - - + Index: linalg-utilities.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/linalg-utilities.lisp,v retrieving revision 1.12 retrieving revision 1.13 diff -u -d -r1.12 -r1.13 --- linalg-utilities.lisp 26 Jan 2006 22:24:06 -0000 1.12 +++ linalg-utilities.lisp 28 Jan 2006 03:10:15 -0000 1.13 @@ -9,6 +9,10 @@ (\$put '\$linalgutilities 2 '\$version) +(defun inform (level pck msg &rest arg) + (if (member level (member (\$get '\$infolevel pck) `(\$debug \$verbose \$silent))) + (apply 'mtell `(,msg ,@arg)))) + (defun \$listp (e &optional (f nil)) (and (op-equalp e 'mlist) (or (eq f nil) (every #'(lambda (s) (eq t (mfuncall f s))) (margs e))))) @@ -36,7 +40,7 @@ (defun array-elem (m i j) (nth j (nth i m))) -(defun \$require_symmetric_matrix (m fun pos) +(defun \$require_symmetric_matrix (m pos fun) (if (not (\$matrixp m)) (merror "The ~:M argument to ~:M must be a matrix" pos fun)) (let ((n (\$matrix_size m))) (if (not (= (\$first n) (\$second n))) @@ -66,11 +70,17 @@ (if (not (and (integerp i) (> i 0))) (merror "The ~:M argument of the function ~:M must be a positive integer" pos fun))) +(defun matrix-map (f mat) + (setq mat (mapcar 'cdr (cdr mat))) + (setq mat (mapcar #'(lambda (s) (mapcar f s)) mat)) + (setq mat (mapcar #'(lambda (s) (cons `(mlist simp) s)) mat)) + `((\$matrix simp) ,@mat)) + ;; Map the lisp function fn over the matrix m. This function is block matrix friendly. -(defun full-matrix-map (m fn) +(defun full-matrix-map (fn m) (if (or (\$listp m) (\$matrixp m)) - (cons (car m) (mapcar #'(lambda (s) (full-matrix-map s fn)) (cdr m))) + (cons (car m) (mapcar #'(lambda (s) (full-matrix-map fn s)) (cdr m))) (funcall fn m))) (defun \$zerofor (mat &optional (fld-name '\$generalring)) @@ -115,7 +125,7 @@ (push '(\$matrix) new-mat))) (defun \$ctranspose (m) - (mfuncall '\$transpose (full-matrix-map m #'(lambda (s) (simplifya `((\$conjugate) ,s) nil))))) + (mfuncall '\$transpose (full-matrix-map #'(lambda (s) (simplifya `((\$conjugate) ,s) nil)) m))) (defun \$zeromatrixp (m) (if (or (\$matrixp m) (\$listp m)) (every '\$zeromatrixp (cdr m)) Index: linalg.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/linalg.mac,v retrieving revision 1.21 retrieving revision 1.22 diff -u -d -r1.21 -r1.22 --- linalg.mac 26 Jan 2006 22:24:06 -0000 1.21 +++ linalg.mac 28 Jan 2006 03:10:15 -0000 1.22 @@ -24,13 +24,16 @@ eval_when([translate, compile], translate_fast_arrays : false); +/* eval_when(translate, declare_translated(ptriangularize_with_proviso,locate_matrix_entry,hankel,toeplitz, nullspace,require_integer,columnspace,rowswap,rowop,require_symbol, mat_fullunblocker,mat_trace,mat_unblocker, column_reduce,good_pivot, hipow_gzeromat_norm))\$ +*/ eval_when([batch,load,translate,compile], + put(infolevel, debug, linalg), load("polynomialp"), load("lu"), load("linalgcholesky"), Index: linalgcholesky.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/linalgcholesky.lisp,v retrieving revision 1.9 retrieving revision 1.10 diff -u -d -r1.9 -r1.10 --- linalgcholesky.lisp 26 Jan 2006 03:12:04 -0000 1.9 +++ linalgcholesky.lisp 28 Jan 2006 03:10:15 -0000 1.10 @@ -8,8 +8,7 @@ ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. (\$put '\$cholesky 1 '\$version) - -(defun \$cholesky (m &optional (fld-name '\$generalring)) + (defun \$cholesky (m &optional (fld-name '\$generalring)) (\$require_nonempty_matrix m "\$first" "\$cholesky") (\$require_selfadjoint_matrix m "\$first" "\$cholesky") @@ -61,7 +60,7 @@ (loop for i from 0 to n do (setq row nil) (loop for j from 0 to n do - (push (full-matrix-map (aref l i j) #'frevert) row)) + (push (full-matrix-map #'frevert (aref l i j)) row)) (setq row (reverse row)) (push '(mlist) row) (push row mat)) Index: lu.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/lu.lisp,v retrieving revision 1.10 retrieving revision 1.11 diff -u -d -r1.10 -r1.11 --- lu.lisp 26 Jan 2006 22:24:06 -0000 1.10 +++ lu.lisp 28 Jan 2006 03:10:15 -0000 1.11 @@ -23,14 +23,6 @@ (setf (nth i p) (nth j p)) (setf (nth j p) x))) -;; Map the lisp function fn over the r by c Maxima matrix m. This function isn't -;; block matrix friendly. - -(defun matrix-map (m r c fn) - (loop for i from 1 to r do - (loop for j from 1 to c do - (setmatelem m (funcall fn (nth j (nth i m))) i j)))) - ;; Return the i,j entry of the Maxima matrix m. The rows of m have been permuted according ;; to the Maxima list p. @@ -63,7 +55,7 @@ (setq fld (\$require_ring fld "\$second" "\$lu_factor")) (let ((m (copy-tree mm)) (c (\$length mm)) (perm) (cnd) (fn)) - (matrix-map m c c (mring-maxima-to-mring fld)) + (setq m (matrix-map (mring-maxima-to-mring fld) m)) (loop for i from 1 to c do (push i perm)) (setq perm (reverse perm)) @@ -134,10 +126,10 @@ (multiple-value-setq (lb ub) (mat-cond-by-lu m perm c (mring-coerce-to-lisp-float fld))) (setq lb (\$limit (mul lb cnd))) (setq ub (\$limit (mul ub cnd))) - (matrix-map m c c (mring-mring-to-maxima fld)) + (setq m (matrix-map (mring-mring-to-maxima fld) m)) `((mlist) ,m ,perm ,(mring-name fld) ,lb ,ub)) (t - (matrix-map m c c (mring-mring-to-maxima fld)) + (setq m (matrix-map (mring-mring-to-maxima fld) m)) `((mlist) ,m ,perm ,(mring-name fld)))))) ;; The first argument should be a matrix in packed LU form. The Maxima list perm @@ -203,13 +195,13 @@ (setq r (\$first n)) (setq c (\$second n)) - (matrix-map mat r c (mring-maxima-to-mring fld)) + (setq mat (matrix-map (mring-maxima-to-mring fld) mat)) ;(displa `((mequal) mat ,mat)) (setq b (copy-tree b1)) (setq c (\$second (\$matrix_size mat))) (setq cc (\$second (\$matrix_size b))) - (matrix-map b r cc (mring-maxima-to-mring fld)) + (setq b (matrix-map (mring-maxima-to-mring fld) b)) (setq bb (copy-tree b)) (loop for i from 1 to r do @@ -232,7 +224,7 @@ (setq acc (funcall fsub acc (funcall fmult (m-elem mat perm i j) (m-elem bb id-perm j q))))) (setmatelem bb (funcall fdiv acc (m-elem mat perm i i)) i q))) - (matrix-map bb r cc (mring-mring-to-maxima fld)) + (setq bb (matrix-map (mring-mring-to-maxima fld) bb)) bb)) (defun \$invert_by_lu (m &optional (fld '\$generalring)) Index: matrixexp.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/matrixexp.lisp,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- matrixexp.lisp 29 Aug 2005 14:03:53 -0000 1.1 +++ matrixexp.lisp 28 Jan 2006 03:10:15 -0000 1.2 @@ -28,27 +28,13 @@ (\$put '\$matrixexp 1 '\$version) -(defun my-matrix-map (f mat) - (setq mat (mapcar 'cdr (cdr mat))) - (setq mat (mapcar #'(lambda (s) (mapcar f s)) mat)) - (setq mat (mapcar #'(lambda (s) (cons `(mlist simp) s)) mat)) - `((\$matrix simp) ,@mat)) - -;; Signal an error if 'mat' is not a square matrix; otherwise, return true. - -(defun require-square-matrix (mat pos fun-name) - (if (not (\$matrixp mat)) - (merror "The ~:M argument to ~:M must be a matrix" pos fun-name)) - (if (not (= (\$length (\$args mat)) (\$length (\$first (\$args mat))))) - (merror "The ~:M argument to ~:M must be a square matrix" pos fun-name))) - ;; When mat is a square matrix, return exp(mat * x). The second ;; argument is optional and it defaults to 1. (defun \$matrixexp (mat &optional (x 1)) (let ((sp) (d) (p) (id) (n (\$length (\$args mat))) (f)) (\$ratvars) - (require-square-matrix mat "\$first" "\$matrixexp") + (\$require_square_matrix mat "\$first" "\$matrixexp") (setq mat (\$spectral_rep mat)) (setq sp (\$first mat)) (setq p (\$second mat)) @@ -77,7 +63,7 @@ (let ((z (gensym)) (expr) (var) (sp) (d) (p) (di) (n (\$length (\$args mat))) (f 0)) - (require-square-matrix mat "\$second" "\$matrixexp") + (\$require_square_matrix mat "\$second" "\$matrixexp") (setq expr (require-lambda lamexpr 1 "\$first" "\$matrixfun")) (setq var (nth 0 (nth 0 expr))) (setq expr (nth 1 expr)) @@ -122,7 +108,7 @@ (defun \$spectral_rep (mat) - (require-square-matrix mat "\$first" "\$spectral_rep") + (\$require_square_matrix mat "\$first" "\$spectral_rep") (let ((\$gcd '\$spmod) (\$algebraic t) (\$resultant '\$subres) (ord) (zi) (\$ratfac nil) (z (gensym)) (res) (m) (n (\$length (\$args mat))) (p) (p1) (p2) (sp) (proj)) @@ -158,7 +144,7 @@ (dotimes (i m) (setq zi (nth i sp)) (setq ord (nth (+ i 1) \$multiplicities)) - (push (my-matrix-map #'(lambda (e) (rational-residue e z zi p ord)) res) + (push (matrix-map #'(lambda (e) (rational-residue e z zi p ord)) res) proj)) (setq proj (nreverse proj)) Index: mring.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/share/linearalgebra/mring.lisp,v retrieving revision 1.9 retrieving revision 1.10 diff -u -d -r1.9 -r1.10 --- mring.lisp 11 Jan 2006 13:00:31 -0000 1.9 +++ mring.lisp 28 Jan 2006 03:10:15 -0000 1.10 @@ -214,8 +214,8 @@ :sub #'(lambda (a b) (\$rectform (sub a b))) :negate #'(lambda (a) (mult -1 a)) :psqrt #'(lambda (s) (if (mlsp s 0) nil (take '(%sqrt) s))) - :add-id #'(lambda () 0) - :mult-id #'(lambda () 1) + :add-id #'(lambda () bigfloatzero) + :mult-id #'(lambda () bigfloatone) :fzerop #'(lambda (s) (like s bigfloatzero)) :adjoint #'cl:identity :mring-to-maxima #'(lambda (s) s) ```