XProfan

WindowTitle "Matrixinversion"
Font 2:randomize:cls rnd(8^8)
' Q: https://www.rhirte.de/vb/gleichsys.htm#mat
' Für XProfan adaptiert von P. Specht 2012-04
' Demoware, keine wie immer geartete Gewähr!
Var n&=12'= Zeilen, Spalten (Testmatrix-Größe)
Declare A![n&,n&],erg$
@MatrixAufruf n&
' Testmatrix mit Zufallszahlen belegen.
' Anschließend:
' Ob die Matrixinversion richtig rechnet, kann man einfach
' testen, denn die Inverse der Inversen muß wieder die
' Ausgangsmatrix ergeben. Im Beispiel wird deshalb die
' größte absolute Abweichung ausgegeben.

Proc @MatrixAufruf : Parameters n&

    Declare i&,j&,Max!,erg$,B![n&,n&],k!

    Whileloop n&:i&=&Loop

        Whileloop n&:j&=&Loop

            A![i&,j&]=(Rnd()-0.5)*1000000
            B![i&,j&]=A![i&,j&]

        Endwhile

    Endwhile

    print "TESTMATRIX:" : @Show A![],n&
    MatInv A![],N&
    print "INVERSE:" : @Show A![],n&
    MatInv A![],N&
    print "INVERSE RÜCKINVERTIERT:" : @Show A![],n&
    ' Fehlerbestimmung und -Ausgabe
    Max! = -1

    Whileloop n&:i&=&Loop

        Whileloop n&:j&=&Loop

            If Abs(A![i&,j&] - B![i&,j&]) > Max!

                Max! = Abs(A![i&,j&] - B![i&,j&])

            endif

            'erg$ = erg$ + Format$("%e",A![i&,j&] - B![i&,j&]) + " "

        endwhile

        'erg$ = erg$ + chr$(10)+chr$(13)

    endwhile

    erg$ = erg$ + "\n Größter Fehler: "+format$("%e",Max!)
    print "DIFFERENZ:" :print erg$
    waitinput
    Clear B![]

Endproc

' Eigentliche Inversion

Proc MatInv :parameters Mat![],N&

    Declare Hlp1&[n&],Hlp2&[n&],Hlp3&[n&]
    Declare Max!,T!,i&,j&,k&,ix&,iy&

    Whileloop n&:i&=&Loop

        Hlp3&[i&]=0

    Endwhile

    Whileloop n&:i&=&Loop

        ' Suche das größte Element
        Max! = 0

        Whileloop n&:j&=&Loop

            If Hlp3&[j&]<>1

                Whileloop n&:k&=&Loop

                    If (Hlp3&[k&]<>1) AND (Max! <= Abs(Mat![j&,k&]))

                        iy& = k&
                        ix& = j&
                        Max! = Abs(Mat![j&,k&])

                    EndIf

                endwhile

            EndIf

        endwhile

        Hlp3&[iy&] = Hlp3&[iy&] + 1
        'Pivotisierung

        If ix&<>iy&

            Whileloop n&:j&=&Loop

                t!=Mat![ix&,j&]
                Mat![ix&,j&]=Mat![iy&,j&]
                Mat![iy&,j&]=t!

            Endwhile

        EndIf

        Hlp1&[i&] = ix&
        Hlp2&[i&] = iy&
        'Kontrolle auf mögliches Verschwinden eines Hauptachsenwertes

        If Abs(Mat![iy&,iy&]) < 10^-300

            Print "Abbruch, Inversion nicht möglich!"
            Waitinput :End

        Else

            T! = Mat![iy&,iy&]
            Mat![iy&,iy&] = 1

            Whileloop n&:j&=&Loop

                Mat![iy&,j&] = Mat![iy&,j&] / T!

            EndWhile

            Whileloop n&:j&=&Loop

                If j&<>iy&

                    T! = Mat![j&,iy&]
                    Mat![j&,iy&] = 0

                    Whileloop n&:k&=&Loop

                        Mat![j&,k&] = Mat![j&,k&]- Mat![iy&,k&] * T!

                    endwhile

                EndIf

            endwhile

        EndIf

    endwhile

    'Rücktausch

    Whileloop n&:i&=&Loop

        j& = N& + 1 - i&

        If Hlp1&[j&]<>Hlp2&[j&]

            ix& = Hlp1&[j&]
            iy& = Hlp2&[j&]

            Whileloop n&:k&=&Loop

                T!=Mat![k&,ix&]
                Mat![k&,ix&]=Mat![k&,iy&]
                Mat![k&,iy&]=T!

            endwhile

        EndIf

    endwhile

    'Hilfsspeicher freigeben
    Clear Hlp1&[],Hlp2&[],Hlp3&[]
    'zur Ausgabe ...

EndProc

' Anzeigen der Matrix

Proc Show :parameters A![],n&

    declare i&,j&

    Whileloop n&:i&=&Loop

        Whileloop n&:j&=&Loop

            erg$ = erg$ + Format$("%e",A![i&,j&])+" "

        Endwhile

        erg$ = erg$+chr$(10)+chr$(13)

    Endwhile

    print erg$
    waitinput 1000
    erg$=""

EndProc

' Eine wichtige Anwendung der Matrizeninversion ist die
' Lösung von linearen Gleichungssystemen. Dieses
' Lösungsverfahren hat ja den immensen Vorteil, sofern man
' die invertierte Matrix kennt, besonders elegant zu sein.
' Denn wenn Ax = b das Gleichungssystem in vektorieller
' Schreibweise beschreibt, dann ist x = inv(A)*A*x = inv(A)*b
' bereits die Lösung. Deshalb:

Proc InvMat : parameters n&,a![],x![]

    Declare i%,j%
    MatInv a![],n&

    WhileLoop n&:i&=&Loop

        x!(i&)=0

        Whileloop n&:j&=&Loop

            x!(i&)=x!(i&)+a![i&,j&] * a![j&,n&+1]'<<< Rechte Seite d.LGS

        Endwhile

    Endwhile

Endproc