@@ -746,35 +746,50 @@ PMMatrix >> productWithVector: aVector [
746746
747747{ #category : #' as yet unclassified' }
748748PMMatrix >> qrFactorization [
749- |identMat q r hh colSize i |
750- self numberOfRows < self numberOfColumns ifTrue: [ self error: ' numberOfRows<numberOfColumns' ].
751- r := PMMatrix rows: (rows deepCopy).
749+
750+ | identMat q r hh colSize i |
751+ self numberOfRows < self numberOfColumns ifTrue: [
752+ self error: ' numberOfRows<numberOfColumns' ].
753+ r := PMMatrix rows: rows deepCopy.
752754 colSize := self numberOfRows.
753- q := PMSymmetricMatrix identity: colSize.
755+ q := PMSymmetricMatrix identity: colSize.
754756 identMat := q deepCopy.
755- 1 to: self numberOfColumns do: [:col |
757+ 1 to: self numberOfColumns do: [ :col |
756758 hh := ((r columnAt: col) copyFrom: col to: colSize) householder.
757- i := (PMVector new : col- 1withAll: 0 ) , (hh at: 2 ).
758- q := q* (identMat - ((hh at: 1 )* i tensorProduct: i ))." not really necessary, should be simplified"
759- i := PMMatrix rows: ( (r rows allButFirst: (col - 1 )) collect: [:aRow | aRow allButFirst: (col - 1 )] ).
760- i := i - ((hh at: 2 ) tensorProduct: ( (hh at: 1 )* (hh at: 2 )* i ) ) .
761- i rows withIndexDo: [ :aRow :index |
762- aRow withIndexDo: [ :n :c | r rowAt: (col + index - 1 ) columnAt: (col + c - 1 ) put: ((n equalsTo: 0 ) ifTrue: [0 ] ifFalse: [n] ) ] ] .
763- " col <colSize ifTrue: [i :=(hh at: 2) copyFrom: 2 to: colSize -col +1. i withIndexDo: [:n :index| r rowAt: col columnAt: index put: n ] ]" " and this part is not correct, dont uncomment before the bug is corrected! useful if q is not explicitely necessary" ].
764- " r rows allButFirst withIndexDo: [:aRow :ri|1 to: (ri min: self numberOfColumns ) do: [:ci|aRow at: ci put:0 ] ] " " not necessary with equalsTo:0"
765- i := 0 .
766- [(r rowAt: colSize) allSatisfy: [:n | n= 0 ] ]whileTrue: [i := i+ 1 .colSize := colSize - 1 ].
767- i> 0 ifTrue: [ r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
768- i := q numberOfColumns - i.
769- q := PMMatrix rows: ( q rows collect: [:row | row copyFrom: 1 to: i]) ].
770- ^ Array with: q with: r
759+ i := (PMVector new : col - 1 withAll: 0 ) , (hh at: 2 ).
760+ q := q * (identMat - ((hh at: 1 ) * i tensorProduct: i)). " not really necessary, should be simplified"
761+ i := PMMatrix rows:
762+ ((r rows allButFirst: col - 1 ) collect: [ :aRow |
763+ aRow allButFirst: col - 1 ]).
764+ i := i - ((hh at: 2 ) tensorProduct: (hh at: 1 ) * (hh at: 2 ) * i).
765+ i rows withIndexDo: [ :aRow :index |
766+ aRow withIndexDo: [ :n :c |
767+ r
768+ rowAt: col + index - 1
769+ columnAt: col + c - 1
770+ put: ((n closeTo: 0 )
771+ ifTrue: [ 0 ]
772+ ifFalse: [ n ]) ] ]
773+ " col <colSize ifTrue: [i :=(hh at: 2) copyFrom: 2 to: colSize -col +1. i withIndexDo: [:n :index| r rowAt: col columnAt: index put: n ] ]" " and this part is not correct, dont uncomment before the bug is corrected! useful if q is not explicitely necessary" ].
774+ " r rows allButFirst withIndexDo: [:aRow :ri|1 to: (ri min: self numberOfColumns ) do: [:ci|aRow at: ci put:0 ] ] " " not necessary with equalsTo:0"
775+ i := 0 .
776+ [ (r rowAt: colSize) allSatisfy: [ :n | n = 0 ] ] whileTrue: [
777+ i := i + 1 .
778+ colSize := colSize - 1 ].
779+ i > 0 ifTrue: [
780+ r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
781+ i := q numberOfColumns - i.
782+ q := PMMatrix rows:
783+ (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
784+ ^ Array with: q with: r
771785]
772786
773787{ #category : #' as yet unclassified' }
774788PMMatrix >> qrFactorizationWithPivoting [
789+
775790 | identMat q r hh colSize i lengthArray rank mx pivot |
776- self numberOfRows < self numberOfColumns
777- ifTrue: [ self error: ' numberOfRows<numberOfColumns' ].
791+ self numberOfRows < self numberOfColumns ifTrue: [
792+ self error: ' numberOfRows<numberOfColumns' ].
778793 lengthArray := self columnsCollect: [ :col | col * col ].
779794 mx := lengthArray indexOf: lengthArray max.
780795 pivot := Array new : lengthArray size.
@@ -791,41 +806,41 @@ PMMatrix >> qrFactorizationWithPivoting [
791806 hh := ((r columnAt: rank) copyFrom: rank to: colSize) householder.
792807 i := (PMVector new : rank - 1 withAll: 0 ) , (hh at: 2 ).
793808 q := q * (identMat - ((hh at: 1 ) * i tensorProduct: i)).
794- i := PMMatrix rows: ((r rows allButFirst: rank - 1 ) collect: [ :aRow | aRow allButFirst: rank - 1 ]).
809+ i := PMMatrix rows:
810+ ((r rows allButFirst: rank - 1 ) collect: [ :aRow |
811+ aRow allButFirst: rank - 1 ]).
795812 i := i - ((hh at: 2 ) tensorProduct: (hh at: 1 ) * (hh at: 2 ) * i).
796- i rows
797- withIndexDo: [ :aRow :index |
798- aRow
799- withIndexDo: [ :n :c |
800- r
801- rowAt: rank + index - 1
802- columnAt: rank + c - 1
803- put:
804- ((n equalsTo: 0 )
805- ifTrue: [ 0 ]
806- ifFalse: [ n ]) ] ].
807- rank + 1 to: lengthArray size do: [ :ind | lengthArray at: ind put: (lengthArray at: ind) - (r rowAt: rank columnAt: ind) squared ].
813+ i rows withIndexDo: [ :aRow :index |
814+ aRow withIndexDo: [ :n :c |
815+ r
816+ rowAt: rank + index - 1
817+ columnAt: rank + c - 1
818+ put: ((n closeTo: 0 )
819+ ifTrue: [ 0 ]
820+ ifFalse: [ n ]) ] ].
821+ rank + 1 to: lengthArray size do: [ :ind |
822+ lengthArray
823+ at: ind
824+ put: (lengthArray at: ind) - (r rowAt: rank columnAt: ind) squared ].
808825 rank < lengthArray size
809826 ifTrue: [
810827 mx := (lengthArray copyFrom: rank + 1 to: lengthArray size) max.
811- (mx equalsTo: 0 )
812- ifTrue: [ mx := 0 ].
828+ (mx closeTo: 0 ) ifTrue: [ mx := 0 ].
813829 mx := mx > 0
814- ifTrue: [ lengthArray indexOf: mx startingAt: rank + 1 ]
815- ifFalse: [ 0 ] ]
830+ ifTrue: [ lengthArray indexOf: mx startingAt: rank + 1 ]
831+ ifFalse: [ 0 ] ]
816832 ifFalse: [ mx := 0 ].
817833 mx > 0 ] whileTrue.
818834 i := 0 .
819- [ (r rowAt: colSize) allSatisfy: [ :n | n = 0 ] ]
820- whileTrue: [
821- i := i + 1 .
822- colSize := colSize - 1 ].
823- i > 0
824- ifTrue: [
825- r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
826- i := q numberOfColumns - i.
827- pivot := pivot copyFrom: 1 to: i.
828- q := PMMatrix rows: (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
835+ [ (r rowAt: colSize) allSatisfy: [ :n | n = 0 ] ] whileTrue: [
836+ i := i + 1 .
837+ colSize := colSize - 1 ].
838+ i > 0 ifTrue: [
839+ r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
840+ i := q numberOfColumns - i.
841+ pivot := pivot copyFrom: 1 to: i.
842+ q := PMMatrix rows:
843+ (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
829844 ^ Array with: q with: r with: pivot
830845]
831846
0 commit comments