DhbIterativeProcess subclass: #DhbJacobiTransformation instanceVariableNames: 'lowerRows transform ' classVariableNames: '' poolDictionaries: ''! DhbIterativeProcess subclass: #DhbLargestEigenValueFinder instanceVariableNames: 'matrix eigenvector transposeEigenvector ' classVariableNames: '' poolDictionaries: ''! SubApplication subclass: #DhbEigenValuesAndVectors instanceVariableNames: '' classVariableNames: '' poolDictionaries: ''! !DhbJacobiTransformation class publicMethods ! matrix: aSymmetricMatrix " (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " ^super new initialize: aSymmetricMatrix! new "Prevent using this message to create instances. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " ^self error: 'Illegal creation message for this class'! ! !DhbJacobiTransformation publicMethods ! evaluateIteration "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | indices | indices := self largestOffDiagonalIndices. self transformAt: ( indices at: 1) and: ( indices at: 2). ^precision! finalizeIterations "Transfer the eigenValues into a vector and set this as the result. eigen values and transform matrix are sorted using a bubble sort. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | n | n := 0. result := lowerRows collect: [:each | n := n + 1. each at: n]. self sortEigenValues! printOn: aStream "Append to the argument aStream, a sequence of characters that describes the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | first | first := true. lowerRows do: [ :each | first ifTrue: [ first := false] ifFalse:[ aStream cr]. each printOn: aStream. ].! transform "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " ^DhbMatrix rows: transform! ! !DhbJacobiTransformation privateMethods ! exchangeAt: anInteger "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | temp n | n := anInteger + 1. temp := result at: n. result at: n put: ( result at: anInteger). result at: anInteger put: temp. transform do: [ :each | temp := each at: n. each at: n put: ( each at: anInteger). each at: anInteger put: temp. ].! initialize: aSymmetricMatrix "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | n m | n := aSymmetricMatrix numberOfRows. lowerRows := Array new: n. transform := Array new: n. 1 to: n do: [ :k | lowerRows at: k put: ( ( aSymmetricMatrix rowAt: k) copyFrom: 1 to: k). transform at: k put: ( ( Array new: n) atAllPut: 0; at: k put: 1; yourself). ]. ^self! largestOffDiagonalIndices "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | n m abs | n := 2. m := 1. precision := ( ( lowerRows at: n) at: m) abs. 1 to: lowerRows size do: [ :i | 1 to: ( i - 1) do: [ :j | abs := ( ( lowerRows at: i) at: j) abs. abs > precision ifTrue: [ n := i. m := j. precision := abs. ]. ]. ]. ^Array with: m with: n! sortEigenValues "Private - Use a bubble sort. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | n bound m | n := lowerRows size. bound := n. [ bound = 0 ] whileFalse: [ m := 0. 1 to: bound - 1 do: [ :j | ( result at: j) abs > ( result at: j + 1) abs ifFalse:[ self exchangeAt: j. m := j. ]. ]. bound := m. ].! transformAt: anInteger1 and: anInteger2 "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 1/6/99 " | d t s c tau apq app aqq arp arq | apq := ( lowerRows at: anInteger2) at: anInteger1. apq = 0 ifTrue: [ ^nil]. app := ( lowerRows at: anInteger1) at: anInteger1. aqq := ( lowerRows at: anInteger2) at: anInteger2. d := aqq - app. arp := d * 0.5 / apq. t := arp > 0 ifTrue: [ 1 / ( ( arp squared + 1) sqrt + arp)] ifFalse:[ 1 / ( arp - ( arp squared + 1) sqrt)]. c := 1 / ( t squared + 1) sqrt. s := t * c. tau := s / ( 1 + c). 1 to: ( anInteger1 - 1) do: [ :r | arp := ( lowerRows at: anInteger1) at: r. arq := ( lowerRows at: anInteger2) at: r. ( lowerRows at: anInteger1) at: r put: ( arp - ( s * (tau * arp + arq))). ( lowerRows at: anInteger2) at: r put: ( arq + ( s * (arp - (tau * arq)))). ]. ( anInteger1 + 1) to: ( anInteger2 - 1) do: [ :r | arp := ( lowerRows at: r) at: anInteger1. arq := ( lowerRows at: anInteger2) at: r. ( lowerRows at: r) at: anInteger1 put: ( arp - ( s * (tau * arp + arq))). ( lowerRows at: anInteger2) at: r put: ( arq + ( s * (arp - (tau * arq)))). ]. ( anInteger2 + 1) to: lowerRows size do: [ :r | arp := ( lowerRows at: r) at: anInteger1. arq := ( lowerRows at: r) at: anInteger2. ( lowerRows at: r) at: anInteger1 put: ( arp - ( s * (tau * arp + arq))). ( lowerRows at: r) at: anInteger2 put: ( arq + ( s * (arp - (tau * arq)))). ]. 1 to: lowerRows size do: [ :r | arp := ( transform at: r) at: anInteger1. arq := ( transform at: r) at: anInteger2. ( transform at: r) at: anInteger1 put: ( arp - ( s * (tau * arp + arq))). ( transform at: r) at: anInteger2 put: ( arq + ( s * (arp - (tau * arq)))). ]. ( lowerRows at: anInteger1) at: anInteger1 put: ( app - (t * apq)). ( lowerRows at: anInteger2) at: anInteger2 put: ( aqq + (t * apq)). ( lowerRows at: anInteger2) at: anInteger1 put: 0.! ! !DhbLargestEigenValueFinder class publicMethods ! matrix: aMatrix "Create a new instance of the receiver for a given matrix and default precision. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 11/2/99 " ^( self new) initialize: aMatrix; yourself! matrix: aMatrix precision: aNumber "Create a new instance of the receiver for a given matrix and desired precision. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 11/2/99 " ^( self new) initialize: aMatrix; desiredPrecision: aNumber; yourself! ! !DhbLargestEigenValueFinder class privateMethods ! defaultMaximumIterations "Private - Answers the default maximum number of iterations for newly created instances. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^100! ! !DhbLargestEigenValueFinder publicMethods ! eigenvalue "Answer the eigen value found by the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 11/2/99 " ^result! eigenvector "Answer the normalized eigen vector found by the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 11/2/99 " ^eigenvector * (1 / eigenvector norm)! evaluateIteration "Iterate the product of the matrix of the eigen vector and the transpose. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " | oldEigenvalue | oldEigenvalue := result. transposeEigenvector := transposeEigenvector * matrix. transposeEigenvector := transposeEigenvector * (1 / (transposeEigenvector at: 1)). eigenvector := matrix * eigenvector. result := eigenvector at: 1. eigenvector := eigenvector * (1 / result). ^oldEigenvalue isNil ifTrue: [2 * desiredPrecision] ifFalse: [(result - oldEigenvalue) abs]! initialize: aMatrix "Defines the matrix for the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 11/2/99 " matrix := aMatrix.! initializeIterations "Initialize the iterations (subclasses must write their own method and call this one last). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " eigenvector := DhbVector new: matrix numberOfRows. eigenvector atAllPut: 1.0. transposeEigenvector := DhbVector new: eigenvector size. transposeEigenvector atAllPut: 1.0! nextLargestEigenValueFinder "Return an eigen value finder for the same eigen values of the receiver except the one found. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 11/2/99 " | norm | norm := 1 / (eigenvector * transposeEigenvector). ^self class new: matrix * ((DhbSymmetricMatrix identity: eigenvector size) - (eigenvector * norm tensorProduct: transposeEigenvector)) precision: desiredPrecision! ! DhbJacobiTransformation initializeAfterLoad! DhbLargestEigenValueFinder initializeAfterLoad! DhbEigenValuesAndVectors initializeAfterLoad!