The C Users' Group Library 1994 August/wc-cdrom-cusersgrouplibrary-1994-08.iso/vol_100/188_01/matfac.doc

Text File | 1987-09-30 | 10KB | 182 lines

.lh 9
.cw 11
             Full Matrix Factorization Functions

     Matri°á factorizatioεá i≤á thσ generiπ terφá fo≥á method≤ ì
whicΦ spli⌠ ß matri° int∩ ß produc⌠ oµ tw∩ o≥ morσ factor≤á t∩ ì
facilitatσá solutioε oµ problem≤ sucΦ a≤ solutioε oµá simultaì
neou≤ linea≥ algebraiπ equations«á  Gaussiaε elimination¼ likσ ì
thσ method≤ showε here¼ factor≤ ß matri° int∩ ß produc⌠ oµ tw∩ ì
triangula≥á matrices¼á bu⌠á onσ oµ thσ triangula≥ matrice≤á i≤ ì
discardeΣ immediately«á  Factorizatioε method≤ arσ alway≤ morσ ì
efficien⌠ thaε truσ inversioε methods╗ wheneve≥ thσ produc⌠ oµ ì
aε inversσ witΦ anothe≥ matri° i≤ desired¼á i⌠ caε bσ replaceΣ ì
b∙á thσá forwarΣ anΣ bacδ substitutioεá usinτá thσá triangula≥ ì
factors.

     Amonτá thσ operation≤ possiblσ bu⌠ no⌠ showε arσá forminτ ì
thσá inversσ b∙ operatinτ oε aε identit∙ matrix¼á anΣ checkinτ ì
thσá factorizatioεá b∙á operatinτ oε ß cop∙á oµá thσá origina∞ ì
matrix¼á t∩ seσ iµ aε identit∙ matri° results«á  Thσ usua∞ usσ ì
i≤á operatioεá oε ß singlσ righ⌠ hanΣ sidσ vecto≥ t∩á solvσá ß ì
singlσ systeφ oµ equations«á  Iε man∙ applications¼á iteratioε ì
i≤á used¼á wherσ ne≈ righ⌠ hanΣ sidσ vector≤ ma∙ bσá developeΣ ì
later¼ fo≥ solutioε witΦ thσ factor≤ obtaineΣ previously.

     Thσá methoΣá showε iε EQD.├ i≤ Crou⌠á factorizatioεá witΦ ì
scaleΣá partia∞á pivoting«á  Equivalen⌠ method≤ arσá showεá iε ì
severa∞á textbooks¼á s∩ wσ wil∞ concentratσ oεá implementatioε ì
question≤ whicΦ arσ harde≥ t∩ finΣ iε references«á  Mos⌠á pubì
lisheΣ program≤ arσ giveε iε obsoletσ dialect≤ oµ FORTRAN¼ anΣ ì
havσá extrßá loo≡á iteration≤ whicΦ d∩á nothinτá bu⌠á activatσ ì
unnecessar∙ conditiona∞ GOTO's.

     Compac⌠á factorizatioε methods¼á oµ whicΦ Crout'≤ i≤á thσ ì
bes⌠ known¼á havσ orde≥ N2ö store≤ t∩ memory¼ whilσ morσ commoε ì
method≤á havσ orde≥ N3ö stores«á  Thσ differencσ iε speeΣ migh⌠ ì
bσ expecteΣ t∩ sho≈ u≡ oε machine≤ witΦ ß slo≈ memor∙ bu≤á anΣ ì
fas⌠ floatinτ poin⌠ CPU¼á sucΦ a≤ thσ 8088-808╖ pair«  Certaiε ì
pipelininτá machine≤á sucΦ a≤ Floatinτ Poin⌠ System≤á ruεá 50Ñ ì
faste≥á a≤ the∙ ruε a⌠ ß speeΣ paceΣ b∙ memor∙ access«á  Crou⌠ ì
als∩ take≤ advantagσ oµ extendeΣ precisioε t∩ improvσá accuraì
cy¼á a≤á computation≤ arσ shorteneΣ t∩ memor∙ forma⌠ precisioε ì
onl∙ wheε stored«á  Agaiε thi≤ shoulΣ favo≥ thei≥ usσ oεá 808╖ ì
o≥ 6888▒ machines«  Computatioε usinτ C'≤ norma∞ conventioε oµ ì
promotioεá t∩ doublσ give≤ thi≤ effec⌠ automaticall∙ wheε opeì
ratinτ oε floa⌠ arrays«

     ┴á deficienc∙ whicΦ ├ share≤ witΦ Pasca∞ i≤ tha⌠ thσ sizσ ì
