$SPAD/src/input testpackage.input

Mark Clements December 3, 2016

Abstract

Contents

2

— testpackage.spad —

)abbrev package TESTP TestPackage
VR ==> Vector(R)
MR ==> Matrix(R)
NNI ==> NonNegativeInteger
TestPackage(F:Join(Field,RadicalCategory,TranscendentalFunctionCategory),_
        R:ComplexCategory(F)): with
transposeH: MR -> MR
unitVector: NNI -> VR
norm: VR -> F
conjugate: VR -> VR
"/": (R,MR) -> MR
"/": (R,VR) -> VR
"/": (VR,R) -> VR
"**": (VR,NNI) -> VR
solveUpperTriangular: (MR,VR) -> VR
householder: VR -> MR
qr: MR -> Record(q:MR,r:MR)
lsqr: (MR,VR) -> VR
polyfit: (VR,VR,NNI) -> VR
== add
import VectorFunctions2(R,F)
transposeH(m:MR):MR == transpose(map(conjugate,m))
conjugate(v:VR):VR == map((vi:R):R +-> conjugate(vi)$R,v)$VR
unitVector(dim) ==
  out := new(dim,0@R)$VR
  out(1) := 1@R
  out
a:R / m:MR == map((mij:R):R +-> a/mij,m)$MR
a:R / v:VR == map((vi:R):R +-> a/vi,v)$VR
v:VR / a:R == map((vi:R):R +-> vi/a,v)$VR
v:VR ** n:NonNegativeInteger == map((vi:R):R +-> vi**n,v)$VR
norm(v:Vector(R)):F == reduce("+",map(norm,v))
solveUpperTriangular(r,b) ==
  n := ncols r
  x := new(n,0@R)$VR
  for k in n..1 by -1 repeat
    index := min(n,k+1)
    x(k) := (b(k)-reduce("+",subMatrix(r,k,k,index,n)*x.(index..n)))/r(k,k)
  x
householder(a) ==
  m := #a
  alpha : R := -a(1)*inv(sqrt norm a(1))*sqrt(norm(a))
  e : VR := unitVector(m)
  u : VR := a + alpha*e
  v : VR := u/complex(sqrt norm(u),0)
  w : R := dot(conjugate a,v)/dot(conjugate v,a)
  scalarMatrix(m,1) - (1+w)*transpose(outerProduct(v,conjugate v))

    qr(a) ==
      m := nrows a
      n := ncols a
      qm := scalarMatrix(m,1)
      rm : MR := copy a
      for i in 1..(if m=n then n-1 else n) repeat
        x := column(subMatrix(rm,i,m,i,n),1)
        h := scalarMatrix(m,1)
        setsubMatrix!(h,i,i,householder x)
        qm := qm*h
        rm := h*rm
      [qm,rm]

———-

— * —

)set break resume
)spool testpackage.output
)set message test on
)set message auto off
)clear all
)sys cp $AXIOM/../../src/input/testpackage.input.pamphlet .
)lisp (tangle "testpackage.input.pamphlet" "testpackage.spad" "testpackage.spad")
)co testpackage.spad
--S 1 of 5
m := matrix [[12, -51, 4+%i], [6, 167, -68], [-4, 24, -41]]
--R
--R        +12   - 51  4 + %i+
--R        |                 |
--R   (1)  | 6   167    - 68 |
--R        |                 |
--R        +- 4   24    - 41 +
--R                                               Type: Matrix(Complex(Integer))
--E 1
--S 2 of 5
out := qr m
--R
--R            + 6     69      58+    +                  6   +
--R            | -    ---   - ---|    |14   21    - 14 + - %i|
--R            | 7    175     175|    |                  7   |
--R            |                 |    |                      |
--R            | 3     158    6  |    |                 69   |
--R   (2)  [q= | -   - ---   --- |,r= |0   - 175  70 + --- %i|]
--R            | 7     175   175 |    |                175   |
--R            |                 |    |                      |
--R            |  2     6     33 |    |                 58   |
--R            |- -  - --   - -- |    |0     0    35 - --- %i|

--R            +  7    35     35 +    +                175   +
--IType: Record(q: Matrix(Complex(Expression(Integer))),...
--E 2
--S 3 of 5
transpose(map(conjugate,out.q))*out.q
--R
--R        +1  0  0+
--R        |       |
--R   (3)  |0  1  0|
--R        |       |
--R        +0  0  1+
--R                                   Type: Matrix(Complex(Expression(Integer)))
--E 3
--S 4 of 5
out.q*transpose(map(conjugate,out.q))
--R
--R        +1  0  0+
--R        |       |
--R   (4)  |0  1  0|
--R        |       |
--R        +0  0  1+
--R                                   Type: Matrix(Complex(Expression(Integer)))
--E 4
--S 5 of 5
out.q * out.r - m
--R
--R        +0  0  0+
--R        |       |
--R   (5)  |0  0  0|
--R        |       |
--R        +0  0  0+
--R                                   Type: Matrix(Complex(Expression(Integer)))
--E 5
)spool
)lisp (bye)

———-

References^{[1] nothing}

6