XProfan

Italia

p.specht

Der Householder-Algorithmus
=====================
Bei kritischen Datenkonstellationen versagen übliche Lösungsverfahren per Gleichungssysteme. Das Gauss'sche Eliminationsverfahren etwa produziert bei grande Werte-Spannweiten innerhalb der Matrizenspalten erhebliche Fehler, die zu absoluter Unbrauchbarkeit der Ergebnisse führen können, ohne dass der Anwender es gleich merken würde. Für solche Fälle hatte der US-Amerikanische Mathematiker Alston Householder 1958 die Idee, die Nicht-Diagonalwerte durch geschickt gewählte Spiegelungen zu Null zu machen. Besonders effizient gelingt das bei Diagonalsymmetrischen Matrizen, wie sie in der Physik oft vorkommen. (Ansonsten kann man jede beliebige Matrix durch Multiplikation mit ihrer Transponierten zu einer Diagonalsymmetrischen Matrix machen und das Resultat dann LU-Zerlegen).

In dieser Demo per private Bildungszwecke wurde Fortran's Arrayindex-Basis '1' beibehalten! In XProfan eine ziemliche Stolperfalle!

Kompilieren Markieren Separieren

WindowTitle "Inversion einer rellen, diagonalsymetrischen Matrix mittels Householder-Algorithmus"
'Lösung per multidimensionaler Spiegelung, Algorithmus von Alston Scott Householder, USA 1958
'Orig.F90-Source: https://jean-pierre.moreau.pagesperso-orange.fr/Fortran/householder_f90.txt
'(D)Demo-Migration nach XProfan11.2a 2014-10 by P.Specht, Wien;  OHNE JEDWEDE GEWÄHR!
'Verfahrensdetails siehe https://de.wikipedia.org/wiki/Householdertransformation
'***************************************************************
'{  Inversion of a real square matrix by Householder's method  *
'
'* SAMPLE RUN [im Original nur single precision, Anm. P.Specht]*
'* 4                                                           *
'* 1.0       0.5        0.333333  0.25                         *
'* 0.5       0.333333   0.25      0.2                          *
'* 0.333333  0.25       0.2       0.166667                     *
'* 0.25      0.2        0.166667  0.142857                     *
'
'* Output data:
'* Input square matrix
'*  0.100000E+01   0.500000E+00   0.333333E+00   0.250000E+00  *
'*  0.500000E+00   0.333333E+00   0.250000E+00   0.200000E+00  *
'*  0.333333E+00   0.250000E+00   0.200000E+00   0.166667E+00  *
'*  0.250000E+00   0.200000E+00   0.166667E+00   0.142857E+00  *
'*                                                             *
'* Transformation # 1                                          *
'* -0.119315E+01  -0.670493E+00  -0.474932E+00  -0.369836E+00  *
'*  0.000000E+00   0.664812E-01   0.657297E-01   0.586884E-01  *
'*  0.000000E+00   0.720990E-01   0.771533E-01   0.724594E-01  *
'*  0.000000E+00   0.665741E-01   0.745319E-01   0.722012E-01  *
'* Transformation # 2                                          *
'* -0.119315E+01  -0.670493E+00  -0.474932E+00  -0.369836E+00  *
'*  0.000000E+00  -0.118533E+00  -0.125656E+00  -0.117542E+00  *
'*  0.000000E+00   0.000000E+00   0.257162E-02   0.378339E-02  *
'*  0.000000E+00   0.000000E+00   0.566533E-02   0.878781E-02  *
'* Transformation # 3                                          *
'* -0.119315E+01  -0.670493E+00  -0.474932E+00  -0.369836E+00  *
'*  0.000000E+00  -0.118533E+00  -0.125656E+00  -0.117542E+00  *
'*  0.000000E+00   0.000000E+00  -0.622167E-02  -0.956580E-02  *
'*  0.000000E+00   0.000000E+00   0.000000E+00   0.187201E-03  *
'
'* Inverted matrix:
'*  0.160327E+02  -0.120367E+03   0.240887E+03  -0.140579E+03  *
'* -0.120367E+03   0.120414E+04  -0.271001E+04   0.168653E+04  *
'*  0.240887E+03  -0.271001E+04   0.650422E+04  -0.421582E+04  *
'* -0.140579E+03   0.168653E+04  -0.421582E+04   0.281033E+04  *
'
'* Determinante: -0.164722E-06
'
'* AP=A*A^-1 matrix:                                           *
'*  0.100000E+01   0.341061E-12  -0.682121E-12   0.341061E-12  *
'*  0.710543E-14   0.100000E+01   0.000000E+00   0.000000E+00  *
'* -0.710543E-14   0.113687E-12   0.100000E+01   0.113687E-12  *
'* -0.355271E-14   0.284217E-13   0.000000E+00   0.100000E+01  *
'* ----------------------------------------------------------- *
'* Ref.: "Méthodes de calcul numérique - Tome 2 - By Claude    *
'*        Nowakowski, PSI Editions, 1984".                     *
'*                  Fortran 90 Version By J-P Moreau, Paris.   *
'*                            (www.jpmoreau.fr)                *
'*       XProfan11.2a Migration by P.Specht, Vienna            *
'}
WindowStyle 24
Window 0,0-%maxx,%maxy-40
Font 2
set("decimals",17)
Declare A![25,50]
Print "\n\n               DER HOUSHOLDER-ALGORITHMUS\n               =========================="
Householder(0,A![])' Bei Angabe n = 0 Dimensionen wird das Testbeispiel berechnet
Waitinput
END

