diff --git a/books/bookvol10.5.pamphlet b/books/bookvol10.5.pamphlet index 58a0132..682a139 100644 --- a/books/bookvol10.5.pamphlet +++ b/books/bookvol10.5.pamphlet @@ -1650,37 +1650,37 @@ a:=[[3.+4.*%i,-4.+5.*%i,5.+6.*%i,7.-8.*%i,-9.-2.*%i]] --E 150 --S 151 of 155 -icamax(5,a,1) -- should be 4 +icamax(5,a,1) -- should be 3 --R ---R (3) 4 +--R (3) 3 --R Type: PositiveInteger --E 151 --S 152 of 155 -icamax(0,a,1) -- should be 1 +icamax(0,a,1) -- should be 0 --R ---R (4) 1 ---R Type: PositiveInteger +--R (4) 0 +--R Type: NonNegativeInteger --E 152 --S 153 of 155 -icamax(5,a,-1) -- should be 1 +icamax(5,a,-1) -- should be 0 --R ---R (5) 1 ---R Type: PositiveInteger +--R (5) 0 +--R Type: NonNegativeInteger --E 153 --S 154 of 155 -icamax(3,a,1) -- should be 3 +icamax(3,a,1) -- should be 2 --R ---R (6) 3 +--R (6) 2 --R Type: PositiveInteger --E 154 --S 155 of 155 -icamax(3,a,2) -- should be 2 +icamax(3,a,2) -- should be 1 --R ---R (7) 2 +--R (7) 1 --R Type: PositiveInteger --E 155 @@ -7706,6 +7706,9 @@ PRIMARR(COMPLEX(FLOAT)) is an array where each array element AE is \end{verbatim} so we need (caar AE) and (cadr AE) to get the real and imag parts. +NOTE: Array indexing in Axiom is 0-based but Fortran is 1-based so +the results are displaced by 1. + \begin{chunk}{BLAS 1 icamax} (defun icamaxSpad (n zx incx) ; Tim Daly May 13, 2012 @@ -7719,6 +7722,7 @@ so we need (caar AE) and (cadr AE) to get the real and imag parts. (icamax n vec incx))) (defun icamax (n cx incx) +; Tim Daly May 13, 2012 (declare (type (simple-array (complex single-float) (*)) cx) (type fixnum incx n)) (labels ( @@ -7727,13 +7731,13 @@ so we need (caar AE) and (cadr AE) to get the real and imag parts. (+ (the single-float (abs (the single-float (realpart zdum)))) (the single-float (abs (the single-float (imagpart zdum)))))))) (declare (ftype (function (complex single-float) single-float) cabs1)) - (let ((ix 0) (smax 0.0f0) (icamax 1) (limit (length cx))) + (let ((ix 0) (smax 0.0f0) (icamax 0) (limit (length cx))) (declare (type (single-float) smax) (type fixnum icamax ix limit)) (when (and (> n 1) (> incx 0)) (setf smax (the single-float (cabs1 (svref cx 0)))) (setf ix (the fixnum (+ ix incx))) - (do ((i 2 (+ i 1))) - ((or (>= i limit) (> i n))) + (do ((i 1 (+ i 1))) + ((or (>= ix limit) (>= i n))) (when (> (cabs1 (the (complex single-float) (svref cx ix))) smax) (setf icamax i) @@ -7741,6 +7745,7 @@ so we need (caar AE) and (cadr AE) to get the real and imag parts. (setf ix (the fixnum (+ ix incx))))) icamax))) + \end{chunk} \begin{chunk}{BLAS 1 icamax lisp test} @@ -7875,6 +7880,54 @@ c \end{chunk} +\begin{chunk}{idamaxEX example} + program idamaxEX +* Tim Daly May 15, 2012 +* unit tests for BLAS idamax + double precision a(5) + integer n + a = (/ 3.0, 4.0, -3.0, 5.0, -1.0 /) + write(6,100)a(1),a(2),a(3),a(4),a(5) + 100 format("a=(/ (/",f6.3," ",f6.3," ",f6.3," ",f6.3," ",f6.3," /)") + n=idamax(5,a,1) + write(6,200)n + 200 format("should be 4",/,"n=",i3) + n=idamax(3,a,1) + write(6,201)n + 201 format("should be 2",/,"n=",i3) + n=idamax(0,a,1) + write(6,202)n + 202 format("should be 0",/,"n=",i3) + n=idamax(-5,a,1) + write(6,203)n + 203 format("should be 0",/,"n=",i3) + n=idamax(5,a,-1) + write(6,204)n + 204 format("should be 0",/,"n=",i3) + n=idamax(5,a,2) + write(6,205)n + 205 format("should be 1",/,"n=",i3) + stop + end +\end{chunk} + +\begin{verbatim} +gcc -o idamaxEX idamaxEX.f -lgfortran idamax.o && ./idamaxEX +a=(/ (/ 3.000 4.000 -3.000 5.000 -1.000 /) +should be 4 +n= 4 +should be 2 +n= 2 +should be 0 +n= 0 +should be 0 +n= 0 +should be 0 +n= 0 +should be 1 +n= 1 +\end{verbatim} + \begin{chunk}{BLAS 1 idamax} (defun idamax (n dx incx) (declare (type (simple-array double-float (*)) dx) diff --git a/changelog b/changelog index 857a77d..92f8683 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,5 @@ +20120515 tpd src/axiom-website/patches.html 20120515.01.tpd.patch +20120515 tpd books/bookvol10.5 BLAS1 icamax 0-based fix 20120514 tpd src/axiom-website/patches.html 20120514.01.tpd.patch 20120514 tpd books/bookvol10.5 BLAS1 icamax 20120508 tpd src/axiom-website/patches.html 20120508.01.tpd.patch diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index eac4a1a..3f9a0f2 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -3914,5 +3914,7 @@ books/bookvol10.5 BLAS1 dzasum
books/bookvol10.5 BLAS1 dznrm2
20120514.01.tpd.patch books/bookvol10.5 BLAS1 icamax
+20120515.01.tpd.patch +books/bookvol10.5 BLAS1 icamax 0-based fix