login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

)abbrev package JACDIAG JacobiDiagonalisation
++ Author: Kurt Pagani
++ Date Created: Sat Jun 16 23:37:27 CEST 2018
++ References:
++  Jacobi, C.G.J. (1846). "Über ein leichtes Verfahren, die in der Theorie
++    der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen".
++  Crelle's Journal 30, pages 51–94.
++  https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
++ Description:
++ The Jacobi eigenvalue algorithm is an iterative method for the calculation
++ of the eigenvalues and eigenvectors of a real symmetric matrix. It is named
++ after Carl Gustav Jacob Jacobi, who first proposed the method in 1846, but
++ only became widely used in the 1950s with the advent of computers.
++ Notes
++  * The tolerance and maximum number of iterations may be changed by the
++    function 'setParams'. Defaults are set in the implementation part.
++  * For the case in which one dimension is an independent subspace, the
++    maximum number of iterations may be reached. If 'checkResult' gives
++    a value smaller than 'eps' then 'jacobi' regularly terminates, otherwise
++    an error message is issued.
++ Usage and Example(s)
++  * See test_jacdiag.input
++  * M:=matrix [[1,1],[1,1]], R:=jacobi(M)
++    R.ev .... eigenvalues
++    R.EV .... eigenvectors as columns of matrix R.EV
++    checkResult(R,M) .... Sum(norm(M*v-l*v))
++    p:=showParams() ..... [maxit= 10000,tol= 0.1 E -8]
++    p.maxit:=20000
++    setParams(p) ........ [maxit= 20000,tol= 0.1 E -8]
++ IDEA: [jedit]return eps, iter in RES ?
++
JacobiDiagonalisation() : Exports == Implementation where
R  ==> Float
PI ==> PositiveInteger
NN ==> NonNegativeInteger
IF ==> R
VIF ==> Vector R
MIF ==> Matrix R
RES ==> Record(ev:VIF,EV:MIF)
PAR ==> Record(maxit:PI,tol:R)
Exports ==  with
jacobi : MIF -> RES
checkResult : (RES,MIF) -> R
showParams : () -> PAR
setParams : PAR -> PAR
eps:R:=0.000000001
maxIter:PI:=10000
showParams():PAR == [maxIter,eps]
setParams(p:PAR):PAR ==
maxIter:=p.maxit
eps:=p.tol
showParams()
maxInd(k:PI,n:NN,S:MIF):PI ==
m:PI:=k+1
for i in k+2..n repeat
if abs(S(k,i)) > abs(S(k,m)) then m:=i::PI
return m
jacobi(M:MIF):RES ==
not square? M => error "not square"
not symmetric? M => error "not symmetric"
S:MIF:=copy M
n:NN:=nrows S
i:PI; k:PI; l:PI; m:PI
state:Integer
s:IF;c:IF;t:IF;p:IF;y:IF;d:IF;r:IF
ind:List(PI):=[1 for i in 1..n]
changed:List Boolean:=[false for i in 1..n]
e:VIF:=new(n,0$IF)$VIF
E:MIF:=new(n,n,0$IF)$MIF
E:=diagonalMatrix [1$IF for i in 1..n] A:IF; B:IF count:Integer:=0 state:=n for k in 1..n repeat ind.k:=maxInd(k::PI,n,S) e.k:=S(k,k) changed.k:=true while (state ~= 0) and (count <= maxIter) repeat m:=1 for k in 2..n-1 repeat if abs(S(k,ind.k)) > abs(S(m,ind.m)) then m:=k::PI -- k:=m l:=ind.m p:=S(k,l) -- y:=(e.l - e.k)/2 d:R:=abs(y)+sqrt(p^2+y^2) r:R:=sqrt(p^2+d^2) c:R:=d/r s:R:=p/r t:R:=p^2/d if y < 0$R then
s:=-s
t:=-t
S(k,l):=0$IF -- y:IF:=e.k e.k:=y-t if changed.k and abs(t)<= eps then changed.k:=false state:=state-1 else if (not changed.k) and abs(t)> eps then changed.k:=true state:=state+1 -- y:IF:=e.l e.l:=y+t if changed.l and abs(t)<= eps then changed.l:=false state:=state-1 else if (not changed.l) and abs(t)> eps then changed.l:=true state:=state+1 -- for i in 1..k-1 repeat A:=S(i,k); B:=S(i,l) S(i,k):=c*A-s*B S(i,l):=s*A+c*B for i in k+1..l-1 repeat A:=S(k,i); B:=S(i,l) S(k,i):=c*A-s*B S(i,l):=s*A+c*B for i in l+1..n repeat A:=S(k,i); B:=S(l,i) S(k,i):=c*A-s*B S(l,i):=s*A+c*B -- for i in 1..n repeat A:=E(i,k); B:=E(i,l) E(i,k):=c*A-s*B E(i,l):=s*A+c*B -- ind.k := maxInd(k,n,S) ind.l := maxInd(l,n,S) -- count:=count+1 --output([count,state]) -- if count > maxIter then if checkResult([e,E]$RES,M) > eps then
error("Maximum iteration count reached")
return [e,E]$RES checkResult(r:RES,M:MIF):R == n:=nrows M s:R:=0$R
for i in 1..n repeat
v:VIF:=column(r.EV,i)
d:VIF:=M*v - r.ev.i * v
s:=s+sqrt(dot(d,d))
return s
   Compiling FriCAS source code from file
