diff --git a/books/bookvol10.5.pamphlet b/books/bookvol10.5.pamphlet index 6943118..fb0225b 100644 --- a/books/bookvol10.5.pamphlet +++ b/books/bookvol10.5.pamphlet @@ -507,26 +507,44 @@ BlasLevelOne() : Exports == Implementation where ++X dasum(3,dx,2) daxpy: (SI, DF, DX, SI, DX,SI) -> DX - ++ daxpy(n,da,dx,incx,dy,incy) computes a y = a*x + y + ++ daxpy(n,da,x,incx,y,incy) computes a y = a*x + y ++ for each of the chosen elements of the vectors x and y ++ and a constant multiplier a - ++ Note that the vector b is modified with the results. + ++ Note that the vector y is modified with the results. ++ - ++X dx:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] - ++X dy:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] - ++X daxpy(6,2.0,dx,1,dy,1) - ++X dy - ++X dm:PRIMARR(DFLOAT):=[[1.0,2.0,3.0]] - ++X dn:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] - ++X daxpy(3,-2.0,dm,1,dn,2) - ++X dn + ++X x:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] + ++X y:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] + ++X daxpy(6,2.0,x,1,y,1) + ++X y + ++X m:PRIMARR(DFLOAT):=[[1.0,2.0,3.0]] + ++X n:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] + ++X daxpy(3,-2.0,m,1,n,2) + ++X n + + dcopy: (SI, DX, SI, DX,SI) -> DX + ++ dcopy(n,x,incx,y,incy) copies y from x + ++ for each of the chosen elements of the vectors x and y + ++ Note that the vector y is modified with the results. + ++ + ++X x:PRIMARR(DFLOAT):=[[1.0,2.0,3.0,4.0,5.0,6.0]] + ++X y:PRIMARR(DFLOAT):=[[0.0,0.0,0.0,0.0,0.0,0.0]] + ++X dcopy(6,x,1,y,1) + ++X y + ++X m:PRIMARR(DFLOAT):=[[1.0,2.0,3.0]] + ++X n:PRIMARR(DFLOAT):=[[0.0,0.0,0.0,0.0,0.0,0.0]] + ++X dcopy(3,m,1,n,2) + ++X n Implementation == add - dcabs1(z:CDF):DF == DCABS1(z)$Lisp - dasum(n:SI,dx:DX,incx:SI):DF == DASUM(n,dx,incx)$Lisp + dcabs1(z:CDF):DF == + DCABS1(z)$Lisp + dasum(n:SI,dx:DX,incx:SI):DF == + DASUM(n,dx,incx)$Lisp daxpy(n:SI,da:DF,dx:DX,incx:SI,dy:DX,incy:SI):DX == DAXPY(n,da,dx,incx,dy,incy)$Lisp + dcopy(n:SI,dx:DX,incx:SI,dy:DX,incy:SI):DX == + DCOPY(n,dx,incx,dy,incy)$Lisp @ <>= @@ -1526,42 +1544,41 @@ RETURN VALUES <>= (defun daxpy (n da dx incx dy incy) - (declare (type (array double-float (*)) dy dx) - (type (double-float) da) (type fixnum incy incx n)) + (declare (type (simple-vector double-float) dx dy) + (type double-float da) (type fixnum incy incx n)) (let ((maxx (length dx)) (maxy (length dy))) (declare (type fixnum maxx maxy)) - (when (and (> n 0) (/= da 0.0) - (> incx 0) (< (* (1- n) incx) maxx) - (> incy 0) (< (* (1- n) incy) maxy)) - (if (and (= incx 1) (= incy 1)) - ; unit increments - (dotimes (i n) - (declare (type fixnum i)) - ; (format t "dy(~s)[~s] = dy(~s)[~s] + ( da[~s] * dx(~s)[~s] )~%" - ; i (+ (svref dy i) (* da (svref dx i))) - ; i (svref dy i) - ; da i (svref dx i)) - (setf (the double-float (svref dy i)) - (+ (the double-float (svref dy i)) - (* (the double-float da) - (the double-float (svref dx i)))))) - ; non-unit increments - (let ((ix 0) (iy 0)) - (declare (type fixnum i ix iy)) - (when (< incx 0) (setq ix (* (1+ (- n)) incx))) - (when (< incy 0) (setq ix (* (1+ (- n)) incy))) + (when (and (> n 0) (/= da 0.0) + (> incx 0) (< (* (1- n) incx) maxx) + (> incy 0) (< (* (1- n) incy) maxy)) + (if (and (= incx 1) (= incy 1)) + ; unit increments (dotimes (i n) ; (format t "dy(~s)[~s] = dy(~s)[~s] + ( da[~s] * dx(~s)[~s] )~%" - ; iy (+ (svref dy iy) (* da (svref dx ix))) - ; iy (svref dy iy) - ; da ix (svref dx ix)) - (setf (the double-float (svref dy iy)) - (+ (the double-float (svref dy iy)) + ; i (+ (svref dy i) (* da (svref dx i))) + ; i (svref dy i) + ; da i (svref dx i)) + (setf (the double-float (svref dy i)) + (+ (the double-float (svref dy i)) (* (the double-float da) - (the double-float (svref dx ix))))) - (setq ix (+ ix incx)) - (setq iy (+ iy incy))))))) - dy) + (the double-float (svref dx i)))))) + ; non-unit increments + (let ((ix 0) (iy 0)) + (declare (type fixnum ix iy)) + (when (< incx 0) (setq ix (* (1+ (- n)) incx))) + (when (< incy 0) (setq ix (* (1+ (- n)) incy))) + (dotimes (i n) + ; (format t "dy(~s)[~s] = dy(~s)[~s] + ( da[~s] * dx(~s)[~s] )~%" + ; iy (+ (svref dy iy) (* da (svref dx ix))) + ; iy (svref dy iy) + ; da ix (svref dx ix)) + (setf (the double-float (svref dy iy)) + (+ (the double-float (svref dy iy)) + (* (the double-float da) + (the double-float (svref dx ix))))) + (setq ix (+ ix incx)) + (setq iy (+ iy incy))))))) + dy) @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -1577,6 +1594,190 @@ RETURN VALUES )set message auto off )clear all +--S 1 of 23 +a:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]] +--R +--R +--R (1) [1.,2.,3.,4.,5.,6.,7.] +--R Type: PrimitiveArray DoubleFloat +--E 1 + +--S 2 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (2) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 2 + +--S 3 of 23 +dcopy(3,a,1,b,1) +--R +--R +--R (3) [1.,2.,3.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 3 + +--S 4 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (4) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 4 + +--S 5 of 23 +dcopy(7,a,1,b,1) +--R +--R +--R (5) [1.,2.,3.,4.,5.,6.,7.] +--R Type: PrimitiveArray DoubleFloat +--E 5 + +--S 6 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (6) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 6 + +--S 7 of 23 +dcopy(8,a,1,b,1) +--R +--R +--R (7) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 7 + +--S 8 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (8) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 8 + +--S 9 of 23 +dcopy(3,a,3,b,3) +--R +--R +--R (9) [1.,0.,0.,4.,0.,0.,7.] +--R Type: PrimitiveArray DoubleFloat +--E 9 + +--S 10 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (10) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 10 + +--S 11 of 23 +dcopy(4,a,2,b,2) +--R +--R +--R (11) [1.,0.,3.,0.,5.,0.,7.] +--R Type: PrimitiveArray DoubleFloat +--E 11 + +--S 12 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (12) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 12 + +--S 13 of 23 +dcopy(5,a,2,b,2) +--R +--R +--R (13) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 13 + +--S 14 of 23 +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] +--R +--R +--R (14) [0.,0.,0.,0.,0.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 14 + +--S 15 of 23 +dcopy(3,a,2,b,2) +--R +--R +--R (15) [1.,0.,3.,0.,5.,0.,0.] +--R Type: PrimitiveArray DoubleFloat +--E 15 + +--S 16 of 23 +a:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0]] +--R +--R +--R (16) [1.,2.,3.] +--R Type: PrimitiveArray DoubleFloat +--E 16 + +--S 17 of 23 +b:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0]] +--R +--R +--R (17) [1.,2.,3.,4.,5.] +--R Type: PrimitiveArray DoubleFloat +--E 17 + +--S 18 of 23 +dcopy(3,a,1,b,1) +--R +--R +--R (18) [1.,2.,3.,4.,5.] +--R Type: PrimitiveArray DoubleFloat +--E 18 + +--S 19 of 23 +b:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0]] +--R +--R +--R (19) [1.,2.,3.,4.,5.] +--R Type: PrimitiveArray DoubleFloat +--E 19 + +--S 20 of 23 +dcopy(3,a,1,b,2) +--R +--R +--R (20) [1.,2.,2.,4.,3.] +--R Type: PrimitiveArray DoubleFloat +--E 20 + +--S 21 of 23 +a:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0]] +--R +--R +--R (21) [1.,2.,3.,4.,5.] +--R Type: PrimitiveArray DoubleFloat +--E 21 + +--S 22 of 23 +b:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0]] +--R +--R +--R (22) [1.,2.,3.] +--R Type: PrimitiveArray DoubleFloat +--E 22 + +--S 23 of 23 +dcopy(5,a,1,b,1) +--R +--R +--R (23) [1.,2.,3.] +--R Type: PrimitiveArray DoubleFloat +--E 23 + )spool )lisp (bye) @ @@ -1585,6 +1786,45 @@ RETURN VALUES dcopy examples ==================================================================== +Assume we have two arrays which we initialize to these values: +a:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]] + +b:PRIMARR(DFLOAT):=[[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]] + +Note that after each call to dcopy the b array is modified. +The changed value is shown. We reset it after each bcopy but we +do not show that here. + +dcopy(3,a,1,b,1) ==> [1.,2.,3.,0.,0.,0.,0.] + +dcopy(7,a,1,b,1) ==> [1.,2.,3.,4.,5.,6.,7.] + +dcopy(8,a,1,b,1) ==> [0.,0.,0.,0.,0.,0.,0.] + +dcopy(3,a,3,b,3) ==> [1.,0.,0.,4.,0.,0.,7.] + +dcopy(4,a,2,b,2) ==> [1.,0.,3.,0.,5.,0.,7.] + +dcopy(5,a,2,b,2) ==> [0.,0.,0.,0.,0.,0.,0.] + +dcopy(3,a,2,b,2) ==> [1.,0.,3.,0.,5.,0.,0.] + +The arrays can be of different lengths: + +a:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0]] + +b:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0]] + +dcopy(3,a,1,b,1) ==> [1.,2.,3.,4.,5.] + +dcopy(3,a,1,b,2) ==> [1.,2.,2.,4.,3.] + +a:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0, 4.0, 5.0]] + +b:PRIMARR(DFLOAT):=[[ 1.0, 2.0, 3.0]] + +dcopy(5,a,1,b,1) ==> [1.,2.,3.] + ==================================================================== Man Page Details ==================================================================== @@ -1623,122 +1863,40 @@ ARGUMENTS incx INTEGER. (input) Increment between elements of x. - If incx = 0, the results will be unpredictable. + If incx = 0, the y vector is unchanged y DOUBLE PRECISION, (output) array of dimension (n-1) * |incy| + 1, result vector. incy INTEGER. (input) - Increment between elements of y. If incy = 0, the results will - be unpredictable. - -NOTES - When working backward (incx < 0 or incy < 0), each routine starts at - the end of the vector and moves backward, as follows: - - x(1-incx * (n-1)), x(1-incx * (n-2)), ..., x(1) - - y(1-incy * (n-1)), y(1-incy * (n-2)), ..., y(1) + Increment between elements of y. + If incy = 0, the y vector is unchanged @ <>= (defun dcopy (n dx incx dy incy) - (declare (type (array double-float (*)) dy dx) - (type fixnum incy incx n)) - (f2cl-lib:with-multi-array-data - ((dx double-float dx-%data% dx-%offset%) - (dy double-float dy-%data% dy-%offset%)) - (prog ((i 0) (ix 0) (iy 0) (m 0) (mp1 0)) - (declare (type fixnum mp1 m iy ix i)) - (if (<= n 0) (go end_label)) - (if (and (= incx 1) (= incy 1)) (go label20)) - (setf ix 1) - (setf iy 1) - (if (< incx 0) - (setf ix - (f2cl-lib:int-add - (f2cl-lib:int-mul (f2cl-lib:int-sub 1 n) incx) - 1))) - (if (< incy 0) - (setf iy - (f2cl-lib:int-add - (f2cl-lib:int-mul (f2cl-lib:int-sub 1 n) incy) - 1))) - (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) - ((> i n) nil) - (tagbody - (setf (f2cl-lib:fref dy-%data% (iy) ((1 *)) dy-%offset%) - (f2cl-lib:fref dx-%data% (ix) ((1 *)) dx-%offset%)) - (setf ix (f2cl-lib:int-add ix incx)) - (setf iy (f2cl-lib:int-add iy incy)))) - (go end_label) - label20 - (setf m (mod n 7)) - (if (= m 0) (go label40)) - (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) - ((> i m) nil) - (tagbody - (setf (f2cl-lib:fref dy-%data% (i) ((1 *)) dy-%offset%) - (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%)))) - (if (< n 7) (go end_label)) - label40 - (setf mp1 (f2cl-lib:int-add m 1)) - (f2cl-lib:fdo (i mp1 (f2cl-lib:int-add i 7)) - ((> i n) nil) - (tagbody - (setf (f2cl-lib:fref dy-%data% (i) ((1 *)) dy-%offset%) - (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%)) - (setf (f2cl-lib:fref dy-%data% - ((f2cl-lib:int-add i 1)) - ((1 *)) - dy-%offset%) - (f2cl-lib:fref dx-%data% - ((f2cl-lib:int-add i 1)) - ((1 *)) - dx-%offset%)) - (setf (f2cl-lib:fref dy-%data% - ((f2cl-lib:int-add i 2)) - ((1 *)) - dy-%offset%) - (f2cl-lib:fref dx-%data% - ((f2cl-lib:int-add i 2)) - ((1 *)) - dx-%offset%)) - (setf (f2cl-lib:fref dy-%data% - ((f2cl-lib:int-add i 3)) - ((1 *)) - dy-%offset%) - (f2cl-lib:fref dx-%data% - ((f2cl-lib:int-add i 3)) - ((1 *)) - dx-%offset%)) - (setf (f2cl-lib:fref dy-%data% - ((f2cl-lib:int-add i 4)) - ((1 *)) - dy-%offset%) - (f2cl-lib:fref dx-%data% - ((f2cl-lib:int-add i 4)) - ((1 *)) - dx-%offset%)) - (setf (f2cl-lib:fref dy-%data% - ((f2cl-lib:int-add i 5)) - ((1 *)) - dy-%offset%) - (f2cl-lib:fref dx-%data% - ((f2cl-lib:int-add i 5)) - ((1 *)) - dx-%offset%)) - (setf (f2cl-lib:fref dy-%data% - ((f2cl-lib:int-add i 6)) - ((1 *)) - dy-%offset%) - (f2cl-lib:fref dx-%data% - ((f2cl-lib:int-add i 6)) - ((1 *)) - dx-%offset%)))) - end_label - (return (values nil nil nil nil nil))))) + (declare (type (simple-vector double-float) dy dx) + (type fixnum incy incx n)) + (let ((maxx (length dx)) (maxy (length dy))) + (declare (type fixnum maxx maxy)) + (when (and (> n 0) + (> incx 0) (< (* (1- n) incx) maxx) + (> incy 0) (< (* (1- n) incy) maxy)) + (if (and (= incx 1) (= incy 1)) + ; unit increments + (dotimes (i n) + (setf (the double-float (svref dy i)) (the double-float (svref dx i)))) + ; non-unit increments + (let ((ix 0) (iy 0)) + (declare (type fixnum ix iy)) + (when (< incx 0) (setq ix (* (1+ (- n)) incx))) + (when (< incy 0) (setq ix (* (1+ (- n)) incy))) + (dotimes (i n) + (setf (the double-float (svref dy iy)) (the double-float (svref dx ix))) + (setq ix (+ ix incx)) + (setq iy (+ iy incy))))))) + dy) @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -75692,12 +75850,12 @@ ARGUMENTS <> <> <> +<> @ <>= <> <> -<> <> <> <> diff --git a/changelog b/changelog index b010f81..da0e041 100644 --- a/changelog +++ b/changelog @@ -1,6 +1,9 @@ +20100418 tpd src/axiom-website/patches.html 20100418.01.tpd.patch +20100418 tpd src/algebra/Makefile add BLAS1 dcopy help and regression test +20100418 tpd books/bookvol10.5 add BLAS1 dcopy 20100417 tpd src/axiom-website/patches.html 20100417.01.tpd.patch -20100417 tpd src/algebra/Makefile add BLAS daxpy help and regression test -20100417 tpd books/bookvol10.5 add BLAS daxpy +20100417 tpd src/algebra/Makefile add BLAS1 daxpy help and regression test +20100417 tpd books/bookvol10.5 add BLAS1 daxpy 20100407 tpd src/axiom-website/patches.html 20100407.01.tpd.patch 20100407 tpd src/algebra/Makefile add dasum regression and help files 20100407 tpd books/bookvol10.5 add BLAS1 dasum function diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet index 557b1dd..484a13e 100644 --- a/src/algebra/Makefile.pamphlet +++ b/src/algebra/Makefile.pamphlet @@ -16561,7 +16561,8 @@ SPADHELP=\ ${HELP}/ZeroDimensionalSolvePackage.help \ ${HELP}/dcabs1.help \ ${HELP}/dasum.help \ - ${HELP}/daxpy.help + ${HELP}/daxpy.help \ + ${HELP}/dcopy.help @ The algebra files contain input chunks in regress format. @@ -16698,7 +16699,8 @@ REGRESS= \ ZeroDimensionalSolvePackage.regress \ dcabs1.regress \ dasum.regress \ - daxpy.regress + daxpy.regress \ + dcopy.regress # these requires graphics # TwoDimensionalViewport @@ -18179,6 +18181,14 @@ ${HELP}/daxpy.help: ${BOOKS}/bookvol10.5.pamphlet >${INPUT}/daxpy.input @echo "daxpy" >>${HELPFILE} +${HELP}/dcopy.help: ${BOOKS}/bookvol10.5.pamphlet + @echo 9130 create dcopy.help from ${BOOKS}/bookvol10.5.pamphlet + @${TANGLE} -R"dcopy.help" ${BOOKS}/bookvol10.5.pamphlet \ + >${HELP}/dcopy.help + @${TANGLE} -R"dcopy.input" ${BOOKS}/bookvol10.5.pamphlet \ + >${INPUT}/dcopy.input + @echo "dcopy" >>${HELPFILE} + @ \section{The Makefile} diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index a29ef4c..b25a835 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -2626,6 +2626,8 @@ src/axiom-website/documentation.html fix typo
20100407.01.tpd.patch books/bookvol10.5 add BLAS1 dasum function
20100417.01.tpd.patch -books/bookvol10.5 add BLAS daxpy
+books/bookvol10.5 add BLAS1 daxpy
+20100418.01.tpd.patch +books/bookvol10.5 add BLAS1 dcopy