> Tim, what's the plan for stuff like this? Is the Lisp level up to > providing the performance for these algorithms that the (presumably) > optimized code already implemented can provide, or at least close > enough so that no one will worry about it? Typically code for symbolic > algebra systems has been a different animal that that for numerically > optimized systems - is Axiom going to attempt to merge the two worlds?
well, i don't plan to rewrite EVERYTHING in lisp, despite what i told Bill, although in principle you can generate machine-instruction equivalent routines from lisp. Before we gave the system to NAG we worked with them to generate an Axiom-Fortran interface which still exists in the code. I'd like to keep it compatible with the NAG library but i don't have access to that anymore. NAG, IBM, and several hundred other places still use the original (or equivalent) fortran code although both NAG and IBMs are known to be much better than the regular BLAS code. I started putting together replacement routines (as pamphlets) for the numerics. I've attached an example from the BLAS library (aka BLAS DROTG, NAG F06AAF, IBM ESSL DROTG) for double-precision rotations in the Givens Plane. For reference see Lawson, C.L., Hanson, R.J., Kincaid, D.R., and Krogh, F.T. "Basic Linear Algebra Subprograms for Fortran Usage" ACM Trans. on Math. Soft. (TOMS) Vol 5, pp308-325 (1979) When I get around to it I'll write and inline the test case code so the testing can be automated from the makefile. Note that the default chunk is a makefile which is useful for standalone testing purposes. you can extract with: tangle -t8 drotg.pamphlet >Makefile.drotg and then extract the code, latex it, and (eventually) run test cases with make -f Makefile.drotg The Makefile.pamphlet in the directory will be used to construct the proper files to hook into Axiom. I don't yet have this fully implemented for all of BLAS but it's coming. I believe Camm is the BLAS maintainer but i'm not certain. t ========================================================================= \documentclass{article} \begin{document} \title{subroutine drotg(da,db,c,s)} \author{Jack Dongarra} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{Purpose} This generates a real Givens plane rotation with parameters $c$ and $s$ such that given a real $a$ and $b$: $$ \left( \begin{array}{cc} c & s\\ -s & c \end{array} \right) \left( \begin{array}{c} da\\ db \end{array} \right) = \left( \begin{array}{c} f\\ 0 \end{array} \right) $$ \noindent The routine computes $c$, $s$, and $f$ as follows: $$ f = \sigma{}\sqrt{da^2 + db^2} $$ $$ c=\left\{ \begin{array}{ccc} da/f & {\rm if\ }& f \ne 0\\ 1 & {\rm if\ }& f = 0 \end{array} \right. s=\left\{ \begin{array}{ccc} db/f & {\rm if\ }& f \ne 0\\ 0 & {\rm if\ }& f = 0 \end{array} \right. $$ $$ {\rm where\ \ \ \ } \sigma=\left\{ \begin{array}{ccc} {\rm sign\ }da & {\rm if\ }& \vert{}da\vert{} > \vert{}db\vert{}\\ {\rm sign\ }db & {\rm if\ }& \vert{}da\vert{} \le \vert{}db\vert{}\\ \end{array} \right. $$ The routine also computes the value of $z$ defined as $$ z=\left\{ \begin{array}{ccl} s & {\rm if\ }& \vert{}s\vert{} < c {\rm\ or\ }c = 0\\ 1/c & {\rm if\ }& 0 < \vert{}c\vert{} \le s \end{array} \right. $$ This enables $c$ and $s$ to be reconstructed from the single value $z$ as $$ c=\left\{ \begin{array}{ccc} \sqrt{1-z^2} & {\rm if\ }& \vert{}z\vert{} \le 1\\ 1/z & {\rm if\ }& \vert{}z\vert{} > 1 \end{array} \right. s=\left\{ \begin{array}{ccc} z & {\rm if\ }& \vert{}z\vert{} \le 1\\ \sqrt{1-c^2} & {\rm if\ }& \vert{}z\vert{} > 1 \end{array} \right. $$ To apply the plane rotation to a pair of real vectors call DROT. \section{Specification} \begin{verbatim} subroutine drotg(double precision da, double precision db, double precision c, double precision s) \end{verbatim} \section{Parameters} \begin{list}{} \item {\bf da} (Double Precision) \begin{list}{} \item {\sl Entry}: first element of the vector which determines the rotation \item {\sl Exit}: the value {f} \end{list} \item {\bf db} (Double Precision) \begin{list}{} \item {\sl Entry}: second element of the vector which determines the rotation \item {\sl Exit}: the value {z} \end{list} \item {\bf c} (Double Precision) \begin{list}{} \item {\sl Entry}: no value \item {\sl Exit}: cosine of the rotation \end{list} \item {\bf s} (Double Precision) \begin{list}{} \item {\sl Entry}: no value \item {\sl Exit}: sine of the rotation \end{list} \end{list} \section{Error Indicators and Warnings} None. \section{Examples} \subsection{Example 1: $f = 0$} \begin{verbatim} call drotg( 0.0, 0,0, c, s) ==> da = 0.0 db = 0.0 c = 1.0 s = 0.0 \end{verbatim} \subsection{Example 2: $c = 0$} \begin{verbatim} call drotg( 0.0, 2.0, c, s) ==> da = 2.0 db = 1.0 c = 0.0 s = 1.0 \end{verbatim} \subsection{Example 3: $\vert{}b\vert{} > \vert{}a\vert{}$} \begin{verbatim} call drotg( 6.0, -8.0, c, s) ==> da = -10.0 db = -1.666... c = -0.5 s = 0.8 \end{verbatim} \subsection{Example 4: $\vert{}a\vert{} > \vert{}b\vert{}$} \begin{verbatim} call drotg( 8.0, 6.0, c, s) ==> da = 10.0 db = 0.6 c = 0.8 s = 0.6 \end{verbatim} <<fortran>>= subroutine drotg(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c double precision da,db,c,s,roe,scale,r,z c roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .ne. 0.0d0 ) go to 10 c = 1.0d0 s = 0.0d0 r = 0.0d0 z = 0.0d0 go to 20 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r z = 1.0d0 if( dabs(da) .gt. dabs(db) ) z = s if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c 20 da = r db = z return end @ \section{Makefile} <<*>>= TANGLE=/usr/local/bin/NOTANGLE WEAVE=/usr/local/bin/NOWEAVE LATEX=/usr/bin/latex all: code doc run code: drotg.pamphlet ${TANGLE} -Rfortran drotg.pamphlet >drotg.f doc: ${WEAVE} -t8 -delay drotg.pamphlet >drotg.tex ${LATEX} drotg.tex 2>/dev/null 1>/dev/null ${LATEX} drotg.tex 2>/dev/null 1>/dev/null run: @echo done remake: ${TANGLE} -t8 drotg.pamphlet >Makefile.drotg @ \eject \begin{thebibliography}{99} \bibitem{1} Lawson, C.L., Hanson, R.J., Kincaid, D.R., and Krogh, F.T. ``Basic Linear Algebra Subprograms for Fortran Usage'' ACM Trans. on Math. Soft. (TOMS) Vol 5, pp308--325 (1979) \bibitem{2} Numerical Algorthms Fortran Library routine F06AAF\\ {\sl http://www.nag.co.uk/numeric/fl/manual/pdf/F06/f06aaf.pdf} \bibitem{3} IBM ESSL Documentation routine DROTG\\ {\sl http://csit1cwe.fsu.edu/extra\_link/essl/essl250.html\#HDRHSROTG} \end{thebibliography} \end{document} _______________________________________________ Axiom-developer mailing list Axiom-developer@nongnu.org http://lists.nongnu.org/mailman/listinfo/axiom-developer