using old system compiler.
JACDIAG abbreviates package JacobiDiagonalisation
------------------------------------------------------------------------
initializing NRLIB JACDIAG for JacobiDiagonalisation
compiling into NRLIB JACDIAG
compiling exported showParams : () -> Record(maxit: PositiveInteger,tol: Float)
Time: 0 SEC.
compiling exported setParams : Record(maxit: PositiveInteger,tol: Float) -> Record(maxit: PositiveInteger,tol: Float)
Time: 0 SEC.
compiling local maxInd : (PositiveInteger,NonNegativeInteger,Matrix Float) -> PositiveInteger
Time: 0.03 SEC.
compiling exported jacobi : Matrix Float -> Record(ev: Vector Float,EV: Matrix Float)
Time: 0.09 SEC.
compiling exported checkResult : (Record(ev: Vector Float,EV: Matrix Float),Matrix Float) -> Float
Time: 0.02 SEC.
(time taken in buildFunctor:  0)
;;;     ***       |JacobiDiagonalisation| REDEFINED
;;;     ***       |JacobiDiagonalisation| REDEFINED
Time: 0 SEC.
Warnings:
[1] setParams:  maxit has no value
[2] setParams:  tol has no value
[3] maxInd:  m has no value
[4] jacobi:  m has no value
[5] jacobi:  t has no value
[6] jacobi:  state has no value
[7] jacobi:  s has no value
[8] checkResult:  EV has no value
[9] checkResult:  ev has no value
Cumulative Statistics for Constructor JacobiDiagonalisation
Time: 0.14 seconds
finalizing NRLIB JACDIAG
Processing JacobiDiagonalisation for Browser database:
--------constructor---------
--->-->JacobiDiagonalisation((jacobi ((Record (: ev (Vector (Float))) (: EV (Matrix (Float)))) (Matrix (Float))))): Not documented!!!!
--->-->JacobiDiagonalisation((checkResult ((Float) (Record (: ev (Vector (Float))) (: EV (Matrix (Float)))) (Matrix (Float))))): Not documented!!!!
--->-->JacobiDiagonalisation((showParams ((Record (: maxit (PositiveInteger)) (: tol (Float)))))): Not documented!!!!
--->-->JacobiDiagonalisation((setParams ((Record (: maxit (PositiveInteger)) (: tol (Float))) (Record (: maxit (PositiveInteger)) (: tol (Float)))))): Not documented!!!!
; compiling file "/var/aw/var/LatexWiki/JACDIAG.NRLIB/JACDIAG.lsp" (written 16 JUL 2018 08:09:39 PM):
; /var/aw/var/LatexWiki/JACDIAG.NRLIB/JACDIAG.fasl written
; compilation finished in 0:00:00.415
------------------------------------------------------------------------
JacobiDiagonalisation is now explicitly exposed in frame initial
JacobiDiagonalisation will be automatically loaded when needed from
/var/aw/var/LatexWiki/JACDIAG.NRLIB/JACDIAG

Test flavours

fricas
-- Unittest .....: Package JacobiDiagonalisation
Author .......: Kurt Pagani -- Date .........: Sun Jun 17 02:14:10 CEST 2018 -- Expected .....: total tests: 16 -- Reference ....: Gregory R.T, Karney D.L., A collection of matrices for testing computational algorithms, -- ............... Krieger (1978), ISBN 0882756494
fricas
)set break resume