Proc Householder :parameters n&,A![]

    declare A0![25,25], V![25],D!,B![25],AI![25,25],AP![25,25],X![25],i&,j&,k&

    if n&=0'Read matrix to be inverted

        n&=4
        A![1,1]=1.0 :A![1,2]=0.5   :A![1,3]=0.333333333333333333:A![1,4]=0.25
        A![2,1]=0.5 :A![2,2]=0.333333333333333333 :A![2,3]=0.25 :A![2,4]=0.2
        A![3,1]=0.333333333333333333:A![3,2]=0.25 :A![3,3]=0.2  :A![3,4]=0.166666666666666666
        A![4,1]=0.25 :A![4,2]=0.2  :A![4,3]=0.166666666666666666:A![4,4]=0.142857142857142857

    endif

    Print "\n ist per diagonalsymetrische Ausgangsmatrizen vorgesehen,   BEISPIEL:\n"

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            If j& = n&

                print tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&])

            Else

                print tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&]),

            endif

        endwhile

    endwhile

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            A![i&,j&+n&] = 0
            A0![i&,j&] = A![i&,j&]

        endwhile

        A![i&,i&+n&] = 1

    endwhile

    'Transform A into triangular matrix
    S1000(n&,A![],V![])
    'N linear systems to solve

    whileloop 1,n&:k&=&Loop

        whileloop 1,n&:i&=&Loop

            B![i&] = A![i&,K&+n&]

        endwhile

        'Solve triangular system
        S2000(n&,A![],B![],X![])

        whileloop 1,n&:i&=&Loop

            AI![i&,K&]=X![i&]

        endwhile

    endwhile'k loop

    print "\n\n Invertierte Matrix:",
    ' Inverted matrix:'

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            If j&=n&

                print tab(j&*27);if(round(AI![i&,j&],8)<0,""," ");format$("%g",round(AI![i&,j&],8))

            Else

                print tab(j&*27);if(round(AI![i&,j&],8)<0,""," ");format$("%g",round(AI![i&,j&],8)),

            endif

        endwhile

    endwhile

    'Calculate determinant
    S4000(n&,A![],D!)
    print "\n\n Determinante = ";format$("%g",D!)
    'Check AP=A*A^-1=Unity Matrix
    print "\n\n A*A^-1=AP=Einheismatrix:",
    ' AP=A*A^-1 matrix:'
    S3000(n&,A0![],AI![],AP![])

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            If j&=n&

                print tab(j&*27);if(round(AP![i&,j&],12)<0,""," ");format$("%g",round(AP![i&,j&],12))

            Else

                print tab(j&*27);if(round(AP![i&,j&],12)<0,""," ");format$("%g",round(AP![i&,j&],12)),

            endif

        endwhile

    endwhile

EndProc'of main program

Proc S1000:parameters n&,A![],V![]

    declare S!,XA!,XB!,XG!,XL!
    print

    whileloop 1,n&-1:k&=&Loop

        S! = 0

        whileloop k&,n&:i&=&Loop

            S!=S!+A![i&,K&]*A![i&,K&]

        endwhile

        XL!= ( (A![K&,K&]<0)-(A![K&,K&]>=0) ) *sqrt(S!)
        XA!= XL!*(XL!-A![K&,K&])
        V![K&]=A![K&,K&]-XL!

        whileloop k&+1,n&:i&=&Loop

            V![i&]=A![i&,K&]

        endwhile

        whileloop k&+1,2*n&:j&=&Loop

            XB!=0

            whileloop k&,n&:i&=&Loop

                XB!=XB!+V![i&]*A![i&,j&]

            endwhile

            XG!=XB!/XA!

            whileloop k&,n&:i&=&Loop

                A![i&,j&]=A![i&,j&]-XG!*V![i&]

            endwhile

        endwhile

        A![K&,K&]=XL!

        whileloop k&+1,n&:i&=&Loop

            A![i&,K&]=0

        endwhile

        print "\n ";format$("0",K&);". Transformation:",

        whileloop 1,n&:i&=&Loop

            whileloop 1,n&:j&=&Loop

                If j&=n&

                    print tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&])

                Else

                    print tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&]),

                endif

            endwhile

        endwhile

    endwhile'k loop

    return

Endproc

Proc S2000 :parameters n&,A![],B![],X![]

    declare i&,j&,S!

    whileloop n&,1,-1:i&=&Loop

        j&=n&
        S!=0
        g20:
        case i&=j&:goto "g40"
        S!=S!-A![i&,j&]*X![j&]
        j&=j&-1
        goto "g20"
        g40:
        X![i&] = (B![i&] + S!) / A![i&,i&]

    endwhile

    return

Endproc

Proc S3000 :parameters n&,A0![],AI![],AP![]

    declare S!

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            S!=0

            whileloop 1,n&:k&=&Loop

                S!=S!+A0![i&,k&]*AI![k&,j&]

            endwhile

            AP![i&,j&]=S!

        endwhile

    endwhile

    return

Endproc

Proc S4000 :parameters n&,A![],D!

    D! = 1

    whileloop 1,n&:i&=&Loop

        D!=D!*A![i&,i&]

    endwhile

    return

Endproc