SubApplication subclass: #DhbMatrixInversion
    instanceVariableNames: ''
    classVariableNames: ''
    poolDictionaries: ''!



!DhbLUPDecomposition publicMethods !

inverseMatrixComponents
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 30/3/99 "
	| n inverseRows column |
	permutation isNil
		ifTrue: [ self protectedDecomposition].
	permutation = 0
		ifTrue: [ ^nil].	"Singular matrix has no inverse"
	n := rows size.
	inverseRows :=( 1 to: n) asVector collect: [ :j | DhbVector new: n].
	1 to: n do:
		[ :j |
		  column := self solve: ( ( Array new: rows size) atAllPut: 0; at: j put: 1; yourself).
		  1 to: n do: [ :i | ( inverseRows at: i) at: j put: ( column at: i)].
		].
	^inverseRows
	
		
	! !

!DhbMatrix class publicMethods !

lupCRLCriticalDimension
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 6/5/99 "
	^40! !

!DhbMatrix publicMethods !

inverse
	"Answer the inverse of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "

	^self isSquare 
		ifTrue: [self lupInverse]
		ifFalse: [self squared inverse * self transpose]!

lupInverse
	^self class rows: self lupDecomposition inverseMatrixComponents! !

!DhbMatrix privateMethods !

largestPowerOf2SmallerThan: anInteger
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved
		 Initial code: 21/3/99 "
	| m m2|
	m := 2.
	[ m2 := m * 2.
	  m2 < anInteger] whileTrue:[ m := m2].
	^m! !

!DhbSymmetricMatrix class publicMethods !

join: anArrayOfMatrices
		"Inverse of the split operation.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	| rows n |
	rows := OrderedCollection new.
	n := 0.
	( anArrayOfMatrices at: 1) rowsDo:
		[ :each |
		  n := n + 1.
		  rows add: each, ( ( anArrayOfMatrices at: 3) columnAt: n).
		].
	n := 0.
	( anArrayOfMatrices at: 2) rowsDo:
		[ :each |
		  n := n + 1.
		  rows add: ( ( anArrayOfMatrices at: 3) rowAt: n), each.
		].
	^self rows: rows !

lupCRLCriticalDimension
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 6/5/99 "
	^36! !

!DhbSymmetricMatrix publicMethods !

crlInverse
	| matrices b1 cb1ct cb1 |
	matrices := self split.
	b1 := (matrices at: 1) inverse.
	cb1 := (matrices at: 3) * b1.
	cb1ct := (cb1 productWithTransposeMatrix: (matrices at: 3)) 
				asSymmetricMatrix.
	matrices at: 3 put: (matrices at: 2) * cb1.
	matrices at: 2 put: ((matrices at: 2) accumulateNegated: cb1ct) inverse.
	matrices at: 1 put: ( b1 accumulate: (cb1 transposeProductWithMatrix: (matrices at: 3))).
	(matrices at: 3) negate.
	^self class join: matrices!

inverse
	"Answer the inverse of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "

	^(rows size < self class lupCRLCriticalDimension or: [lupDecomposition notNil]) 
			ifTrue: [self lupInverse]
			ifFalse: [self crlInverse]! !

!DhbSymmetricMatrix privateMethods !

split
		"Private - 
		 Answers an array of 3 matrices split from the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	| n b d c |
	n := self largestPowerOf2SmallerThan: rows size.
	^Array with: ( self class rows: ( ( 1 to: n) asVector collect: [ :k | ( rows at: k) copyFrom: 1 to: n]))
			  with:( self class rows: ( ( (n+1) to: rows size) asVector collect: [ :k | ( rows at: k) copyFrom: (n+1) to: rows size]))
			  with: ( self class superclass rows: ( ( (n+1) to: rows size) asVector collect: [ :k | ( rows at: k) copyFrom: 1 to: n]))! !

DhbMatrixInversion initializeAfterLoad!