fricas
)expose UnittestCount UnittestAux Unittest
UnittestCount is now explicitly exposed in frame initial
UnittestAux is now explicitly exposed in frame initial
Unittest is now explicitly exposed in frame initial
--)co jacdiag
fricas
)expose JACDIAG
JacobiDiagonalisation is already explicitly exposed in frame initial
-- Test Suite
testsuite "Package JacobiDiagonalisation"
All user variables and function definitions have been cleared.
WARNING: string for testsuite should have less than 15 characters!
Type: Void
fricas
-- Cases
testcase "jacobi"
All user variables and function definitions have been cleared.
Type: Void
fricas
-- https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
S:= matrix [[4,-30,60,-35], [-30,300,-675,420],
[60,-675,1620,-1050],[-35,420,-1050,700]]
 (1)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (2)
Type: Matrix(Float)
fricas
R:=jacobi(M)
 (3)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
ev1:= vector [0.1666428611_7189046264, 37.1014913651_27658188,
2585.2538109289_223144, 1.4780548447_781369118]
 (4)
Type: Vector(Float)
fricas
EV1:= matrix [[0.7926082911_6404284196, -0.1791862905_3191911245, _
0.0291933231_7921032382_1, -0.5820756994_9722239026], _
[0.4519231209_0048280218, 0.7419177906_0266858684, _
-0.3287120558_2291241902, 0.3705021850_6710175866], _
[0.3224163985_8196511663, -0.1002281368_8300231266, _
0.7914111458_4119456574, 0.5095786345_0180583406], _
[0.2521611696_891868645, -0.6382825282_3465850686, _
-0.5145527499_457719896, 0.5140482722_2216914813]]
 (5)
Type: Matrix(Float)
fricas
eps:=0.0000000001
 (6)
Type: Float
fricas
floatNull(v:Vector Float):Boolean == sqrt dot(v,v) <= eps
Function declaration floatNull : Vector(Float) -> Boolean has been
added to workspace.
Type: Void
fricas
vector members(map(abs,R.EV-EV1))
 (7)
Type: Vector(Float)
fricas
testEquals("floatNull(R.ev-ev1)", "true")
fricas
Compiling function floatNull with type Vector(Float) -> Boolean
Type: Void
fricas
testEquals("floatNull(vector members(map(abs,R.EV-EV1)))", "true")
Type: Void
fricas
-- Gregory R.T, Karney D.L, A collection of matrices for testing
-- computational algorithms, Krieger, 1978
-- CHAPTER IV, Example 4.1
S:= matrix [[5,4,1,1], [4,5,1,1], [1,1,4,2],[1,1,2,4]]
 (8)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (9)
Type: Matrix(Float)
fricas
p:=showParams()
 (10)
Type: Record(maxit: PositiveInteger?,tol: Float)
fricas
p.maxit:=20000
 (11)
Type: PositiveInteger?
fricas
p.tol:=0.00000000001
 (12)
Type: Float
fricas
setParams(p)
 (13)
Type: Record(maxit: PositiveInteger?,tol: Float)
fricas
R1:=jacobi(M)
 ?}}, \: \right. \ \ \displaystyle \left.{ \begin{array}{@{}l} \displaystyle EV = \ \ \displaystyle {\left[ \begin{array}{cccc} {0.7071067811 \ 865475244}&{0.6324555320 \ 336758664}&{0.0}& -{0.3162277660 \ 168379332} \ -{0.7071067811 \ 865475244}&{0.6324555320 \ 336758664}&{0.0}& -{0.3162277660 \ 168379332} \ {0.0}&{0.3162277660 \ 168379332}&{0.7071067811 \ 865475244}&{0.6 324555320 \ 336758664} \ {0.0}&{0.3162277660 \ 168379332}& -{0.7071067811 \ 86547524 4}&{0.6324555320 \ 336758664} " class="equation" src="images/3192297519243230694-16.0px.png" align="bottom" Style="vertical-align:text-bottom" width="793" height="155"/> (14)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R1,M)<eps","true")
Type: Void
fricas
-- Example 4.2
S:= matrix [[6,4,4,1], [4,6,1,4], [4,1,6,4],[1,4,4,6]]
 (15)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (16)
Type: Matrix(Float)
fricas
R2:=jacobi(M)
 (17)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R2,M)<eps","true")
Type: Void
fricas
-- Example 4.3
S:= matrix [[2,1,3,4], [1,-3,1,5], [3,1,6,-2],[4,5,-2,-1]]
 (18)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (19)
Type: Matrix(Float)
fricas
R3:=jacobi(M)
 (20)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R3,M)<eps","true")
