XProfan

English

p.specht

The Householder-Algorithmus
=====================
with censorial Datenkonstellationen failure übliche Lösungsverfahren for Gleichungssysteme. the Gauss'sche Eliminationsverfahren about produced with large values-Spannweiten inside the Matrizenspalten vastly Error, The To absolute Unbrauchbarkeit the Results lead can, without that the users it same remember would. for such Cases having the US-Amerikanische Mathematiker Alston Householder 1958 The idea, The not-Diagonalwerte through skillful chosen Spiegelungen To zero To make. particularly efficient succeed the with Diagonalsymmetrischen Matrizen, How tappt im dunkeln in the Physics often vorkommen. (otherwise can each arbitrary Matrix through Multiplikation with yours Transponierten to a Diagonalsymmetrischen Matrix make and the result then LU-decompose).

in this demonstration for private Bildungszwecke watts Fortran's Arrayindex-Base '1' to maintain! In XProfan a ziemliche Stolperfalle!

Compile Mark Separation

Window Title "Inversion of/ one rellen, diagonalsymetrischen Matrix through Householder-Algorithmus"
'Solution by multidimensionaler reflection, Algorithmus of Alston Scott Householder, USA 1958
'Orig.F90-Source: https://jean-pierre.moreau.pagesperso-orange.fr/Fortran/householder_f90.txt
'(D)demonstration-Migration to XProfan11.2a 2014-10 by P.woodpecker, Wien;  OHNE JEDWEDE GEWÄHR!
'Verfahrensdetails see https://de.wikipedia.org/wiki/Householdertransformation
'***************************************************************
'{  Inversion of a real square matrix by Householder's method  *
'
'* SAMPLE RUN [in the Original only single precision, Anm. P.woodpecker]*
'* 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.woodpecker, Vienna            *
'}
Window Style 24
Window 0,0-%maxx,%maxy-40
Font 2
set("decimals",17)
Declare A![25,50]
Print "\n\n               THE HOUSHOLDER-ALGORITHMUS\n               =========================="
Householder(0,A![])' with indicated n = 0 Dimensionen becomes the Testbeispiel accounts
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 is for 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,""," ");stature$("%g",A![i&,j&])

            Else

                print tab(j&*27);if(A![i&,j&]<0,""," ");stature$("%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,""," ");stature$("%g",round(AI![i&,j&],8))

            Else

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

            endif

        endwhile

    endwhile

    'Calculate determinant
    s4000(n&,A![],D!)
    print "\n\n Determinante = ";stature$("%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,""," ");stature$("%g",round(AP![i&,j&],12))

            Else

                print tab(j&*27);if(round(AP![i&,j&],12)<0,""," ");stature$("%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 ";stature$("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,""," ");stature$("%g",A![i&,j&])

                Else

                    print tab(j&*27);if(A![i&,j&]<0,""," ");stature$("%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