[Axiom-developer] Re: Axiom + High Energy Physics

[root]

(this is a resend) .... i just clipped and tested what i sent out
and i managed to accidently delete the \usepackage line without
which the pamphlet stuff is broken. sorry about the flying fingers. 
this one is correct..... --t


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


> 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}
\usepackage{axiom}
\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}