Type: Void
fricas
-- Example 4.4
S:= matrix [[5,7,6,5], [7,10,8,7], [6,8,10,9],[5,7,9,10]]
 (21)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (22)
Type: Matrix(Float)
fricas
R4:=jacobi(M)
 (23)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R4,M)<eps","true")
Type: Void
fricas
-- Example 4.5
eps:=0.00000001  -- eps from above too small
 (24)
Type: Float
fricas
M:= matrix [[0.81321,-0.00013,0.00014,0.00011,0.00021], _
[-0.00013,0.93125,0.23567,0.41235,0.41632], _
[0.00014,0.23567,0.18765,0.50632,0.30697],  _
[0.00011,0.41235,0.50632,0.27605,0.46322],  _
[0.00021,0.41632,0.30697,0.46322,0.41931]]
 (25)
Type: Matrix(Float)
fricas
R5:=jacobi(M)
 (26)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R5,M)<eps","true")
Type: Void
fricas
-- Example 4.6
eps:=0.0000001  -- eps from above too small
 (27)
Type: Float
fricas
S:= matrix [[5,4,3,2,1],[4,6,0,4,3],[3,0,7,6,5],[2,4,6,8,7],[1,3,5,7,9]]
 (28)
Type: Matrix(NonNegativeInteger?)
fricas
M:=map(s+->s::Float, S)
 (29)
Type: Matrix(Float)
fricas
R6:=jacobi(M)
 (30)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R6,M)<eps","true")
Type: Void
fricas
-- Example 4.7
S:= matrix [[10,1,2,3,4],[1,9,-1,2,-3],[2,-1,7,3,-5],[3,2,3,12,-1],_
[4,-3,-5,-1,15]]
 (31)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (32)
Type: Matrix(Float)
fricas
R7:=jacobi(M)
 (33)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R7,M)<eps","true")
Type: Void
fricas
-- Example 4.8
S:= matrix [[5,1,-2,0,-2,5],[1,6,-3,2,0,6],[-2,-3,8,-5,-6,0],[0,2,-5,5,1,-2],_
[-2,0,-6,1,6,-3],[5,6,0,-2,-3,8]]
 (34)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (35)
Type: Matrix(Float)
fricas
R8:=jacobi(M)
 (36)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R8,M)<eps","true")
Type: Void
fricas
-- Example 4.9
S:= matrix [[1,2,3,0,1,2],[2,4,5,-1,0,3],[3,5,6,-2,-3,0],[0,-1,-2,1,2,3],_
[1,0,-3,2,4,5],[2,3,0,3,5,6]]
 (37)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (38)
Type: Matrix(Float)
fricas
R9:=jacobi(M)
 (39)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R9,M)<eps","true")
Type: Void
fricas
-- Example 4.10
eps:=0.00001  -- eps from above too small
 (40)
Type: Float
fricas
S:= matrix [[611,196,-192,407,-8,-52,-49,29], _
[196,899,113,-192,-71,-43,-8,-44], _
[-192,113,899,196,61,49,8,52], _
[407,-192,196,611,8,44,59,-23],_
[-8,-71,61,8,411,-599,208,208], _
[-52,-43,49,44,-599,411,208,208], _
[-49,-8,8,59,208,208,99,-911], _
[29,-44,52,-23,208,208,-911,99]]
 (41)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float, S)
 (42)
Type: Matrix(Float)
fricas
R10:=jacobi(M)
 (43)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R10,M)<eps","true")
Type: Void
fricas
-- Example 4.11
eps:=0.00001  -- eps from above too small
 (44)
Type: Float
fricas
S:= diagonalMatrix [5,6,6,6,6,6,6,6,6,6,5]
 (45)
Type: Matrix(Integer)
fricas
for i in 2..11 repeat S(i-1,i):=S(i,i-1):=3
Type: Void
fricas
for i in 3..11 repeat S(i-2,i):=S(i,i-2):=1
Type: Void
fricas
for i in 4..11 repeat S(i-3,i):=S(i,i-3):=1
Type: Void
fricas
S(1,2):=S(2,1):=S(10,11):=S(11,10):=2
 (46)
Type: PositiveInteger?
fricas
M:=map(s+->s::Float, S)
 (47)
Type: Matrix(Float)
fricas
R11:=jacobi(M)
 (48)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R11,M)<eps","true")
