Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, N.Y. 13730 PR Pan-Reif May 21, 1986 Adam Fritz 133 Main Street Afton, New York 13730 BYTÅ magazinå foò April¬ 198¶ presentó aî articlå entitled¬ "Thå Inversioî oæ Largå Matrices¢ bù Thomaó E® Phipps¬ Jr.¬ REFERENCE [1]® Thió articlå presentó á BASIÃ prograí foò aî algorithí bù Victoò Paî anä Johî Reiæ tï inverô matrices® Dr. Phipps§ prograí amendó onå versioî oæ theiò algorithí tï applù 'quartic§ estimatioî tï updatå thå inverså estimate. Thió writeuð discusseó conversioî anä applicatioî oæ thaô program. CONVERSION ========== Thå BASIÃ prograí oæ Referencå [1Ý ió converteä tï Microsofô BASIC¬ Microsofô FORTRAN¬ TURBÏ Pascal¬ anä Ecosofô Ecï C® Seå DISTRIBUTION. Thå FORTRAN¬ Pascal¬ anä Ã versionó conform¬ loosely¬ tï thå LINPACË stylå followeä bù Referencå [2]®  Wherå  possible¬  loopó havå beeî replaceä bù correspondinç BLAÓ routinå calls. Thå conversionó uså á residuaì matriø ratheò thaî 'erroò' matriø witè suitablå adaptatioî oæ terminology. See DISCUSSION. Phipps useó á 'quartic§ modificatioî foò residuaì matriø extrapolatioî tï acceleratå convergence® Thió modificatioî alsï reduceó stabilitù becauså thå regioî oæ convergencå ió con tracted® Thå versionó distributeä witè thió writeuð uså firsô ordeò expressionó foò extrapolatioî factors® Seå DISCUSSION. DISCUSSION ========== Reference [1] uses an 'error§ matrix. Huh¿ Iî fact¬ residualó - noô erroró - arå beinç defineä witè respecô tï thå identitù matrix. Bù erroò wå wanô tï kno÷ thå differencå betweeî 'true§ anä computeä values. Thió distinctioî ió importanô becauså iô ió noô truå thaô smalleò residualó meaî smalleò errors® Second¬ wå wanô tï qualifù 'true§ valueó iî termó oæ theiò numer icaì representation. Thaô is¬ probleí Á ió numericallù encodeä A§ Š(1/³ <¾ .33...3© anä somå algorithí leadó tï solutioî B§ whicè maù oò maù noô be¬ iî somå sense¬ best. Whaô doeó B§ havå tï dï witè Â - thå 'true§ solution¿ Foò physicaì problems¬ aó opposeä tï mathematicaì problems¬ onlù A§ maù bå knowable® Seå Noteó (3© and (4). Paragrapè fouò says that:    "Thå probleí isº Giveî A¬ whaô ió A^-1? Iæ thå dimensioî î oæ    Á ió small¬ thió probleí ió easy.¢ É wonder. Supposå wå trù á particularlù simplå ³ ø 3¬ viz.; _ _ | | | 1 1/2 1/3 | | 1/2 1/3 1/4 | | 1/3 1/4 1/5 | |_ _| wherå a[i,jÝ ½ 1/(i+j-1)® Thå followinç tablå waó produceä bù Phipps' program; ITER. E(1,1) E(1,N) B(1,1) B(1,N) 0 .595041 -.156198 .297521 .0991736 1 .213926 -.300395 .590827 .17486 2 .481107 -.183538 .459955 .0225266 3 .0702667 -.293482 .991716 -.182955 4 .134694 -.111926 1.68666 -1.18852 5 .0223539 .0580739 2.78269 -2.49344 6 .0519888 .0383779 2.63106 -2.21623 7 .044227 .0421093 2.73269 -1.93643 8 .0476556 .0382407 2.91608 -.852009 9 .0394996 .0354136 3.60656 2.5629 10 .0284853 .0234768 5.30509 11.2484 11 8.97074E-03 7.82442E-03 7.79567 23.8781 12 6.19888E-04 5.23567E-04 8.91858 29.5865 13 3.8147E-06 1.90735E-06 8.99976 29.9988 14 0 0 8.99999 29.9999 15 0 0 8.99999 29.9999 Run number 1 Duration = No. iter.= 15 Matrix dim.= 3 Max. error= .000003 Another iteration (y,n)? n The maximum error was calculated by multiplying A^-1*A. As a check on the true error in the approximated inverse matrix the maximum error in A*A^-1 will now be determined. The max. error found this way is 4.69208E-04 Ok That wasn't too bad. Now try a similar 4 by 4; ITER. E(1,1) E(1,N) B(1,1) B(1,N) 0 .672 -.101669 .2304 .0576 Š 1 .253518 -.227796 .540098 .115097 2 .559764 -.113231 .411087 -.0264034 3 -.0250019 -.225365 1.1246 -.233489 4 -.165854 -.0644244 2.11955 -1.12598 5 .2051 .206097 2.56083 -1.8227 6 .319758 .399902 3.22895 -2.37613 7 18.6404 26.701 77.4394 -108.269 Overflow ... 8 6.89423E+16 9.91363E+16 2.76255E+17 -4.10944E+17 Overflow ... 9 1.41784E+37 4.41557E+37 1.70141E+38 -1.70141E+38 Overflow ... Oops¡ Maybå ³ ió smalì anä ´ ió transfinite. Yoõ know¬ Hottentoô ordinaì numbers¬ viz.¬ 1¬ 2¬ 3¬ Aleph0. Thió example¬ oæ course¬ ió Hilbert'ó matrix. Iô haó brokeî morå thaî onå algo rithm. Iô ió easù tï remembeò and¬ witè somå attention¬ singlå precisioî computeò algorithmó caî geô gooä resultó foò orderó uð tï anä includinç 5¡ Doublå precisioî doesn'ô helð thå author'ó program® Thå tablå belo÷ showó á Hilberô ´ bù ´ ruî througè thå author'ó prograí witè everythinç defineä doublå anä tabulateä figureó truncateä tï ninå digits; ITER. E(1,1) E(1,N) B(1,1) B(1,N) 0 .672000001 -.101668571 .230399997 .057599999 1 .253517594 -.227796094 .540097556 .115097049 2 .559763986 -.113231208 .411087080 -.026403395 3 -.025001561 .225365214 1.12459632 -.233489548 4 -.165853433 -.064424267 2.11955069 -1.12597738 5 .205099223 .206096680 2.56082681 -1.82269827 6 .319755948 .399899500 3.22894220 -2.37612685 7 18.6400009 26.7004894 77.4380076 -108.266784 Overflow ... 8 6.8931D+16 9.9121D+16 2.7621D+17 -4.1088D+17 Overflow ... 9 1.4178D+37 4.4155D+37 1.7014D+38 -1.7014D+38 Overflow ... Thå probleí ió Phipps' 'quartic§ modification® Hå modifieä thå algorithí tï acceleratå convergence® Stabilitù waó sacrificeä becauså thå convergencå regioî waó contracted® Iæ thå prograí ió modifieä tï uså Paî anä Reif'ó unarù extrapolatioî formulation; 1+e[i,i] insteaä oæ Phipps§ quartiã; 1+(1+(1+(1+e[i,i])*e[i,i])*e[i,i])*e[i,i] Š theî foò Hilberô orderó ´ anä µ thå singlå precisioî versioî produceó thå followinç tabulations; Order 4 ------- ITER. E(1,1) E(1,N) B(1,1) B(1,N) 0 .672 -.101669 .2304 .0576 1 .406052 -.181823 .427215 .0941397 2 .363117 -.185065 .501535 .0567919 ... 25 7.62939E-06 3.8147E-06 16.0004 -140.01 ... 51 7.62939E-06 3.8147E-06 16.0004 -140.01 STUCK! Order 5 ------- ITER. E(1,1) E(1,N) B(1,1) B(1,N) 0 .719271 -.0737348 .191806 .0383611 1 .461549 -.139085 .377856 .0652478 2 .388011 -.147169 .475874 .0353428 ... 31 -1.52588E-05 -7.62939E-06 25.0139 630.75 ... 51 6.10352E-05 2.28882E-05 24.9914 629.529 STUCK! Thå caseó sho÷ stability¬ buô thå convergencå criterion¬ viz.¬ 3e-6¬ ió noô achievablå witè singlå precision.. TURBÏ Pascaì useó á fivå bytå floatinç poinô encodinç foò versionó withouô numeriã coprocessors® Therå arå 4° significanô bits® Cf. MANTISSA.PAS® Thå quartiã formulatioî is¬ oæ course¬ stilì unstable® (Iô waó unstablå foò doublå precisioî BASIC¬ 5¶ bits.) Thå unarù formulatioî produceó thå followinç tabulatioî foò ordeò 6; Pan-Reif Test Program. Dimension = 6 Iterations = 100 Residual = 1.e-6 Original Matrix (by column): Š 1.000000E+00 5.000000E-01 3.333333E-01 2.500000E-01 2.000000E-01 5.000000E-01 3.333333E-01 2.500000E-01 2.000000E-01 1.666667E-01 3.333333E-01 2.500000E-01 2.000000E-01 1.666667E-01 1.428571E-01 2.500000E-01 2.000000E-01 1.666667E-01 1.428571E-01 1.250000E-01 2.000000E-01 1.666667E-01 1.428571E-01 1.250000E-01 1.111111E-01 1.666667E-01 1.428571E-01 1.250000E-01 1.111111E-01 1.000000E-01 1.666667E-01 1.428571E-01 1.250000E-01 1.111111E-01 1.000000E-01 9.090909E-02 Iterations: 44 Inverse (by column): 3.600002E+01 -6.30000E+02 3.360003E+03 -7.56001E+03 7.560007E+03 -6.30000E+02 1.470001E+04 -8.82001E+04 2.116802E+05 -2.20500E+05 3.360000E+03 -8.82000E+04 5.644800E+05 -1.41120E+06 1.512000E+06 -7.56000E+03 2.116800E+05 -1.41120E+06 3.628800E+06 -3.96900E+06 7.560002E+03 -2.20500E+05 1.512000E+06 -3.96900E+06 4.410001E+06 -2.77200E+03 8.316003E+04 -5.82120E+05 1.552320E+06 -1.74636E+06 -2.77200E+03 8.316007E+04 -5.82120E+05 1.552320E+06 -1.74636E+06 6.985442E+05 Residuals (by column): -1.08979E-05 1.642644E-03 -1.64233E-02 5.220413E-02 -6.53734E-02 9.408221E-06 7.874668E-04 -9.55081E-03 3.239799E-02 -4.19393E-02 1.658825E-05 4.084706E-04 -6.29961E-03 2.271676E-02 -3.02548E-02 1.904229E-05 2.159923E-04 -4.51177E-03 1.720142E-02 -2.34582E-02 1.960201E-05 1.085103E-04 -3.42065E-03 1.371288E-02 -1.90763E-02 1.934567E-05 4.445761E-05 -2.70283E-03 1.133776E-02 -1.60398E-02 2.813005E-02 1.841569E-02 1.350582E-02 1.061213E-02 8.723021E-03 7.399201E-03 End of Test. Á routinå nameä MANTISSA.ª ió useä tï estimatå thå numbeò oæ significanô bits® Foò BASIÃ anä FORTRAÎ thå resulô ió 2µ bits® Foò TURBÏ Pascaì thå resulô ió 4° bits® Foò Ecï Ã thå resulô ió 2µ bits® Notå thaô Ecosoft'ó compileò detecteä thå constanô Šfoldinç opportunitù oò kepô MANTISSA'ó intermediatå resultó iî á floatinç poinô pseudï registeò anä claimeä 5· bits® Thió probablù meanó thaô theiò arithmetiã ió donå doublå precision® Defeat inç their optimizeò produces; Significant Bits: 25. Belo÷ arå listeä runtimeó computeä froí [2]¬ [3]¬ anä [5Ý foò inversioî oæ á 10° bù 10° singlå precisioî matriø usinç onå calì eacè tï LINPACK'ó SGEFÁ anä SGEDÉ writteî iî FORTRAN¬ where: SGEFÁ  - Singlå precisioî GEneraì matriø FActorization;  SGEDÉ  - Singlå  precisioî  GEneraì Determinanô anä Inverse»  anä cBLAÓ  - codeä (assembly) Basiã Lineaò Algebrá Subprograms. Machine Software Time (seconds) ------- -------- -------------- IBM PC Microsoft 3.1 3568. SBC-200 4MHz Z80 Microsoft 3.2 1785. CO16 8MHz/8087 Microsoft 3.2 138. Apollo DN460/660 AEGIS 8.0 (cBLAS) 16.2 VAX 11/780 FPA VMS 4.1 (cBLAS) 5.88 NEC SX-2 F77 (cBLAS) .048 Runtimeó scalå witè cubå oæ dimension® Á PÃ runninç Micro sofô FORTRAN¬ invertinç á 6° ø 6° matrix¬ anä usinç LINPACË algorithmó woulä requirå 77± seconds. Therå woulä bå adequatå timå betweeî 77± secondó anä 15-2´ houró lateò tï consideò thå author'ó conclusion;    "Newton'ó iteratioî makeó thå besô oæ á baä joâ anä doeó aó    welì aó anù methoä coulä do." DISTRIBUTION ============ PR.DOC This writeup and discussion PR10.BAS Reference 1, text format PR11.BAS Microsoft BASIC, ASCII format PR11/T.LBR TURBO Pascal (tm) version PRDR.PAS Pan-Reif driver HILGEN.PAS Hilbert test matrix generator RANGEN.PAS Random test matrix generator PR.PAS Pan-Reif routine BLAS.PAS Basic linear algebra subroutines [2] OUT.PAS Vector and matrix output routine [2] MANTISSA.PAS Compute number significant bits [2] PR11/F.LBR Microsoft FORTRAN version. PRDR.FOR Pan-Reif driver HILGEN.FOR Hilbert test matrix generator RANGEN.FOR Random test matrix generator Š PR.FOR Pan-Reif routine BLAS.FOR Basic linear algebra subroutines [2] OUT.FOR Vector and matrix output routine [2] MANTISSA.PAS Compute number significant bits [2] PR11/C.LBR Ecosoft Eco-C (tm) version PRDR.C Pan-Reif driver HILGEN.C Hilbert test matrix generator RANGEN.C Random test matrix generator PR.C Pan-Reif routine BLAS.C Basic linear algebra subroutines [2] OUT.C Vector and matrix output routine [2] ABS.C Float absolute value routine MANTISSA.C Compute number significant bits [2] NOTES ===== (1) Everyone¬ excepô MACSYMÁ freaks¬ knowó thaô matriø compu        tationó involvå floatinç poinô numberó and¬ therefore¬        non-commutativå arithmetic® Á carefullù crafteä sequen        tiaì algorithí wilì alwayó beaô á paralleì algorithí iî        accuracù and¬ sometimes¬ speed® Amdahì haó argueä foò        yearó thaô anù algorithm'ó sequentiaì componenô limitó        speeä performancå gaiî availablå froí 'paralyzing'. (2) Thå thirä paragrapè makeó á claií thaô inversioî ió RAÍ        intensive® Maybe® Referencå [3Ý reportó á familù oæ        algorithmó foò dealinç witè determinants¬ lineaò equa        tions¬ anä inverses®        LINPACË routineó depend¬ mostly¬ oî LÕ factorizatioî baseä        oî Gaussiaî eliminatioî witè partiaì pivoting® Operationó        arå organizeä sï thaô foò languageó likå FORTRAÎ thaô        storå matriceó bù column¬ successivå operandó arå storeä        contiguously® Whaô thió meanó ió thaô onlù onå oò twï        columnó arå accesseä concurrently® Reductioî iî pagå        thrash-inç foò virtuaì storagå machineó ió obviouó anä        reaì RAÍ memorù ió similarlù reduced® Pascal¬ TURBÏ Pas        caì anyway¬ storeó matriceó bù row® LINPACË routineó arå        easilù adapteä - somå oæ theí - tï microcomputeò environ        ments. [2]    (3) Thå author'ó prograí defineó residualó witè respecô tï thå        matriø producô BA® Thå solutioî foò Hilberô ordeò ³ usinç        Eco-Ã is;         Inverså (bù column):         8.999969E° -3.599984E± 2.999984E±         -3.599989E± 1.919994E² -1.799994E²         2.999982E± -1.799990E² 1.799991E² Š        Residualsº (bù column):         3.576279E-· 0.000000E° 1.966958E-·         -3.298129E-¶ 0.000000E° -1.752378E-¶         -5.165643E-· 0.000000E° 2.503395E-¶        Residualó foò matriø producô AÂ computeä foò thió solutioî        are;         Residualsº (bù column):         3.278255E-µ -1.902183E-´ 1.825889E-´         2.241135E-µ -1.296997E-´ 1.258850E-´         1.717210E-µ -9.940863E-µ 9.787083E-µ        Notå thaô thå residualó arå consistentlù differenô bù aó        mucè aó twï tï threå orderó oæ magnitude® Iæ thå author'ó        prograí ió changeä tï minimizå residualó foò producô AÂ        (bù reorderinç botè residuaì anä ne÷ inverså computation©        thå solutioî is;         Inverså (bù column):         8.999982E° -3.599989E± 2.999986E±         -3.599991E± 1.919994E² -1.799993E²         2.999992E± -1.799994E² 1.799993E²         Residualsº (bù column):         5.722046E-¶ 5.245209E-¶ 4.583598E-¶         -2.618632E-µ -2.670288E-µ -2.311468E-µ         1.601379E-µ 1.907349E-µ 1.698732E-µ        wherå thå residualó arå computeä oî BA® Thå solutionó arå        similaò buô arå noô identical® Iæ thå residualó arå com        puteä foò AÂ witè solutioî oî AÂ thå residualó are;         Residualsº (bù column):         -5.960464E-· -3.298130E-¶ 7.549993E-·         -4.768372E-· -3.814697E-¶ 0.000000E°         -3.755090E-· -3.278257E-¶ 1.728535E-¶        Notå thaô residualó oæ AÂ anä BÁ foò solutioî oî AÂ arå oî        thå  samå  ordeò  oæ magnitude®  Thå poinô waó  madå  thaô        residualó arå noô errors® Clearly¬ residualó arå noô resi        duals, either! (4) Thå authoò reportó that¬ "Anù real-numbeò matriø caî bå        expresseä iî thió fractional-elemenô forí bù dividinç alì        elements¬ etc.¢ Sucè aî operatioî transformó probleí A§        intï A'§ becauså roundofæ erroò appearó unlesó purelù        exponentiaì scalinç ió applied® Thió probleí haó solutioî        B'§ whoså relatioî tï B§ anä Â theî becomeó aî issue. Š   Therå ió á furtheò probleí dealinç witè erroò becauså run        timeó oæ tenó oæ houró arå reported® Iô woulä bå wiså tï        perforí somå posô ruî validitù checë tï sho÷ thaô thå        probleí solveä waó thaô formulated® Thió algorithm'ó        convergencå maù havå tï dï witè RAÍ biô erroò rate® Onå        ió remindeä oæ aî ancienô storù abouô carryinç bitó iî á        sieve. (5) Thå articlå discusseó iteration¬ initiaì estimates¬        refinement¬ anä prograí outline® Speakinç oæ refinement¬        referencå [4Ý applieó iteratioî iî tandeí witè modifieä        Gaussiaî elimination® Barrodalå anä Stuarô introducå thå        strikinç innovatioî oæ applyinç thå concepô oæ residuaì tï        Gaussiaî eliminatioî foò solutioî oæ lineaò systems®        MODGÅ selectó rowó foò pivotinç baseä oî righô hanä sidå        'residuals§ anä theî ro÷ pivotó oî maximuí absolutå        coefficienô vectoò elemenô value¬ pluó somå detailó foò        thå reader® MODGÅ ió followeä bù REFINÅ whicè applieó        iteration® Iæ yoõ like¬ MODGÅ constructó thå initiaì        estimatå anä REFINÅ iterateó tï reducå residuals® Thå        poinô is¬ iteration¬ Newtoniaî anä otherwiså haó beeî iî        uså foò á lonç timå foò thió clasó oæ problems® (6) REFINE¬ likå thå author'ó program¬ doesn'ô quiô whilå it'ó        ahead® Wheî residualó begiî tï grow¬ iteratioî shoulä        terminate® Notå thaô maximuí absolutå residuaì valuå doeó        noô changå monotonicallù witè iteration® Thå exiô        criterioî woulä havå tï bå somå norí defineä oî thå        residual matrix.    (7) Thå authoò conjectureó thaô "accuracù ..® [isÝ fairlù        independenô oæ n.¢ Thió ió onlù partlù true® Thå authoò        generateó randoí vectoró uniformlù distributeä iî á hyper        cubå oæ dimensioî î witè sidå 2® Singularitù ariseó froí        colinearitù amonç twï oò morå vectors®   Iî effect¬ onå ió asking» whaô ió thå probabilitù thaô twï        oò morå vectoró havå theiò uniô vectoró iî á penciì thaô        ió computationallù singular® Uniformlù distributeä        matriceó arå usefuì foò speeä benchmarking® Foò measure        menô oæ accuracù anä convergencå morå contriveä caseó arå        necessary.    (8) The author claims;        "Aî ilì conditioneä matriø cannoô normallù bå spotteä bù        prioò inspection."   Iæ thió means» caî someonå starå aô á matriø anä saù        whetheò hió prograí wilì takå á hike¿ Probablù not® Johî        Voî Neumann¬ maybe® Thió writer¬ no® I'lì beô Dongarrá        anä thå resô oæ thå lineaò systemó mafiá arå annoyeä thaô        theù spenô alì thaô taxpayeò moneù anä peoplå stilì don'ô        believå yoõ caî avoiä ruî timå erroò messages® Thå pur        poså oæ LINPACË routineó suffixeä witè CÏ ió tï computå á Š       COnditioî estimatå foò á matrix® Foò example¬ SGECÏ com        puteó thå conditioî - thå reciprocaì condition¬ really¬        nameä RCONÄ - foò Singlå precision¬ GEneraì matrices.   Foò example¬ usinç DIDÒ froí [2]¬ RCONÄ anä inverses¬ bù        column¬ foò Hilberô orderó 3¬ 4¬ anä µ arå showî below®        Effectivå RCONÄ foò ordeò ¶ vanishes. Ordeò ³ RCONDº .14688422E-0² ------- 8.99999E+00 -3.59999E+01 2.99999E+01 -3.60000E+01 1.92000E+02 -1.80000E+02 3.00000E+01 -1.80000E+02 1.80000E+02 Order 4 RCOND: .46460722E-04 ------- 1.60000E+01 -1.20000E+02 2.40000E+02 -1.40000E+02 -1.20000E+02 1.20000E+03 -2.70000E+03 1.68000E+03 2.40000E+02 -2.70000E+03 6.48002E+03 -4.20001E+03 -1.40000E+02 1.68000E+03 -4.20001E+03 2.80001E+03 Order 5 RCOND: .14412633E-05 ------- 2.49957E+01 -2.99925E+02 1.04968E+03 -1.39953E+03 6.29770E+02 -2.99917E+02 4.79852E+03 -1.88937E+04 2.68707E+04 -1.25955E+04 1.04963E+03 -1.88933E+04 7.93519E+04 -1.17558E+05 5.66797E+04 -1.39942E+03 2.68696E+04 -1.17556E+05 1.79135E+05 -8.81684E+04 6.29710E+02 -1.25948E+04 5.66782E+04 -8.81675E+04 4.40842E+04 REFERENCE ========= [1] The Inversion of Large Matrices Thomas E. Phipps, Jr. BYTE Magazine Volume , Number 4 April, 1986, Pages 181-188 [2] SIG/M Whetstone Benchmarks and LINPACK Conversions Adam Fritz Volume 244 September 20, 1985 [3] LINPACK User's Guide J.J. Dongarra, C.B. Moler, J.R. Bunch, & G.W. Stewart Society for Industrial and Applied Mathematics Fourth Printing, 1984 [4] Collected Algorithms from ACM Algorithm 576 Š A FORTRAN Program for Solving Ax = b I. Barrodale & G.F. Stuart ACM Transactions on Mathematical Software Volume 7, Number 3 September, 1981 Pages 391-397 [5] Performance of Various Computers Using Standard Linear Equations Software in A FORTRAN Environment Jack J. Dongarra Technical Memorandum No. 23 October 11, 1985 Mathematics and Computer Science Division Argonne National Laboratory NOTICES ======= TURBO Pascal (tm) Borland International Eco-C (tm) Ecosoft, Inc. CP/M (tm) Digital Research, Inc. MS-DOS (tm) Microsoft, Inc. Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, N.Y. 13730