[RELEASE] - Matrix Computations

[ctl3d32 0]

Hi Folks!

Here is an AutoIT script that performs some matrix computations.

Reference:

FORTRAN 90 for Scientists and Engineers by Brian Hahn,

March 1994 in paperback,

368 pages,

ISBN 0 340 60034 9 ?6.95

Enjoy,

ctl3d32

expandcollapse popup

#include <Array.au3>

Local $A   [3][3] = [[2,2,2],[3,2,2],[3,2,3]]
Local $Aa  [2][2] = [[1,2],[3,4]]
Local $Bb  [2][2] = [[5,6],[0,-1]]
Local $b   [3][1] = [[0],[1],[1]]

_ArrayDisplay(_Inv($A), "Matrix Inversion")

_ArrayDisplay(_sMul($A,2), "Matrix * scalar")

_ArrayDisplay(_mMul($Aa,$Bb), "Matrix * Matrix")

_ArrayDisplay(_SolveLE($A,$b), "LE Solution")

_ArrayDisplay(_mAdd($Aa,$Bb), "Matrix + Matrix")

_ArrayDisplay(_mSub($Aa,$Bb), "Matrix - Matrix")

Func _Inv($M)
    ; Matrix inversion by Gauss reduction
    
    ; Reference:
    ;      FORTRAN 90 for Scientists and Engineers by Brian Hahn,
    ;      March 1994 in paperback,
    ;      368 pages,
    ;      ISBN 0 340 60034 9   ?6.95
    
    Local $N = UBound($M,1) ;Number of Equations
    Local $A [$N][2*$N]     ;Augmented Matrix
    Local $Inv [$N][$N]     ;Inverted Matrix
    Local $TempRow[2*$N]    ;Temporary Row
    local $PivElt, $TarElt, $PivRow, $TarRow, $K, $i, $j
    
    #region - Filling Matrix A
    For $i = 0 To $N-1
        For $j = 0 To $N-1
            $A [$i][$j] = $M [$i][$j]
        Next
    Next
    For $i = 0 To $N-1
        $A [$i][$N+$i] = 1
    Next
    #endregion - Filling Matrix A
    
    #region - Inverting Matrix $M
    For $PivRow = 0 To $N-1                             ;Process every row
        $PivElt = $A [$PivRow][$PivRow]                 ;Choose pivot element
        If ($PivElt == 0) Then
            $K = $PivRow + 1                            ;Run down rows to find a non-zero pivot
            While($PivElt == 0 And $K <= $N)
                $PivElt = $A[$K][$PivRow]               ;Try next row
                $K = $K + 1
            WEnd
            If ($PivElt == 0) Then
                SetError(1)                             ;Couldn't find a non-zero pivot: solution rubbish
                Return
            Else                                        ;Non-zero pivot in row K, so swop rows PivRow and K
                For $i = 0 To 2*$N-1
                    $TempRow [$PivRow][$i] = $A [$PivRow][$i]
                    $K = $K - 1                         ;Adjust for overcount
                Next
                For $i = 0 To 2*$N-1
                    $A [$PivRow][$i] = $A [$K][$i]
                Next
                For $i = 0 To 2*$N-1
                    $A [$K][$i] = $TempRow [$K][$i]
                Next
            EndIf
        EndIf
        For $i = 0 To 2*$N-1
            $A [$PivRow][$i] = $A [$PivRow][$i]/$PivElt ;Divide whole row
        Next
        #region - Now replace all other rows by target row minus pivot row times element in target row and pivot column 
        For $TarRow = 0 To $N-1
            If ($TarRow <> $PivRow) Then
                $TarElt = $A [$TarRow][$PivRow]
                For $i = 0 To 2*$N-1
                    $A [$TarRow][$i] = $A [$TarRow][$i] - $A [$PivRow][$i]*$TarElt
                Next
            EndIf
        Next
        #endregion - Now replace all other rows by target row minus pivot row times element in target row and pivot colum
    Next
    #region - finally extract the inverse from columns N+1 to 2N
    For $i = 0 To $N-1
        For $j = 0 To $N-1
            $Inv [$i][$j] = $A [$i][$N+$j]
        Next
    Next
    #endregion - finally extract the inverse from columns N+1 to 2N
    Return $Inv
    #endregion - Inverting Matrix $M