Type: Void
fricas
-- Example 4.12
M:=matrix [[0.2500,0.06675,0.04000,0.02475,0.07050,0.06375,0.06925,0.02050,_
0.03600, -0.01025,-0.00175,0.02750,0.02300,0.00200], _
[0.06675,0.25000,0.10400,0.07475,0.03625,0.11675,0.11050,0.06225,0.05100, _
0.03250,0.02400,0.03600,0.06350,0.05300], _
[0.04000,0.10400,0.25000,0.14575,0.03725,0.07175,0.07800,0.12200,0.11275, _
0.09375,0.10175,0.09600,0.14300,0.11550], _
[0.02475,0.07475,0.14575,0.25000,0.05375,0.07000,0.05225,0.12800,0.12475, _
0.10550,0.13000,0.14575,0.13975,0.13375], _
[0.07050,0.03625,0.03725,0.05375,0.25000,0.04575,0.05750,0.05700,0.05050, _
0.01475,0.04500,0.07150,0.05300,0.01600], _
[0.06375,0.11675,0.07175,0.07000,0.04575,0.25000,0.08625,0.08800,0.07150, _
0.04850,0.03200,0.04475,0.03300,0.04500], _
[0.06925,0.11050,0.07800,0.05225,0.05750,0.08625,0.25000,0.08725,0.06950, _
0.03725,0.04025,0.04300,0.04075,0.01450], _
[0.02050,0.06225,0.12200,0.12800,0.05700,0.08800,0.08725,0.25000,0.14100, _
0.13275,0.15550,0.13050,0.11825,0.09125], _
[0.03600,0.05100,0.11275,0.12475,0.05050,0.07150,0.06950,0.14100,0.25000, _
0.07425,0.10750,0.09175,0.10725,0.08225], _
[-0.01025,0.03250,0.09375,0.10550,0.01475,0.04850,0.03725,0.13275,0.07425, _
0.25000,0.15500,0.09625,0.09950,0.09425], _
[-0.00175,0.02400,0.10175,0.13000,0.04500,0.03200,0.04025,0.15550,0.10750, _
0.15500,0.25000,0.13350,0.14850,0.13050], _
[0.02750,0.03600,0.09600,0.14575,0.07150,0.04475,0.04300,0.13050,0.09175, _
0.09625,0.13350,0.25000,0.11100,0.10075], _
[0.02300,0.06350,0.14300,0.13975,0.05300,0.03300,0.04075,0.11825,0.10725, _
0.09950,0.14850,0.11100,0.25000,0.14325], _
[0.00200,0.05300,0.11550,0.13375,0.01600,0.04500,0.01450,0.09125,0.08225, _
0.09425,0.13050,0.10075,0.14325,0.25000]]
 (49)
Type: Matrix(Float)
fricas
R12:=jacobi(M)
 (50)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R12,M)<eps","true")
Type: Void
fricas
-- Example 4.13 (Hilbert Matrix)
H(n) == matrix [[1/(i+j-1) for i in 1..n] for j in 1..n]
Type: Void
fricas
M:=H(20);
fricas
Compiling function H with type PositiveInteger -> Matrix(Fraction(
Integer))
Type: Matrix(Fraction(Integer))
fricas
R13:=jacobi(M)
 (51)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R13,M)<eps","true")
Type: Void
fricas
-- Example 4.23
M:=matrix [[1 for i in 1..50] for j in 1..50];
Type: Matrix(Integer)
fricas
R23:=jacobi(M);
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
testEquals("checkResult(R23,M)<eps","true")
Type: Void
fricas
-- Results/Statistics
fricas
)set output algebra on

fricas
)version
Value = "FriCAS 1.3.4 compiled at Wed Jun 27 14:53:01 UTC 2018"
fricas
)lisp (lisp-implementation-type)
Value = "SBCL"
fricas
)lisp (lisp-implementation-version)
Value = "1.1.1"
statistics()
=============================================================================
General WARNINGS:
* do not use ')clear completely' before having used 'statistics()'
It clears the statistics without warning!
* do not forget to pass the arguments of the testXxxx functions as Strings!
Otherwise, the test will fail and statistics() will not notice!
* testLibraryError does not prevent FriCAS from aborting the current block.
Thus, if a block contains other test functions, they will not be executed
and statistics() will not notice!
=============================================================================
Testsuite: Package JacobiDiagonalisation
failed (total): 0 (1)
=============================================================================
testsuite | testcases: failed (total) | tests: failed (total)
Package JacobiDiagonalisation    0     (1)               0    (16)
=============================================================================
File summary.
unexpected failures: 0
expected failures: 0
unexpected passes: 0
total tests: 16
Type: Void

 Subject:   Be Bold !! ( 14 subscribers )