!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! COPYRIGHT NOTE: If you use results generated by this code, ! please refer to the paper: S.\ Bhattacharya, J.\ B{\l}awzdziewicz, ! and E.\ Wajnryb, Physica A {\bf 356}, 294 (2005), where our ! numerical method is described. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Fortran 95 routine ! ! WALL2_poiseu2vel(Utrans,Orot,Z,HH) ! ! to evaluate translational and angular velocity ! of a sphere at a position Z in Poiseuille flow ! in parallel-wall channel of width HH ! under creeping flow conditions ! ! INPUT: ! ! Z - distance of particle center from the lower wall ! (normalized by particle diameter d) ! ! HH - wall separation (normalized by particle diameter d) ! ! OUTPUT: ! ! Utrans - translational velocity of the particle ! (normalized by fluid velocity in channel center Ufluid) ! ! Orot - rotational velocity of the particle ! (normalized by Ufluid/d) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! NOTE: THE ROUTINE WALL2_poiseu2vel REQUIRES THAT DATA FILE ! rig_wall_tt.cf IS IN THE DIRECTORY FROM WHICH THE ROUTINE IS ! CALLED ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PROGRAM MAIN IMPLICIT NONE REAL*8 Utrans,Orot REAL*8 Z,HH 1 WRITE(*,*) 'WRITE Z and HH' READ(*,*) Z,HH CALL WALL2_poiseu2vel(Utrans,Orot,Z,HH) WRITE(*,*) 'Utrans, Orot: ',Utrans,Orot WRITE(*,*) GOTO 1 STOP END ********************************************** ***************************************** MODULE SIZE ! NN,LLG INTEGER NN,LLG END MODULE SIZE ***************************************** MODULE MULT ! LL INTEGER LL END MODULE MULT ***************************************** MODULE MULT_FULL_polym ! LL2,LLG1 INTEGER LL2,LLG1 END MODULE MULT_FULL_polym ***************************************** MODULE MULT_WALL_polym ! LC,LC2,LCG1 INTEGER LC,LC2,LCG1 END MODULE MULT_WALL_polym C*********************************************************** MODULE Z_COEF REAL*8, ALLOCATABLE :: Z0(:),Z2(:),Z1(:),ZB(:) END MODULE Z_COEF ***************************************** MODULE S_COEF_polym REAL*8, ALLOCATABLE :: S11P(:,:,:),So00P(:,:,:),S00P(:,:,:), * S01P(:,:,:), S02P(:,:,:),S10P(:,:,:), * S20P(:,:,:) END MODULE S_COEF_polym ***************************************** MODULE FACT REAL*8, ALLOCATABLE :: SIL(:) END MODULE FACT ***************************************** MODULE C_COEF_polym REAL*8, ALLOCATABLE :: C01(:,:),C02(:,:),C03(:,:),C04(:,:), * C11(:,:),C12(:,:),C13(:,:),C14(:,:), * C21(:,:),C22(:,:),C23(:,:),C24(:,:), * C25(:,:),C26(:,:),C27(:,:) END MODULE C_COEF_polym ***************************************** MODULE RIG_WALL_COEF INTEGER MM PARAMETER(MM=450) REAL*8 XA(0:MM) REAL*8 YA(0:MM) REAL*8 YB(0:MM) REAL*8 XC(0:MM) REAL*8 YC(0:MM) DATA XA/ * 9.5238095238095238D-01, * -5.1190476190476190D-02, * 1.7356150793650794D-01, * -1.3887028769841270D-01, * -4.3468656994047619D-02, * 1.7670570979042659D-01, * -1.0226409116299515D-01, * -6.0164611842356572D-02, * 1.2978188051748528D-01, * -1.0280710836409261D-01, * 2.0401779394999985D-02, * 7.7241563477720071D-02, * -1.0534623431182061D-01, * 5.6064930574563372D-02, * 7.0993228769989266D-03, * -4.9203622467965860D-02, * 4.9148669182127192D-02, * -1.7208325753720004D-02, * -9.7105767827494033D-03, * 1.5947261217746953D-02, * -2.7757377462879601D-03, * -1.3517840132789262D-02, * 1.6751323984348835D-02, * -5.9990489156973835D-03, * -1.1177751592713057D-02, * 2.1665779275267472D-02, * -2.0381498878674745D-02, * 1.0685535254547363D-02, * 4.6512765490163544D-04, * -6.1792886523887267D-03, * 6.3479092570419252D-03, * -3.8002471187771852D-03, * 2.5416792660147241D-03, * -4.0132851559881031D-03, * 6.4874990527691414D-03, * -7.7673713524265220D-03, * 6.1521256266634450D-03, * -2.5526871716271322D-03, * -1.3257640899816405D-03, * 3.7374260476572417D-03, * -4.2462597577914645D-03, * 3.6332562711491642D-03, * -2.8474940773272364D-03, * 2.5986313352081607D-03, * -2.6387554726100231D-03, * 2.5718884489997532D-03, * -1.9702386609148141D-03, * 9.5474585457224278D-04, * 9.2852076097682363D-05, * -7.9996969859048812D-04, * 1.0110791900069201D-03, * -9.4210731747863647D-04, * 8.1521261523556961D-04, * -8.6285826595829178D-04, * 1.0328907473366725D-03, * -1.2328294196438381D-03, * 1.3144880077576746D-03, * -1.2569440977051753D-03, * 1.1167533392820420D-03, * -9.6158321102220864D-04, * 8.4184261029768135D-04, * -7.1036903621594301D-04, * 5.3851480365063560D-04, * -2.7544530806411168D-04, * -4.1710238142883719D-05, * 3.7649171294785923D-04, * -6.6347317732891396D-04, * 8.7838467752370194D-04, * -1.0174654659438292D-03, * 1.0911180550118885D-03, * -1.1155717982950519D-03, * 1.0837772391183746D-03, * -1.0014955602289274D-03, * 8.5983234476680634D-04, * -6.8082711114994940D-04, * 4.7581660990321763D-04, * -2.7382931191245693D-04, * 8.3511663966472414D-05, * 8.4308487402709365D-05, * -2.3128139231235363D-04, * 3.5675237769373846D-04, * -4.5854650219102089D-04, * 5.3605170163175884D-04, * -5.8450600656460830D-04, * 6.0839989228000980D-04, * -6.0691968118627369D-04, * 5.8833472835652320D-04, * -5.5217311829981788D-04, * 5.0412879942969494D-04, * -4.4265216033271572D-04, * 3.7179902600148679D-04, * -2.9265526447540346D-04, * 2.0961978357285721D-04, * -1.2618631425429635D-04, * 4.5093733871756066D-05, * 3.0266405369553636D-05, * -9.9534179198418776D-05, * 1.5955283846062932D-04, * -2.1106755035037456D-04, * 2.5116938021767555D-04, * -2.8152382857228225D-04, * 3.0048133462262794D-04, * -3.1060545612640070D-04, * 3.1150999936435050D-04, * -3.0563355668694883D-04, * 2.9313435785145245D-04, * -2.7565913856305918D-04, * 2.5376849290676552D-04, * -2.2842253432912317D-04, * 2.0059470791835156D-04, * -1.7063847408444569D-04, * 1.3968270584319173D-04, * -1.0755989104357157D-04, * 7.5439165809406779D-05, * -4.2968754340404261D-05, * 1.1379025918223755D-05, * 1.9605482423974600D-05, * -4.8828013195681455D-05, * 7.6446328216006944D-05, * -1.0155256913219249D-04, * 1.2416742191533851D-04, * -1.4369270181995866D-04, * 1.5999432554210893D-04, * -1.7280962794179930D-04, * 1.8190657980808575D-04, * -1.8737332206897948D-04, * 1.8895071798224684D-04, * -1.8702776235108129D-04, * 1.8136029778154745D-04, * -1.7255451021425158D-04, * 1.6042983125599939D-04, * -1.4573872970435586D-04, * 1.2841391812428139D-04, * -1.0928124447663964D-04, * 8.8403559786672991D-05, * -6.6598232311465952D-05, * 4.4045878067551530D-05, * -2.1484140715440892D-05, * -8.1717978092818428D-07, * 2.2256098015444599D-05, * -4.2518598156324211D-05, * 6.1180678133428123D-05, * -7.7935860030809552D-05, * 9.2557368198589551D-05, * -1.0478829043684166D-04, * 1.1459097667962487D-04, * -1.2178539178099348D-04, * 1.2649490054405102D-04, * -1.2863033807688422D-04, * 1.2843526989117940D-04, * -1.2591212911589904D-04, * 1.2137713604312790D-04, * -1.1491448607576457D-04, * 1.0686676377143118D-04, * -9.7384586908680333D-05, * 8.6797926368758288D-05, * -7.5305828930683911D-05, * 6.3196014379708248D-05, * -5.0696768817859077D-05, * 3.8034031291874197D-05, * -2.5446992981242184D-05, * 1.3090547145490357D-05, * -1.1988478216000729D-06, * -1.0144652404871586D-05, * 2.0723865628612409D-05, * -3.0520438085760335D-05, * 3.9346097514432780D-05, * -4.7236565475268989D-05, * 5.4037725198188023D-05, * -5.9825890599477191D-05, * 6.4484485074712287D-05, * -6.8116660420672403D-05, * 7.0644399943520956D-05, * -7.2185225119943908D-05, * 7.2698719634075598D-05, * -7.2306473581468860D-05, * 7.1003044367500646D-05, * -6.8906345554265934D-05, * 6.6041885033310362D-05, * -6.2518641583877643D-05, * 5.8387846328817343D-05, * -5.3746394374925486D-05, * 4.8665053319635560D-05, * -4.3227158404427227D-05, * 3.7516249797657265D-05, * -3.1601900060713226D-05, * 2.5573533001934159D-05, * -1.9487666684669063D-05, * 1.3433048277102020D-05, * -7.4544580025401003D-06, * 1.6341524516373497D-06, * 3.9931431915433267D-06, * -9.3564080297686762D-06, * 1.4429134879375826D-05, * -1.9155027061609469D-05, * 2.3513988734102015D-05, * -2.7466617857864256D-05, * 3.0997597518737768D-05, * -3.4085318902349694D-05, * 3.6717889723001532D-05, * -3.8891258749609697D-05, * 4.0595948606747158D-05, * -4.1844296081256835D-05, * 4.2628598477215648D-05, * -4.2975707262360569D-05, * 4.2879407444825549D-05, * -4.2378722141747315D-05, * 4.1468938189073049D-05, * -4.0198653221633490D-05, * 3.8564887628747442D-05, * -3.6623141506432741D-05, * 3.4372511403500508D-05, * -3.1872788573118286D-05, * 2.9125483635364399D-05, * -2.6192223318517702D-05, * 2.3077157301358421D-05, * -1.9841509393317357D-05, * 1.6492094135443730D-05, * -1.3087749963418428D-05, * 9.6377301994179023D-06, * -6.1967871050560015D-06, * 2.7761216015458556D-06, * 5.7499991122350134D-07, * -3.8441599830363341D-06, * 6.9886608419491837D-06, * -9.9957823592281771D-06, * 1.2830149368408573D-05, * -1.5479756666119765D-05, * 1.7916971161232452D-05, * -2.0131534679575608D-05, * 2.2103644604401611D-05, * -2.3825748367451667D-05, * 2.5285646222225899D-05, * -2.6479299056527586D-05, * 2.7401594286065755D-05, * -2.8052604008949961D-05, * 2.8433539462845426D-05, * -2.8548934168896329D-05, * 2.8405363477722848D-05, * -2.8011910613719901D-05, * 2.7379417889836614D-05, * -2.6521351793205395D-05, * 2.5451649609534155D-05, * -2.4187765381845108D-05, * 2.2745546536048837D-05, * -2.1145848852287864D-05, * 1.9405289560019396D-05, * -1.7547398100963161D-05, * 1.5588515006104954D-05, * -1.3554024795075417D-05, * 1.1459077653048569D-05, * -9.3300555517772685D-06, * 7.1801650357870434D-06, * -5.0359367160533232D-06, * 2.9080563620130500D-06, * -8.2240504753025196D-07, * -1.2132484032560851D-06, * 3.1743857665509307D-06, * -5.0563752710276717D-06, * 6.8366641926348681D-06, * -8.5138147931179516D-06, * 1.0067714684454892D-05, * -1.1500033135142539D-05, * 1.2793436345311317D-05, * -1.3952492985053901D-05, * 1.4962848457412728D-05, * -1.5831668584378367D-05, * 1.6547646652007621D-05, * -1.7120175529368557D-05, * 1.7540944410843945D-05, * -1.7821161441534867D-05, * 1.7955354637229681D-05, * -1.7956118343070191D-05, * 1.7820575633853717D-05, * -1.7562282231660846D-05, * 1.7180645819832760D-05, * -1.6689780970689589D-05, * 1.6091023130164893D-05, * -1.5398679580441309D-05, * 1.4615629803119671D-05, * -1.3756054002370808D-05, * 1.2823982952037291D-05, * -1.1833202270488916D-05, * 1.0788507847757263D-05, * -9.7030776902641276D-06, * 8.5821061166531367D-06, * -7.4380039514375457D-06, * 6.2760268543315656D-06, * -5.1077087335200449D-06, * 3.9380664582732913D-06, * -2.7776916105766331D-06, * 1.6311039038712993D-06, * -5.0792482272984464D-07, * -5.8803709843752087D-07, * 1.6481277958492347D-06, * -2.6694188884606706D-06, * 3.6441977252255222D-06, * -4.5705407957827071D-06, * 5.4416341057826453D-06, * -6.2566414177939368D-06, * 7.0095936215188999D-06, * -7.7007846304559008D-06, * 8.3250301200426764D-06, * -8.8837609033745209D-06, * 9.3725146072938901D-06, * -9.7938333670412770D-06, * 1.0143913719916085D-05, * -1.0426355286017902D-05, * 1.0637951847940580D-05, * -1.0783282634971262D-05, * 1.0859679184663076D-05, * -1.0872602703121114D-05, * 1.0819865378094035D-05, * -1.0707697282552442D-05, * 1.0534336322984564D-05, * -1.0306657064233838D-05, * 1.0023269934725998D-05, * -9.6915624160744400D-06, * 9.3104654762979780D-06, * -8.8877446011520664D-06, * 8.4226000302106949D-06, * -7.9230406117338176D-06, * 7.3884849248773369D-06, * -6.8270540930657995D-06, * 6.2383342131668565D-06, * -5.6304340672929052D-06, * 5.0030565449740201D-06, * -4.3641824281016580D-06, * 3.7135809602468189D-06, * -3.0590002755755351D-06, * 2.4002261797419013D-06, * -1.7446821212809072D-06, * 1.0921218417360804D-06, * -4.4956581883914742D-07, * -1.8331106573651572D-07, * 7.9995517489664136D-07, * -1.4008147036924987D-06, * 1.9798519923188912D-06, * -2.5376792886247538D-06, * 3.0688094281477110D-06, * -3.5740555111043246D-06, * 4.0485012474713137D-06, * -4.4931927246172610D-06, * 4.9037914700389922D-06, * -5.2816030175884440D-06, * 5.6228614408610974D-06, * -5.9291522603403506D-06, * 6.1972653878727727D-06, * -6.4290803966905015D-06, * 6.6219167005241381D-06, * -6.7779553026517810D-06, * 6.8950100428495164D-06, * -6.9755636504746302D-06, * 7.0178828687874523D-06, * -7.0247456961080480D-06, * 6.9948242515017531D-06, * -6.9311788838031094D-06, * 6.8328354303689952D-06, * -6.7031175807107278D-06, * 6.5413508644067778D-06, * -6.3510969366191999D-06, * 6.1319255208221508D-06, * -5.8876057179739511D-06, * 5.6178971294784538D-06, * -5.3267403819301296D-06, * 5.0140304419784174D-06, * -4.6838432368804453D-06, * 4.3361549142405498D-06, * -3.9751343227882453D-06, * 3.6007907180181730D-06, * -3.2173477708852921D-06, * 2.8247936828797687D-06, * -2.4273658215110127D-06, * 2.0249995821495212D-06, * -1.6218927315999274D-06, * 1.2178918237459578D-06, * -8.1711900556130910D-07, * 4.1929306577006217D-07, * -2.8428316982312489D-08, * -3.5590920125904210D-07, * 7.2985852405553631D-07, * -1.0940126525111726D-06, * 1.4447106427723065D-06, * -1.7827189473434927D-06, * 2.1045817813599399D-06, * -2.4112724383302782D-06, * 2.6995694960202286D-06, * -2.9706338847677153D-06, * 3.2214799482981913D-06, * -3.4535030365953574D-06, * 3.6639584252151663D-06, * -3.8544414784824231D-06, * 4.0224536041594515D-06, * -4.1698043323474116D-06, * 4.2942859164058593D-06, * -4.3977584760562776D-06, * 4.4784066062716450D-06, * -4.5382115687354863D-06, * 4.5756409255322494D-06, * -4.5926515350930762D-06, * 4.5882034100500729D-06, * -4.5642493545751706D-06, * 4.5197443458245845D-06, * -4.4572434693017528D-06, * 4.3755563887292700D-06, * -4.2771367554300104D-06, * 4.1614974268732739D-06, * -4.0307340605504180D-06, * 3.8848871673730193D-06, * -3.7254126128839394D-06, * 3.5530943592666756D-06, * -3.3693701424209017D-06, * 3.1746101656740207D-06, * -2.9708574551335947D-06, * 2.7581718776018704D-06, * -2.5388995011191691D-06, * 2.3126332171300616D-06, * -2.0823797611341706D-06, * 1.8472815011040878D-06, * -1.6100294837825499D-06, * 1.3700979276271842D-06, * -1.1310377076963919D-06, * 8.8901929186427732D-07, * -6.5273877826295355D-07, * 4.1382745118168761D-07, * -1.8301165728680858D-07, * -4.8078272950490442D-08, * 2.6981277891346056D-07, * -4.9229052487533637D-07, * 7.0078832203130170D-07, * -9.1071647684298568D-07, * 1.1039002695898686D-06, * -1.2974894182633235D-06, * 1.4748384977820659D-06, * -1.6469467251106648D-06, * 1.8048558659197624D-06, * -1.9571909175162168D-06, * 2.0925796387521033D-06, * -2.2277588118430895D-06, * 2.3381019065765265D-06, * -2.4535387071181527D-06, * 2.5561313943034940D-06, * -2.6294058685356340D-06, * 2.7244717055903965D-06, * -2.7497066359395556D-06, * 2.8405863397306857D-06, * -2.7982433682974968D-06, * 2.9141061748765552D-06, * -2.8579375755580532D-06, * 2.9194873371055231D-06, * -2.8512042498196344D-06/ DATA YA/ * 8.2933333333333333D-01, * 1.1450000000000000D-01, * 7.8184027777777778D-02, * -1.1057703993055556D-01, * 1.0469165039062500D-02, * 9.6937299431694878D-02, * -8.4414303590380956D-02, * -5.6632727767739977D-03, * 6.4604008373218002D-02, * -7.2066827995023966D-02, * 3.4400919207390160D-02, * 2.8814493053670935D-02, * -5.7951299112550897D-02, * 3.7818258857280833D-02, * -2.8448501251530938D-03, * -2.3102846404265437D-02, * 2.3528460018514020D-02, * -4.3035078032105014D-03, * -1.0656532945644665D-02, * 1.0642637228062229D-02, * 2.1581740583442461D-03, * -1.4830443315316760D-02, * 1.6129826812978278D-02, * -6.6567937126443173D-03, * -6.3344597246591213D-03, * 1.3126961196616027D-02, * -1.1230960882338794D-02, * 4.6625390980100508D-03, * 8.5403438730850489D-04, * -1.5917783612250220D-03, * -6.9143086074647723D-04, * 2.6397568610607282D-03, * -1.6847749453564665D-03, * -1.6881308830944241D-03, * 4.8438963362524243D-03, * -5.7939601546069945D-03, * 4.1238835834772206D-03, * -1.5134612173358018D-03, * -3.6804827546977780D-04, * 7.4177418459994729D-04, * -2.0034202174837875D-04, * -1.5251381119354697D-04, * -2.2679508164038070D-04, * 1.1601476142174437D-03, * -1.8748989055010818D-03, * 1.9330395098230725D-03, * -1.3564172906566997D-03, * 6.5532712637277912D-04, * -2.7240756836841390D-04, * 3.0087617292313123D-04, * -5.1547193144298510D-04, * 5.6948774061114589D-04, * -3.6050935793379798D-04, * -3.2323938853241288D-05, * 3.6575938437461777D-04, * -5.3098181306471145D-04, * 5.4091927124391765D-04, * -5.2998319703744083D-04, * 5.9226958573290352D-04, * -7.1726283890530817D-04, * 8.2592380703368731D-04, * -8.2249209112966131D-04, * 7.0136797274331208D-04, * -5.0174475246458743D-04, * 3.0631258734293069D-04, * -1.4623786855780323D-04, * 2.2058999108432494D-05, * 9.9057426582421449D-05, * -2.3233915110560631D-04, * 3.6527812654829408D-04, * -4.7219662546146732D-04, * 5.2661533823802540D-04, * -5.3364871828558942D-04, * 5.0513447245376129D-04, * -4.6746850710901190D-04, * 4.2453190983349648D-04, * -3.7828334503205035D-04, * 3.1723668795577598D-04, * -2.4246546013314357D-04, * 1.5818060613710245D-04, * -7.5810567175748433D-05, * 2.8385375199850398D-06, * 5.9898664190896720D-05, * -1.1399701784567389D-04, * 1.6396203020637879D-04, * -2.0731916247658088D-04, * 2.4287165951879294D-04, * -2.6510461974965171D-04, * 2.7495015138631007D-04, * -2.7261112825933565D-04, * 2.6249838594832221D-04, * -2.4610560196487065D-04, * 2.2487679741902188D-04, * -1.9877709407895707D-04, * 1.6800165520006106D-04, * -1.3454848750983914D-04, * 9.9671602654391703D-05, * -6.6284284809533768D-05, * 3.4227095839847832D-05, * -4.9409296062724337D-06, * -2.2970735532788267D-05, * 4.8186664670311279D-05, * -7.1414366773706695D-05, * 9.1065949393601089D-05, * -1.0769220187706756D-04, * 1.2059182105515542D-04, * -1.3056105353527003D-04, * 1.3743567001334779D-04, * -1.4148467741465329D-04, * 1.4255426186539893D-04, * -1.4053843332914077D-04, * 1.3570947357993842D-04, * -1.2809260314983849D-04, * 1.1834353734198644D-04, * -1.0643842122175826D-04, * 9.3000743542878084D-05, * -7.7926367012726372D-05, * 6.1840227002161269D-05, * -4.4819737619271955D-05, * 2.7518203115789541D-05, * -1.0158172568374961D-05, * -6.7811343617928660D-06, * 2.3039902842330890D-05, * -3.8339944054552044D-05, * 5.2327235996433270D-05, * -6.4858920616530454D-05, * 7.5518358073312167D-05, * -8.4343734439729912D-05, * 9.0965999780270057D-05, * -9.5608331708407706D-05, * 9.7970201125214964D-05, * -9.8386131340307565D-05, * 9.6630644111156399D-05, * -9.3115184880716614D-05, * 8.7727367826899149D-05, * -8.0928233939260582D-05, * 7.2712907186125891D-05, * -6.3530771065509738D-05, * 5.3450859142450857D-05, * -4.2869744877799374D-05, * 3.1914096620563594D-05, * -2.0913178878597597D-05, * 1.0029248987914598D-05, * 4.9287396224743236D-07, * -1.0488445957065251D-05, * 1.9803341930212461D-05, * -2.8290532975561540D-05, * 3.5873861098742812D-05, * -4.2436168896152893D-05, * 4.7964850503012226D-05, * -5.2385874388742845D-05, * 5.5735682364622558D-05, * -5.7989003086165253D-05, * 5.9211540528559577D-05, * -5.9424419717652693D-05, * 5.8704740271450789D-05, * -5.7116575831339524D-05, * 5.4736040739774226D-05, * -5.1664287028929266D-05, * 4.7966276627092813D-05, * -4.3770871465947337D-05, * 3.9125044940674901D-05, * -3.4176098303186616D-05, * 2.8951034224701445D-05, * -2.3607004985194391D-05, * 1.8151524651878518D-05, * -1.2742675571201839D-05, * 7.3702386960681866D-06, * -2.1852790044061641D-06, * -2.8368342378219988D-06, * 7.5582482693801343D-06, * -1.2014223057438631D-05, * 1.6085080032504231D-05, * -1.9813399910840916D-05, * 2.3101102753155720D-05, * -2.5995064643020912D-05, * 2.8420267400462518D-05, * -3.0425299741024430D-05, * 3.1958093258098706D-05, * -3.3067169387321214D-05, * 3.3722002590295273D-05, * -3.3969896623846837D-05, * 3.3799240689598642D-05, * -3.3255436868087260D-05, * 3.2342349204234625D-05, * -3.1103306355513684D-05, * 2.9553836827751426D-05, * -2.7735397755966170D-05, * 2.5671237376162422D-05, * -2.3401316190924954D-05, * 2.0952798424927497D-05, * -1.8364628440569169D-05, * 1.5664499863265481D-05, * -1.2890834151115824D-05, * 1.0069003641058315D-05, * -7.2372774706582303D-06, * 4.4164478496768161D-06, * -1.6448193586329136D-06, * -1.0630199496490311D-06, * 3.6687503596031375D-06, * -6.1650102229764321D-06, * 8.5137220927674854D-06, * -1.0715246854368898D-05, * 1.2732225705795982D-05, * -1.4572778427704040D-05, * 1.6200899718275681D-05, * -1.7632121100135022D-05, * 1.8832518016041130D-05, * -1.9824388599826851D-05, * 2.0576630477698457D-05, * -2.1117435477199776D-05, * 2.1419190888844738D-05, * -2.1514943501505371D-05, * 2.1381090468872432D-05, * -2.1054392153164388D-05, * 2.0515575420666622D-05, * -1.9803926121861488D-05, * 1.8904585079731053D-05, * -1.7858177901871562D-05, * 1.6654097839368693D-05, * -1.5333174636826049D-05, * 1.3888661945541559D-05, * -1.2360550891008691D-05, * 1.0745366050600826D-05, * -9.0813446177560752D-06, * 7.3675421283194973D-06, * -5.6396816932944041D-06, * 3.8985160518542740D-06, * -2.1766685327678951D-06, * 4.7572038001060497D-07, * 1.1752074393937537D-06, * -2.7745532068734712D-06, * 4.2969187511999869D-06, * -5.7415398787075432D-06, * 7.0867890463596519D-06, * -8.3333677871361587D-06, * 9.4633119511287640D-06, * -1.0479318986197167D-05, * 1.1366849344785543D-05, * -1.2130971807546423D-05, * 1.2760228189179806D-05, * -1.3262273494585288D-05, * 1.3628312200605141D-05, * -1.3868644259910007D-05, * 1.3976671708311944D-05, * -1.3965254669531503D-05, * 1.3829508559061649D-05, * -1.3584644407422319D-05, * 1.3227012193110347D-05, * -1.2773863629181092D-05, * 1.2222330547516437D-05, * -1.1591320229979342D-05, * 1.0878336170009107D-05, * -1.0103507959593514D-05, * 9.2643549301021592D-06, * -8.3817734803486897D-06, * 7.4530047174000440D-06, * -6.4992578693122895D-06, * 5.5172667000288834D-06, * -4.5281214282603570D-06, * 3.5278837156939953D-06, * -2.5371318833102400D-06, * 1.5511513712435236D-06, * -5.8966735102957510D-07, * -3.5286060489182862D-07, * 1.2578415903267361D-06, * -2.1316449623146446D-06, * 2.9570305408858365D-06, * -3.7411619535231984D-06, * 4.4683005958965249D-06, * -5.1463411313388535D-06, * 5.7611338381013366D-06, * -6.3212214364371660D-06, * 6.8140713326681106D-06, * -7.2487796089496554D-06, * 7.6144062372168993D-06, * -7.9205011221239882D-06, * 8.1576468530210679D-06, * -8.3357482843283031D-06, * 8.4468033423612479D-06, * -8.5009777051132669D-06, * 8.4915482506967712D-06, * -8.4288547962197552D-06, * 8.3072957203201782D-06, * -8.1373071912454888D-06, * 7.9142381508530498D-06, * -7.6485526508725299D-06, * 7.3363727050901530D-06, * -6.9881307377138786D-06, * 6.6005439523715887D-06, * -6.1839616928559100D-06, * 5.7355234090706149D-06, * -5.2654507746340119D-06, * 4.7711419516024215D-06, * -4.2626519314468867D-06, * 3.7374870807574028D-06, * -3.2055010632564563D-06, * 2.6641737581193852D-06, * -2.1231262029704257D-06, * 1.5796949152940000D-06, * -1.0432396009036050D-06, * 5.1085562508873506D-07, * 8.3852106666225366D-09, * -5.1770680985321263D-07, * 1.0083498827063380D-06, * -1.4839160629333655D-06, * 1.9359811220546484D-06, * -2.3685637853165155D-06, * 2.7735950731029629D-06, * -3.1555320775066523D-06, * 3.5066784378322502D-06, * -3.8319341951626383D-06, * 4.1239896176735147D-06, * -4.3881772745223077D-06, * 4.6175841461105192D-06, * -4.8179518122943729D-06, * 4.9827695897573246D-06, * -5.1181535002792161D-06, * 5.2179959119635716D-06, * -5.2887437785818857D-06, * 5.3246879984807020D-06, * -5.3325561150057837D-06, * 5.3070276221717660D-06, * -5.2550555262789391D-06, * 5.1716925592449926D-06, * -5.0640592176003833D-06, * 4.9275608541641044D-06, * -4.7694264246615343D-06, * 4.5853883612220220D-06, * -4.3827255966621739D-06, * 4.1574676703316183D-06, * -3.9168869650784691D-06, * 3.6572763604568356D-06, * -3.3858483108328643D-06, * 3.0991222279431709D-06, * -2.8042013853087085D-06, * 2.4977927232489310D-06, * -2.1868459736861435D-06, * 1.8682153141319560D-06, * -1.5486580285735561D-06, * 1.2251348924645337D-06, * -9.0417766394260003D-07, * 5.8281367119231817D-07, * -2.6732209905971170D-07, * -4.5241708025853451D-08, * 3.4887210679416923D-07, * -6.4652190420274280D-07, * 9.3247397310903621D-07, * -1.2097222898379437D-06, * 1.4728474246840546D-06, * -1.7249137117820813D-06, * 1.9608031042646170D-06, * -2.1836750918353719D-06, * 2.3887115771261092D-06, * -2.5791871837230939D-06, * 2.7505776530550106D-06, * -2.9062876141945756D-06, * 3.0420763506347710D-06, * -3.1614881204481823D-06, * 3.2605517947182592D-06, * -3.3429556400419283D-06, * 3.4049810561793748D-06, * -3.4504595982949337D-06, * 3.4759055986889595D-06, * -3.4852883570985913D-06, * 3.4753335757071735D-06, * -3.4501383237592122D-06, * 3.4066167084423665D-06, * -3.3489796561089324D-06, * 3.2743058603526479D-06, * -3.1869028300417813D-06, * 3.0839896966934547D-06, * -2.9699505502974026D-06, * 2.8421209728505368D-06, * -2.7049395808198881D-06, * 2.5558350308236257D-06, * -2.3992770474873493D-06, * 2.2327650014419519D-06, * -2.0607756277058589D-06, * 1.8808580169968691D-06, * -1.6974717965450066D-06, * 1.5081964450847527D-06, * -1.3174509596069395D-06, * 1.1228277729897111D-06, * -9.2868288349603689D-07, * 7.3260631894543965D-07, * -5.3887035264872685D-07, * 3.4504943995340296D-07, * -1.5531345482738765D-07, * -3.2789635705307792D-08, * 2.1520865428472598D-07, * -3.9443078068034541D-07, * 5.6653721563949908D-07, * -7.3405673034755446D-07, * 8.9321381682093808D-07, * -1.0465831557329601D-06, * 1.1905406054565580D-06, * -1.3277090724052371D-06, * 1.4546210259185094D-06, * -1.5739480679428582D-06, * 1.6823817705091909D-06, * -1.7826412334116479D-06, * 1.8715770090754406D-06, * -1.9519530256232473D-06, * 2.0207762700340773D-06, * -2.0808515630377459D-06, * 2.1293375887217611D-06, * -2.1690750679229454D-06, * 2.1973677159293721D-06, * -2.2170863120651643D-06, * 2.2256712927856066D-06, * -2.2260170062546609D-06, * 2.2156909520017641D-06, * -2.1976040999012486D-06, * 2.1694403055074501D-06, * -2.1341199326672362D-06, * 2.0894317313541431D-06, * -2.0382981119598770D-06, * 1.9786007857335790D-06, * -1.9132568859751689D-06, * 1.8402289466130794D-06, * -1.7624216493027050D-06, * 1.6778662513411920D-06, * -1.5894480642685323D-06, * 1.4952552288160353D-06, * -1.3981471131589398D-06, * 1.2962573540971915D-06, * -1.1924132205679038D-06, * 1.0847830756573039D-06, * -9.7615640693885293D-07, * 8.6472628800472264D-07, * -7.5323925870583653D-07, * 6.3990394952054487D-07, * -5.2741933121268763D-07, * 4.1400138055345445D-07, * -3.0229744089694493D-07, * 1.9052362279439036D-07, * -8.1272155366000661D-08, * -2.7246900367985701D-08, * 1.3249934447078706D-07, * -2.3628631020594793D-07, * 3.3613397621108348D-07, * -4.3385878627822361D-07, * 5.2704909322354282D-07, * -6.1753976610143546D-07, * 7.0298258988451872D-07, * -7.8523301357277020D-07, * 8.6200699430507315D-07, * -9.3518180913829562D-07, * 1.0025378318569673D-06, * -1.0659745722596427D-06/ DATA YB/ * 3.4400000000000000D-01, * 2.8000000000000000D-02, * 4.2666666666666667D-02, * 3.8000000000000000D-02, * -1.5470000000000000D-01, * 9.8845833333333333D-02, * 6.4693638392857143D-02, * -1.5161317661830357D-01, * 1.0030225457085503D-01, * 1.4464305305480957D-02, * -1.0586455375172875D-01, * 1.0860190060881503D-01, * -3.0035581703494804D-02, * -5.1865354391069639D-02, * 7.7595361374427159D-02, * -4.6567280970084801D-02, * -6.5130242677305672D-03, * 3.5818029104909360D-02, * -2.6235245900191135D-02, * -2.0973445864916419D-03, * 2.0760017731973319D-02, * -1.4783355388867999D-02, * -6.9263196909097121D-03, * 2.4180140314812979D-02, * -2.4042774776681733D-02, * 8.3717363423324704D-03, * 9.0844877373187935D-03, * -1.6881874146075945D-02, * 1.2799380838097530D-02, * -3.6735207528461191D-03, * -1.9332211228965751D-03, * 9.8380596786565263D-04, * 3.5383179211796536D-03, * -6.0998723213829736D-03, * 3.9137141785963466D-03, * 1.4245898768749853D-03, * -6.0737976482944627D-03, * 7.1303552410808283D-03, * -4.8859561832736350D-03, * 1.6494488390042989D-03, * 1.5228800746136123D-04, * 2.6379902235052066D-05, * -1.0895767435606099D-03, * 1.5029032611225051D-03, * -6.5230828551534134D-04, * -8.5514175421028706D-04, * 1.9637935489304898D-03, * -2.0432957310123992D-03, * 1.3444298149348353D-03, * -5.9182134108211271D-04, * 3.6165783804787020D-04, * -6.8261542295831656D-04, * 1.1266687534059483D-03, * -1.2732295864410117D-03, * 9.8151213074288225D-04, * -4.8508565918349199D-04, * 7.0229342016979338D-05, * 9.8843338786337969D-05, * -1.2258781476941329D-04, * 1.7715537940170311D-04, * -3.8432111542678374D-04, * 6.9287406638899473D-04, * -9.6215482685026843D-04, * 1.0724633285429212D-03, * -1.0107490310802640D-03, * 8.5656148183947472D-04, * -7.0034372117423028D-04, * 5.8169908961642128D-04, * -4.7584007785053466D-04, * 3.3826942897811062D-04, * -1.5353773275122099D-04, * -5.4151363407508753D-05, * 2.3546452586479855D-04, * -3.6268263213064517D-04, * 4.3215561854952479D-04, * -4.7225359925206396D-04, * 5.0071077491085095D-04, * -5.2810307126720791D-04, * 5.3879997165512184D-04, * -5.2429072508634864D-04, * 4.7746113147592690D-04, * -4.1156501972987601D-04, * 3.3737617936914973D-04, * -2.6664562414888755D-04, * 1.9760773520137086D-04, * -1.2774219561801659D-04, * 5.3111364069808147D-05, * 2.2418117761486844D-05, * -9.3058490127275079D-05, * 1.5206775220102640D-04, * -1.9778944912593785D-04, * 2.3125645364529598D-04, * -2.5588711191138775D-04, * 2.7255134527919336D-04, * -2.8156394791930599D-04, * 2.8091798741786317D-04, * -2.7153805696720584D-04, * 2.5388712730440023D-04, * -2.3181002697390863D-04, * 2.0616302161897061D-04, * -1.7949351580032003D-04, * 1.5041439622973747D-04, * -1.2044603366672626D-04, * 8.8484335700220169D-05, * -5.7074250855695914D-05, * 2.5811716813551258D-05, * 3.1212831378025934D-06, * -3.0860530459592881D-05, * 5.6116406459764796D-05, * -7.9983705527276291D-05, * 1.0098500559302100D-04, * -1.1960093772391646D-04, * 1.3447443247980434D-04, * -1.4623882261067433D-04, * 1.5417664364042460D-04, * -1.5907843130015723D-04, * 1.6039759334197207D-04, * -1.5871250617336641D-04, * 1.5352908452742843D-04, * -1.4549308148183567D-04, * 1.3441168281984183D-04, * -1.2113378963441093D-04, * 1.0559250506158825D-04, * -8.8636216845056818D-05, * 7.0152910591523912D-05, * -5.1022037300404419D-05, * 3.1172122896047846D-05, * -1.1586623458569225D-05, * -7.8127742264476840D-06, * 2.6063103148636928D-05, * -4.3369113573024255D-05, * 5.8854697335080722D-05, * -7.2833567033513208D-05, * 8.4499586523322491D-05, * -9.4259703701302294D-05, * 1.0142518906144890D-04, * -1.0652723602532573D-04, * 1.0902876285263386D-04, * -1.0955584249312887D-04, * 1.0769191524705983D-04, * -1.0410549523545201D-04, * 9.8481433501074033D-05, * -9.1510799954605334D-05, * 8.2970720041135862D-05, * -7.3547387909929108D-05, * 6.3075321254199464D-05, * -5.2198843054047043D-05, * 4.0778316780576580D-05, * -2.9401305349487740D-05, * 1.7938512361604202D-05, * -6.9154033960524673D-06, * -3.8059605194935538D-06, * 1.3769824733501799D-05, * -2.3139587426235182D-05, * 3.1526714465770109D-05, * -3.9124840951119910D-05, * 4.5602349678179444D-05, * -5.1184761738136830D-05, * 5.5588971587402965D-05, * -5.9073424276714289D-05, * 6.1394146127086218D-05, * -6.2838260938790675D-05, * 6.3189615966786362D-05, * -6.2757658915071891D-05, * 6.1345600499903273D-05, * -5.9279809079567441D-05, * 5.6376835351230610D-05, * -5.2973498636133322D-05, * 4.8894086086953830D-05, * -4.4478666720221248D-05, * 3.9555349237625895D-05, * -3.4461339733338948D-05, * 2.9026425804978193D-05, * -2.3579379483676993D-05, * 1.7949918005861141D-05, * -1.2453101933878252D-05, * 6.9173707995691686D-06, * -1.6400410906347600D-06, * -3.5520830165088489D-06, * 8.3821835371205464D-06, * -1.3025194262528906D-05, * 1.7226611286525217D-05, * -2.1163086433169281D-05, * 2.4603069500992625D-05, * -2.7724574440183311D-05, * 3.0318490898125594D-05, * -3.2563715583155982D-05, * 3.4272292617305536D-05, * -3.5623562062583836D-05, * 3.6448813652283632D-05, * -3.6927338420384948D-05, * 3.6907264287929438D-05, * -3.6567332483952423D-05, * 3.5769899901081927D-05, * -3.4692740374200804D-05, * 3.3209821645391444D-05, * -3.1497562143912526D-05, * 2.9438996760003930D-05, * -2.7208787761893229D-05, * 2.4696687830740092D-05, * -2.2075196459453580D-05, * 1.9238704812912784D-05, * -1.6357092863577645D-05, * 1.3327576794477666D-05, * -1.0316864582476197D-05, * 7.2234522250248573D-06, * -4.2102336271307681D-06, * 1.1757146407802695D-06, * 1.7217418927952805D-06, * -4.5846336001565605D-06, * 7.2598989191592421D-06, * -9.8518078618850156D-06, * 1.2213369943121936D-05, * -1.4451149531347898D-05, * 1.6424933034479340D-05, * -1.8243859873760788D-05, * 1.9775079698687070D-05, * -2.1130356739336351D-05, * 2.2184617338587232D-05, * -2.3052078352631394D-05, * 2.3615640287166140D-05, * -2.3991597243567400D-05, * 2.4070784464603234D-05, * -2.3971024925964178D-05, * 2.3590810270257545D-05, * -2.3048806477592677D-05, * 2.2250658743958569D-05, * -2.1315094191624169D-05, * 2.0154216322583601D-05, * -1.8885982285553034D-05, * 1.7428110893141435D-05, * -1.5896955028710343D-05, * 1.4214908650754536D-05, * -1.2495925323404714D-05, * 1.0666090143065741D-05, * -8.8362390679320823D-06, * 6.9351725156414336D-06, * -5.0699971208197838D-06, * 3.1713072086607493D-06, * -1.3419955304932869D-06, * -4.8637831480514732D-07, * 2.2154809973191154D-06, * -3.9134940134421737D-06, * 5.4868465836225638D-06, * -7.0040258977260557D-06, * 8.3762889773317768D-06, * -9.6728418575451684D-06, * 1.0809707622241573D-05, * -1.1857072003904072D-05, * 1.2735551147071872D-05, * -1.3516426019082073D-05, * 1.4124640815093744D-05, * -1.4632555441879321D-05, * 1.4969106346564432D-05, * -1.5207602495663466D-05, * 1.5280587719900838D-05, * -1.5262097987506183D-05, * 1.5087871666154388D-05, * -1.4832380542604135D-05, * 1.4434136148723420D-05, * -1.3967708977633646D-05, * 1.3373970797932666D-05, * -1.2727229839850644D-05, * 1.1970327605528348D-05, * -1.1176945219981915D-05, * 1.0291536554486750D-05, * -9.3868040713570609D-06, * 8.4084971567687988D-06, * -7.4280152119576828D-06, * 6.3921309882191903D-06, * -5.3706540262772652D-06, * 4.3111543683032744D-06, * -3.2816095727760150D-06, * 2.2302059799097745D-06, * -1.2228956920870052D-06, * 2.0834261473350782D-07, * 7.4967016653277544D-07, * -1.7021017337602196D-06, * 2.5874326246452715D-06, * -3.4563120890013196D-06, * 4.2495326272220391D-06, * -5.0174725382812542D-06, * 5.7032422889862039D-06, * -6.3569647726289517D-06, * 6.9240368305084930D-06, * -7.4543227409484232D-06, * 7.8954417177083249D-06, * -8.2969814972814828D-06, * 8.6086925962709558D-06, * -8.8798569827016560D-06, * 9.0622436207448857D-06, * -9.2047909836914639D-06, * 9.2611569122703660D-06, * -9.2798923996530969D-06, * 9.2164026355313571D-06, * -9.1188026601643725D-06, * 8.9440959155403429D-06, * -8.7399099367812528D-06, * 8.4646937350325046D-06, * -8.1655338581130703D-06, * 7.8021721790257749D-06, * -7.4210998305296259D-06, * 6.9832019524540892D-06, * -6.5343197827031860D-06, * 6.0363379559731194D-06, * -5.5343941393875378D-06, * 4.9912368030040609D-06, * -4.4512480211749035D-06, * 3.8779144271778528D-06, * -3.3148129910684416D-06, * 2.7260542949412747D-06, * -2.1543640631045786D-06, * 1.5643751494486447D-06, * -9.9792010884904877D-07, * 4.2006563435414488D-07, * 1.2828578144954716D-07, * -6.8170843761355059D-07, * 1.2002610139887423D-06, * -1.7182199973018942D-06, * 2.1966299036873654D-06, * -2.6695172242836803D-06, * 3.0989299789061495D-06, * -3.5186796140196597D-06, * 3.8918279086975399D-06, * -4.2519940815516808D-06, * 4.5632562502661827D-06, * -4.8590529999570710D-06, * 5.1044735870440790D-06, * -5.3327773876714302D-06, * 5.5100518808093951D-06, * -5.6693696441172042D-06, * 5.7777959742651554D-06, * -5.8682012701143240D-06, * 5.9086011384858557D-06, * -5.9316418862829342D-06, * 5.9062553535914082D-06, * -5.8648365729234794D-06, * 5.7771936374521141D-06, * -5.6754390976688660D-06, * 5.5302121973454227D-06, * -5.3733089756651957D-06, * 5.1761504780325406D-06, * -4.9701805272537536D-06, * 4.7275493238061891D-06, * -4.4793121679061271D-06, * 4.1982934698991463D-06, * -3.9151240929232604D-06, * 3.6032464386635171D-06, * -3.2928323134682268D-06, * 2.9578856462504210D-06, * -2.6280866687483951D-06, * 2.2779451338008938D-06, * -1.9366199891095606D-06, * 1.5790728308766014D-06, * -1.2339150240476668D-06, * 8.7650864589825607D-07, * -5.3489508604978068D-07, * 1.8478896775883953D-07, * 1.4635639435064935D-07, * -4.8251763451232018D-07, * 7.9684304272103295D-07, * -1.1130395976443886D-06, * 1.4048659856470465D-06, * -1.6957934464612702D-06, * 1.9601950565939423D-06, * -2.2213347482592139D-06, * 2.4541992282438331D-06, * -2.6818678950516039D-06, * 2.8799357159575715D-06, * -3.0713146413972508D-06, * 3.2321981554388340D-06, * -3.3853422643026042D-06, * 3.5075251705964740D-06, * -3.6213530927788669D-06, * 3.7041714910086520D-06, * -3.7784379549097009D-06, * 3.8220445303250094D-06, * -3.8572967892265748D-06, * 3.8626099785934594D-06, * -3.8601302479287204D-06, * 3.8287704781436971D-06, * -3.7905065696123034D-06, * 3.7247218343364495D-06, * -3.6532083112636991D-06, * 3.5557914543676664D-06, * -3.4540637011080636D-06, * 3.3282638097763323D-06, * -3.1997674083964404D-06, * 3.0491976866974509D-06, * -2.8976954310374269D-06, * 2.7262398320751590D-06, * -2.5557185889768430D-06, * 2.3674393376512546D-06, * -2.1820187956375685D-06, * 1.9810667431748329D-06, * -1.7849118792883693D-06, * 1.5754414045719504D-06, * -1.3726802602058156D-06, * 1.1587701814831071D-06, * -9.5341827783304478D-07, * 7.3899997149845896D-07, * -5.3489242446434118D-07, * 3.2368607480534831D-07, * -1.2441819689464129D-07, * -8.0122169362778756D-08, * 2.7124525811845504D-07, * -4.6597755824709604D-07, * 6.4598204009003610D-07, * -8.2811259130749615D-07, * 9.9438774314769828D-07, * -1.1614984573534462D-06, * 1.3118192957787654D-06, * -1.4618858773831635D-06, * 1.5944270073566040D-06, * -1.7258285853632168D-06, * 1.8391697558539986D-06, * -1.9506905188797842D-06, * 2.0438145113720023D-06, * -2.1346380224091696D-06, * 2.2069215918807961D-06, * -2.2766185392057293D-06, * 2.3278171930685149D-06, * -2.3763273864521875D-06, * 2.4065548270747517D-06, * -2.4341642766088062D-06, * 2.4438673499049261D-06, * -2.4511812710139056D-06, * 2.4411112599981903D-06, * -2.4290238355056873D-06, * 2.4002049167145823D-06, * -2.3698666182256253D-06, * 2.3235622635343284D-06, * -2.2763454929305743D-06, * 2.2140235525588732D-06, * -2.1514873130023505D-06, * 2.0747844602273733D-06, * -1.9986387074602344D-06, * 1.9093248643450292D-06, * -1.8213951257285484D-06, * 1.7213384243582203D-06, * -1.6235312071912660D-06, * 1.5146639745026892D-06, * -1.4089334185606596D-06, * 1.2932196065585708D-06, * -1.1815357700413957D-06, * 1.0609401883973271D-06, * -9.4525929280114757D-07, * 8.2171893861313853D-07, * -7.0395583738270580D-07, * 5.7935355803279976D-07, * -4.6135663688065497D-07/ DATA XC/ * 1.5000000000000000D+00, * -2.5000000000000000D-01, * -8.3333333333333333D-02, * 8.3333333333333333D-02, * -2.5000000000000000D-02, * -1.6666666666666667D-02, * 3.7202380952380952D-03, * -8.9285714285714286D-03, * 4.7743055555555556D-03, * -3.6024305555555556D-03, * 1.3139204545454545D-03, * -8.5819128787878788D-04, * -5.1958133012820513D-04, * 1.8243475274725275D-04, * -9.1610863095238095D-04, * 3.8859049479166667D-04, * -7.3960248161764706D-04, * 2.8862049377042484D-04, * -4.5872292323419225D-04, * 1.3188813862047697D-04, * -2.3489452543712798D-04, * 9.2291728758708739D-06, * -9.4072620859259202D-05, * -6.1795331429744112D-05, * -1.9372304280598958D-05, * -9.0631888462946965D-05, * 1.2482163573262359D-05, * -9.2516382219930174D-05, * 2.0188210633000717D-05, * -8.0470929200621857D-05, * 1.6256654134360693D-05, * -6.3367668659456315D-05, * 8.1525790984883453D-06, * -4.6346754743759114D-05, * -1.4321588981552284D-07, * -3.1872967162006904D-05, * -6.8329445694016667D-06, * -2.0745243096543982D-05, * -1.1375072067856409D-05, * -1.2853785851919570D-05, * -1.3878858851789917D-05, * -7.6769118292931412D-06, * -1.4735024535702190D-05, * -4.5740007331399651D-06, * -1.4406518706515143D-05, * -2.9388752456938436D-06, * -1.3321627343044042D-05, * -2.2663516268149185D-06, * -1.1828848098409159D-05, * -2.1694106070422280D-06, * -1.0186558416514643D-05, * -2.3712489967138187D-06, * -8.5703843892047465D-06, * -2.6867826527376113D-06, * -7.0879904955530762D-06, * -3.0016850027999601D-06, * -5.7955023298229911D-06, * -3.2529932953964178D-06, * -4.7126105819518746D-06, * -3.4129390311560785D-06, * -3.8351096536068256D-06, * -3.4763627127104309D-06, * -3.1445867953274529D-06, * -3.4514343519799726D-06, * -2.6154719265328675D-06, * -3.3531417305725839D-06, * -2.2198707506399539D-06, * -3.1989537339969685D-06, * -1.9306564544852060D-06, * -3.0061117071933445D-06, * -1.7232651018572899D-06, * -2.7900881151033900D-06, * -1.5765731817439904D-06, * -2.5638468441384344D-06, * -1.4731590908883593D-06, * -2.3376275814146932D-06, * -1.3991774848455102D-06, * -2.1190514435149442D-06, * -1.3440126905230122D-06, * -1.9134049394077030D-06, * -1.2998266484871192D-06, * -1.7240054265459364D-06, * -1.2610777447610126D-06, * -1.5525854546752844D-06, * -1.2240579248803742D-06, * -1.3996580399100217D-06, * -1.1864748231422857D-06, * -1.2648420962176947D-06, * -1.1470914767032236D-06, * -1.1471388075454853D-06, * -1.1054269237433427D-06, * -1.0451571528820220D-06, * -1.0615152769832697D-06, * -9.5729128101001251D-07, * -1.0157176646253370D-06, * -8.8185487527751804D-07, * -9.6857993527455093D-07, * -8.1717872457679041D-07, * -9.2072863551311475D-07, * -7.6167791445146884D-07, * -8.7279805997655628D-07, * -7.1389471938613082D-07, * -8.2538184496980256D-07, * -6.7252265410104636D-07, * -7.7900342609489841D-07, * -6.3641638900592534D-07, * -7.3410057737591757D-07, * -6.0459145638444036D-07, * -6.9102011395384808D-07, * -5.7621693226497938D-07, * -6.5001962769024342D-07, * -5.5060360907445426D-07, * -6.1127381366185704D-07, * -5.2718959201575184D-07, * -5.7488352950028140D-07, * -5.0552476142304289D-07, * -5.4088621213339348D-07, * -4.8525514044302298D-07, * -5.0926666630196989D-07, * -4.6610788438296436D-07, * -4.7996754738889920D-07, * -4.4787735508370770D-07, * -4.5289909953003837D-07, * -4.3041255027855705D-07, * -4.2794789029574168D-07, * -4.1360601396823146D-07, * -4.0498441614351475D-07, * -3.9738425008413894D-07, * -3.8386954782481638D-07, * -3.8169958991035488D-07, * -3.6445985013599408D-07, * -3.6652341683194657D-07, * -3.4661185268364758D-07, * -3.5184062405580626D-07, * -3.3018537332733885D-07, * -3.3764516714887495D-07, * -3.1504600822737862D-07, * -3.2393656963962651D-07, * -3.0106690558756812D-07, * -3.1071724343142599D-07, * -2.8812993707616886D-07, * -2.9799049397238926D-07, * -2.7612637370990722D-07, * -2.8575909118146931D-07, * -2.6495716333596472D-07, * -2.7402429967868077D-07, * -2.5453289595492654D-07, * -2.6278527491802908D-07, * -2.4477353186135226D-07, * -2.5203874468028553D-07, * -2.3560795653942291D-07, * -2.4177890758434017D-07, * -2.2697341583592807D-07, * -2.3199749153049563D-07, * -2.1881487538432363D-07, * -2.2268392513733099D-07, * -2.1108433970027574D-07, * -2.1382558421139263D-07, * -2.0374015885727821D-07, * -2.0540808309892139D-07, * -1.9674634416933435D-07, * -1.9741558745967362D-07, * -1.9007190880658811D-07, * -1.8983113065232451D-07, * -1.8369024467534122D-07, * -1.8263692062340961D-07, * -1.7757854311886682D-07, * -1.7581462804944016D-07, * -1.7171726394694966D-07, * -1.6934564959817202D-07, * -1.6608965488710954D-07, * -1.6321134265031496D-07, * -1.6068132167951692D-07, * -1.5739322975148737D-07, * -1.5547984762715441D-07, * -1.5187317253241943D-07, * -1.5047446038661908D-07, * -1.4663351592108556D-07, * -1.4565574307490746D-07, * -1.4165720424246699D-07, * -1.4101538631353572D-07, * -1.3692787132008801D-07, * -1.3654597758131755D-07, * -1.3242990701006942D-07, * -1.3224082415643738D-07, * -1.2814850275720800D-07, * -1.2809380595940030D-07, * -1.2406967880051353D-07, * -1.2409925472952977D-07, * -1.2018029560344396D-07, * -1.2025185615297465D-07, * -1.1646805196699325D-07, * -1.1654657178858863D-07, * -1.1292147212223454D-07, * -1.1297857789266597D-07, * -1.0952988390922156D-07, * -1.0954321851105652D-07, * -1.0628338994412838D-07, * -1.0623597047737109D-07, * -1.0317283346605371D-07, * -1.0305241822109372D-07, * -1.0018976034648236D-07, * -9.9988236543799000D-08, * -9.7326378543431607D-08, * -9.7039179761399905D-08, * -9.4575516092632027D-08, * -9.4201075832868931D-08, * -9.1930578552213454D-08, * -9.1469824299717658D-08, * -8.9385506656765769D-08, * -8.8841397045075188D-08, * -8.6934734791997956D-08, * -8.6311841046510650D-08, * -8.4573150772591081D-08, * -8.3877282443320366D-08, * -8.2296057292842038D-08, * -8.1533931367714647D-08, * -8.0099135321609721D-08, * -7.9278087101297707D-08, * -7.7978409628966224D-08, * -7.7106143214691097D-08, * -7.5930216560745668D-08, * -7.5014592430439758D-08, * -7.3951174117720153D-08, * -7.3000031018819540D-08, * -7.2038154347238563D-08, * -7.1059162594179937D-08, * -7.0188258015674040D-08, * -6.9188801227305295D-08, * -6.8398791498814021D-08, * -6.7385873828195121D-08, * -6.6667245803325677D-08, * -6.5647421784802526D-08, * -6.4991277614644086D-08, * -6.3970601867693835D-08, * -6.3368692254144343D-08, * -6.2352686429273421D-08, * -6.1797428420439073D-08, * -6.0791062940033853D-08, * -6.0275544585339958D-08, * -5.9283232914028153D-08, * -5.8801206913768488D-08, * -5.7826810282115668D-08, * -5.7372678578108289D-08, * -5.6419519275119202D-08, * -5.5988310340644517D-08, * -5.5059191880383935D-08, * -5.4646532282389241D-08, * -5.3743764934811702D-08, * -5.3345846562363280D-08, * -5.2471276915654930D-08, * -5.2084821097970111D-08, * -5.1239864487531243D-08, * -5.0862084064183349D-08, * -5.0047758860547568D-08, * -4.9676319116649699D-08, * -4.8893282010339343D-08, * -4.8526261251298487D-08, * -4.7774842806432328D-08, * -4.7410693220497277D-08, * -4.6690933090781493D-08, * -4.6328442433082258D-08, * -4.5640123743761498D-08, * -4.5278378272630326D-08, * -4.4621060770377520D-08, * -4.4259409775058572D-08, * -4.3632461435111174D-08, * -4.3270483612987034D-08, * -4.2673110469672002D-08, * -4.2310582340249247D-08, * -4.1741856374031877D-08, * -4.1378722855463075D-08, * -4.0837607827505194D-08, * -4.0473955048673295D-08, * -3.9959330223318290D-08, * -3.9595360599748213D-08, * -3.9106042337094096D-08, * -3.8742051901463031D-08, * -3.8276813136962395D-08, * -3.7913171084046150D-08, * -3.7470758740586109D-08, * -3.7107889121418532D-08, * -3.6687039522259333D-08, * -3.6325405002441209D-08, * -3.5924857371369606D-08, * -3.5564944953224489D-08, * -3.5183453101909043D-08, * -3.4825761698968171D-08, * -3.4462104011348841D-08, * -3.4107133755919516D-08, * -3.3760121586040859D-08, * -3.3408364745879151D-08, * -3.3076849349359462D-08, * -3.2728782727278716D-08, * -3.2411660848027976D-08, * -3.2067739538221029D-08, * -3.1763957771468430D-08, * -3.1424610148036513D-08, * -3.1133168198552971D-08, * -3.0798792014889651D-08, * -3.0518744965803691D-08, * -3.0189704447786604D-08, * -2.9920164150868398D-08, * -2.9596787972008316D-08, * -2.9336923664978486D-08, * -2.9019503697540064D-08, * -2.8768541948057489D-08, * -2.8457332690504272D-08, * -2.8214556760182684D-08, * -2.7909775347943239D-08, * -2.7674524063195580D-08, * -2.7376350776555535D-08, * -2.7148016986399944D-08, * -2.6856596176176168D-08, * -2.6634624870468772D-08, * -2.6350066228917123D-08, * -2.6133952383896034D-08, * -2.5856332494961105D-08, * -2.5645618706567583D-08, * -2.5374982816035960D-08, * -2.5169256775282048D-08, * -2.4905620727597794D-08, * -2.4704512586320967D-08, * -2.4447864880724068D-08, * -2.4251044550443114D-08, * -2.4001348474669560D-08, * -2.3808522895956743D-08, * -2.3565718700973282D-08, * -2.3376629115801708D-08, * -2.3140636199927398D-08, * -2.2955055454848270D-08, * -2.2725774530133752D-08, * -2.2543504433888348D-08, * -2.2320819651782818D-08, * -2.2141688407056103D-08, * -2.1925469424196407D-08, * -2.1749329149666465D-08, * -2.1539433118081470D-08, * -2.1366157473701331D-08, * -2.1162430942849608D-08, * -2.0991912868402824D-08, * -2.0794193589266841D-08, * -2.0626343163650560D-08, * -2.0434461787611940D-08, * -2.0269204214004967D-08, * -2.0082985881440040D-08, * -1.9920259601491197D-08, * -1.9739525416971946D-08, * -1.9579280355377977D-08, * -1.9403848748058944D-08, * -1.9246044687373134D-08, * -1.9075732656608325D-08, * -1.8920337740812633D-08, * -1.8754961988296315D-08, * -1.8601951352563185D-08, * -1.8441329303342698D-08, * -1.8290683826490271D-08, * -1.8134634542075080D-08, * -1.7986339717464220D-08, * -1.7834684704970217D-08, * -1.7688729624987314D-08, * -1.7541293546825009D-08, * -1.7397669995625320D-08, * -1.7254281284680291D-08, * -1.7112982933517984D-08, * -1.6973474319096174D-08, * -1.6834496018325357D-08, * -1.6698704968358081D-08, * -1.6562042130041000D-08, * -1.6429811215177430D-08, * -1.6295459280169746D-08, * -1.6166636465439770D-08, * -1.6034590448827191D-08, * -1.5909029318545826D-08, * -1.5779283427371209D-08, * -1.5656843348886838D-08, * -1.5529390666222858D-08, * -1.5409936897994666D-08, * -1.5284769127575697D-08, * -1.5168172876908833D-08, * -1.5045280142729232D-08, * -1.4931418578306855D-08, * -1.4810789273814324D-08, * -1.4699545497950466D-08, * -1.4581166179706517D-08, * -1.4472429165008421D-08, * -1.4356284485947554D-08, * -1.4249948980826224D-08, * -1.4136021658516485D-08, * -1.4031988065724076D-08, * -1.3920258881309835D-08, * -1.3818433113416446D-08, * -1.3708880937205795D-08, * -1.3609174252659547D-08, * -1.3501776092600512D-08, * -1.3404104915746716D-08, * -1.3298835985315650D-08, * -1.3203121713485802D-08, * -1.3099955515785639D-08, * -1.3006124316307229D-08, * -1.2905032741440677D-08, * -1.2813015341166145D-08, * -1.2713968774207907D-08, * -1.2623700243916946D-08, * -1.2526667681058239D-08, * -1.2438087216853289D-08, * -1.2343036387530415D-08, * -1.2256087091121484D-08, * -1.2162984584167140D-08, * -1.2077613243729718D-08, * -1.1986424635800585D-08, * -1.1902581508889849D-08, * -1.1813271493626474D-08, * -1.1730910093442534D-08, * -1.1643442610007328D-08, * -1.1562519496130038D-08, * -1.1476857855946463D-08, * -1.1397332430494307D-08, * -1.1313439441174951D-08, * -1.1235273751190643D-08, * -1.1153111836794176D-08, * -1.1076270383519616D-08, * -1.0995801700416857D-08, * -1.0920251255991655D-08, * -1.0841437803749433D-08, * -1.0767147235750060D-08, * -1.0689950962558769D-08, * -1.0616891066688984D-08, * -1.0541273968966012D-08, * -1.0469417310113191D-08, * -1.0395341526010399D-08, * -1.0324662287796187D-08, * -1.0252090184425730D-08, * -1.0182564027302566D-08, * -1.0111458281572171D-08, * -1.0043062209449159D-08, * -9.9733858824660883D-09, * -9.9060981177878703D-09, * -9.8378147228506466D-09, * -9.7716145900008379D-09, * -9.7046881542500855D-09, * -9.6395559711058906D-09, * -9.5739510909507708D-09, * -9.5098680683771460D-09, * -9.4455499588524434D-09, * -9.3824981078920212D-09, * -9.3194326461334610D-09, * -9.2573946926219519D-09, * -9.1955484556743042D-09, * -9.1345077619897345D-09, * -9.0738480591841733D-09/ DATA YC/ * 4.7200000000000000D-01, * -1.3600000000000000D-01, * -1.1200000000000000D-01, * 2.2316666666666667D-01, * -7.3600000000000000D-02, * -6.2400000000000000D-02, * 1.0215476190476190D-01, * -4.7714285714285714D-02, * -3.0947916666666667D-02, * 6.8111675347222222D-02, * -3.9411425781250000D-02, * -1.8253071178089489D-02, * 5.6571238640027168D-02, * -4.6609217993505708D-02, * 6.6346292387161936D-03, * 2.9116390173913290D-02, * -3.6209065064597834D-02, * 1.7082143957057448D-02, * 7.2223386813957597D-03, * -1.6829550212806958D-02, * 9.0308070502676157D-03, * 5.3125869191405562D-03, * -1.1961801174688067D-02, * 6.6254327977171207D-03, * 4.9209579093679157D-03, * -1.2557766563381617D-02, * 1.1221342023251912D-02, * -3.3058868572595320D-03, * -4.3781511286035770D-03, * 6.9378942454524954D-03, * -4.2668979836015938D-03, * 2.6465507967097246D-04, * 1.4770745346071534D-03, * -6.4876376275806497D-06, * -2.5435181163049142D-03, * 3.6728741033678549D-03, * -2.3199942773761088D-03, * -3.3218980738680673D-04, * 2.4010107921624552D-03, * -2.6997127838461383D-03, * 1.6107792375962346D-03, * -3.1978392960669817D-04, * -1.4462864800190784D-04, * -2.5714426547225104D-04, * 9.0366080756162712D-04, * -1.0613081862189289D-03, * 5.9777064256183724D-04, * 1.4368562658462104D-04, * -6.0775656522722038D-04, * 6.0165188281673950D-04, * -2.7144370397287921D-04, * 8.5155721169875297D-06, * -1.6439045383982801D-05, * 2.7537985585513706D-04, * -5.3929022410230438D-04, * 6.3587995404547473D-04, * -5.2279463148163987D-04, * 3.3883728529370889D-04, * -2.0715390621295863D-04, * 1.8762269032536883D-04, * -2.0503971254850025D-04, * 1.8045621043936247D-04, * -5.8940519057838148D-05, * -1.1033071557612852D-04, * 2.6639026156049494D-04, * -3.4455089067637721D-04, * 3.5848138996003487D-04, * -3.3493459632205689D-04, * 3.2327758866596014D-04, * -3.2051885809497983D-04, * 3.1876721731771198D-04, * -2.8495011491961848D-04, * 2.2367294042628626D-04, * -1.4107283489882928D-04, * 6.8450984714095897D-05, * -1.0697058599387422D-05, * -2.6536960943398819D-05, * 6.1603193134024803D-05, * -9.5798723182375631D-05, * 1.3508243318683469D-04, * -1.6553733994519131D-04, * 1.8615360206642408D-04, * -1.8969031500893858D-04, * 1.8623971789173978D-04, * -1.7597455776753032D-04, * 1.6710819924613960D-04, * -1.5328602284743927D-04, * 1.3705156651218424D-04, * -1.1230880829190950D-04, * 8.5731676090384513D-05, * -5.5923510080196872D-05, * 3.0443826502259130D-05, * -5.7291038137283416D-06, * -1.4365085064233328D-05, * 3.4973661522325679D-05, * -5.1929992822539663D-05, * 6.8370644318674676D-05, * -7.9367323143204268D-05, * 8.8311242877927294D-05, * -9.2067994111711346D-05, * 9.5032434822759751D-05, * -9.4523258931762346D-05, * 9.3957892563391577D-05, * -9.0029693845736955D-05, * 8.5632565369465131D-05, * -7.8025557136697927D-05, * 7.0483721026399589D-05, * -6.0741048558321154D-05, * 5.1743935360352776D-05, * -4.0989633879045605D-05, * 3.1009969358956206D-05, * -1.9467463402090508D-05, * 9.0141907826644552D-06, * 2.3540698701376435D-06, * -1.2068486418547503D-05, * 2.2186714142104905D-05, * -3.0431245911807810D-05, * 3.8862500961161882D-05, * -4.5269950454951726D-05, * 5.1555708095830356D-05, * -5.5582554878339440D-05, * 5.9239826295566304D-05, * -6.0603879113242124D-05, * 6.1586742678943573D-05, * -6.0405043779777482D-05, * 5.8894704412044139D-05, * -5.5367430797640363D-05, * 5.1613758279034810D-05, * -4.6098937146643610D-05, * 4.0600788648012361D-05, * -3.3711453531294057D-05, * 2.7140109558764976D-05, * -1.9549246415857849D-05, * 1.2562104528696746D-05, * -4.9125868815033165D-06, * -1.8491507669988706D-06, * 8.9323261738307704D-06, * -1.4882318972410812D-05, * 2.0887807414564620D-05, * -2.5602166966663558D-05, * 3.0199128245651781D-05, * -3.3424912172676970D-05, * 3.6435130777413034D-05, * -3.8062975668282837D-05, * 3.9454124516484501D-05, * -3.9525334694418640D-05, * 3.9410041448127110D-05, * -3.8092378184912795D-05, * 3.6683156413239050D-05, * -3.4220362775234927D-05, * 3.1788354379800357D-05, * -2.8471708285321236D-05, * 2.5324244070303637D-05, * -2.1466967718730786D-05, * 1.7915686280047773D-05, * -1.3818873178524654D-05, * 1.0150688100907209D-05, * -6.0833219917725053D-06, * 2.5484110737538267D-06, * 1.2611824376403452D-06, * -4.4578959625590546D-06, * 7.8309560089204244D-06, * -1.0536772956905408D-05, * 1.3346970707661510D-05, * -1.5459841422425887D-05, * 1.7629058235227153D-05, * -1.9092888001657933D-05, * 2.0587204099827716D-05, * -2.1388126404741792D-05, * 2.2213308631503426D-05, * -2.2373736707215305D-05, * 2.2568136650874861D-05, * -2.2139412733381035D-05, * 2.1767083026392473D-05, * -2.0823203020527720D-05, * 1.9967792324524959D-05, * -1.8599050767582485D-05, * 1.7357067408540357D-05, * -1.5663269303262007D-05, * 1.4137623723362232D-05, * -1.2222219613174739D-05, * 1.0516872846243478D-05, * -8.4819539822332596D-06, * 6.6971134120912273D-06, * -4.6388803751578484D-06, * 2.8668958536040440D-06, * -8.7214084509234568D-07, * -8.0540998219700290D-07, * 2.6616947619732297D-06, * -4.1759842884295163D-06, * 5.8322239138527139D-06, * -7.1284091791836209D-06, * 8.5372660003782969D-06, * -9.5748819029134785D-06, * 1.0703285738757720D-05, * -1.1455982922812851D-05, * 1.2284661564373073D-05, * -1.2739551341060177D-05, * 1.3262212826431420D-05, * -1.3418877928645660D-05, * 1.3641086453692390D-05, * -1.3510345295768775D-05, * 1.3448223826224468D-05, * -1.3050779429342174D-05, * 1.2729618524759661D-05, * -1.2094609275792830D-05, * 1.1547364624877339D-05, * -1.0710804203303536D-05, * 9.9764949530084980D-06, * -8.9796309024329404D-06, * 8.1016651522528682D-06, * -6.9892667690250671D-06, * 6.0136947319747133D-06, * -4.8322770389790563D-06, * 3.8059990124700243D-06, * -2.6020315861256727D-06, * 1.5710225729374221D-06, * -3.8918714764938743D-07, * -6.0319989120806344D-07, * 1.7216407909531038D-06, * -2.6362253102094234D-06, * 3.6548746075971186D-06, * -4.4579756221764938D-06, * 5.3464902055195090D-06, * -6.0109666796824807D-06, * 6.7459355380831763D-06, * -7.2518967214864909D-06, * 7.8173737038249186D-06, * -8.1525329605102430D-06, * 8.5402084505331727D-06, * -8.6998774807637719D-06, * 8.9088990782522834D-06, * -8.8956444761987508D-06, * 8.9321184462592334D-06, * -8.7551204499319011D-06, * 8.6313409836875282D-06, * -8.3055082410918161D-06, * 8.0389743460070626D-06, * -7.5838792034685149D-06, * 7.1961664681550028D-06, * -6.6348690785502634D-06, * 6.1504242161509896D-06, * -5.5082522472003163D-06, * 4.9531752926514037D-06, * -4.2565215890650059D-06, * 3.6573945620649294D-06, * -2.9325872687407283D-06, * 2.3153987219905281D-06, * -1.5876879545868273D-06, * 9.7689192680490916D-07, * -2.6958598641107617D-07, * -3.1267732384611541D-07, * 9.7890455506493213D-07, * -1.5134105009755640D-06, * 2.1210306278610715D-06, * -2.5918891526227867D-06, * 3.1268486333309461D-06, * -3.5217450094985069D-06, * 3.9736213483479260D-06, * -4.2838909323984975D-06, * 4.6458913813851200D-06, * -4.8664469569371721D-06, * 5.1352707026269825D-06, * -5.2644052733228976D-06, * 5.4399934081079628D-06, * -5.4790837086481074D-06, * 5.5642828400148302D-06, * -5.5174195642897539D-06, * 5.5175848568138316D-06, * -5.3911547978495994D-06, * 5.3137167716086096D-06, * -5.1159600268227424D-06, * 4.9699769247120209D-06, * -4.7105394303592488D-06, * 4.5062537749029603D-06, * -4.1957520171585079D-06, * 3.9441664103620917D-06, * -3.5937775286095240D-06, * 3.3062611246005868D-06, * -2.9273480790090059D-06, * 2.6152817427335702D-06, * -2.2190599745428558D-06, * 1.8935250950353540D-06, * -1.4907742948973469D-06, * 1.1622876555593314D-06, * -7.6310995491356790D-07, * 4.4140502655962355D-07, * -5.5029153912126644D-08, * -2.5111734940967919D-07, * 6.1648769387587612D-07, * -8.9938622789894970D-07, * 1.2366823950992479D-06, * -1.4898208889766730D-06, * 1.7931808471101548D-06, * -2.0112782149985639D-06, * 2.2760830856489322D-06, * -2.4551094504773397D-06, * 2.6779885509182654D-06, * -2.8151557338413330D-06, * 2.9939636642066399D-06, * -3.0876894613461062D-06, * 3.2214587198112625D-06, * -3.2713084165799760D-06, * 3.3601809453569366D-06, * -3.3667894325682067D-06, * 3.4119304095991562D-06, * -3.3769081762694265D-06, * 3.3804052790728855D-06, * -3.3062314676700500D-06, * 3.2709827467124652D-06, * -3.1608883645318090D-06, * 3.0904817684585056D-06, * -2.9483260473258055D-06, * 2.8469135344824410D-06, * -2.6770563134076716D-06, * 2.5492253570139198D-06, * -2.3563982229956140D-06, * 2.2070433659841856D-06, * -1.9962221227220877D-06, * 1.8304190596019048D-06, * -1.6066999007048155D-06, * 1.4295843574004248D-06, * -1.1980659004870506D-06, * 1.0147193506590811D-06, * -7.8039244417508695D-07, * 5.9573644553695755D-07, * -3.6338339544942316D-07, * 1.8208405705436667D-07, * 4.3811358654923739D-08, * -2.1742755140477999D-07, * 4.3275820027971250D-07, * -5.9477626588591796D-07, * 7.9587712051163951D-07, * -9.4285251903762870D-07, * 1.1265529550454999D-06, * -1.2555572182823719D-06, * 1.4192199251416351D-06, * -1.5278728862802534D-06, * 1.6694194297568264D-06, * -1.7559082164142414D-06, * 1.8738315394265510D-06, * -1.9369167325609752D-06, * 2.0302811095924411D-06, * -2.0692906890502266D-06, * 2.1377198563258234D-06, * -2.1525317498091425D-06, * 2.1961861176634411D-06, * -2.1872003425123203D-06, * 2.2067443566320520D-06, * -2.1748458919235503D-06, * 2.1714067454387124D-06, * -2.1179203945289214D-06, * 2.0930394022574137D-06, * -2.0196780020074222D-06, * 1.9752560304875273D-06, * -1.8840634319661995D-06, * 1.8223018331902294D-06, * -1.7155921186196676D-06, * 1.6389306406639579D-06, * -1.5192250519313482D-06, * 1.4302781952796975D-06, * -1.3002412300054248D-06, * 1.2017344840597395D-06, * -1.0641105660649297D-06, * 9.5881789654856439D-07, * -8.1636994941284314D-07, * 7.0705381531822510D-07, * -5.6250496222429055D-07, * 4.5186002275356826D-07, * -3.0783950556355396D-07, * 1.9844103619998866D-07, * -5.7435295191808042D-08, * -4.8306828466531523D-08, * 1.8399714124635946D-07, * -2.8387820830028111D-07, * 4.1217468956713764D-07, * -5.0422816804346872D-07, * 6.2330419221138620D-07, * -7.0582768159421093D-07, * 8.1413078061422650D-07, * -8.8570477784694097D-07, * 9.8197183780878940D-07, * -1.0414713709357292D-06, * 1.1247367848164764D-06, * -1.1713361351709300D-06, * 1.2409332101001310D-06, * -1.2741040940772127D-06, * 1.3296601519843084D-06, * -1.3491638634338323D-06, * 1.3905895928735762D-06, * -1.3964637144081396D-06, * 1.4239374266105891D-06, * -1.4164778010299042D-06, * 1.4304253136714822D-06, * -1.4101640246236637D-06, * 1.4112349423111588D-06, * -1.3789150865504903D-06, * 1.3679562682934156D-06, * -1.3245043116034164D-06, * 1.3025313141522373D-06, * -1.2490278110557519D-06, * 1.2171950750850317D-06, * -1.1548445013074571D-06, * 1.1144150077010837D-06, * -1.0445154068524887D-06, * 9.9683047581508093D-07, * -9.2074356103063131D-07, * 8.6719340059387635D-07, * -7.8631568113950727D-07, * 7.2831121707251157D-07, * -6.4404664236551291D-07, * 5.8299308165090120D-07, * -4.9672761378583376D-07, * 4.3400011159494443D-07, * -3.4707855507171352D-07, * 2.8400027317506864D-07, * -1.9770560953858427D-07, * 1.3552837458030817D-07, * -5.1063772166317441D-08, * -9.0485329171626414D-09, * 9.0574936309349525D-08, * -1.4756018376763088D-07, * 2.2514769342019229D-07, * -2.7805652352873536D-07, * 3.5082152238869002D-07, * -3.9882561094144452D-07, * 4.6600849392616066D-07, * -5.0840617454477740D-07, * 5.6937570310101832D-07, * -6.0559508338271450D-07, * 6.5985031573376915D-07, * -6.8945005729491105D-07, * 7.3662003706719416D-07, * -7.5928799314453127D-07, * 7.9912939915137918D-07, * -8.1467932014477702D-07, * 8.4707229362763449D-07, * -8.5543882146275503D-07, * 8.8038119474227483D-07, * -8.8161337192153831D-07, * 8.9921352493502257D-07, * -8.9346704438334281D-07, * 9.0393561373631633D-07, * -8.9146403178686043D-07, * 8.9510469146720371D-07, * -8.7624981929808703D-07, * 8.7344934379631114D-07, * -8.4863102300222159D-07, * 8.3984883290277030D-07, * -8.0955428928635624D-07, * 7.9531166702840931D-07/ END MODULE RIG_WALL_COEF C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE WALL2_poiseu2vel(Utrans,Orot,Z,HH) USE SIZE ! NN,LLG USE MULT ! LL IMPLICIT NONE REAL*8 Utrans,Orot REAL*8 Z,HH INTEGER NG REAL*8 KUU REAL*8 A(2,2,0:1),CONF(1),UX(1),OY(1) CHARACTER*1 LUB_WALL REAL*8, ALLOCATABLE :: G(:,:) REAL*8 FX(1),TY(1) INTEGER M NN=1 LL=40 NG=50 KUU=7.D0 LUB_WALL ='Y' ! 'Y' , 'y' or 'N' , 'n' , wall lubrication CALL INIT_LL_NN_polym CALL ZINIT CALL GINIT_polym CALL CINIT_polym CONF(1)=Z CALL LUB2WALL_polym(A,CONF,HH,LUB_WALL) ALLOCATE( G(LLG,LLG) ) M=1 CALL G2WALL_DEF_polym(G,CONF,HH,KUU,NG,.FALSE.,M) CALL MATREV(G,LLG,'M') CALL EXTR_G_polym(A(1,1,M),G) CALL MATREV(A(1,1,M),2*NN,'M') CALL FORCES_POISEU_polym(FX,TY,G,CONF,HH) CALL VELOCITIES_polym(UX,OY,A,FX,TY) Utrans=ux(1) Orot= OY(1) DEALLOCATE(G) RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE INIT_LL_NN_polym USE SIZE ! NN,LLG USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 USE MULT_WALL_polym ! LC,LC2,LCG1 USE Z_COEF USE S_COEF_polym USE FACT USE C_COEF_polym IMPLICIT NONE INTEGER LL_OLD,NN_OLD DATA LL_OLD,NN_OLD /-1,-1/ SAVE LL_OLD,NN_OLD IF( LL_OLD.NE.LL) THEN ! --------------------- START OF NEW LL C WRITE(*,*) 'NEW LL=',LL LL2 =2*LL LLG1=3*1 *LL LLG =3*NN*LL LC =LL+2 LC2 =LL+LC LCG1=3*1 *LC LL_OLD=LL ENDIF !---------------------------------------- END OF NEW LL IF( NN_OLD.NE.NN) THEN ! --------------------- START OF NEW NN C WRITE(*,*) 'NEW NN=',NN NN_OLD=NN LLG=3*NN*LL ENDIF RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE ZINIT USE MULT ! LL USE Z_COEF IMPLICIT NONE INTEGER L REAL*8 A0,A2,A1,B,D,ZZ INTEGER :: LL_OLD=-1 SAVE LL_OLD IF( LL_OLD.EQ.LL) RETURN C WRITE(*,*) 'NEW ZINIT' LL_OLD=LL IF(ALLOCATED(Z0)) DEALLOCATE( Z0,Z1,Z2,ZB ) ALLOCATE( Z0(LL),Z1(LL),Z2(LL),ZB(LL) ) ZZ=0.D0 DO 1 L=1,LL C ---------------------------------------------------------------------- A0=2.D0*(2*L+1)*L*(2*L-1)*(2*L+1)*(1-ZZ)/(2.D0**(2*L)*(L+1)* *(1+2*(L-1)*ZZ)) C ---------------------------------------------------------------------- A1=.5D0*(2*L+1)*L*(L+1)*(1-(L+2)*ZZ)/((1+(L-1)*ZZ)*2.D0**(2*L)) C ---------------------------------------------------------------------- A2=.25D0*(2*L+1)*(2*L-1)*(2*L+1)*(2*L+3)*(1-3*ZZ)/ * ((1+2*(L-1)*ZZ)*2.D0**(2*L)) C ---------------------------------------------------------------------- B =.25D0*(2*L+1)*(L+1)*(2*L+1)*(2*L+1)*(2*L+3)*(2*L+1)*(1-5*ZZ)/ * (2.D0**(2*L)*8*L*(1+2*(L-1)*ZZ)) C ----------------------------------------------------------------------- Z1(L)= 1.D0/A1 D=1.D0/(A0*B-A2*A2) Z0(L)= D*B Z2(L)=-D*A2 ZB(L)= D*A0 1 CONTINUE RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE GINIT_polym USE MULT ! LL USE MULT_WALL_polym ! LC,LC2,LCG1 USE S_COEF_polym USE FACT IMPLICIT NONE REAL*8 S11,So00,S00,S01,S02,S10,S20, SM REAL*8 NLM INTEGER L1,M1,L,I,M INTEGER :: LL_OLD=-1 SAVE LL_OLD C LOCAL FUNCTIONS C ---------------------------------------------------------------------- S11(L1,M1,L,M)=(-1)**(L1+M1+1)*SIL(L+L1-M+M1)/ *(SIL(L-M)*SIL(L1+M1)*(L+1)*(2*L+1)*(L1+1)) C ---------------------------------------------------------------------- C comment; SM=S11(L1,M1,L,M), SM is used in the following definitions C --------------------------------------------------------------------- So00(L1,M1,L,M)=.5D0*(L+1)*(L1+1)*SM/(2*L+2*L1-1) C --------------------------------------------------------------------- S00(L1,M1,L,M)=SM*(L+1)*(L1+1)*(L*L1-2*L-2*L1+1)* * (-L*L1*(L+L1)+2*M1*M1*L*L+2*M*M*L1*L1+(4*M*M1+1)*L*L1- * M1*(2*M+M1)*L-M*(2*M1+M)*L1+M*M1)/ * (1.D0*L*L1*(2*L-1)*(2*L1-1)*(2*L+2*L1-1)*(L+L1-M+M1)* * (L+L1-M+M1-1)) C --------------------------------------------------------------------- S01(L1,M1,L,M)=SM*(M1*L+M*L1)*(L1+1)/(1.D0*L*L1*(L+L1-M+M1)) C --------------------------------------------------------------------- S02(L1,M1,L,M)=-SM*L*(L1+1)/((2*L+1)*(2*L+3)) C ---------------------------------------------------------------------- S10(L1,M1,L,M)=-SM*(M1*L+M*L1)*(L+1)/(1.D0*L*L1*(L+L1-M+M1)) C ---------------------------------------------------------------------- S20(L1,M1,L,M)=-SM*(L+1)*L1/((2*L1+1)*(2*L1+3)) C ---------------------------------------------------------------------- IF( LL_OLD.EQ.LL) RETURN C WRITE(*,*) 'NEW GINIT' LL_OLD=LL IF(ALLOCATED(S11P)) THEN DEALLOCATE( S11P,So00P,S00P, * S01P, S02P,S10P, * S20P ) DEALLOCATE( SIL ) ENDIF ALLOCATE( S11P(LL,LC,0:1),So00P(LL,LC,0:1),S00P(LL,LC,0:1), * S01P(LL,LC,0:1), S02P(LL,LC,0:1),S10P(LL,LC,0:1), * S20P(LL,LC,0:1) ) ALLOCATE( SIL(0:2*LC2) ) SIL(0)=1.D0 DO I=1,2*LC2 SIL(I)=I*SIL(I-1) !!!!!!!!!!!!! SILNIA CCC SIL(I)=I*SIL(I-1)/750 ENDDO DO M=0,1 M1=M DO 20 L =1,LC DO 21 L1=1,LL SM=NLM(L1,M1,SIL,LC2)*S11(L1,M1,L,M)/NLM(L,M,SIL,LC2) C TRANSITION 1 <- 1 in SIGMA --------------------------------------------- S11P(L1,L,M)=SM C TRANSITION 0 <- 0 in SIGMA, PART 0 -------------------------------------- So00P(L1,L,M)=So00(L1,M1,L,M) C TRANSITION 0 <- 0 in SIGMA, PART 2 -------------------------------------- S00P(L1,L,M)=S00(L1,M1,L,M) C TRANSITION 0 <- 1 and 1 <- 0 in SIGMA ----------------------------------- S01P(L1,L,M)=S01(L1,M1,L,M) S10P(L1,L,M)=S10(L1,M1,L,M) C TRANSITION 0 <- 2 and 2 <- 0 in SIGMA ----------------------------------- S02P(L1,L,M)=S02(L1,M1,L,M) S20P(L1,L,M)=S20(L1,M1,L,M) 21 CONTINUE 20 CONTINUE ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** FUNCTION NLM(L,M,SIL,LC2) IMPLICIT NONE REAL*8 NLM INTEGER L,M,LC2 REAL*8 SIL(0:2*LC2) NLM=SQRT( SIL(L+M)/( (2.D0*L+1.D0)*SIL(L-M) ) ) RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE CINIT_polym USE MULT ! LL USE MULT_WALL_polym ! LC,LC2,LCG1 USE FACT USE C_COEF_polym IMPLICIT NONE REAL*8 A,NLM INTEGER L,M INTEGER :: LL_OLD=-1 SAVE LL_OLD IF( LL_OLD.EQ.LL) RETURN C WRITE(*,*) 'NEW CINIT' LL_OLD=LL IF(ALLOCATED(C01)) THEN DEALLOCATE( C01,C02,C03,C04, * C11,C12,C13,C14, * C21,C22,C23,C24, * C25,C26,C27 ) ENDIF ALLOCATE( C01(LC,0:1),C02(LC,0:1),C03(LC,0:1),C04(LC,0:1), * C11(LC,0:1),C12(LC,0:1),C13(LC,0:1),C14(LC,0:1), * C21(LC,0:1),C22(LC,0:1),C23(LC,0:1),C24(LC,0:1), * C25(LC,0:1),C26(LC,0:1),C27(LC,0:1) ) DO M=0,1 DO 1 L=1,LC C----------------------------------------------------------------- IF(L-2.GE.1) THEN A=-(-1)**(L+M)*4.D0*(L-2)*(L+M)*(L+M-1)/ * ( 1.D0*(L-1)*(2*L-1)*(2*L-3) ) A=A*NLM(L-2,M,SIL,LC)/NLM(L,M,SIL,LC) C01(L,M)=A ENDIF IF(L-1.GE.1.AND.L-1.LE.LL) THEN A=-(-1)**(L+M)*2.D0*(L+M) A=A*NLM(L-1,M,SIL,LC)/NLM(L,M,SIL,LC) C02(L,M)=A ENDIF IF(L-1.GE.1.AND.L-1.LE.LL) THEN A=(-1)**(L+M)*4.D0*M*(L+M)/( 1.D0*L*(L-1) ) A=A*NLM(L-1,M,SIL,LC)/NLM(L,M,SIL,LC) C03(L,M)=A ENDIF IF(L.LE.LL) THEN A=-(-1)**(L+M)*( 1.D0 + 2.D0*(L-2)*(L+M)*(L-M)/( 1.D0*L*(2*L-1)) ) C04(L,M)=A ENDIF C----------------------------------------------------------------- IF(L-1.GE.1.AND.L-1.LE.LL) THEN A=-(-1)**(L+M+1)*4.D0*(L-1)*M*(L+M)/( 1.D0*L*(2*L-1)*(2*L+1) ) A=A*NLM(L-1,M,SIL,LC)/NLM(L,M,SIL,LC) C11(L,M)=A ENDIF IF(L.LE.LL) THEN A=-(-1)**(L+M+1)*2.D0*M C12(L,M)=A ENDIF IF(L.LE.LL) THEN A=-(-1)**(L+M+1)*( 1.D0 - 4.D0*M*M/(1.D0*L*(L+1)) ) C13(L,M)=A ENDIF IF(L+1.LE.LL) THEN A=-(-1)**(L+M+1)*2.D0*(L-1)*M*(L-M+1)/( 1.D0*(L+1)*(2*L+1) ) A=A*NLM(L+1,M,SIL,LC)/NLM(L,M,SIL,LC) C14(L,M)=A ENDIF C----------------------------------------------------------------- IF(L-1.GE.1.AND.L-1.LE.LL) THEN A= * (-1)**(L+M+2)*2.D0*(L-1)*(L+1)*(2*L+3)*(L+M)/( 1.D0*L*L*(2*L-1) ) A=A*NLM(L-1,M,SIL,LC)/NLM(L,M,SIL,LC) C21(L,M)=A ENDIF IF(L.LE.LL) THEN A=(-1)**(L+M+2)*1.D0*(L+1)*(2*L+1)*(2*L+3)/L C22(L,M)=A ENDIF IF(L.LE.LL) THEN A=-(-1)**(L+M+2)*2.D0*M*(2*L+1)*(2*L+3)/( 1.D0*L*L ) C23(L,M)=A ENDIF IF(L.LE.LL) THEN A=-(-1)**(L+M+2)*(1.D0 - 2.D0*(L+3)*(L-M+1)*(L+M+1)/ * (1.D0*(L+1)*(2*L+3)) ) C24(L,M)=A ENDIF IF(L+1.LE.LL) THEN A=(-1)**(L+M+2)*4.D0*(L+2)*(L-M+1) A=A*NLM(L+1,M,SIL,LC)/NLM(L,M,SIL,LC) C25(L,M)=A ENDIF IF(L+1.LE.LL) THEN A=-(-1)**(L+M+2)*2.D0*M*(L+3)*(2*L+1)*(L-M+1)/ * ( 1.D0*L*(L+1)*(L+2) ) A=A*NLM(L+1,M,SIL,LC)/NLM(L,M,SIL,LC) C26(L,M)=A ENDIF IF(L+2.LE.LL) THEN A=(-1)**(L+M+2)*1.D0*(L+3)*(2*L+1)*(L-M+1)*(L-M+2)/ * (1.D0*(L+2)*(2*L+3)) A=A*NLM(L+2,M,SIL,LC)/NLM(L,M,SIL,LC) C27(L,M)=A ENDIF C----------------------------------------------------------------- 1 CONTINUE ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** C 'l' lubrication SUBROUTINE LUB2WALL_polym(A,CONF,HH,LUB_WALL) USE SIZE ! NN,LLG IMPLICIT NONE REAL*8 A(2*NN,2*NN,0:1),CONF(NN),HH CHARACTER*1 LUB_WALL REAL*8 Z REAL*8 A1(2,2,0:1) INTEGER I A=0.D0 IF(LUB_WALL.EQ.'Y'.OR.LUB_WALL.EQ.'y') THEN ! WALL LUBRICATION DO I=1,NN Z=CONF(I) CALL LUB_L_II(A1,Z) ! LOWER WALL LUBRICATION A( (I-1)*2+1:I*2 , (I-1)*2+1:I*2 , 0:1 )= * A( (I-1)*2+1:I*2 , (I-1)*2+1:I*2 , 0:1 )+A1 CALL LUB_U_II(A1,Z,HH) ! UPPER WALL LUBRICATION A( (I-1)*2+1:I*2 , (I-1)*2+1:I*2 , 0:1 )= * A( (I-1)*2+1:I*2 , (I-1)*2+1:I*2 , 0:1 )+A1 ENDDO ENDIF ! END OF WALL LUBRICATION RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE G2WALL_DEF_polym(G,CONF,HH,KUU,NG,SYMM,M) USE SIZE ! NN,LLG USE MULT_FULL_polym ! LL2,LLG1 IMPLICIT NONE REAL*8 G(LLG,LLG),CONF(NN),HH REAL*8 KUU INTEGER NG LOGICAL SYMM INTEGER M REAL*8 R,Z1,Z2 REAL*8 G1(LLG1,LLG1),GL1(LLG1,LLG1) REAL*8 GU1(LLG1,LLG1),GF1(LLG1,LLG1) INTEGER I,J CALL Z_II_polym(G1) DO I=1,NN Z1=CONF(I) Z2=CONF(I) CALL G_LOWER_WALL_IJ_polym(GL1,Z1,Z2,M) CALL G_UPPER_WALL_IJ_polym(GU1,Z1,Z2,HH,M) CALL G2WALL_FOURIER_IJ(GF1,Z1,Z2,HH,KUU,NG,M) G( (I-1)*LLG1+1:I*LLG1 , (I-1)*LLG1+1:I*LLG1 )=G1+GL1+GU1+GF1 ENDDO DO I=1,NN-1 DO J=I+1,NN Z1=CONF(I) Z2=CONF(J) R=CONF(I)-CONF(J) CALL G_LOWER_WALL_IJ_polym(GL1,Z1,Z2,M) CALL G_UPPER_WALL_IJ_polym(GU1,Z1,Z2,HH,M) CALL G2WALL_FOURIER_IJ(GF1,Z1,Z2,HH,KUU,NG,M) CALL GDIR_IJ_polym(G1,R,M) G( (I-1)*LLG1+1:I*LLG1 , (J-1)*LLG1+1:J*LLG1 )=GL1+GU1+GF1+G1 ENDDO ENDDO ! END OF GRAND MOBILITY IF(.NOT.SYMM) RETURN DO I=1,LLG-1 ! SYMMETRIZATION DO J=I+1,LLG G(J,I)=G(I,J) ENDDO ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE Z_II_polym(G1) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 USE Z_COEF IMPLICIT NONE REAL*8 G1(LLG1,LLG1) INTEGER B,L,H H(B,L)=B*LL + L G1=0.D0 DO L=1,LL G1(H(0,L),H(0,L))=Z0(L) G1(H(0,L),H(2,L))=Z2(L) G1(H(2,L),H(0,L))=Z2(L) G1(H(2,L),H(2,L))=ZB(L) G1(H(1,L),H(1,L))=Z1(L) ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE GDIR_IJ_polym(G1,R,M) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 USE S_COEF_polym IMPLICIT NONE REAL*8 G1(LLG1,LLG1),R INTEGER M REAL*8 YC1(0:LL2),YC2(0:LL2) INTEGER L1,L,B,H INTEGER HLC,HRC H(B,L)=B*LL+L C H(I,B,L)=NN*LL*B + (I-1)*LL + L C CLEARING G ************************************************ G1=0.D0 CALL YLMR_polym(YC1,YC2,R,LL2) DO 20 L =1,LL DO 20 L1=1,LL C TRANSITION 1 <- 1 in SIGMA --------------------------------------------- HLC=H(1,L1) HRC=H(1,L ) G1(HLC,HRC)=S11P(L1,L,M)*YC1(L+L1) C TRANSITION 0 <- 0 in SIGMA, PART o ------------------------------------ HLC=H(0,L1) HRC=H(0,L ) G1(HLC,HRC)=So00P(L1,L,M)*YC2(L+L1) C TRANSITION 0 <- 0 in SIGMA, part 2 ************************************ G1(HLC,HRC)=G1(HLC,HRC)+S00P(L1,L,M)*YC1(L+L1-2) C TRANSITION 0 <- 1 in SIGMA ******************************************** HLC=H(0,L1) HRC=H(1,L ) G1(HLC,HRC)= S01P(L1,L,M)*YC1(L+L1-1) C TRANSITION 1 <- 0 in SIGMA ******************************************** HLC=H(1,L1) HRC=H(0,L ) G1(HLC,HRC)=S10P(L1,L,M)*YC1(L+L1-1) C TRANSITION 0 <- 2 in SIGMA ******************************************** HLC=H(0,L1) HRC=H(2,L ) G1(HLC,HRC)=S02P(L1,L,M)*YC1(L+L1) C TRANSITION 2 <- 0 in SIGMA ******************************************** HLC=H(2,L1) HRC=H(0,L ) G1(HLC,HRC)=S20P(L1,L,M)*YC1(L+L1) C END OF THE L,M,L1,M1 LOOPS 20 20 CONTINUE RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE YLMR_polym(YC1,YC2,R,LL2) IMPLICIT NONE INTEGER LL2 REAL*8 R,YC1(0:LL2),YC2(0:LL2) INTEGER L,RU IF(R.GT.0) THEN RU=1 ELSE RU=-1 ENDIF R=ABS(R) DO L=0,LL2 YC1(L) = RU**L/( R**(L+1) ) YC2(L) = RU**L/( R**(L-1) ) ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE G_LOWER_WALL_IJ_polym(G1,Z1,Z2,M) USE MULT_FULL_polym ! LL2,LLG1 USE MULT_WALL_polym ! LC,LC2,LCG1 IMPLICIT NONE REAL*8 G1(LLG1,LLG1),Z1,Z2 INTEGER M REAL*8 R REAL*8 GP1(LLG1,LCG1) R=Z1+Z2 CALL GP_IJ_polym(GP1,R,M) CALL GPXC_IJ_polym(G1,GP1,Z2,M) RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE G_UPPER_WALL_IJ_polym(G1,Z1,Z2,HH,M) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 USE MULT_WALL_polym ! LC,LC2,LCG1 IMPLICIT NONE REAL*8 G1(LLG1,LLG1),Z1,Z2,HH INTEGER M REAL*8 GP1(LLG1,LCG1) REAL*8 R INTEGER B,L,H,B1,L1,B2,L2 H(B,L) =B*LL + L R=2.D0*HH-Z1-Z2 CALL GP_IJ_polym(GP1,R,M) CALL GPXC_IJ_polym(G1,GP1,HH-Z2,M) DO B1=0,2 DO L1=1,LL DO B2=0,2 DO L2=1,LL G1(H(B1,L1),H(B2,L2))= * (-1)**(B1+L1+B2+L2)*G1(H(B1,L1),H(B2,L2)) ENDDO ENDDO ENDDO ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE GP_IJ_polym(GP1,R,M) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 USE MULT_WALL_polym ! LC,LC2,LCG1 USE S_COEF_polym IMPLICIT NONE REAL*8 GP1(LLG1,LCG1),R INTEGER M REAL*8 YC1(0:LC2),YC2(0:LC2) INTEGER L1,L,B,H,HC INTEGER HLC,HRC H(B,L) =B*LL + L HC(B,L)=B*LC + L C H(I,B,L)=NN*LL*B + (I-1)*LL + L C CLEARING GP ********************************************** GP1=0.D0 CALL YLMR_polym(YC1,YC2,R,LC2) DO 20 L =1,LC DO 20 L1=1,LL C TRANSITION 1 <- 1 in SIGMA --------------------------------------------- HLC= H(1,L1) HRC=HC(1,L ) GP1(HLC,HRC)=S11P(L1,L,M)*YC1(L+L1) C TRANSITION 0 <- 0 in SIGMA, PART o ------------------------------------ HLC= H(0,L1) HRC=HC(0,L ) GP1(HLC,HRC)=So00P(L1,L,M)*YC2(L+L1) C TRANSITION 0 <- 0 in SIGMA, part 2 ************************************ GP1(HLC,HRC)=GP1(HLC,HRC)+S00P(L1,L,M)*YC1(L+L1-2) C TRANSITION 0 <- 1 in SIGMA ******************************************** HLC= H(0,L1) HRC=HC(1,L ) GP1(HLC,HRC)= S01P(L1,L,M)*YC1(L+L1-1) C TRANSITION 1 <- 0 in SIGMA ******************************************** HLC= H(1,L1) HRC=HC(0,L ) GP1(HLC,HRC)=S10P(L1,L,M)*YC1(L+L1-1) C TRANSITION 0 <- 2 in SIGMA ******************************************** HLC= H(0,L1) HRC=HC(2,L ) GP1(HLC,HRC)=S02P(L1,L,M)*YC1(L+L1) C TRANSITION 2 <- 0 in SIGMA ******************************************** HLC= H(2,L1) HRC=HC(0,L ) GP1(HLC,HRC)=S20P(L1,L,M)*YC1(L+L1) C END OF THE L,M,L1,M1 LOOPS 20 20 CONTINUE RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE GPXC_IJ_polym(G1,GP1,Z,M) ! Z=Z2 USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 USE MULT_WALL_polym ! LC,LC2,LCG1 USE C_COEF_polym IMPLICIT NONE REAL*8 G1(LLG1,LLG1),GP1(LLG1,LCG1) INTEGER M INTEGER B,L,H,HC,HGL REAL*8 Z H(B,L) =B*LL + L HC(B,L)=B*LC + L C CLEARING G1 ********************************************** G1=0.D0 DO 1 HGL=1,LLG1 DO 1 L=1,LC C----------------------------------------------------------------- IF(L-2.GE.1) THEN G1(HGL,H(2,L-2))= * G1(HGL,H(2,L-2))+GP1(HGL,HC(0,L))*C01(L,M) ENDIF IF(L-1.GE.1.AND.L-1.LE.LL) THEN G1(HGL,H(0,L-1))= * G1(HGL,H(0,L-1))+GP1(HGL,HC(0,L))*C02(L,M)*Z ENDIF IF(L-1.GE.1.AND.L-1.LE.LL) THEN G1(HGL,H(1,L-1))= * G1(HGL,H(1,L-1))+GP1(HGL,HC(0,L))*C03(L,M) ENDIF IF(L.LE.LL) THEN G1(HGL,H(0,L ))= * G1(HGL,H(0,L ))+GP1(HGL,HC(0,L))*C04(L,M) ENDIF C----------------------------------------------------------------- IF(L-1.GE.1.AND.L-1.LE.LL) THEN G1(HGL,H(2,L-1))= * G1(HGL,H(2,L-1))+GP1(HGL,HC(1,L))*C11(L,M) ENDIF IF(L.LE.LL) THEN G1(HGL,H(0,L ))= * G1(HGL,H(0,L ))+GP1(HGL,HC(1,L))*C12(L,M)*Z ENDIF IF(L.LE.LL) THEN G1(HGL,H(1,L ))= * G1(HGL,H(1,L ))+GP1(HGL,HC(1,L))*C13(L,M) ENDIF IF(L+1.LE.LL) THEN G1(HGL,H(0,L+1))= * G1(HGL,H(0,L+1))+GP1(HGL,HC(1,L))*C14(L,M) ENDIF C----------------------------------------------------------------- IF(L-1.GE.1.AND.L-1.LE.LL) THEN G1(HGL,H(2,L-1))= * G1(HGL,H(2,L-1))+GP1(HGL,HC(2,L))*C21(L,M)*Z ENDIF IF(L.LE.LL) THEN G1(HGL,H(0,L ))= * G1(HGL,H(0,L ))+GP1(HGL,HC(2,L))*C22(L,M)*Z*Z ENDIF IF(L.LE.LL) THEN G1(HGL,H(1,L ))= * G1(HGL,H(1,L ))+GP1(HGL,HC(2,L))*C23(L,M)*Z ENDIF IF(L.LE.LL) THEN G1(HGL,H(2,L ))= * G1(HGL,H(2,L ))+GP1(HGL,HC(2,L))*C24(L,M) ENDIF IF(L+1.LE.LL) THEN G1(HGL,H(0,L+1))= * G1(HGL,H(0,L+1))+GP1(HGL,HC(2,L))*C25(L,M)*Z ENDIF IF(L+1.LE.LL) THEN G1(HGL,H(1,L+1))= * G1(HGL,H(1,L+1))+GP1(HGL,HC(2,L))*C26(L,M) ENDIF IF(L+2.LE.LL) THEN G1(HGL,H(0,L+2))= * G1(HGL,H(0,L+2))+GP1(HGL,HC(2,L))*C27(L,M) ENDIF C----------------------------------------------------------------- 1 CONTINUE RETURN END C*********************************************************** C*********************************************************** C*********************************************************** C NEW, February 02, 2007 SUBROUTINE G2WALL_FOURIER_IJ(G1,Z1,Z2,HH,KUU,NG,M) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLG1 IMPLICIT NONE REAL*8 G1(LLG1,LLG1),Z1,Z2,HH,KUU INTEGER NG,M REAL*8 UNIT,KU,BU,AL REAL*8 SZSKJ(6,6,0:2*LL+2) REAL*8 SZSI(6,6,0:2*LL+2) REAL*8 GR(6,0:2*LL+2,LLG1) REAL*8 T(6,3,LL) REAL*8 H INTEGER I,B1,B2,C1,C2,L1,L2 INTEGER HG,L,B,j REAL*8 dx,xm,xr,w(5),x(5) SAVE w,x DATA w/ 0.29552422471475287017389299465133830D0, * 0.26926671930999635509122692156946934D0, * 0.21908636251598204399553493422816310D0, * 0.14945134915058059314577633965769722D0, * 6.6671344308688137593568809893331896D-0002 / DATA x/ 0.14887433898163121088482600112971997D0, * 0.43339539412924719079926594316578413D0, * 0.67940956829902440623432736511487354D0, * 0.86506336668898451073209668842349306D0, * 0.97390652851717172007796401208445200D0 / HG(B,L)=B*LL + L UNIT=(2*LL+2)/MAX( Z1+HH+(HH-Z2) , (HH-Z1)+HH+Z2 ) KU=KUU*UNIT if(4.D0*KU*HH.GT.709.D0) KU=709.D0/(4.D0*HH) ! PATCH SZSI=0.D0 H=KU/NG DO I=0,NG-1 ! GAUSS MAIN LOOP START AL=I*H BU=AL+H xm=0.5*(BU+AL) xr=0.5*(BU-AL) do j=1,5 dx=xr*x(j) CALL SZSKJ_DEF(SZSKJ,Z1,Z2,HH,xm+dx) SZSI=SZSI+w(J)*SZSKJ enddo do j=1,5 dx=xr*x(j) CALL SZSKJ_DEF(SZSKJ,Z1,Z2,HH,xm-dx) SZSI=SZSI+w(J)*SZSKJ enddo ENDDO ! GAUSS MAIN LOOP END SZSI=xr*SZSI ! GAUSS SCAILING TO THE RANGE OF INTEGRATION *------------------ DO L=1,LL CALL T_DEF(T(1,1,L),L,M) ENDDO G1=0.D0 GR=0.D0 DO B2=0,2 DO L2=1,LL C DO M2=-L2,L2 DO C1=1,6 DO L=1+L2+B2-2,LL+L2+B2 ! INSTAED OF C DO L=0,2*LL+2 C DO M=0,LL+ABS(M2) ! INSTAED OF C DO M=0,2*LL DO C2=1,6 GR(C1,L,HG(B2,L2))=GR(C1,L,HG(B2,L2))+ * SZSI(C1,C2,L)*T(C2,B2+1,L2) ENDDO ENDDO ENDDO ENDDO ENDDO DO L1=1,LL DO B1=0,2 DO L2=1,LL DO B2=0,2 DO C1=1,6 G1(HG(B1,L1),HG(B2,L2))=G1(HG(B1,L1),HG(B2,L2))+ * T(C1,B1+1,L1)*GR(C1,L1+B1+L2+B2-2,HG(B2,L2)) ENDDO ENDDO ENDDO ENDDO ENDDO RETURN END ********************************************************************** C NEW, November 28, 2005 SUBROUTINE SZSKJ_DEF(SZSKJ,Z1,Z2,HH,K) USE MULT ! LL IMPLICIT NONE REAL*8 SZSKJ(6,6,0:2*LL+2),Z1,Z2,HH,K REAL*8 SZS(6,6),KP(0:2*LL+2) INTEGER B1,B2,L CALL SZS_DEF(SZS,Z1,Z2,HH,K) DO L=0,2*LL+2 KP(L)=K**L ENDDO DO B1=1,6 DO B2=1,6 DO L=0,2*LL+2 SZSKJ(B1,B2,L)=-SZS(B1,B2)*KP(L) ! MINUS INCLUDED !!! ENDDO ENDDO ENDDO RETURN END ********************************************************************** SUBROUTINE SZS_DEF(SZS,Z1,Z2,H,K) IMPLICIT NONE REAL*8 SZS(6,6),Z1,Z2,H,K REAL*8 EZ1,EZ2,E2H,D1,D2,D3 EZ1=EXP(k*z1) EZ2=EXP(k*z2) E2H=EXP(2.D0*k*H) CALL D1_DEF(D1,k*H) D2=-D1 D3=-1+E2H SZS(1,1)=(-1.D0+E2H*(1.D0+4.D0*H**2*k**2))/(D1*Ez1*Ez2) SZS(1,2)=0.D0 SZS(1,3)=(-2.D0*k*(-z2+E2H*(H+z2+4*H**2*k**2*z2)))/(D1*Ez1*Ez2) SZS(1,4)=-((Ez2*(-1+E2H*(1+4*H**2*k**2-4*H*k**2*z2)))/(D1*Ez1)) SZS(1,5)=0.D0 SZS(1,6)=(2.D0*E2H*Ez2*H*k)/(D1*Ez1) SZS(2,1)=0.D0 SZS(2,2)=1.D0/(D3*Ez1*Ez2) SZS(2,3)=0.D0 SZS(2,4)=0.D0 SZS(2,5)=-(Ez2/(D3*Ez1)) SZS(2,6)=0.D0 SZS(3,1)=(-2.D0*k*(-z1+E2H*(H+z1+4*H**2*k**2*z1)))/(D1*Ez1*Ez2) SZS(3,2)=0.D0 SZS(3,3)=(-1.D0-4.D0*k**2*z1*z2+E2H*(1.D0+4*k**2*z1*z2+16.D0* * H**2*k**4*z1*z2+4.D0*H*k**2*(z1+z2)))/(D1*Ez1*Ez2) SZS(3,4)=(2.D0*Ez2*k*(-z1+z2+E2H*(H+z1+4*H**2*k**2*z1-z2-4*H* * k**2*z1*z2)))/(D1*Ez1) SZS(3,5)=0.D0 SZS(3,6)=-((Ez2*(-1.D0+E2H*(1.D0+4*H*k**2*z1)))/(D1*Ez1)) SZS(4,1)=-((Ez1*(-1.D0+E2H*(1.D0+4*H**2*k**2-4.D0*H*k**2*z1)))/ * (D1*Ez2)) SZS(4,2)=0.D0 SZS(4,3)=(2.D0*Ez1*k*(z1-z2+E2H*(H-z1+z2+4*H**2*k**2*z2- * 4.D0*H*k**2*z1*z2)))/(D1*Ez2) SZS(4,4)=(Ez1*Ez2*(-1.D0-4.D0*H**2*k**2-4.D0*k**2*z1*z2+4.D0*H* * k**2*(z1+z2)+E2H*(1.D0+16.D0*H**4*k**4+4.D0*k**2*z1*z2-8.D0*H* * k**2*(z1+z2)-16.D0*H**3*k**4*(z1+z2)+4.D0*H**2*k**2*(3.D0+4.D0* * k**2*z1*z2))))/(D1*E2H) SZS(4,5)=0.D0 SZS(4,6)=(2.D0*Ez1*Ez2*k*((-1.D0+2.D0*E2H)*H+4.D0*E2H*H**3* * k**2+z1-E2H*z1-4.D0*E2H*H**2*k**2*z1))/(D2*E2H) SZS(5,1)=0.D0 SZS(5,2)=-(Ez1/(D3*Ez2)) SZS(5,3)=0.D0 SZS(5,4)=0.D0 SZS(5,5)=(Ez1*Ez2)/(D3*E2H) SZS(5,6)=0.D0 SZS(6,1)=(2*E2H*Ez1*H*k)/(D1*Ez2) SZS(6,2)=0.D0 SZS(6,3)=-((Ez1*(-1.D0+E2H*(1.D0+4.D0*H*k**2*z2)))/(D1*Ez2)) SZS(6,4)=(2.D0*Ez1*Ez2*k*((-1.D0+2.D0*E2H)*H+4.D0*E2H*H**3* * k**2+z2-E2H*z2-4.D0*E2H*H**2*k**2*z2))/(D2*E2H) SZS(6,5)=0.D0 SZS(6,6)=(Ez1*Ez2*(1.D0-1.D0/E2H+4*H**2*k**2))/D1 RETURN END ********************************************************************** SUBROUTINE D1_DEF(D1,K) IMPLICIT NONE REAL*8 D1,K REAL*8 C(4:25),S INTEGER N INTEGER INIC DATA INIC/1/ SAVE INIC,C IF(INIC.EQ.1) THEN INIC=0 C( 4)=1.3333333333333333D0 C( 5)=2.6666666666666667D0 C( 6)=2.8444444444444444D0 C( 7)=2.1333333333333333D0 C( 8)=1.2571428571428571D0 C( 9)=0.61798941798941799D0 C(10)=0.26299823633156966D0 C(11)=0.099329805996472663D0 C(12)=0.033879536101758324D0 C(13)=0.010569183902517236D0 C(14)=0.0030445744731459017D0 C(15)=0.0008157956094464031D0 C(16)=0.00020451897700574949D0 C(17)=0.00004819949793494767D0 C(18)=0.000010720826296439258D0 C(19)=2.2581883770803783D-6 C(20)=4.5176956081668844D-7 C(21)=8.606538313937843D-8 C(22)=1.5649669482100519D-8 C(23)=2.7218179330265624D-9 C(24)=4.5364881486794064D-10 C(25)=7.2584904977070422D-11 ENDIF IF(K.GT.0.5D0) THEN D1=1.D0+EXP(4.D0*K)-2.D0*EXP(2.D0*K)*(1.D0+2.D0*K**2) ELSE S=0.D0 DO N=25,4,-1 S=C(N)+K*S ENDDO S=S*K**4 D1=S ENDIF RETURN END ********************************************************************** C NEW, February 14, 2006 SUBROUTINE T_DEF(TCS,L,M) USE FACT ! SIL IMPLICIT NONE INTEGER L,M REAL*8 TCS(6,3) REAL*8 A,B,C,N2,N3 C N1= 1.D0 N2=-1.D0/(L+1.D0) N3= 1.D0*L/( (L+1.D0)*(2*L+1.D0)*(2*L+3.D0) ) A=1.D0/SQRT(4.D0*SIL(L-M)*SIL(L+M)*(2*L+1)) B=2.D0*A*M/L C=A*(L*(2*L**2-2.D0*L-1.D0)-2.D0*M**2*(L-2.D0))/(L*(2.D0*L-1.D0)) TCS(1,1)= (-1)**(L+M)*C TCS(1,2)=-(-1)**(L+M)*2.D0*B*N2 TCS(1,3)= (-1)**(L+M)*4.D0*A*N3 TCS(2,1)= (-1)**(L+M)*B TCS(2,2)=-(-1)**(L+M)*2.D0*A*N2 TCS(2,3)= 0.D0 TCS(3,1)= (-1)**(L+M)*A TCS(3,2)= 0.D0 TCS(3,3)= 0.D0 TCS(4,1)= A TCS(4,2)= 0.D0 TCS(4,3)= 0.D0 TCS(5,1)= B TCS(5,2)= 2.D0*A*N2 TCS(5,3)= 0.D0 TCS(6,1)= C TCS(6,2)= 2.D0*B*N2 TCS(6,3)= 4.D0*A*N3 RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE SOLVE_polym(G,A,NN,LLG,LL) IMPLICIT NONE INTEGER NN,LLG,LL REAL*8 G(LLG,LLG),A(2*NN,2*NN) REAL*8 GB(LLG,2*NN) CHARACTER*1 UPLO PARAMETER (UPLO = 'U') INTEGER I,J,B,B1,B2,L,H,HT, INFO REAL*8 P PARAMETER(P=4.D0/3.D0) H(I,B,L)=(3*(I-1)+B)*LL + L HT(I,B)=2*(I-1)+B+1 C H(I,B,L)=NN*LLF*B + (I-1)*LLF + L GB=0.D0 DO I=1,NN DO B=0,1 GB(H(I,B,1),HT(I,B))=1.D0 ENDDO ENDDO CALL DPOTRF(UPLO,LLG,G,LLG,INFO) * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, the leading minor of order k is not * positive definite, and the factorization could not be * completed. IF(INFO.GT.0) THEN WRITE(*,*) 'k=',INFO WRITE(*,*) 'The leading minor of order k is not' WRITE(*,*) 'positive definite, and the factorization could not be' WRITE(*,*) 'completed.' ENDIF CALL DPOTRS(UPLO,LLG,2*NN,G,LLG,GB,LLG,INFO) DO I=1,NN DO B1=0,1 DO J=1,NN DO B2=0,1 A(HT(I,B1),HT(J,B2))=A(HT(I,B1),HT(J,B2))+ * P*GB(H(I,B1,1),HT(J,B2)) ENDDO ENDDO ENDDO ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE EXTR_G_polym(A,G) USE SIZE ! NN,LLG USE MULT ! LL IMPLICIT NONE REAL*8 A(2*NN,2*NN),G(LLG,LLG) INTEGER I,J,B,B1,B2,L,H,HT REAL*8 P PARAMETER(P=4.D0/3.D0) H(I,B,L)=(3*(I-1)+B)*LL + L HT(I,B)=2*(I-1)+B+1 DO I=1,NN DO J=1,NN DO B1=0,1 DO B2=0,1 A(HT(I,B1),HT(J,B2))=A(HT(I,B1),HT(J,B2))+ * P*G(H(I,B1,1),H(J,B2,1)) ENDDO ENDDO ENDDO ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE FORCES_POISEU_polym(FX,TY,G,CONF,HH) USE SIZE ! NN,LLG USE MULT ! LL IMPLICIT NONE REAL*8 FX(NN),TY(NN),G(LLG,LLG),CONF(NN),HH REAL*8 FX1(NN),TY1(NN),FX2(NN),TY2(NN) INTEGER I,J,B,L,H REAL*8 PTT,PTD,PT3,SQ5,SQ2,P115,SQ23,P2SQ15 SAVE PTT,PTD,PT3,SQ5,SQ2,P115,SQ23,P2SQ15 C PARAMETER(PTT = 4.D0/3.D0) C PARAMETER(PTD = 0.4216370213557839D0) C PARAMETER(PT3 = 0.1840174824912945D0) C PARAMETER(SQ5 = 0.7071067811865476D0) C PARAMETER(SQ2 = 1.414213562373095D0) C PARAMETER(P115 = 1.D0/15.D0) C PARAMETER(SQ23 = 0.4714045207910317D0) C PARAMETER(P2SQ15 = 0.5163977794943222D0) LOGICAL :: FIX=.TRUE. SAVE FIX H(I,B,L)=(3*(I-1)+B)*LL + L IF(FIX) THEN FIX=.FALSE. PTT=4.D0/3.D0 PTD=SQRT(8.D0/45.D0) PT3=SQRT(32.D0/945.D0) SQ5=SQRT(0.5D0) SQ2=SQRT(2.D0) P115=1.D0/15.D0 SQ23=SQRT(2.D0)/3.D0 P2SQ15=2.D0/SQRT(15.D0) ENDIF FX1=0.D0 TY1=0.D0 FX2=0.D0 TY2=0.D0 DO I=1,NN DO J=1,NN FX1(I)=FX1(I) * -PTT*G(H(I,0,1),H(J,0,1))*(-CONF(J)) * +PTT*G(H(I,0,1),H(J,1,1))*(-0.5D0) * -PTD*G(H(I,0,1),H(J,0,2))*(-SQ5) TY1(I)=TY1(I) * +PTT*G(H(I,1,1),H(J,0,1))*(-CONF(J)) * -PTT*G(H(I,1,1),H(J,1,1))*(-0.5D0) * +PTD*G(H(I,1,1),H(J,0,2))*(-SQ5) FX2(I)=FX2(I) * -PTT*G(H(I,0,1),H(J,0,1))*(-CONF(J)**2) * +PTT*G(H(I,0,1),H(J,1,1))*(-CONF(J)) * -PTD*G(H(I,0,1),H(J,0,2))*(-SQ2*CONF(J)) * -PTT*G(H(I,0,1),H(J,2,1))*(-P115) * +PTD*G(H(I,0,1),H(J,1,2))*(-SQ23) * -PT3*G(H(I,0,1),H(J,0,3))*(-P2SQ15) TY2(I)=TY2(I) * +PTT*G(H(I,1,1),H(J,0,1))*(-CONF(J)**2) * -PTT*G(H(I,1,1),H(J,1,1))*(-CONF(J)) * +PTD*G(H(I,1,1),H(J,0,2))*(-SQ2*CONF(J)) * +PTT*G(H(I,1,1),H(J,2,1))*(-P115) * -PTD*G(H(I,1,1),H(J,1,2))*(-SQ23) * +PT3*G(H(I,1,1),H(J,0,3))*(-P2SQ15) ENDDO ENDDO FX=FX1*(4.D0/HH)+FX2*(-4.D0/HH**2) TY=TY1*(4.D0/HH)+TY2*(-4.D0/HH**2) RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE VELOCITIES_polym(UX,OY,A,FX,TY) USE SIZE ! NN,LLG IMPLICIT NONE REAL*8 UX(NN),OY(NN),A(2*NN,2*NN,0:1),FX(NN),TY(NN) INTEGER I,J UX=0.D0 OY=0.D0 DO I=1,NN DO J=1,NN UX(I)=UX(I)+A((I-1)*2+1,(J-1)*2+1,1)*FX(J)- * A((I-1)*2+1,(J-1)*2+2,1)*TY(J) OY(I)=OY(I)-A((I-1)*2+2,(J-1)*2+1,1)*FX(J)+ * A((I-1)*2+2,(J-1)*2+2,1)*TY(J) ENDDO ENDDO RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE MATREV(A,NN,FLAG) IMPLICIT NONE INTEGER NN REAL*8 A(NN,NN) CHARACTER*1 FLAG INTEGER INFO,I,J CHARACTER*1 UPLO PARAMETER (UPLO = 'U') IF(FLAG.EQ.'F'.OR.FLAG.EQ.'f') RETURN CALL DPOTRF(UPLO,NN,A,NN,INFO) * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, the leading minor of order k is not * positive definite, and the factorization could not be * completed. IF(INFO.GT.0) THEN WRITE(*,*) 'k=',INFO WRITE(*,*) 'The leading minor of order k is not' WRITE(*,*) 'positive definite, and the factorization could not be' WRITE(*,*) 'completed.' ENDIF CALL DPOTRI(UPLO,NN,A,NN,INFO) DO 334 I=2,NN DO 334 J=1,I-1 334 A(I,J)=A(J,I) RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE NORMA_polym(A,NN,FLAG,NORM) IMPLICIT NONE INTEGER NN REAL*8 A(2*NN,2*NN) INTEGER HT,I,J,B CHARACTER*1 FLAG,NORM HT(I,B)=2*(I-1)+B+1 IF(NORM.EQ.'N'.OR.NORM.EQ.'n') RETURN IF(FLAG.EQ.'F'.OR.FLAG.EQ.'f') THEN DO 1 I=1,NN DO 1 J=1,NN A(HT(I,0),HT(J,0))=A(HT(I,0),HT(J,0))/3.D0 1 CONTINUE ELSE DO 2 I=1,NN DO 2 J=1,NN A(HT(I,0),HT(J,0))=A(HT(I,0),HT(J,0))*3.D0 2 CONTINUE ENDIF RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE LUB_L_II(A1,Z) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLF,LLM,LLG1 IMPLICIT NONE REAL*8 A1(2,2,0:1),Z REAL*8 G1(LLG1,LLG1),GL1(LLG1,LLG1),AL1(2,2,0:1) INTEGER M REAL*8 XA1F,YA1F,YB1F,XC1F,YC1F REAL*8 XA1,YA1,YB1,XC1,YC1, X C CLEARING MATRIX AE1 **************************************** A1=0.D0 X=1.D0/Z XA1=XA1F(X)*3.D0 YA1=YA1F(X)*3.D0 YB1=YB1F(X) XC1=XC1F(X) YC1=YC1F(X) A1(1,1,0)=+XA1 A1(2,2,0)=+XC1 A1(1,1,1)=+YA1 A1(1,2,1)=+YB1 A1(2,1,1)=+YB1 A1(2,2,1)=+YC1 AL1=0.D0 DO M=0,1 CALL Z_II_polym(G1) CALL G_LOWER_WALL_IJ_polym(GL1,Z,Z,M) G1=G1+GL1 CALL SOLVE_polym(G1,AL1(1,1,M),1,LLG1,LL) ENDDO A1=A1-AL1 RETURN END C*********************************************************** C*********************************************************** C*********************************************************** SUBROUTINE LUB_U_II(A1,Z,HH) USE MULT ! LL USE MULT_FULL_polym ! LL2,LLF,LLM,LLG1 IMPLICIT NONE REAL*8 A1(2,2,0:1),Z,HH REAL*8 G1(LLG1,LLG1),GU1(LLG1,LLG1),AU1(2,2,0:1) INTEGER M REAL*8 XA1F,YA1F,YB1F,XC1F,YC1F REAL*8 XA1,YA1,YB1,XC1,YC1, X C CLEARING MATRIX AE1 **************************************** A1=0.D0 X=1.D0/(HH-Z) XA1=XA1F(X)*3.D0 YA1=YA1F(X)*3.D0 YB1=YB1F(X) XC1=XC1F(X) YC1=YC1F(X) A1(1,1,0)=+XA1 A1(2,2,0)=+XC1 A1(1,1,1)=+YA1 A1(1,2,1)=-YB1 A1(2,1,1)=-YB1 A1(2,2,1)=+YC1 AU1=0.D0 DO M=0,1 CALL Z_II_polym(G1) CALL G_UPPER_WALL_IJ_polym(GU1,Z,Z,HH,M) G1=G1+GU1 CALL SOLVE_polym(G1,AU1(1,1,M),1,LLG1,LL) ENDDO A1=A1-AU1 RETURN END C*********************************************************** C*********************************************************** C*********************************************************** FUNCTION XA1F(X) USE RIG_WALL_COEF IMPLICIT NONE REAL*8 XA1F,X REAL*8 T,S,A INTEGER I T=0.5D0*X A=t/(1.0D0 -t) - 0.2D0*log(1.0D0 - t) * -(1.0D0/21.0D0)*((1.0D0 - t)/t)*log(1.0D0 - t) S=0.D0 DO 1 I=MM-3,0,-1 1 S=S*T+XA(I) S=S+T**(MM-2)*( 0.875D0*XA(MM-2)+0.5D0*XA(MM-1)*T+ * 0.125D0*XA(MM)*T**2 ) XA1F=A+S RETURN END C--------------------------------------- FUNCTION YA1F(X) USE RIG_WALL_COEF IMPLICIT NONE REAL*8 YA1F,X REAL*8 T,S,A INTEGER I T=0.5D0*X A=-(8.0D0/15.0D0)*log(1.0D0 - t) * -(64.0D0/375.0D0)*((1.0D0 - t)/t)*log(1.0D0 - t) S=0.D0 DO 1 I=MM-3,0,-1 1 S=S*T+YA(I) S=S+T**(MM-2)*( 0.875D0*YA(MM-2)+0.5D0*YA(MM-1)*T+ * 0.125D0*YA(MM)*T**2 ) YA1F=A+S RETURN END C--------------------------------------- FUNCTION YB1F(X) USE RIG_WALL_COEF IMPLICIT NONE REAL*8 YB1F,X REAL*8 T,S,A INTEGER I T=0.5D0*X A=0.1D0*log(1.0D0 - t) * +(43.0D0/250.0D0)*((1.0D0 - t)/t)*log(1.0D0 - t) A=2.D0*A S=0.D0 DO 1 I=MM-3,0,-1 1 S=S*T+YB(I) S=S+T**(MM-2)*( 0.875D0*YB(MM-2)+0.5D0*YB(MM-1)*T+ * 0.125D0*YB(MM)*T**2 ) YB1F=-(A+S) RETURN END C--------------------------------------- FUNCTION XC1F(X) USE RIG_WALL_COEF IMPLICIT NONE REAL*8 XC1F,X REAL*8 T,S,A INTEGER I T=0.5D0*X A=0.5D0*((1.0D0 - t)/t)*log(1.0D0 - t) S=0.D0 DO 1 I=MM-3,0,-1 1 S=S*T+XC(I) S=S+T**(MM-2)*( 0.875D0*XC(MM-2)+0.5D0*XC(MM-1)*T+ * 0.125D0*XC(MM)*T**2 ) XC1F=A+S RETURN END C--------------------------------------- FUNCTION YC1F(X) USE RIG_WALL_COEF IMPLICIT NONE REAL*8 YC1F,X REAL*8 T,S,A INTEGER I T=0.5D0*X A=-0.4D0*log(1.0D0 - t) * -(66.0D0/125.0D0)*((1.0D0 - t)/t)*log(1.0D0 - t) S=0.D0 DO 1 I=MM-3,0,-1 1 S=S*T+YC(I) S=S+T**(MM-2)*( 0.875D0*YC(MM-2)+0.5D0*YC(MM-1)*T+ * 0.125D0*YC(MM)*T**2 ) YC1F=A+S RETURN END C*********************************************************** C*********************************************************** C*********************************************************** C Lapack Subroutines * ************************************************************************ * * File of the REAL*8 Level-3 BLAS. * ========================================== * * SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, * $ BETA, C, LDC ) * * SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * See: * * Dongarra J. J., Du Croz J. J., Duff I. and Hammarling S. * A set of Level 3 Basic Linear Algebra Subprograms. Technical * Memorandum No.88 (Revision 1), Mathematics and Computer Science * Division, Argonne National Laboratory, 9700 South Cass Avenue, * Argonne, Illinois 60439. * * ************************************************************************ * SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC REAL*8 ALPHA, BETA * .. Array Arguments .. REAL*8 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - REAL*8. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL*8 array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - REAL*8 array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - REAL*8. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - REAL*8 array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB REAL*8 TEMP * .. Parameters .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of DGEMM . * END SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO ) * * -- LAPACK auxiliary routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DLAUUM computes the product U * U' or L' * L, where the triangular * factor U or L is stored in the upper or lower triangular part of * the array A. * * If UPLO = 'U' or 'u' then the upper triangle of the result is stored, * overwriting the factor U in A. * If UPLO = 'L' or 'l' then the lower triangle of the result is stored, * overwriting the factor L in A. * * This is the blocked form of the algorithm, calling Level 3 BLAS. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the triangular factor stored in the array A * is upper or lower triangular: * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The order of the triangular factor U or L. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * On entry, the triangular factor U or L. * On exit, if UPLO = 'U', the upper triangle of A is * overwritten with the upper triangle of the product U * U'; * if UPLO = 'L', the lower triangle of A is overwritten with * the lower triangle of the product L' * L. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. REAL*8 ONE PARAMETER ( ONE = 1.0Q+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER I, IB, NB * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. EXTERNAL DGEMM, DLAUU2, DSYRK, DTRMM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAUUM', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * NB = ILAENV( 1, 'DLAUUM', UPLO, N, -1, -1, -1 ) * IF( NB.LE.1 .OR. NB.GE.N ) THEN * * Use unblocked code * CALL DLAUU2( UPLO, N, A, LDA, INFO ) ELSE * * Use blocked code * IF( UPPER ) THEN * * Compute the product U * U'. * DO 10 I = 1, N, NB IB = MIN( NB, N-I+1 ) CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', $ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ), $ LDA ) CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) IF( I+IB.LE.N ) THEN CALL DGEMM( 'No transpose', 'Transpose', I-1, IB, $ N-I-IB+1, ONE, A( 1, I+IB ), LDA, $ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA ) CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1, $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), $ LDA ) END IF 10 CONTINUE ELSE * * Compute the product L' * L. * DO 20 I = 1, N, NB IB = MIN( NB, N-I+1 ) CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', IB, $ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA ) CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) IF( I+IB.LE.N ) THEN CALL DGEMM( 'Transpose', 'No transpose', IB, I-1, $ N-I-IB+1, ONE, A( I+IB, I ), LDA, $ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA ) CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, ONE, $ A( I+IB, I ), LDA, ONE, A( I, I ), LDA ) END IF 20 CONTINUE END IF END IF * RETURN * * End of DLAUUM * END SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DPOTF2 computes the Cholesky factorization of a real symmetric * positive definite matrix A. * * The factorization has the form * A = U' * U , if UPLO = 'U', or * A = L * L', if UPLO = 'L', * where U is an upper triangular matrix and L is lower triangular. * * This is the unblocked version of the algorithm, calling Level 2 BLAS. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the upper or lower triangular part of the * symmetric matrix A is stored. * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * On entry, the symmetric matrix A. If UPLO = 'U', the leading * n by n upper triangular part of A contains the upper * triangular part of the matrix A, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading n by n lower triangular part of A contains the lower * triangular part of the matrix A, and the strictly upper * triangular part of A is not referenced. * * On exit, if INFO = 0, the factor U or L from the Cholesky * factorization A = U'*U or A = L*L'. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, the leading minor of order k is not * positive definite, and the factorization could not be * completed. * * ===================================================================== * * .. Parameters .. REAL*8 ONE, ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER J REAL*8 AJJ * .. * .. External Functions .. LOGICAL LSAME REAL*8 DDOT EXTERNAL LSAME, DDOT * .. * .. External Subroutines .. EXTERNAL DGEMV, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPOTF2', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * IF( UPPER ) THEN * * Compute the Cholesky factorization A = U'*U. * DO 10 J = 1, N * * Compute U(J,J) and test for non-positive-definiteness. * AJJ = A( J, J ) - DDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 ) IF( AJJ.LE.ZERO ) THEN A( J, J ) = AJJ GO TO 30 END IF AJJ = SQRT( AJJ ) A( J, J ) = AJJ * * Compute elements J+1:N of row J. * IF( J.LT.N ) THEN CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), $ LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA ) CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) END IF 10 CONTINUE ELSE * * Compute the Cholesky factorization A = L*L'. * DO 20 J = 1, N * * Compute L(J,J) and test for non-positive-definiteness. * AJJ = A( J, J ) - DDOT( J-1, A( J, 1 ), LDA, A( J, 1 ), $ LDA ) IF( AJJ.LE.ZERO ) THEN A( J, J ) = AJJ GO TO 30 END IF AJJ = SQRT( AJJ ) A( J, J ) = AJJ * * Compute elements J+1:N of column J. * IF( J.LT.N ) THEN CALL DGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ), $ LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 ) CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) END IF 20 CONTINUE END IF GO TO 40 * 30 CONTINUE INFO = J * 40 CONTINUE RETURN * * End of DPOTF2 * END SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DPOTRF computes the Cholesky factorization of a real symmetric * positive definite matrix A. * * The factorization has the form * A = U' * U , if UPLO = 'U', or * A = L * L', if UPLO = 'L', * where U is an upper triangular matrix and L is lower triangular. * * This is the block version of the algorithm, calling Level 3 BLAS. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the upper or lower triangular part of the * symmetric matrix A is stored. * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * On entry, the symmetric matrix A. If UPLO = 'U', the leading * n by n upper triangular part of A contains the upper * triangular part of the matrix A, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading n by n lower triangular part of A contains the lower * triangular part of the matrix A, and the strictly upper * triangular part of A is not referenced. * * On exit, if INFO = 0, the factor U or L from the Cholesky * factorization A = U'*U or A = L*L'. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, the leading minor of order k is not * positive definite, and the factorization could not be * completed. * * ===================================================================== * * .. Parameters .. REAL*8 ONE PARAMETER ( ONE = 1.0Q+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER J, JB, NB * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. EXTERNAL DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPOTRF', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.N ) THEN * * Use unblocked code. * CALL DPOTF2( UPLO, N, A, LDA, INFO ) ELSE * * Use blocked code. * IF( UPPER ) THEN * * Compute the Cholesky factorization A = U'*U. * DO 10 J = 1, N, NB * * Update and factorize the current diagonal block and test * for non-positive-definiteness. * JB = MIN( NB, N-J+1 ) CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE, $ A( 1, J ), LDA, ONE, A( J, J ), LDA ) CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN * * Compute the current block row. * CALL DGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1, $ J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ), $ LDA, ONE, A( J, J+JB ), LDA ) CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', $ JB, N-J-JB+1, ONE, A( J, J ), LDA, $ A( J, J+JB ), LDA ) END IF 10 CONTINUE * ELSE * * Compute the Cholesky factorization A = L*L'. * DO 20 J = 1, N, NB * * Update and factorize the current diagonal block and test * for non-positive-definiteness. * JB = MIN( NB, N-J+1 ) CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE, $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN * * Compute the current block column. * CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, $ J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ), $ LDA, ONE, A( J+JB, J ), LDA ) CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit', $ N-J-JB+1, JB, ONE, A( J, J ), LDA, $ A( J+JB, J ), LDA ) END IF 20 CONTINUE END IF END IF GO TO 40 * 30 CONTINUE INFO = INFO + J - 1 * 40 CONTINUE RETURN * * End of DPOTRF * END SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DPOTRI computes the inverse of a real symmetric positive definite * matrix A using the Cholesky factorization A = U'*U or A = L*L' * computed by DPOTRF. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the factor stored in A is upper or lower * triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * On entry, the triangular factor U or L from the Cholesky * factorization A = U'*U or A = L*L', as computed by DPOTRF. * On exit, the upper or lower triangle of the (symmetric) * inverse of A, overwriting the input factor U or L. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, the (k,k) element of the factor U or L is * zero, and the inverse could not be computed. * * ===================================================================== * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLAUUM, DTRTRI, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPOTRI', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Invert the triangular Cholesky factor U or L. * CALL DTRTRI( UPLO, 'Non-unit', N, A, LDA, INFO ) IF( INFO.GT.0 ) $ RETURN * * Form inv(U)*inv(U)' or inv(L)'*inv(L). * CALL DLAUUM( UPLO, N, A, LDA, INFO ) * RETURN * * End of DPOTRI * END * ************************************************************************ * SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC REAL*8 ALPHA, BETA * .. Array Arguments .. REAL*8 A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * DSYRK performs one of the symmetric rank k operations * * C := alpha*A*A' + beta*C, * * or * * C := alpha*A'*A + beta*C, * * where alpha and beta are scalars, C is an n by n symmetric matrix * and A is an n by k matrix in the first case and a k by n matrix * in the second case. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array C is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of C * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of C * is to be referenced. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. * * TRANS = 'T' or 't' C := alpha*A'*A + beta*C. * * TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with TRANS = 'N' or 'n', K specifies the number * of columns of the matrix A, and on entry with * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number * of rows of the matrix A. K must be at least zero. * Unchanged on exit. * * ALPHA - REAL*8. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL*8 array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * BETA - REAL*8. * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - REAL*8 array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array C must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA REAL*8 TEMP * .. Parameters .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. $ ( .NOT.LSAME( TRANS, 'T' ) ).AND. $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN INFO = 2 ELSE IF( N .LT.0 )THEN INFO = 3 ELSE IF( K .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYRK ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) )THEN * * Form C := alpha*A*A' + beta*C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A'*A + beta*C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP = ZERO DO 220, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * End of DSYRK . * END * ************************************************************************ * SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB REAL*8 ALPHA * .. Array Arguments .. REAL*8 A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - REAL*8. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - REAL*8 array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - REAL*8 array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA REAL*8 TEMP * .. Parameters .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of DTRSM . * END SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DTRTRI computes the inverse of a real upper or lower triangular * matrix A. * * This is the Level 3 BLAS version of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the matrix A is upper or lower triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * DIAG (input) CHARACTER*1 * Specifies whether or not the matrix A is unit triangular. * = 'N': Non-unit triangular * = 'U': Unit triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * * On entry, the triangular matrix A. If UPLO = 'U', the * leading n by n upper triangular part of the array A contains * the upper triangular matrix, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading n by n lower triangular part of the array A contains * the lower triangular matrix, and the strictly upper * triangular part of A is not referenced. If DIAG = 'U', the * diagonal elements of A are also not referenced and are * assumed to be 1. * * On exit, the (triangular) inverse of the original matrix, in * the same storage format. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * > 0: if INFO = k, A(k,k) is exactly zero. The triangular * matrix is singular and its inverse can not be computed. * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. REAL*8 ONE, ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER J, JB, NB, NN * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) NOUNIT = LSAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DTRTRI', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Check for singularity if non-unit. * IF( NOUNIT ) THEN DO 10 INFO = 1, N IF( A( INFO, INFO ).EQ.ZERO ) $ RETURN 10 CONTINUE INFO = 0 END IF * * Determine the block size for this environment. * NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.N ) THEN * * Use unblocked code * CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) ELSE * * Use blocked code * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix * DO 20 J = 1, N, NB JB = MIN( NB, N-J+1 ) * * Compute rows 1:j-1 of current block column * CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1, $ JB, ONE, A, LDA, A( 1, J ), LDA ) CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1, $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA ) * * Compute inverse of current diagonal block * CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO ) 20 CONTINUE ELSE * * Compute inverse of lower triangular matrix * NN = ( ( N-1 ) / NB )*NB + 1 DO 30 J = NN, 1, -NB JB = MIN( NB, N-J+1 ) IF( J+JB.LE.N ) THEN * * Compute rows j+jb:n of current block column * CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG, $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA, $ A( J+JB, J ), LDA ) CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG, $ N-J-JB+1, JB, -ONE, A( J, J ), LDA, $ A( J+JB, J ), LDA ) END IF * * Compute inverse of current diagonal block * CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO ) 30 CONTINUE END IF END IF * RETURN * * End of DTRTRI * END INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, $ N4 ) * * -- LAPACK auxiliary routine (preliminary version) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 20, 1992 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 * .. * * Purpose * ======= * * ILAENV is called from the LAPACK routines to choose problem-dependent * parameters for the local environment. See ISPEC for a description of * the parameters. * * This version provides a set of parameters which should give good, * but not optimal, performance on many of the currently available * computers. Users are encouraged to modify this subroutine to set * the tuning parameters for their particular machine using the option * and problem size information in the arguments. * * This routine will not function correctly if it is converted to all * lower case. Converting it to all upper case is allowed. * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies the parameter to be returned as the value of * ILAENV. * = 1: the optimal blocksize; if this value is 1, an unblocked * algorithm will give the best performance. * = 2: the minimum block size for which the block routine * should be used; if the usable block size is less than * this value, an unblocked routine should be used. * = 3: the crossover point (in a block routine, for N less * than this value, an unblocked routine should be used) * = 4: the number of shifts, used in the nonsymmetric * eigenvalue routines * = 5: the minimum column dimension for blocking to be used; * rectangular blocks must have dimension at least k by m, * where k is given by ILAENV(2,...) and m by ILAENV(5,...) * = 6: the crossover point for the SVD (when reducing an m by n * matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds * this value, a QR factorization is used first to reduce * the matrix to a triangular form.) * = 7: the number of processors * = 8: the crossover point for the multishift QR and QZ methods * for nonsymmetric eigenvalue problems. * * NAME (input) CHARACTER*(*) * The name of the calling subroutine, in either upper case or * lower case. * * OPTS (input) CHARACTER*(*) * The character options to the subroutine NAME, concatenated * into a single character string. For example, UPLO = 'U', * TRANS = 'T', and DIAG = 'N' for a triangular routine would * be specified as OPTS = 'UTN'. * * N1 (input) INTEGER * N2 (input) INTEGER * N3 (input) INTEGER * N4 (input) INTEGER * Problem dimensions for the subroutine NAME; these may not all * be required. * * (ILAENV) (output) INTEGER * >= 0: the value of the parameter specified by ISPEC * < 0: if ILAENV = -k, the k-th argument had an illegal value. * * Further Details * =============== * * The following conventions have been used when calling ILAENV from the * LAPACK routines: * 1) OPTS is a concatenation of all of the character options to * subroutine NAME, in the same order that they appear in the * argument list for NAME, even if they are not used in determining * the value of the parameter specified by ISPEC. * 2) The problem dimensions N1, N2, N3, N4 are specified in the order * that they appear in the argument list for NAME. N1 is used * first, N2 second, and so on, and unused problem dimensions are * passed a value of -1. * 3) The parameter value returned by ILAENV is checked for validity in * the calling subroutine. For example, ILAENV is used to retrieve * the optimal blocksize for STRTRI as follows: * * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) * IF( NB.LE.1 ) NB = MAX( 1, N ) * * ===================================================================== * * .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER I, IC, IZ, NB, NBMIN, NX * .. * .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL * .. * .. Executable Statements .. * GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC * * Invalid value for ISPEC * ILAENV = -1 RETURN * 100 CONTINUE * * Convert NAME to upper case if the first character is lower case. * ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN * * ASCII character set * IF( IC.GE.97 .AND. IC.LE.122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF * ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN * * EBCDIC character set * IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) $ SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF * ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN * * Prime machines: ASCII+128 * IF( IC.GE.225 .AND. IC.LE.250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF * C1 = SUBNAM( 1:1 ) SNAME = C1.EQ.'S' .OR. C1.EQ.'D' CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' IF( .NOT.( CNAME .OR. SNAME ) ) $ RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) * GO TO ( 110, 200, 300 ) ISPEC * 110 CONTINUE * * ISPEC = 1: block size * * In these examples, separate code is provided for setting NB for * real and complex. We assume that NB will take the same value in * single or REAL*8. * NB = 1 * IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'PO' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NB = 1 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRF' ) THEN NB = 64 ELSE IF( C3.EQ.'TRD' ) THEN NB = 1 ELSE IF( C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( C2.EQ.'GB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'PB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'TR' ) THEN IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN IF( C3.EQ.'EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN * 200 CONTINUE * * ISPEC = 2: minimum block size * NBMIN = 2 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN * 300 CONTINUE * * ISPEC = 3: crossover point * NX = 0 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( SNAME .AND. C3.EQ.'TRD' ) THEN NX = 1 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NX = 1 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN * 400 CONTINUE * * ISPEC = 4: number of shifts (used by xHSEQR) * ILAENV = 6 RETURN * 500 CONTINUE * * ISPEC = 5: minimum column dimension (not used) * ILAENV = 2 RETURN * 600 CONTINUE * * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) * ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN * 700 CONTINUE * * ISPEC = 7: number of processors (not used) * ILAENV = 1 RETURN * 800 CONTINUE * * ISPEC = 8: crossover point for multishift (used by xHSEQR) * ILAENV = 50 RETURN * * End of ILAENV * END LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * .. Executable Statements .. * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END REAL*8 function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c REAL*8 dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0Q0 dtemp = 0.0Q0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end * ************************************************************************ * * File of the REAL*8 Level-2 BLAS. * =========================================== * * SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) * * SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * * SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * * SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * * SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * * SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * * SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * * SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * * SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) * * SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP ) * * SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * * SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) * * See: * * Dongarra J. J., Du Croz J. J., Hammarling S. and Hanson R. J.. * An extended set of Fortran Basic Linear Algebra Subprograms. * * Technical Memoranda Nos. 41 (revision 3) and 81, Mathematics * and Computer Science Division, Argonne National Laboratory, * 9700 South Cass Avenue, Argonne, Illinois 60439, US. * * Or * * NAG Technical Reports TR3/87 and TR4/87, Numerical Algorithms * Group Ltd., NAG Central Office, 256 Banbury Road, Oxford * OX2 7DE, UK, and Numerical Algorithms Group Inc., 1101 31st * Street, Suite 100, Downers Grove, Illinois 60515-1263, USA. * ************************************************************************ * SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. REAL*8 ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. REAL*8 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL*8. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL*8 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - REAL*8 array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - REAL*8. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - REAL*8 array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. Local Scalars .. REAL*8 TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * END SUBROUTINE DLAUU2( UPLO, N, A, LDA, INFO ) * * -- LAPACK auxiliary routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DLAUU2 computes the product U * U' or L' * L, where the triangular * factor U or L is stored in the upper or lower triangular part of * the array A. * * If UPLO = 'U' or 'u' then the upper triangle of the result is stored, * overwriting the factor U in A. * If UPLO = 'L' or 'l' then the lower triangle of the result is stored, * overwriting the factor L in A. * * This is the unblocked form of the algorithm, calling Level 2 BLAS. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the triangular factor stored in the array A * is upper or lower triangular: * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The order of the triangular factor U or L. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * On entry, the triangular factor U or L. * On exit, if UPLO = 'U', the upper triangle of A is * overwritten with the upper triangle of the product U * U'; * if UPLO = 'L', the lower triangle of A is overwritten with * the lower triangle of the product L' * L. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. REAL*8 ONE PARAMETER ( ONE = 1.0Q+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER I REAL*8 AII * .. * .. External Functions .. LOGICAL LSAME REAL*8 DDOT EXTERNAL LSAME, DDOT * .. * .. External Subroutines .. EXTERNAL DGEMV, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAUU2', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * IF( UPPER ) THEN * * Compute the product U * U'. * DO 10 I = 1, N AII = A( I, I ) IF( I.LT.N ) THEN A( I, I ) = DDOT( N-I+1, A( I, I ), LDA, A( I, I ), LDA ) CALL DGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ), $ LDA, A( I, I+1 ), LDA, AII, A( 1, I ), 1 ) ELSE CALL DSCAL( I, AII, A( 1, I ), 1 ) END IF 10 CONTINUE * ELSE * * Compute the product L' * L. * DO 20 I = 1, N AII = A( I, I ) IF( I.LT.N ) THEN A( I, I ) = DDOT( N-I+1, A( I, I ), 1, A( I, I ), 1 ) CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA, $ A( I+1, I ), 1, AII, A( I, 1 ), LDA ) ELSE CALL DSCAL( I, AII, A( I, 1 ), LDA ) END IF 20 CONTINUE END IF * RETURN * * End of DLAUU2 * END subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increment, 8/21/90. c REAL*8 da,dx(1) integer i,incx,ix,m,mp1,n c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n dx(ix) = da*dx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end * ************************************************************************ * SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB REAL*8 ALPHA * .. Array Arguments .. REAL*8 A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRMM performs one of the matrix-matrix operations * * B := alpha*op( A )*B, or B := alpha*B*op( A ), * * where alpha is a scalar, B is an m by n matrix, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) multiplies B from * the left or right as follows: * * SIDE = 'L' or 'l' B := alpha*op( A )*B. * * SIDE = 'R' or 'r' B := alpha*B*op( A ). * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - REAL*8. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - REAL*8 array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - REAL*8 array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the matrix B, and on exit is overwritten by the * transformed matrix. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA REAL*8 TEMP * .. Parameters .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0Q+0, ZERO = 0.0Q+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*A*B. * IF( UPPER )THEN DO 50, J = 1, N DO 40, K = 1, M IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) DO 30, I = 1, K - 1 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 30 CONTINUE IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) B( K, J ) = TEMP END IF 40 CONTINUE 50 CONTINUE ELSE DO 80, J = 1, N DO 70 K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) B( K, J ) = TEMP IF( NOUNIT ) $ B( K, J ) = B( K, J )*A( K, K ) DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * Form B := alpha*B*A'. * IF( UPPER )THEN DO 110, J = 1, N DO 100, I = M, 1, -1 TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE B( I, J ) = ALPHA*TEMP 100 CONTINUE 110 CONTINUE ELSE DO 140, J = 1, N DO 130, I = 1, M TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 120, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 120 CONTINUE B( I, J ) = ALPHA*TEMP 130 CONTINUE 140 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*A. * IF( UPPER )THEN DO 180, J = N, 1, -1 TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = 1, M B( I, J ) = TEMP*B( I, J ) 150 CONTINUE DO 170, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 160, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE ELSE DO 220, J = 1, N TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 190, I = 1, M B( I, J ) = TEMP*B( I, J ) 190 CONTINUE DO 210, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 200, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 200 CONTINUE END IF 210 CONTINUE 220 CONTINUE END IF ELSE * * Form B := alpha*B*A'. * IF( UPPER )THEN DO 260, K = 1, N DO 240, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 230, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 230 CONTINUE END IF 240 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 250, I = 1, M B( I, K ) = TEMP*B( I, K ) 250 CONTINUE END IF 260 CONTINUE ELSE DO 300, K = N, 1, -1 DO 280, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 270, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 270 CONTINUE END IF 280 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE END IF 300 CONTINUE END IF END IF END IF * RETURN * * End of DTRMM . * END SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. REAL*8 A( LDA, * ) * .. * * Purpose * ======= * * DTRTI2 computes the inverse of a real upper or lower triangular * matrix. * * This is the Level 2 BLAS version of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the matrix A is upper or lower triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * DIAG (input) CHARACTER*1 * Specifies whether or not the matrix A is unit triangular. * = 'N': Non-unit triangular * = 'U': Unit triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) REAL*8 array, dimension (LDA,N) * On entry, the triangular matrix A. If UPLO = 'U', the * leading n by n upper triangular part of the array A contains * the upper triangular matrix, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading n by n lower triangular part of the array A contains * the lower triangular matrix, and the strictly upper * triangular part of A is not referenced. If DIAG = 'U', the * diagonal elements of A are also not referenced and are * assumed to be 1. * * On exit, the (triangular) inverse of the original matrix, in * the same storage format. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. REAL*8 ONE PARAMETER ( ONE = 1.0Q+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER J REAL*8 AJJ * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DSCAL, DTRMV, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) NOUNIT = LSAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DTRTI2', -INFO ) RETURN END IF * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix. * DO 10 J = 1, N IF( NOUNIT ) THEN A( J, J ) = ONE / A( J, J ) AJJ = -A( J, J ) ELSE AJJ = -ONE END IF * * Compute elements 1:j-1 of j-th column. * CALL DTRMV( 'Upper', 'No transpose', DIAG, J-1, A, LDA, $ A( 1, J ), 1 ) CALL DSCAL( J-1, AJJ, A( 1, J ), 1 ) 10 CONTINUE ELSE * * Compute inverse of lower triangular matrix. * DO 20 J = N, 1, -1 IF( NOUNIT ) THEN A( J, J ) = ONE / A( J, J ) AJJ = -A( J, J ) ELSE AJJ = -ONE END IF IF( J.LT.N ) THEN * * Compute elements j+1:n of j-th column. * CALL DTRMV( 'Lower', 'No transpose', DIAG, N-J, $ A( J+1, J+1 ), LDA, A( J+1, J ), 1 ) CALL DSCAL( N-J, AJJ, A( J+1, J ), 1 ) END IF 20 CONTINUE END IF * RETURN * * End of DTRTI2 * END * ************************************************************************ * SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. REAL*8 A( LDA, * ), X( * ) * .. * * Purpose * ======= * * DTRMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - REAL*8 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - REAL*8 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL*8 ZERO PARAMETER ( ZERO = 0.0Q+0 ) * .. Local Scalars .. REAL*8 TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LSAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := A*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*A( I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*A( I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A'*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = N, 1, -1 TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 110, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 130, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of DTRMV . * END ************************************************************************ SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. REAL*8 A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DPOTRS solves a system of linear equations A*X = B with a symmetric * positive definite matrix A using the Cholesky factorization A = U'*U * or A = L*L' computed by DPOTRF. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the factor stored in A is upper or lower * triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input) REAL*8 array, dimension (LDA,N) * The triangular factor U or L from the Cholesky factorization * A = U'*U or A = L*L', as computed by DPOTRF. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * B (input/output) REAL*8 array, dimension (LDB,NRHS) * On entry, the right hand side vectors B for the system of * linear equations. * On exit, the solution vectors, X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. REAL*8 ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL UPPER * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( NRHS.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPOTRS', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * IF( UPPER ) THEN * * Solve A*X = B where A = U'*U. * * Solve U'*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, $ ONE, A, LDA, B, LDB ) * * Solve U*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, $ NRHS, ONE, A, LDA, B, LDB ) ELSE * * Solve A*X = B where A = L*L'. * * Solve L*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, $ NRHS, ONE, A, LDA, B, LDB ) * * Solve L'*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', N, NRHS, $ ONE, A, LDA, B, LDB ) END IF * RETURN * * End of DPOTRS * END