diff --git a/books/bookvol10.5.pamphlet b/books/bookvol10.5.pamphlet
index ee2146c..99b5790 100644
--- a/books/bookvol10.5.pamphlet
+++ b/books/bookvol10.5.pamphlet
@@ -348,7 +348,7 @@ For complex symmetric matrices, TRANSx=H is not allowed.
)set message auto off
)clear all
---S 1 of 96
+--S 1 of 104
t1:Complex DoubleFloat := complex(1.0,0)
--R
--R
@@ -356,7 +356,7 @@ t1:Complex DoubleFloat := complex(1.0,0)
--R Type: Complex(DoubleFloat)
--E 1
---S 2 of 96
+--S 2 of 104
dcabs1(t1)
--R
--R
@@ -364,7 +364,7 @@ dcabs1(t1)
--R Type: DoubleFloat
--E 2
---S 3 of 96
+--S 3 of 104
t2:Complex DoubleFloat := complex(1.0,1.0)
--R
--R
@@ -372,7 +372,7 @@ t2:Complex DoubleFloat := complex(1.0,1.0)
--R Type: Complex(DoubleFloat)
--E 3
---S 4 of 96
+--S 4 of 104
dcabs1(t2)
--R
--R
@@ -380,7 +380,7 @@ dcabs1(t2)
--R Type: DoubleFloat
--E 4
---S 5 of 96
+--S 5 of 104
t3:Complex DoubleFloat := complex(1.0,-1.0)
--R
--R
@@ -388,7 +388,7 @@ t3:Complex DoubleFloat := complex(1.0,-1.0)
--R Type: Complex(DoubleFloat)
--E 5
---S 6 of 96
+--S 6 of 104
dcabs1(t3)
--R
--R
@@ -396,7 +396,7 @@ dcabs1(t3)
--R Type: DoubleFloat
--E 6
---S 7 of 96
+--S 7 of 104
t4:Complex DoubleFloat := complex(-1.0,-1.0)
--R
--R
@@ -404,7 +404,7 @@ t4:Complex DoubleFloat := complex(-1.0,-1.0)
--R Type: Complex(DoubleFloat)
--E 7
---S 8 of 96
+--S 8 of 104
dcabs1(t4)
--R
--R
@@ -412,7 +412,7 @@ dcabs1(t4)
--R Type: DoubleFloat
--E 8
---S 9 of 96
+--S 9 of 104
t5:Complex DoubleFloat := complex(-2.0,-2.0)
--R
--R
@@ -420,7 +420,7 @@ t5:Complex DoubleFloat := complex(-2.0,-2.0)
--R Type: Complex(DoubleFloat)
--E 9
---S 10 of 96
+--S 10 of 104
dcabs1(t5)
--R
--R
@@ -430,196 +430,196 @@ dcabs1(t5)
)clear all
---S 11 of 96
+--S 11 of 104
a:PRIMARR(DFLOAT):=[ [1.0,2.0,3.0,4,0,5,0,6,0] ]
--R
--R (1) [1.,2.,3.,4.,0.,5.,0.,6.,0.]
--R Type: PrimitiveArray(DoubleFloat)
--E 11
---S 12 of 96
+--S 12 of 104
dasum(3,a,-1) -- 0.0 neg incx
--R
--R (2) 0.
--R Type: DoubleFloat
--E 12
---S 13 of 96
+--S 13 of 104
dasum(3,a,0) -- 0.0 zero incx
--R
--R (3) 0.
--R Type: DoubleFloat
--E 13
---S 14 of 96
+--S 14 of 104
dasum(-1,a,1) -- 0.0 neg elements
--R
--R (4) 0.
--R Type: DoubleFloat
--E 14
---S 15 of 96
+--S 15 of 104
dasum(0,a,1) -- 0.0 no elements
--R
--R (5) 0.
--R Type: DoubleFloat
--E 15
---S 16 of 96
+--S 16 of 104
dasum(1,a,1) -- 1.0 1.0
--R
--R (6) 1.
--R Type: DoubleFloat
--E 16
---S 17 of 96
+--S 17 of 104
dasum(2,a,1) -- 3.0 1.0+2.0
--R
--R (7) 3.
--R Type: DoubleFloat
--E 17
---S 18 of 96
+--S 18 of 104
dasum(3,a,1) -- 6.0 1.0+2.0+3.0
--R
--R (8) 6.
--R Type: DoubleFloat
--E 18
---S 19 of 96
+--S 19 of 104
dasum(4,a,1) -- 10.0 1.0+2.0+3.0+4.0
--R
--R (9) 10.
--R Type: DoubleFloat
--E 19
---S 20 of 96
+--S 20 of 104
dasum(5,a,1) -- 15.0 1.0+2.0+3.0+4.0+5.0
--R
--R (10) 10.
--R Type: DoubleFloat
--E 20
---S 21 of 96
+--S 21 of 104
dasum(6,a,1) -- 21.0 1.0+2.0+3.0+4.0+5.0+6.0
--R
--R (11) 15.
--R Type: DoubleFloat
--E 21
---S 22 of 96
+--S 22 of 104
dasum(7,a,1) -- 21.0 1.0+2.0+3.0+4.0+5.0+6.0
--R
--R (12) 15.
--R Type: DoubleFloat
--E 22
---S 23 of 96
+--S 23 of 104
dasum(1,a,2) -- 1.0 1.0
--R
--R (13) 1.
--R Type: DoubleFloat
--E 23
---S 24 of 96
+--S 24 of 104
dasum(2,a,2) -- 4.0 1.0+3.0
--R
--R (14) 4.
--R Type: DoubleFloat
--E 24
---S 25 of 96
+--S 25 of 104
dasum(3,a,2) -- 9.0 1.0+3.0+5.0
--R
--R (15) 4.
--R Type: DoubleFloat
--E 25
---S 26 of 96
+--S 26 of 104
dasum(4,a,2) -- 9.0 1.0+3.0+5.0
--R
--R (16) 4.
--R Type: DoubleFloat
--E 26
---S 27 of 96
+--S 27 of 104
dasum(1,a,3) -- 1.0 1.0
--R
--R (17) 1.
--R Type: DoubleFloat
--E 27
---S 28 of 96
+--S 28 of 104
dasum(2,a,3) -- 5.0 1.0+4.0
--R
--R (18) 5.
--R Type: DoubleFloat
--E 28
---S 29 of 96
+--S 29 of 104
dasum(3,a,3) -- 5.0 1.0+4.0
--R
--R (19) 5.
--R Type: DoubleFloat
--E 29
---S 30 of 96
+--S 30 of 104
dasum(1,a,4) -- 1.0 1.0
--R
--R (20) 1.
--R Type: DoubleFloat
--E 30
---S 31 of 96
+--S 31 of 104
dasum(2,a,4) -- 6.0 1.0+5.0
--R
--R (21) 1.
--R Type: DoubleFloat
--E 31
---S 32 of 96
+--S 32 of 104
dasum(3,a,4) -- 6.0 1.0+5.0
--R
--R (22) 1.
--R Type: DoubleFloat
--E 32
---S 33 of 96
+--S 33 of 104
dasum(1,a,5) -- 1.0 1.0
--R
--R (23) 1.
--R Type: DoubleFloat
--E 33
---S 34 of 96
+--S 34 of 104
dasum(2,a,5) -- 7.0 1.0+6.0
--R
--R (24) 6.
--R Type: DoubleFloat
--E 34
---S 35 of 96
+--S 35 of 104
dasum(3,a,5) -- 7.0 1.0+6.0
--R
--R (25) 6.
--R Type: DoubleFloat
--E 35
---S 36 of 96
+--S 36 of 104
dasum(1,a,6) -- 1.0 1.0
--R
--R (26) 1.
--R Type: DoubleFloat
--E 36
---S 37 of 96
+--S 37 of 104
dasum(2,a,6) -- 1.0 1.0
--R
--R (27) 1.
--R Type: DoubleFloat
--E 37
---S 38 of 96
+--S 38 of 104
dasum(1,a,7) -- 1.0 1.0
--R
--R (28) 1.
@@ -628,7 +628,7 @@ dasum(1,a,7) -- 1.0 1.0
)clear all
---S 39 of 96
+--S 39 of 104
a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -636,7 +636,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 39
---S 40 of 96
+--S 40 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -644,7 +644,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 40
---S 41 of 96
+--S 41 of 104
daxpy(3,2.0,a,1,b,1)
--R
--R
@@ -652,7 +652,7 @@ daxpy(3,2.0,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 41
---S 42 of 96
+--S 42 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -660,7 +660,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 42
---S 43 of 96
+--S 43 of 104
daxpy(7,2.0,a,1,b,1)
--R
--R
@@ -668,7 +668,7 @@ daxpy(7,2.0,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 43
---S 44 of 96
+--S 44 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -681,7 +681,7 @@ Note that Axiom properly handles array indexes that are out of bounds.
The BLAS daxpy routine cannot check this condition.
\begin{chunk}{BlasLevelOne.input}
---S 45 of 96
+--S 45 of 104
daxpy(8,2.0,a,1,b,1)
--R
--R
@@ -689,7 +689,7 @@ daxpy(8,2.0,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 45
---S 46 of 96
+--S 46 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -697,7 +697,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 46
---S 47 of 96
+--S 47 of 104
daxpy(3,2.0,a,3,b,3)
--R
--R
@@ -705,7 +705,7 @@ daxpy(3,2.0,a,3,b,3)
--R Type: PrimitiveArray(DoubleFloat)
--E 47
---S 48 of 96
+--S 48 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -713,7 +713,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 48
---S 49 of 96
+--S 49 of 104
daxpy(4,2.0,a,2,b,2)
--R
--R
@@ -721,7 +721,7 @@ daxpy(4,2.0,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 49
---S 50 of 96
+--S 50 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -729,7 +729,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 50
---S 51 of 96
+--S 51 of 104
daxpy(5,2.0,a,2,b,2)
--R
--R
@@ -737,7 +737,7 @@ daxpy(5,2.0,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 51
---S 52 of 96
+--S 52 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -745,7 +745,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 52
---S 53 of 96
+--S 53 of 104
daxpy(3,2.0,a,2,b,2)
--R
--R
@@ -753,7 +753,7 @@ daxpy(3,2.0,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 53
---S 54 of 96
+--S 54 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -761,7 +761,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 54
---S 55 of 96
+--S 55 of 104
daxpy(3,-2.0,a,2,b,2)
--R
--R
@@ -769,7 +769,7 @@ daxpy(3,-2.0,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 55
---S 56 of 96
+--S 56 of 104
a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ]
--R
--R
@@ -777,7 +777,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 56
---S 57 of 96
+--S 57 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R
--R
@@ -785,7 +785,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 57
---S 58 of 96
+--S 58 of 104
daxpy(3,-2.0,a,1,b,2)
--R
--R
@@ -793,7 +793,7 @@ daxpy(3,-2.0,a,1,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 58
---S 59 of 96
+--S 59 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -801,7 +801,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 59
---S 60 of 96
+--S 60 of 104
daxpy(3,0.0,a,1,b,2)
--R
--R
@@ -811,7 +811,7 @@ daxpy(3,0.0,a,1,b,2)
)clear all
---S 61 of 96
+--S 61 of 104
a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R
--R
@@ -819,7 +819,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 61
---S 62 of 96
+--S 62 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -827,7 +827,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 62
---S 63 of 96
+--S 63 of 104
dcopy(3,a,1,b,1)
--R
--R
@@ -835,7 +835,7 @@ dcopy(3,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 63
---S 64 of 96
+--S 64 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -843,7 +843,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 64
---S 65 of 96
+--S 65 of 104
dcopy(7,a,1,b,1)
--R
--R
@@ -851,7 +851,7 @@ dcopy(7,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 65
---S 66 of 96
+--S 66 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -859,7 +859,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 66
---S 67 of 96
+--S 67 of 104
dcopy(8,a,1,b,1)
--R
--R
@@ -867,7 +867,7 @@ dcopy(8,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 67
---S 68 of 96
+--S 68 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -875,7 +875,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 68
---S 69 of 96
+--S 69 of 104
dcopy(3,a,3,b,3)
--R
--R
@@ -883,7 +883,7 @@ dcopy(3,a,3,b,3)
--R Type: PrimitiveArray(DoubleFloat)
--E 69
---S 70 of 96
+--S 70 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -891,7 +891,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 70
---S 71 of 96
+--S 71 of 104
dcopy(4,a,2,b,2)
--R
--R
@@ -899,7 +899,7 @@ dcopy(4,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 71
---S 72 of 96
+--S 72 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -907,7 +907,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 72
---S 73 of 96
+--S 73 of 104
dcopy(5,a,2,b,2)
--R
--R
@@ -915,7 +915,7 @@ dcopy(5,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 73
---S 74 of 96
+--S 74 of 104
b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R
--R
@@ -923,7 +923,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 74
---S 75 of 96
+--S 75 of 104
dcopy(3,a,2,b,2)
--R
--R
@@ -931,7 +931,7 @@ dcopy(3,a,2,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 75
---S 76 of 96
+--S 76 of 104
a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ]
--R
--R
@@ -939,7 +939,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 76
---S 77 of 96
+--S 77 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R
--R
@@ -947,7 +947,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 77
---S 78 of 96
+--S 78 of 104
dcopy(3,a,1,b,1)
--R
--R
@@ -955,7 +955,7 @@ dcopy(3,a,1,b,1)
--R Type: PrimitiveArray(DoubleFloat)
--E 78
---S 79 of 96
+--S 79 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R
--R
@@ -963,7 +963,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 79
---S 80 of 96
+--S 80 of 104
dcopy(3,a,1,b,2)
--R
--R
@@ -971,7 +971,7 @@ dcopy(3,a,1,b,2)
--R Type: PrimitiveArray(DoubleFloat)
--E 80
---S 81 of 96
+--S 81 of 104
a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R
--R
@@ -979,7 +979,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 81
---S 82 of 96
+--S 82 of 104
b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ]
--R
--R
@@ -987,7 +987,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ]
--R Type: PrimitiveArray(DoubleFloat)
--E 82
---S 83 of 96
+--S 83 of 104
dcopy(5,a,1,b,1)
--R
--R
@@ -997,63 +997,63 @@ dcopy(5,a,1,b,1)
)clear all
---S 84 of 96
+--S 84 of 104
a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ]
--R
--R (1) [1.,2.,3.,4.,5.]
--R Type: PrimitiveArray(DoubleFloat)
--E 84
---S 85 of 96
+--S 85 of 104
b:PRIMARR(DFLOAT):=[ [ 5.0, 6.0, 7.0, 8.0, 9.0] ]
--R
--R (2) [5.,6.,7.,8.,9.]
--R Type: PrimitiveArray(DoubleFloat)
--E 85
---S 86 of 96
+--S 86 of 104
ddot(0,a,1,b,1)
--R
--R (3) 0.
--R Type: DoubleFloat
--E 86
---S 87 of 96
+--S 87 of 104
ddot(3,a,1,b,1)
--R
--R (4) 38.
--R Type: DoubleFloat
--E 87
---S 88 of 96
+--S 88 of 104
ddot(3,a,1,b,2)
--R
--R (5) 46.
--R Type: DoubleFloat
--E 88
---S 89 of 96
+--S 89 of 104
ddot(3,a,2,b,1)
--R
--R (6) 58.
--R Type: DoubleFloat
--E 89
---S 90 of 96
+--S 90 of 104
ddot(3,a,1,b,-2)
--R
--R (7) 38.
--R Type: DoubleFloat
--E 90
---S 91 of 96
+--S 91 of 104
ddot(3,a,-2,b,1)
--R
--R (8) 50.
--R Type: DoubleFloat
--E 91
---S 92 of 96
+--S 92 of 104
ddot(3,a,-2,b,-2)
--R
--R (9) 71.
@@ -1062,34 +1062,125 @@ ddot(3,a,-2,b,-2)
)clear all
---S 93 of 96
+--S 93 of 104
a:PRIMARR(DFLOAT):=[ [ 3.0, -4.0, 5.0, -7.0, 9.0] ]
--R
--R (1) [3.,- 4.,5.,- 7.,9.]
--R Type: PrimitiveArray(DoubleFloat)
--E 93
---S 94 of 96
+--S 94 of 104
dnrm2(3,a,1)
--R
--R (2) 7.0710678118654755
--R Type: DoubleFloat
--E 94
---S 95 of 96
+--S 95 of 104
dnrm2(5,a,1)
--R
--R (3) 13.416407864998739
--R Type: DoubleFloat
--E 95
---S 96 of 96
+--S 96 of 104
dnrm2(3,a,2)
--R
---R (4) 10.72380529476361
+--R (4) 10.723805294763608
--R Type: DoubleFloat
--E 96
+)clear all
+--S 97 of 104
+a:MATRIX(DFLOAT):=[[6,5,0],[5,1,4],[0,4,3]]
+--R
+--R +6. 5. 0.+
+--R | |
+--R (1) |5. 1. 4.|
+--R | |
+--R +0. 4. 3.+
+--R Type: Matrix(DoubleFloat)
+--E 97
+
+--S 98 of 104
+t1:=drotg(elt(a,1,1),elt(a,1,2),0.0,0.0)
+--R
+--R (2)
+--R [7.810249675906654, 0.64018439966447993, 0.76822127959737585,
+--R 0.64018439966447993]
+--R Type: PrimitiveArray(DoubleFloat)
+--E 98
+
+--S 99 of 104
+g1:MATRIX(DFLOAT):=[[elt(t1,2), elt(t1,3),0.0],_
+ [-elt(t1,3),elt(t1,2),0.0],_
+ [0.0, 0.0, 1.0]]
+--R
+--R
+--R + 0.76822127959737585 0.64018439966447993 0.+
+--R | |
+--R (3) |- 0.64018439966447993 0.76822127959737585 0.|
+--R | |
+--R + 0. 0. 1.+
+--R Type: Matrix(DoubleFloat)
+--E 99
+
+--S 100 of 104
+t2:=g1*a
+--R
+--R + 7.810249675906654 4.4812907976513596 2.5607375986579197+
+--R | |
+--R (4) |- 4.4408920985006262E-16 - 2.4327007187250236 3.0728851183895034|
+--R | |
+--R + 0. 4. 3. +
+--R Type: Matrix(DoubleFloat)
+--E 100
+
+--S 101 of 104
+t3:=drotg(elt(t2,2,2),elt(a,3,2),0.0,0.0)
+--R
+--R (5)
+--R [4.6816698716254272, - 1.924474241977076, - 0.51962243930719854,
+--R 0.85439599751428896]
+--R Type: PrimitiveArray(DoubleFloat)
+--E 101
+
+--S 102 of 104
+g2:MATRIX(DFLOAT):=[[1.0, 0.0, 0.0],_
+ [0.0, elt(t3,2),elt(t3,3)],_
+ [0.0,-elt(t3,3),elt(t3,2)]]
+--R
+--R
+--R +1. 0. 0. +
+--R | |
+--R (6) |0. - 0.51962243930719854 0.85439599751428896 |
+--R | |
+--R +0. - 0.85439599751428896 - 0.51962243930719854+
+--R Type: Matrix(DoubleFloat)
+--E 102
+
+--S 103 of 104
+g2*g1*a
+--R
+--R + 7.810249675906654 4.4812907976513596 2.5607375986579197 +
+--R | |
+--R (7) | 2.2204460492503131E-16 4.6816698716254272 0.96644793161452336 |
+--R | |
+--R +- 4.4408920985006262E-16 0. - 4.1843280638948093+
+--R Type: Matrix(DoubleFloat)
+--E 103
+
+--S 104 of 104
+q:=transpose(g1)*transpose(g2)
+--R
+--R +0.76822127959737585 0.33265417936007158 0.54697098874441952 +
+--R | |
+--R (8) |0.64018439966447993 - 0.39918501523208583 - 0.65636518649330344|
+--R | |
+--R + 0. 0.85439599751428896 - 0.51962243930719854+
+--R Type: Matrix(DoubleFloat)
+--E 104
+
)spool
)lisp (bye)
\end{chunk}
@@ -1221,6 +1312,13 @@ BlasLevelOne() : Exports == Implementation where
++X dnrm2(5,a,1) -- 13.416407864998739 = sqrt(180.0)
++X dnrm2(3,a,2) -- 10.72380529476361 = sqrt(115.0)
+ drotg: (DF, DF, DF, DF) -> DX
+ ++ drotg computes a 2D plane Givens rotation spanned by two
+ ++ coordinate axes.
+ ++
+ ++X a:MATRIX(DFLOAT):=[[6,5,0],[5,1,4],[0,4,3]]
+ ++X drotg(elt(a,1,1),elt(a,1,2),0.0D0,0.0D0)
+
Implementation == add
dcabs1(z:CDF):DF ==
@@ -1235,6 +1333,8 @@ BlasLevelOne() : Exports == Implementation where
DDOT(n,dx,incx,dy,incy)$Lisp
dnrm2(n:SI,dx:DX,incx:SI):DF ==
DNRM2(n,dx,incx)$Lisp
+ drotg(a:DF,b:DF,c:DF,s:DF):DX ==
+ DROTG(a,b,c,s)$Lisp
\end{chunk}
\begin{chunk}{BLAS1.dotabb}
@@ -1441,6 +1541,16 @@ c NEW
end
\end{chunk}
+\begin{verbatim}
+gcc -o dcabs1EX dcabs1EX.f -lgfortran dcabs1.o && ./dcabs1EX
+
+ a=( 2.100, 2.100)
+ b=( 310.000, 4100.000)
+ a+b=( 312.100, 4102.100)
+dcabs1(c)=( 4414.200, 0.000)
+
+\end{verbatim}
+
\begin{chunk}{BLAS 1 dcabs1}
(declaim (ftype (function (cons) double-float) dcabs1))
(defun dcabs1 (z)
@@ -2196,6 +2306,41 @@ c
\end{chunk}
+\begin{verbatim}
+gcc -o dasumEX dasumEX.f -lgfortran dasum.o && ./dasumEX
+
+a(1)= 1.000 a(2)= 2.000 a(3)= 3.000
+a(4)= 4.000 a(5)= 5.000 a(6)= 6.000
+d= 0.000 should be 0.0, negative index
+d= 0.000 should be 0.0, zero increment
+d= 0.000 should be 0.0, negative elements
+d= 0.000 should be 0.0, no elements
+d= 1.000 should be 1.0
+d= 3.000 should be 3.0 = 1.0+2.0
+d= 6.000 should be 6.0 = 1.0+2.0+3.0
+d=10.000 should be 10.0 = 1.0+2.0+3.0+4.0
+d=15.000 should be 15.0 = 1.0+2.0+3.0+4.0+5.0
+d=21.000 should be 21.0 = 1.0+2.0+3.0+4.0+5.0+6.0
+d=21.000 should be 21.0 = 1.0+2.0+3.0+4.0+5.0+6.0
+d= 1.000 should be 1.0 = 1.0
+d= 4.000 should be 4.0 = 1.0+3.0
+d= 9.000 should be 9.0 = 1.0+3.0+5.0
+d= 9.000 should be 9.0 = 1.0+3.0+5.0
+d= 1.000 should be 1.0 = 1.0
+d= 5.000 should be 5.0 = 1.0+4.0
+d= 5.000 should be 5.0 = 1.0+4.0
+d= 1.000 should be 1.0 = 1.0
+d= 6.000 should be 6.0 = 1.0+5.0
+d= 6.000 should be 6.0 = 1.0+5.0
+d= 1.000 should be 1.0 = 1.0
+d= 7.000 should be 7.0 = 1.0+6.0
+d= 7.000 should be 7.0 = 1.0+6.0
+d= 1.000 should be 1.0 = 1.0
+d= 1.000 should be 1.0 = 1.0
+d= 1.000 should be 1.0 = 1.0
+
+\end{verbatim}
+
\begin{chunk}{BLAS 1 dasum}
(declaim (ftype (function (fixnum (simple-array double-float (*)) fixnum)
double-float) dasum))
@@ -2697,6 +2842,44 @@ gcc -o daxpyEx daxpyEX.f -lgfortran daxpy.o
end
\end{chunk}
+\begin{verbatim}
+gcc -o daxpyEX daxpyEX.f -lgfortran daxpy.o && ./daxpyEX
+
+a(1)= 1.000 a(2)= 2.000 a(3)= 3.000
+a(4)= 4.000 a(5)= 5.000 a(6)= 6.000
+b(1)= 1.000 b(2)= 2.000 b(3)= 3.000
+b(4)= 4.000 b(5)= 5.000 b(6)= 6.000
+
+t200 is (/ 3.0, 6.0, 9.0, 4.0, 5.0, 6.0, 7.0 /)
+b(1)= 3.000 b(2)= 6.000 b(3)= 9.000
+b(4)= 4.000 b(5)= 5.000 b(6)= 6.000 b(7)= 7.000
+
+t300 is (/ 3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0 /)
+b(1)= 3.000 b(2)= 6.000 b(3)= 9.000
+b(4)=12.000 b(5)=15.000 b(6)=18.000 b(7)=21.000
+
+t302 is (/ 3.0, 2.0, 3.0, 12.0, 5.0, 6.0, 21.0 /)
+b(1)= 3.000 b(2)= 2.000 b(3)= 3.000
+b(4)=12.000 b(5)= 5.000 b(6)= 6.000 b(7)=21.000
+
+t303 is (/ 3.0, 2.0, 9.0, 4.0, 15.0, 6.0, 21.0 /)
+b(1)= 3.000 b(2)= 2.000 b(3)= 9.000
+b(4)= 4.000 b(5)=15.000 b(6)= 6.000 b(7)=21.000
+
+t305 is (/ 3.0, 2.0, 9.0, 4.0, 15.0, 6.0, 7.0 /)
+b(1)= 3.000 b(2)= 2.000 b(3)= 9.000
+b(4)= 4.000 b(5)=15.000 b(6)= 6.000 b(7)= 7.000
+
+t306 is (/ -1.0, 2.0, -1.0, 4.0, -1.0, 6.0, 7.0 /)
+b(1)=-1.000 b(2)= 2.000 b(3)=-1.000
+b(4)= 4.000 b(5)=-1.000 b(6)= 6.000 b(7)= 7.000
+
+t307 is (/ 2.0, 2.0, 7.0, 4.0, 11.0, 6.0, 7.0 /)
+b(1)= 3.000 b(2)= 2.000 b(3)= 7.000
+b(4)= 4.000 b(5)=11.000 b(6)= 6.000 b(7)= 7.000
+
+\end{verbatim}
+
\begin{chunk}{BLAS 1 daxpy}
(declaim (ftype (function (fixnum double-float
(simple-array double-float (*)) fixnum
@@ -3186,6 +3369,40 @@ gcc -o dcopyEx dcopyEX.f -lgfortran dcopy.o
end
\end{chunk}
+\begin{verbatim}
+gcc -o dcopyEX dcopyEX.f -lgfortran dcopy.o && ./dcopyEX
+
+a(1)= 1.000 a(2)= 2.000 a(3)= 3.000
+a(4)= 4.000 a(5)= 5.000 a(6)= 6.000
+b(1)= 0.000 b(2)= 0.000 b(3)= 0.000
+b(4)= 0.000 b(5)= 0.000 b(6)= 0.000
+
+t200 is (/ 1.0 2.0 3.0 0.0 0.0 0.0 0.0 /)
+b(1)= 1.000 b(2)= 2.000 b(3)= 3.000
+b(4)= 0.000 b(5)= 0.000 b(6)= 0.000 b(7)= 0.000
+
+t300 is (/ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 /)
+b(1)= 1.000 b(2)= 2.000 b(3)= 3.000
+b(4)= 4.000 b(5)= 5.000 b(6)= 6.000 b(7)= 7.000
+
+t302 is (/ 1.0, 0.0, 0.0, 4.0, 0.0, 0.0, 7.0 /)
+b(1)= 1.000 b(2)= 0.000 b(3)= 0.000
+b(4)= 4.000 b(5)= 0.000 b(6)= 0.000 b(7)= 7.000
+
+t303 is (/ 1.0, 0.0, 3.0, 0.0, 5.0, 0.0, 7.0 /)
+b(1)= 1.000 b(2)= 0.000 b(3)= 3.000
+b(4)= 0.000 b(5)= 5.000 b(6)= 0.000 b(7)= 7.000
+
+t305 is (/ 1.0, 2.0, 3.0, 4.0, 5.0, 0.0, 0.0 /)
+b(1)= 1.000 b(2)= 2.000 b(3)= 3.000
+b(4)= 4.000 b(5)= 5.000 b(6)= 0.000 b(7)= 0.000
+
+t306 is (/ 1.0, 2.0, 2.0, 4.0, 3.0, 0.0, 0.0 /)
+b(1)= 1.000 b(2)= 2.000 b(3)= 2.000
+b(4)= 4.000 b(5)= 3.000 b(6)= 0.000 b(7)= 0.000
+
+\end{verbatim}
+
\begin{chunk}{BLAS 1 dcopy}
(declaim (ftype (function (fixnum (simple-array double-float (*)) fixnum
(simple-array double-float (*)) fixnum)
@@ -3502,6 +3719,32 @@ gcc -o ddotEx ddotEX.f -lgfortran ddot.o
end
\end{chunk}
+\begin{verbatim}
+a=(/ 1.000 2.000 3.000 4.000 5.000/)
+a=(/ 5.000 6.000 7.000 8.000 9.000/)
+
+t200 is 0.0
+c= 0.000
+
+t201 is 28.0
+c=38.000
+
+t202 is 46.0
+c=46.000
+
+t203 is 58.0
+c=58.000
+
+t204 is 38.0
+c=38.000
+
+t205 is 50.0
+c=50.000
+
+t206 is 71.0
+c=71.000
+\end{verbatim}
+
\begin{chunk}{BLAS 1 ddot}
(declaim (ftype (function (fixnum (simple-array double-float (*)) fixnum
(simple-array double-float (*)) fixnum)
@@ -3754,6 +3997,20 @@ gcc -o dnrm2Ex dnrm2EX.f -lgfortran dnrm2.o
end
\end{chunk}
+\begin{verbatim}
+gcc -o dnrm2EX dnrm2EX.f -lgfortran dnrm2.o && ./dnrm2EX
+a=(/ 3.000 -4.000 5.000/)
+
+t200 is sqrt(50.0)=7.071
+c= 7.071
+
+t201 is sqrt(180.0)=13.416
+c=13.416
+
+t202 is sqrt(115.0)=10.724
+c= 10.724
+\end{verbatim}
+
\begin{chunk}{BLAS 1 dnrm2}
(declaim (ftype (function (fixnum (simple-array double-float (*)) fixnum)
double-float) dnrm2))
@@ -3827,6 +4084,94 @@ gcc -o dnrm2Ex dnrm2EX.f -lgfortran dnrm2.o
drotg examples
====================================================================
+A Givens rotation is a rotation in the plane spanned by two
+coordinate axes, named after Wallace Givens. [REF-Wiki3]
+
+A Givens rotation is represented by a matrix of the form
+
+ +- -+
+ | 1 ... 0 ... 0 ... 0 |
+ | . . . . . . . |
+ | . . . . . . . |
+ | . . . . . . . |
+ | 0 ... c ... -s ... 0 |
+ | . . . . . . . |
+G(i,j,theta) = | . . . . . . . |
+ | . . . . . . . |
+ | 0 ... s ... c ... 0 |
+ | . . . . . . . |
+ | . . . . . . . |
+ | . . . . . . . |
+ | 0 ... 0 ... 0 ... 1 |
+ +- -+
+
+where c = cos(theta) and s = sin(theta) appear at the intersections
+of the ith, jth rows and columns. The non-zero elements of a Givens
+matrix are given by:
+
+g = 1 for k != i,j
+ kk
+g = c
+ ii
+g = c (sign of sine switches for j>i)
+ jj
+g = -s
+ ji
+g = s for i > j
+ ij
+
+The product G(i,j,theta)*x represents the counterclockwise rotation
+of the vector in the (i,j) plane of theta radians.
+
+The main use of Givens rotations in numerical linear algebra is to
+introduce zeros in vectors or matrices. This effect can be employed
+for computing the QR decomposition of a matrix. One advantage over
+the Householder tranformations is that they can easily be parallelized,
+and another is that often for very sparse matrices they have a lower
+operation count.
+
+When a Givens rotation matrix, G(i,k,theta) multiplies another
+matrix A from the left, GA, only rows i and j of A are affected.
+Thus we restrict attention to the following problem. Given a and b,
+find c=cos(theta) and s=sin(theta) such that
+
++- -+ +- -+ +- -+
+| c -s | | a | | r |
+| | | | = | |
+| s c | | b | | 0 |
++- -+ +- -+ +- -+
+
+The solution is
+
+r = sqrt(a^2+b^2)
+c = a/r
+s = -b/r
+
+However, the computation for r may overflow or underflow. An alternative
+formulation avoiding this problem [REF-GC96, p5.1.8] is implemented as the
+hypot function in many programming languages.
+
+As Anderson [REF-And00] discovered in improving LAPACK, a previously
+overlooked numerical consideration is continuity. To achieve this
+we require r to be positive.
+
+ if (b = 0) then { c = copysign(1,a); s=0; r=abs(a) }
+ else if (a = 0) then { c = 0; s = copysign(1,b); r = abs(b) }
+ else if (abs(b) > abs(a)) then
+ t = a/b
+ u = copysign(sqrt(1+t*t),b)
+ s = 1/u
+ c = s*t
+ r = b*u
+ else
+ t = b/a
+ u = copysign(sqrt(1+t*t),a)
+ c = 1/u
+ s = c*t
+ r = a*u
+
+copysign can be implemented as x*sgn(y) using the sign function.
+
====================================================================
Man Page Details
====================================================================
@@ -3860,7 +4205,7 @@ ARGUMENTS
a (input and output) DOUBLE PRECISION
First vector component. On input, the first component of the
- vector to be rotated. On output, a is overwritten by by r, the
+ vector to be rotated. On output, a is overwritten by r, the
first component of the vector in the rotated coordinate system
where:
@@ -3957,39 +4302,294 @@ c
\end{chunk}
+\begin{chunk}{drotg example}
+ program drotgEX
+* Tim Daly May 2, 2012
+* unit tests for BLAS drotg
+ double precision a11,a12,a13,a21,a22,a23,a31,a32,a33
+ double precision b11,b12,b13,b21,b22,b23,b31,b32,b33
+ double precision c11,c12,c13,c21,c22,c23,c31,c32,c33
+ double precision d11,d12,d13,d21,d22,d23,d31,d32,d33
+ double precision e11,e12,e13,e21,e22,e23,e31,e32,e33
+ double precision f11,f12,f13,f21,f22,f23,f31,f32,f33
+ double precision a,b,c,s
+ a11=6.0d0
+ a12=5.0d0
+ a13=0.0d0
+ a21=5.0d0
+ a22=1.0d0
+ a23=4.0d0
+ a31=0.0d0
+ a32=4.0d0
+ a33=3.0d0
+ write(6,10)
+ write(6,20)
+ write(6,30)a11,a12,a13
+ write(6,30)a21,a22,a23
+ write(6,30)a31,a32,a33
+ write(6,20)
+ 10 format("A=");
+ 20 format(" +- -+")
+ 30 format(" | ",f6.3," ",f6.3," ",f6.3," |")
+
+ a=a11
+ b=a21
+ c=0.0d0
+ s=0.0d0
+ write(6,100)a,b,c,s
+ 100 format(/,"a=",f6.3," b=",f6.3," c=",f6.3," s=",f6.3);
+
+ call drotg(a,b,c,s)
+ write(6,100)a,b,c,s
+
+ b11=c
+ b12=s
+ b13=0.0d0
+ b21=-s
+ b22=c
+ b23=0.0d0
+ b31=0.0d0
+ b32=0.0d0
+ b33=1.0d0
+ write(6,11)
+ 11 format(/,"G1=")
+ write(6,20)
+ write(6,30)b11,b12,b13
+ write(6,30)b21,b22,b23
+ write(6,30)b31,b32,b33
+ write(6,20)
+
+ c11=b11*a11+b12*a21+b13*a31
+ c21=b21*a11+b22*a21+b23*a31
+ c31=b31*a11+b32*a21+b33*a31
+ c12=b11*a12+b12*a22+b13*a32
+ c22=b21*a12+b22*a22+b23*a32
+ c32=b31*a12+b32*a22+b33*a32
+ c13=b11*a13+b12*a23+b13*a33
+ c23=b21*a13+b22*a23+b23*a33
+ c33=b31*a13+b32*a23+b33*a33
+ write(6,12)
+ 12 format(/,"G1*A=")
+ write(6,20)
+ write(6,30)c11,c12,c13
+ write(6,30)c21,c22,c23
+ write(6,30)c31,c32,c33
+ write(6,20)
+
+ a=c22
+ b=c32
+ c=0.0d0
+ s=0.0d0
+ write(6,100)a,b,c,s
+
+ call drotg(a,b,c,s)
+ write(6,100)a,b,c,s
+
+ d11=1.0d0
+ d12=0.0d0
+ d13=0.0d0
+ d21=0.0d0
+ d22=c
+ d23=s
+ d31=0.0d0
+ d32=-s
+ d33=c
+ write(6,13)
+ 13 format(/,"G2=")
+ write(6,20)
+ write(6,30)d11,d12,d13
+ write(6,30)d21,d22,d23
+ write(6,30)d31,d32,d33
+ write(6,20)
+
+ e11=d11*c11+d12*c21+d13*c31
+ e21=d21*c11+d22*c21+d23*c31
+ e31=d31*c11+d32*c21+d33*c31
+ e12=d11*c12+d12*c22+d13*c32
+ e22=d21*c12+d22*c22+d23*c32
+ e32=d31*c12+d32*c22+d33*c32
+ e13=d11*c13+d12*c23+d13*c33
+ e23=d21*c13+d22*c23+d23*c33
+ e33=d31*c13+d32*c23+d33*c33
+ write(6,14)
+ 14 format(/,"G2*G1*A=")
+ write(6,20)
+ write(6,30)e11,e12,e13
+ write(6,30)e21,e22,e23
+ write(6,30)e31,e32,e33
+ write(6,20)
+
+ f11=b31*d13 + b21*d12 + b11*d11
+ f12=b31*d23 + b21*d22 + b11*d21
+ f13=b31*d33 + b21*d32 + b11*d31
+ f21=b32*d13 + b22*d12 + b12*d11
+ f22=b32*d23 + b22*d22 + b12*d21
+ f23=b32*d33 + b22*d32 + b12*d31
+ f31=b33*d13 + b23*d12 + b13*d11
+ f32=b33*d23 + b23*d22 + b13*d21
+ f33=b33*d33 + b23*d32 + b13*d31
+ write(6,15)
+ 15 format(/,"Q=transpose(G1)*transpose(G2)")
+ write(6,20)
+ write(6,30)f11,f12,f13
+ write(6,30)f21,f22,f23
+ write(6,30)f31,f32,f33
+ write(6,20)
+
+ stop
+ end
+\end{chunk}
+
+\begin{verbatim}
+gcc -o drotgEX drotgEX.f -lgfortran drotg.o && ./drotgEX
+
+A=
+ +- -+
+ | 6.000 5.000 0.000 |
+ | 5.000 1.000 4.000 |
+ | 0.000 4.000 3.000 |
+ +- -+
+
+a= 6.000 b= 5.000 c= 0.000 s= 0.000
+
+a= 7.810 b= 0.640 c= 0.768 s= 0.640
+
+G1=
+ +- -+
+ | 0.768 0.640 0.000 |
+ | -0.640 0.768 0.000 |
+ | 0.000 0.000 1.000 |
+ +- -+
+
+G1*A=
+ +- -+
+ | 7.810 4.481 2.561 |
+ | -0.000 -2.433 3.073 |
+ | 0.000 4.000 3.000 |
+ +- -+
+
+a=-2.433 b= 4.000 c= 0.000 s= 0.000
+
+a= 4.682 b=-1.924 c=-0.520 s= 0.854
+
+G2=
+ +- -+
+ | 1.000 0.000 0.000 |
+ | 0.000 -0.520 0.854 |
+ | 0.000 -0.854 -0.520 |
+ +- -+
+
+G2*G1*A=
+ +- -+
+ | 7.810 4.481 2.561 |
+ | 0.000 4.682 0.966 |
+ | 0.000 0.000 -4.184 |
+ +- -+
+
+Q=transpose(G1)*transpose(G2)
+ +- -+
+ | 0.768 0.333 0.547 |
+ | 0.640 -0.399 -0.656 |
+ | 0.000 0.854 -0.520 |
+ +- -+
+
+\end{verbatim}
+
\begin{chunk}{BLAS 1 drotg}
(defun drotg (da db c s)
+ ; Tim Daly May 2, 2012
(declare (type (double-float) s c db da))
- (prog ((roe 0.0) (scale 0.0) (r 0.0) (z 0.0))
+ (let ((roe 0.0d0) (scale 0.0d0) (r 0.0d0) (z 0.0d0))
(declare (type (double-float) z r scale roe))
(setf roe db)
(when (> (the double-float (abs da)) (the double-float (abs db)))
(setf roe da))
(setf scale (+ (the double-float (abs da)) (the double-float (abs db))))
- (if (/= scale 0.0) (go label10))
- (setf c 1.0)
- (setf s 0.0)
- (setf r 0.0)
- (setf z 0.0)
- (go label20)
- label10
- (setf r
- (* scale (f2cl-lib:dsqrt (+ (expt (/ da scale) 2) (expt (/ db scale) 2)))))
- (setf r (* (f2cl-lib:dsign 1.0 roe) r))
- (setf c (/ da r))
- (setf s (/ db r))
- (setf z 1.0)
- (when (> (the double-float (abs da)) (the double-float (abs db)))
- (setf z s))
- (if (and (>= (the double-float (abs db)) (the double-float (abs da)))
- (/= c 0.0))
- (setf z (/ 1.0 c)))
- label20
- (setf da r)
- (setf db z)
- (return (values da db c s))))
+ (if (/= scale 0.0d0)
+ (progn
+ (setf r
+ (the double-float (*
+ scale
+ (the double-float (sqrt
+ (the double-float (+
+ (the double-float (*
+ (the double-float (/ da scale))
+ (the double-float (/ da scale))))
+ (the double-float (*
+ (the double-float (/ db scale))
+ (the double-float (/ db scale)))))))))))
+ (setf r (* (float-sign roe 1.0d0) r))
+ (setf c (the double-float (/ da r)))
+ (setf s (the double-float (/ db r)))
+ (setf z 1.0d0)
+ (when (> (the double-float (abs da)) (the double-float (abs db)))
+ (setf z s))
+ (if (and (>= (the double-float (abs db)) (the double-float (abs da)))
+ (/= c 0.0))
+ (setf z (the double-float (/ 1.0d0 c))))
+ (make-array 4 :initial-contents (list r z c s)))
+ (make-array 4 :initial-contents (list 0.0d0 0.0d0 1.0d0 0.0d0)))))
+
+\end{chunk}
+
+\begin{chunk}{BLAS 1 drotg lisp test}
+(load "drotg.lisp")
+(defun m (m i j) (svref (svref m (1- i)) (1- j)))
+(defun matprint (mat)
+ (format t "+- -+~%")
+ (format t "| ~6,3f ~6,3f ~6,3f |~%" (m mat 1 1) (m mat 1 2) (m mat 1 3))
+ (format t "| ~6,3f ~6,3f ~6,3f |~%" (m mat 2 1) (m mat 2 2) (m mat 2 3))
+ (format t "| ~6,3f ~6,3f ~6,3f |~%" (m mat 3 1) (m mat 3 2) (m mat 3 3))
+ (format t "+- -+~%"))
+(defun matmult (a b)
+ (vector
+ (vector (+ (* (m a 1 1) (m b 1 1))
+ (* (m a 1 2) (m b 2 1))
+ (* (m a 1 3) (m b 3 1)))
+ (+ (* (m a 1 1) (m b 1 2))
+ (* (m a 1 2) (m b 2 2))
+ (* (m a 1 3) (m b 3 2)))
+ (+ (* (m a 1 1) (m b 1 3))
+ (* (m a 1 2) (m b 2 3))
+ (* (m a 1 3) (m b 3 3))))
+ (vector (+ (* (m a 2 1) (m b 1 1))
+ (* (m a 2 2) (m b 2 1))
+ (* (m a 2 3) (m b 3 1)))
+ (+ (* (m a 2 1) (m b 1 2))
+ (* (m a 2 2) (m b 2 2))
+ (* (m a 2 3) (m b 3 2)))
+ (+ (* (m a 2 1) (m b 1 3))
+ (* (m a 2 2) (m b 2 3))
+ (* (m a 2 3) (m b 3 3))))
+ (vector (+ (* (m a 3 1) (m b 1 1))
+ (* (m a 3 2) (m b 2 1))
+ (* (m a 3 3) (m b 3 1)))
+ (+ (* (m a 3 1) (m b 1 2))
+ (* (m a 3 2) (m b 2 2))
+ (* (m a 3 3) (m b 3 2)))
+ (+ (* (m a 3 1) (m b 1 3))
+ (* (m a 3 2) (m b 2 3))
+ (* (m a 3 3) (m b 3 3))))))
+(defun transpose (mat)
+ (vector (vector (m mat 1 1) (m mat 1 2) (m mat 1 3))
+ (vector (m mat 2 1) (m mat 2 2) (m mat 2 3))
+ (vector (m mat 3 1) (m mat 3 2) (m mat 3 3))))
+(setq a #(#(6.0d0 5.0d0 0.0d0)
+ #(5.0d0 1.0d0 4.0d0)
+ #(0.0d0 4.0d0 3.0d0)))
+(multiple-value-setq (x y c s) (drotg (m a 1 1) (m a 2 1) 0.0d0 0.0d0))
+(list x y c s)
+(setq g1 (vector (vector c s 0.0d0)
+ (vector (- s) c 0.0d0)
+ (vector 0.0d0 0.0d0 1.0d0)))
+(matprint (setq g1a (matmult g1 a)))
+(multiple-value-setq (xx yy cc ss) (drotg (m g1a 2 2) (m g1a 3 2) 0.0d0 0.0d0))
+(list xx yy cc ss)
+(matprint (setq g2 (vector (vector 1.0d0 0.0d0 0.0d0) (vector 0.0d0 cc ss) (vector 0.0d0 (- ss) cc))))
+(matprint (setq g2 (matmult g1 a)))
\end{chunk}
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{drot BLAS}
%\pagehead{drot}{drot}
@@ -136771,12 +137371,12 @@ Warning: Types of argument 1 in call to ZLARFB do not match.
\getchunk{BLAS 1 dcopy}
\getchunk{BLAS 1 ddot}
\getchunk{BLAS 1 dnrm2}
+\getchunk{BLAS 1 drotg}
\end{chunk}
\begin{chunk}{untested}
\getchunk{BLAS lsame}
\getchunk{BLAS xerbla}
-\getchunk{BLAS 1 drotg}
\getchunk{BLAS 1 drot}
\getchunk{BLAS 1 dscal}
\getchunk{BLAS 1 dswap}
diff --git a/books/bookvolbib.bib b/books/bookvolbib.bib
index cb44193..36bbdd6 100644
--- a/books/bookvolbib.bib
+++ b/books/bookvolbib.bib
@@ -1,6 +1,12 @@
%% Created for Timothy Daly at 2012-03-10 06:07:15 -0500
%% Saved with string encoding Unicode (UTF-8)
+@misc{REF-Wiki3,
+ Date-Added = {2012-04-02 07:24:41 -0500},
+ Date-Modified = {2012-04-02 07:27:05 -0500},
+ Howpublished = {\verb|en.wikipedia.org/wiki/Givens_rotation|},
+ Title = {Givens rotation}}
+
@book{REF-Pea56,
Author = {T. Pearcey},
Date-Added = {2012-03-10 06:06:24 -0500},
@@ -229,6 +235,17 @@ real or complex matrix, with applications to condition estimation},
Volume = {139},
Year = {1969}}
+@techreport{REF-And00,
+ Author = {Edward Anderson},
+ Date-Added = {2012-04-02 05:23:40 -0500},
+ Date-Modified = {2012-04-02 05:24:46 -0500},
+ Institution = {University of Tennessee},
+ Number = {UT-CS-00-454},
+ Title = {Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem},
+ Type = {LAPACK Working Note 150},
+ Month = {Dec 4},
+ Year = {2000}}
+
@techreport{REF-Ris88,
Author = {Robert Risch},
Date-Added = {2012-03-10 05:23:40 -0500},
@@ -736,6 +753,15 @@ real or complex matrix, with applications to condition estimation},
Title = {Matrix Computations},
Year = {1989}}
+@book{REF-GC96,
+ Author = {Gene H. Golub, Charles F. Van Loan},
+ Date-Added = {2012-04-02 04:33:05 -0500},
+ Date-Modified = {2012-04-02 04:34:57 -0500},
+ Keywords = {ISBN 978-0-8018-5414-9},
+ Publisher = {Johns Hopkins University Press},
+ Title = {Matrix Computations},
+ Year = {1996}}
+
@article{REF-Flo63,
Author = {R. W. Floyd},
Date-Added = {2012-03-09 04:31:02 -0500},
diff --git a/books/bookvolbib.pamphlet b/books/bookvolbib.pamphlet
index 3d70300..a722655 100644
--- a/books/bookvolbib.pamphlet
+++ b/books/bookvolbib.pamphlet
@@ -1200,6 +1200,11 @@ ISBN 0-486-61272-4
\bibitem[Alt05]{Alt05}
Altmann, Simon L. Rotations, Quaternions, and Double Groups
Dover Publications, Inc. 2005 ISBN 0-486-44518-6
+\bibitem[And00]{And00}
+Anderson, Edward
+``Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem''
+LAPACK Working Note 150, University of Tennessee, UT-CS-00-454,
+December 4, 2000.
\bibitem[Ba10]{Ba10}
Baker, Martin ``3D World Simulation''
\verb|www.euclideanspace.com|
@@ -1271,6 +1276,10 @@ Drinfeld-Vladut bound'' Invent. Math., vol. 121, 1995, pp. 211--222.
Golub, Gene H. and Van Loan, Charles F.
``Matrix Computations''
Johns Hopkins University Press ISBN 0-8018-3772-3 (1989)
+\bibitem[GL96]{GL96}
+Golub, Gene H. and Van Loan, Charles F.
+``Matrix Computations''
+Johns Hopkins University Press ISBN 978-0-8018-5414-9 (1996)
\bibitem[Ha1896]{Ha1896}
Hathway, Arthur S., "A Primer Of Quaternions" (1896)
\bibitem[Ha95]{Ha95}
@@ -1496,5 +1505,7 @@ Wolfram Research, \verb|mathworld.wolfram.com/Quaternion.html|
\bibitem[Yu76]{Yu76}
D.Y.Y. Yun. ``On square-free decomposition algorithms''
{\sl Proceedings of SYMSAC'76} pages 26-35, 1976
+\bibitem[Wiki3]{Wiki3}
+\verb|en.wikipedia.org/wiki/Givens_rotation|
\end{thebibliography}
\end{document}
diff --git a/changelog b/changelog
index aa6d0c5..abf94b4 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,7 @@
+20120503 tpd src/axiom-website/patches.html 20120503.01.tpd.patch
+20120503 tpd books/bookvolbib.bib BLAS1 drotg references
+20120503 tpd books/bookvolbib BLAS1 drotg references
+20120503 tpd books/bookvol10.5 BLAS1 drotg
20120501 tpd src/axiom-website/patches.html 20120501.01.fmp.patch
20120501 fmp books/bookvol10.3 change HASHEQ to SXHASH
20120501 fmp src/interp/sys-pkg.lisp remove HASHEQ
diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html
index b444ae3..8dd1425 100644
--- a/src/axiom-website/patches.html
+++ b/src/axiom-website/patches.html
@@ -3900,5 +3900,7 @@ books/bookvol10.5 BLAS1 dnrm2
books/bookvol5 reset si::*system-directory* to the null string
20120501.01.fmp.patch
books/bookvol10.3 change HASHEQ to SXHASH
+20120503.01.tpd.patch
+books/bookvol10.5 BLAS1 drotg