EndFunc

Func _sMul($M,$s)
    ; Matrix multiplication by a scalar
    
    Local $p   = UBound($M,1)
    Local $q   = UBound($M,2)
    Local $Mat [$p][$q]
    Local $i, $j    

    For $i = 0 To $p-1
        For $j = 0 To $q-1
            $Mat [$i][$j] = $M [$i][$j]*$s
        Next
    Next
    Return $Mat
EndFunc

Func _mMul($A,$B)
    ; multiplies A (n x m) by B (m x p) and returns product in C (n x p)
    ; Formula: cij = Sum(k=1:m,aik*bik)
    
    ; Reference:
    ;      FORTRAN 90 for Scientists and Engineers by Brian Hahn,
    ;      March 1994 in paperback,
    ;      368 pages,
    ;      ISBN 0 340 60034 9   ?6.95
    
    Local $n = UBound($A,1)
    Local $m = UBound($A,2)
    Local $o = UBound($B,1)
    Local $p = UBound($B,2)
    Local $C [$n][$p]
    Local $s, $i, $j, $k
    
    If ($m <> $o) Then
        SetError(1) ;Matrix dimensions incompatibles
        Return
    Else
        For $i = 0 To $n-1
            For $j = 0 To $p-1
                For $k = 0 To $m-1
                    $s = $s + $A [$i][$k] * $B [$k][$j]
                Next
                $C [$i][$j] = $s
                $s = 0
            Next
        Next    
    EndIf
    Return $C
EndFunc

Func _SolveLE($A,$b)
    ;Solution of a system of linear equations, e.g.: 
    ;   2x + 2y + 2z = 0
    ;   3x + 2y + 2z = 1
    ;   3x + 2y + 3z = 1
    
    ;       2 2 2         0
    ;   A = 3 2 2     b = 1
    ;       3 2 3         1
    
    ; Reference:
    ;      FORTRAN 90 for Scientists and Engineers by Brian Hahn,
    ;      March 1994 in paperback,
    ;      368 pages,
    ;      ISBN 0 340 60034 9   ?6.95
    
    Local $mA = UBound($A,1)
    Local $mb = UBound($b,1)
    Local $nb = UBound($b,2)
    
    If (($mA <> $mb) Or ($nb <> 1)) Then
        SetError(1) ;Matrix dimensions does not match
        Return
    Else
        $Solution = _mMul(_Inv($A),$b)
        Return $Solution
    EndIf
EndFunc

Func _mAdd($A,$B)
    ; Add Matrix A to Matrix B
    
    Local $mA = UBound($A,1)
    Local $nA = UBound($A,2)
    Local $mB = UBound($B,1)
    Local $nB = UBound($B,2)
    
    If (($mA<>$mB) Or ($nA<>$nB)) Then
        SetError(1) ;Matrix dimensions does not match
        Return
    Else
        Local $C [$mA][$nB] 
        For $i = 0 To $mA-1
            For $j = 0 To $nA-1
                $C [$i][$j] = $A [$i][$j] + $B [$i][$j]
            Next
        Next
        Return $C
    EndIf   
EndFunc

Func _mSub($A,$B)
    ; Subtract Matrix B to Matrix A
    
    Local $mA = UBound($A,1)
    Local $nA = UBound($A,2)
    Local $mB = UBound($B,1)
    Local $nB = UBound($B,2)
    
    If (($mA<>$mB) Or ($nA<>$nB)) Then
        SetError(1) ;Matrix dimensions does not match
        Return
    Else
        Local $C [$mA][$nB] 
        For $i = 0 To $mA-1
            For $j = 0 To $nA-1
                $C [$i][$j] = $A [$i][$j] - $B [$i][$j]
            Next
        Next
        Return $C
    EndIf   
EndFunc

Edited June 1, 2010 by ctl3d32

[jchd 1,012]

Thank you for this nice "addition", even if AutoIt isn't particularly efficient doing matrix/tensors operations!

It might prove useful some day.

[ctl3d32 0]

Thanks!

I will soon add more features.

Enjoy