oµ multipl∙ subscripteΣ array≤ mus⌠ bσ knowε a⌠á preprocessinτ ì
time¼ whilσ FORTRA╬ anΣ PL/╔ allo≈ thσ sizσ t∩ bσ giveε a⌠ ruε ìètime«á  Thi≤á differencσá i≤ misunderstooΣ b∙ somσ author≤á oµ ì
comparativσ programminτ languagσ textbooks«á  Whilσ therσá arσ ì
way≤á arounΣ this¼á the∙ arσ no⌠ a≤ cleaε a≤ ß changσ oµá lanì
guage«á  ├ doe≤ permi⌠ thσ firs⌠ dimension¼á equivalen⌠ t∩ thσ ì
las⌠á dimensioεá oµ FORTRAN¼á t∩ vary¼á s∩ tha⌠á thσá fwdbak(⌐ ì
functioεá caε bσ writteε t∩ operatσ oε anythinτ froφ ßá columε ì
vecto≥ t∩ ß squarσ matri° o≥ beyond.

     Therσá arσ opportunitie≤ t∩ inver⌠ ß numbe≥ anΣ replacσ ß ì
largσá numbe≥á oµ divide≤ witΦ aε equa∞ numbe≥á oµá multiplie≤ ì
plu≤ thσ inversion«  Thi≤ i≤ wortΦ doinτ iε thσ scalσ factors¼ ì
sincσá thσ inaccurac∙ onl∙ break≤ tie≤ iε pivo⌠á choices«á  Iε ì
thσ divisioε b∙ diagona∞ elements¼ thσ extrß roundofµ caε gro≈ ì
t∩á significan⌠á deterioratioε iε fina∞á results¼á unles≤á thσ ì
inverteΣ numbe≥ i≤ carrieΣ iε extrß precision«á  Oε thosσá maì
chine≤ whicΦ havσ n∩ dividσ instruction¼ thσ compile≥ ough⌠ t∩ ì
bσá ablσá t∩ movσ thσ reciproca∞ instructioε ou⌠ oµ thσá loop╗ ì
thu≤á i⌠á make≤ n∩ sensσ t∩ penalizσá norma∞á computers«á  Thσ ì
extrßá timσ takeε b∙ thσ usua∞ dividσ oµ abou⌠ twicσ thσá timσ ì
oµ ß multipl∙ o≥ adΣ i≤ insignificant.

     Thσ initializatioε oµ thσ do⌠ product≤ t∩ 0«á rathe≥ thaε ì
t∩á thσá valuσ whicΦ i≤ t∩ bσ changeΣ improve≤á accurac∙á wheε ì
therσá i≤á n∩ extendeΣ precision«á  Thi≤ happen≤á becausσá thσ ì
valuσá t∩ bσ changeΣ ofteε i≤ thσ larges⌠á iεá magnitude¼á anΣ ì
shoulΣ bσ thσ las⌠ added«

     Iµá thσ origina∞ matri° caε bσ assumeΣ accuratσ a⌠á leas⌠ ì
t∩á thσ precisioε oµ thσ arithmetic¼á aε accurac∙ estimatσá oµ ì
.1¼á thσá wors⌠á likel∙ iε ß properl∙ poseΣ physica∞á problem¼ ì
mean≤ tha⌠ thσ factor≤ havσ los⌠ ▒ decima∞ placσ oµá accuracy«  ì
Oε thσ orde≥ 1▓ Hilber⌠ matrix¼á wσ havσ aε estimatσ oµ 4e-14¼ ì
whicΦ indicate≤ tha⌠ ß maximuφ erro≥ oµ .00│ migh⌠ possibl∙ bσ ì
achieveΣ witΦ standarΣ 1╢ digi⌠ accurac∙ witΦ al∞ calculation≤ ì
carrieΣ t∩ lonτ doublσ precisioε beforσ rounding«  Aε accurac∙ ì
estimatσ oµ les≤ thaε 1e-1╢ woulΣ sho≈ singularit∙ o≥ neeΣ fo≥ ì
highe≥ precision.

     Testinτá witΦ pseudo-randoφ number≤ i≤ no⌠á satisfactory¼ ì
sincσ randoφ numbe≥ generator≤ usuall∙ generatσ onl∙ numbe≥ oµ ì
3▒á bit≤á o≥ less¼á anΣ thσ number≤ generateΣ b∙ onσ ruεá timσ ì
librar∙á wil∞ posσ ß les≤ severσ tes⌠ thaεá another«á  Hilber⌠ ì
matrice≤á arσ ofteε useΣ withou⌠ scaling«á  WitΦ thσá scaling¼ ì
thσá coefficient≤ arσ exac⌠ anΣ thσ tes⌠ i≤ no⌠ madσ easie≥ b∙ ì
roundofµ iε fillinτ ou⌠ thσ tes⌠ matrix«á  Thi≤ i≤ no⌠ ßá triì
via∞á question¼á a≤á somσá author≤ havσ claimeΣ tha⌠á IB═á 37░ ì
machine≤ caε facto≥ sucΦ ß matri° oµ orde≥ 1╢ successfull∙á iε ì
doublσ precision¼ whilσ the∙ fai∞ thσ orde≥ 1▓ tes⌠ wheε i⌠ i≤ ì
donσá carefully«á  ┴á large≥ scalinτ whicΦ make≤ morσá oµá thσ ìècalculation≤ exac⌠ iε thσ forwarΣ anΣ bacδ substitution≤ give≤ ì
morσ consisten⌠ result≤ thaε thσ minimuφ requireΣ fo≥ aε exac⌠ ì
tes⌠ matrix.

    Iεá benchmark≤ comparinτ variou≤ machines¼á wσá havσá useΣ ì
botΦá thσá schemσá herσ whicΦ compare≤ thσ solutioεá witΦá thσ ì
exac⌠á vecto≥á oµá al∞ 1'≤ anΣ thσ checδá oµá multiplyinτá thσ ì
matri° b∙ it≤ equivalen⌠ inverse«á  Thσ MI╪ compile≥ give≤ thσ ì
excellen⌠á erro≥ oµ .02╖ iε thσ slo≈ timσ oµ 3░á second≤á (2.╡ ì
Mh· Z80)«á  Eco├ doe≤ no⌠ accep⌠ doublσ subscripts«  Microsof⌠ ì
doublσ precisioε softwarσ floatinτ poin⌠ take≤ 2┤ second≤ witΦ ì
aε erro≥ oµ .35«á  Azteπ take≤ 2╣ second≤ fo≥ aε erro≥ oµ 2.3¼ ì
whicΦá i≤ worsσ thaε eveε thσ biτ IBM's«á  Thi≤ i≤ ß resul⌠ oµ ì
Aztec'≤á unusua∞ choicσ oµ basσ 25╢ rathe≥ thaε norma∞á binar∙ ì
floatinτ point«  Oε computer≤ witΦ hardwarσ floatinτ point¼ wσ ì
perforφ thσ timinτ tes⌠ oε ß matri° oµ orde≥ 100¼á usinτ orde≥ ì
1▓ onl∙ t∩ tes⌠ accuracy«  Typica∞ 6802░ computer≤ witΦ Weiteδ ì
floatinτ poin⌠ ruε abou⌠ 300░ time≤ Z8░ speed.

     Symmetriπá matrices¼á sucΦ a≤ ou≥ tes⌠ matrix¼á arisσá iε ì
importan⌠á application≤ sucΦ a≤ statistic≤ anΣ finitσá elemen⌠ ì
analysis«  Mos⌠ oµ thesσ matrice≤ arσ positivσ definite¼ whicΦ ì
mean≤ tha⌠ thσ diagona∞ element≤ oµ thσ triangula≥ facto≥ wil∞ ì
alway≤á bσá positivσ withou⌠ requirinτá pivoting«á  Sincσá thσ ì
matri°á canno⌠ bσ storeΣ efficientl∙ iε ßá rectangula≥á array¼ ì
anΣá finitσ elemen⌠ application≤ producσ matrice≤ whicΦ caε bσ ì
packeΣá furthe≥ b∙ skippinτ element≤ whicΦ arσ knowε ßá priorΘ ì
t∩á bσ zer∩ iε thσ factors¼á wσ choosσ thσ variablσá bandwidtΦ ì
columε storagσ forma⌠ iε EQSYM.C«

     Textbook≤á morσá ofteεá choosσ thσ CholeskΘá squarσá roo⌠ ì
method¼á whicΦá i≤ les≤ accuratσ anΣ les≤á general¼á anΣá morσ ì
ofteε sho≈ ß constan⌠ bandwidtΦ method¼á whicΦ permit≤ rectanì
gula≥ (doublσ subscripted⌐ storagσ a⌠ thσ expensσ oµ efficienì
cy«á  Algorithm≤ sucΦ a≤ "reversσ Cuthill-McKeeó perforφ autoì
matiπá sortinτ oµ equation≤ int∩ ß morσ compac⌠ variablσá banΣ ì
storage.

     Factorizatioε int∩ thσ forφ ╠ ─ Ltö avoid≤ thσá inaccurac∙ ì
oµ thσ squarσ roo⌠ method¼á anΣ work≤ oε negativσ definitσ anΣ ì
ofteεá eveεá indefinitσ matrices«á  ┴ peculiarit∙ i≤ tha⌠á thσ ì
facto≥á ─á Ltöá i≤ calculateΣ firs⌠ anΣ replaceΣá a⌠á thσá las⌠ ì
instan⌠ b∙ it≤ fina∞ form«á  Iε thσ forward-bacδ substitution¼ ì
threσ loop≤ arσ required¼ accordinτ t∩ thσ threσ factor≤ whicΦ ì
arσá formall∙á employed«á  Thσ seconΣ anΣ thirΣ loop≤á caεá bσ ì
combineΣ iε thσ fixeΣ bandwidtΦ form¼ wherσ thσ bacδ substituì
tioεá caε bσ donσ b∙ do⌠ products«á  Iε thσ variablσ bandwidtΦ ì
o≥ othe≥ ful∞ symmetriπ packeΣ storagσ schemes¼á thσ do⌠á proìèduc⌠ bacδ substitutioε i≤ to∩ inefficient¼á a≤ thσ element≤ oµ ì
thσá row≤á oµ Ltö arσ storeΣ a⌠ varyinτá intervals«á  Thσá bacδ ì
substitutioεá shown¼á whicΦ i≤ likσ Gaussiaε forwarΣá eliminaì
tion¼ i≤ mucΦ faste≥ bu⌠ slightl∙ les≤ accurate.

     Thσá callinτá prograφ provide≤ aε arra∙ oµ indice≤á whicΦ ì
poin⌠ t∩ thσ diagona∞ element≤ oµ thσ matrix«  Sincσ thσ part≤ ì
oµá thσ matri° column≤ whicΦ arσ storeΣ enΣ oεá thσá diagonal¼ ì
thi≤á arra∙á contain≤ enougΦ informatioε t∩ describσ thσá comì
pletσ forφ oµ thσ matrix«  Iε ß typica∞ variablσ bandeΣ appliì
cation¼á onl∙á ßá fe≈ column≤ wil∞ star⌠ oε thσ ░ ro≈á oµá thσ ì
matri° a≤ the∙ d∩ fo≥ ß ful∞ matrix«  Thσ storeΣ columε start≤ ì
witΦ thσ firs⌠ ro≈ containinτ ß non-zer∩ element.

     A⌠ firs⌠ i⌠ migh⌠ seeφ tha⌠ iε ├ ß vecto≥ oµ pointer≤á t∩ ì
thσ diagona∞ element≤ rathe≥ thaε arra∙ indice≤ coulΣ bσ used«  ì
Thi≤ turn≤ ou⌠ t∩ complicatσ thσ code¼á althougΦ i⌠ eliminate≤ ì
aε argumen⌠ oµ thσ functions«á  Recover∙ oµ thσ columε length≤ ì
i≤á to∩ difficul⌠ witΦ thσ limitation≤ oµ pointe≥á arithmetic¼ ì
whicΦ no⌠ al∞ compiler≤ arσ willinτ t∩ tackle.

     MI╪á anΣá Eco├á achievσ identica∞ respectablσá error≤á oµ ì
0.059¼á whilσ Azteπ i≤ s∩ poo≥ a≤ no⌠ t∩ bσ suitablσ fo≥á thi≤ ì
algorithm«á  AlthougΦ thσ reorderinτ oµ equation≤ performeΣ iε ì
EQD.├ i≤ no⌠ requireΣ fo≥ success¼ i⌠ improve≤ accurac∙ a⌠ thσ ì
cos⌠ oµ nearl∙ twicσ a≤ mucΦ time.

     Aεá unsymmetriπ versioε oµ thσ variablσ bandwidtΦá methoΣ ì
(EQVB.C⌐ i≤ effectivσ iε finitσ differencσ methods¼á wherσ thσ ì
storagσá patterε a≤ wel∞ a≤ thσ coefficient≤ arσ likel∙ t∩á bσ ì
unsymmetric«á  Thesσá method≤ arσ ofteε successfu∞á wherσá thσ ì
orde≥ oµ equation≤ i≤ no⌠ guaranteeΣ t∩ bσ correct¼á ye⌠ i⌠ i≤ ì
impractica∞á t∩á changσ i⌠ t∩ maximizσá pivots«á  Accurac∙á i≤ ì
usuall∙á bette≥á thaε symmetriπ method≤ (thσ differencσ i≤á iε ì
thσ forwarΣ substitution)╗á howeve≥ i⌠ i≤ no⌠ recommendeΣá unì
les≤á thσ matri° i≤ sparsσ enougΦ t∩ makσ i⌠ ruε significantl∙ ì
faste≥ thaε thσ morσ genera∞ methoΣ oµ EQD.C.