home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / exampleCode / MP / timer / regular / linpackd.m < prev    next >
Encoding:
Text File  |  1994-08-02  |  37.7 KB  |  1,486 lines
  1. # 1 "linpackd.f"
  2. # 0 "linpackd.f"
  3. C*$* Padding ( DD3, DD2, DD1 )
  4. # 0 "linpackd.f"
  5. C*$* Storage Order ( AA, A, B, DD1, X, DD2, TIME, DD3, IPVT )
  6. *
  7. *PLEASE NOTE THAT netlib HAS MOVED, THE NEW ADDRESS IS netlib@ornl.gov.
  8. *THE OLD ADDRESS, netlib@mcs.anl.gov, WILL BE TURNED OFF SOON.
  9. *
  10. *** from netlib, Fri Jul 27 14:07:10 EDT 1990 ***
  11. # 6 "linpackd.f"
  12.       program MAIN
  13. # 6 "linpackd.f"
  14.       double precision second
  15.       DOUBLEPRECISION AA(200,200), A(201,200), B(200), X(200)
  16. # 8 "linpackd.f"
  17.       DOUBLEPRECISION TIME(8,6), CRAY, OPS, TOTAL, NORMA, NORMX
  18. # 9 "linpackd.f"
  19.       double precision resid,residn,eps,epslon
  20.       INTEGER IPVT(200)
  21. # 11 "linpackd.f"
  22.       INTEGER II1, II2, II3
  23. # 11 "linpackd.f"
  24. # 11 "linpackd.f"
  25. # 11 "linpackd.f"
  26.       DOUBLE PRECISION DD1 (129), DD2 (129), DD3 (665)
  27.       LDA = 201
  28. # 12 "linpackd.f"
  29.       LDAA = 200
  30. c
  31. # 14 "linpackd.f"
  32.       N = 100
  33. # 15 "linpackd.f"
  34.       cray = .056
  35.       WRITE (6, 1)
  36. # 17 "linpackd.f"
  37.     1 format(' Please send the results of this run to:'//
  38.      $       ' Jack J. Dongarra'/
  39.      $       ' Computer Science Department'/
  40.      $       ' University of Tennessee'/
  41.      $       ' Knoxville, Tennessee 37996-1300'//
  42.      $       ' Fax: 615-974-8296'//
  43.      $       ' Internet: dongarra@cs.utk.edu'/)
  44.       OPS = (2.D0 * 1000000) / 3.D0 + 2.D0 * 10000
  45. c
  46. # 26 "linpackd.f"
  47.          call matgen(a,lda,n,b,norma)
  48.          t1 = second()
  49.          call dgefa(a,lda,n,ipvt,info)
  50.          TIME(1,1) = SECOND(  ) - T1
  51. # 30 "linpackd.f"
  52.          t1 = second()
  53.          CALL DGESL (A,LDA,N,IPVT,B,0)
  54. # 32 "linpackd.f"
  55.          TIME(1,2) = SECOND(  ) - T1
  56. # 33 "linpackd.f"
  57.          TOTAL = TIME(1,1) + TIME(1,2)
  58. # 37 "linpackd.f"
  59.          II1 = MOD (N, 4)
  60. # 37 "linpackd.f"
  61.          DO 3 I=1,II1
  62. # 38 "linpackd.f"
  63.             X(I) = B(I)
  64. # 39 "linpackd.f"
  65.     3       CONTINUE
  66. c
  67. c     compute a residual to verify results.
  68. c
  69. # 37 "linpackd.f"
  70. C$DOACROSS IF(N .GT. 166),SHARE(II1,N,X,B),LOCAL(I)
  71. # 37 "linpackd.f"
  72.          DO 4 I=II1+1,N,4
  73. # 38 "linpackd.f"
  74.             x(i) = b(i)
  75.             X(I+1) = B(I+1)
  76. # 38 "linpackd.f"
  77.             X(I+2) = B(I+2)
  78. # 38 "linpackd.f"
  79.             X(I+3) = B(I+3)
  80. # 39 "linpackd.f"
  81.     4    CONTINUE
  82. # 40 "linpackd.f"
  83.          call matgen(a,lda,n,b,norma)
  84. # 41 "linpackd.f"
  85.          II2 = MOD (N, 4)
  86. # 41 "linpackd.f"
  87.          DO 5 I=1,II2
  88. # 42 "linpackd.f"
  89.             B(I) = -B(I)
  90. # 43 "linpackd.f"
  91.     5       CONTINUE
  92. # 41 "linpackd.f"
  93. C$DOACROSS IF(N .GT. 125),SHARE(II2,N,B),LOCAL(I)
  94. # 41 "linpackd.f"
  95.          DO 6 I=II2+1,N,4
  96. # 42 "linpackd.f"
  97.             B(I) = -B(I)
  98. # 42 "linpackd.f"
  99.             B(I+1) = -B(I+1)
  100. # 42 "linpackd.f"
  101.             B(I+2) = -B(I+2)
  102. # 42 "linpackd.f"
  103.             B(I+3) = -B(I+3)
  104. # 43 "linpackd.f"
  105.     6    CONTINUE
  106. # 44 "linpackd.f"
  107.          call dmxpy(n,b,n,lda,x,a)
  108.          RESID = 0.
  109. # 46 "linpackd.f"
  110.          NORMX = 0.
  111. # 47 "linpackd.f"
  112.          II3 = MOD (N, 4)
  113. # 47 "linpackd.f"
  114.          DO 2 I=1,II3
  115. # 48 "linpackd.f"
  116.             RESID = DMAX1 (RESID, DABS (B(I)))
  117. # 49 "linpackd.f"
  118.             NORMX = DMAX1 (NORMX, DABS (X(I)))
  119. # 50 "linpackd.f"
  120.     2    CONTINUE
  121. # 47 "linpackd.f"
  122.          DO 30 I=II3+1,N,4
  123. # 48 "linpackd.f"
  124.             resid = dmax1( resid, dabs(b(i)) )
  125.             normx = dmax1( normx, dabs(x(i)) )
  126. # 48 "linpackd.f"
  127.             RESID = DMAX1 (RESID, DABS (B(I+1)))
  128. # 49 "linpackd.f"
  129.             NORMX = DMAX1 (NORMX, DABS (X(I+1)))
  130. # 48 "linpackd.f"
  131.             RESID = DMAX1 (RESID, DABS (B(I+2)))
  132. # 49 "linpackd.f"
  133.             NORMX = DMAX1 (NORMX, DABS (X(I+2)))
  134. # 48 "linpackd.f"
  135.             RESID = DMAX1 (RESID, DABS (B(I+3)))
  136. # 49 "linpackd.f"
  137.             NORMX = DMAX1 (NORMX, DABS (X(I+3)))
  138. # 50 "linpackd.f"
  139.    30    continue
  140.          EPS = EPSLON( 1.D0 )
  141. # 52 "linpackd.f"
  142.          residn = resid/( n*norma*normx*eps )
  143.          WRITE (6, 40)
  144. # 54 "linpackd.f"
  145.    40    format('     norm. resid      resid           machep',
  146.      $          '         x(1)          x(n)')
  147.          WRITE (6, 50) RESIDN, RESID, EPS, X(1), X(N)
  148. # 57 "linpackd.f"
  149.    50    format(1p5e16.8)
  150. c
  151.          WRITE (6, 60) N
  152. # 60 "linpackd.f"
  153.    60    format(//'    times are reported for matrices of order ',i5)
  154.          WRITE (6, 70)
  155. # 62 "linpackd.f"
  156.    70    format(6x,'dgefa',6x,'dgesl',6x,'total',5x,'mflops',7x,'unit',
  157.      $         6x,'ratio')
  158. c
  159.          TIME(1,3) = TOTAL
  160. # 66 "linpackd.f"
  161.          TIME(1,4) = OPS / (1.0D6 * TOTAL)
  162. # 67 "linpackd.f"
  163.          TIME(1,5) = 2.D0 / TIME(1,4)
  164. # 68 "linpackd.f"
  165.          TIME(1,6) = TOTAL / CRAY
  166. # 69 "linpackd.f"
  167.          WRITE (6, 80) LDA
  168. # 70 "linpackd.f"
  169.    80    format(' times for array with leading dimension of',i4)
  170.          WRITE (6, 110) (TIME(1,I), I=1,6)
  171. c
  172. # 73 "linpackd.f"
  173.          call matgen(a,lda,n,b,norma)
  174.          t1 = second()
  175.          call dgefa(a,lda,n,ipvt,info)
  176.          TIME(2,1) = SECOND(  ) - T1
  177. # 77 "linpackd.f"
  178.          t1 = second()
  179.          CALL DGESL (A,LDA,N,IPVT,B,0)
  180. # 79 "linpackd.f"
  181.          TIME(2,2) = SECOND(  ) - T1
  182. # 80 "linpackd.f"
  183.          TOTAL = TIME(2,1) + TIME(2,2)
  184. # 81 "linpackd.f"
  185.          TIME(2,3) = TOTAL
  186. # 82 "linpackd.f"
  187.          TIME(2,4) = OPS / (1.0D6 * TOTAL)
  188. # 83 "linpackd.f"
  189.          TIME(2,5) = 2.D0 / TIME(2,4)
  190. # 84 "linpackd.f"
  191.          TIME(2,6) = TOTAL / CRAY
  192. c
  193. # 86 "linpackd.f"
  194.          call matgen(a,lda,n,b,norma)
  195.          t1 = second()
  196.          call dgefa(a,lda,n,ipvt,info)
  197.          TIME(3,1) = SECOND(  ) - T1
  198. # 90 "linpackd.f"
  199.          t1 = second()
  200.          CALL DGESL (A,LDA,N,IPVT,B,0)
  201. # 92 "linpackd.f"
  202.          TIME(3,2) = SECOND(  ) - T1
  203. # 93 "linpackd.f"
  204.          TOTAL = TIME(3,1) + TIME(3,2)
  205. # 94 "linpackd.f"
  206.          TIME(3,3) = TOTAL
  207. # 95 "linpackd.f"
  208.          TIME(3,4) = OPS / (1.0D6 * TOTAL)
  209. # 96 "linpackd.f"
  210.          TIME(3,5) = 2.D0 / TIME(3,4)
  211. # 97 "linpackd.f"
  212.          TIME(3,6) = TOTAL / CRAY
  213. c
  214. # 100 "linpackd.f"
  215.          TM2 = 0
  216. # 101 "linpackd.f"
  217.          t1 = second()
  218.          DO 90 I=1,10
  219. # 103 "linpackd.f"
  220.             tm = second()
  221.             call matgen(a,lda,n,b,norma)
  222.             tm2 = tm2 + second() - tm
  223.             call dgefa(a,lda,n,ipvt,info)
  224.    90    continue
  225.          TIME(4,1) = (SECOND(  ) - T1 - TM2) / 10
  226. # 109 "linpackd.f"
  227.          t1 = second()
  228.          DO 100 I=1,10
  229. # 111 "linpackd.f"
  230.             CALL DGESL (A,LDA,N,IPVT,B,0)
  231. # 112 "linpackd.f"
  232.   100    continue
  233.          TIME(4,2) = (SECOND(  ) - T1) / 10
  234. # 114 "linpackd.f"
  235.          TOTAL = TIME(4,1) + TIME(4,2)
  236. # 115 "linpackd.f"
  237.          TIME(4,3) = TOTAL
  238. # 116 "linpackd.f"
  239.          TIME(4,4) = OPS / (1.0D6 * TOTAL)
  240. # 117 "linpackd.f"
  241.          TIME(4,5) = 2.D0 / TIME(4,4)
  242. # 118 "linpackd.f"
  243.          TIME(4,6) = TOTAL / CRAY
  244. c
  245. # 120 "linpackd.f"
  246.          WRITE (6, 110) (TIME(2,I), I=1,6)
  247. # 121 "linpackd.f"
  248.          WRITE (6, 110) (TIME(3,I), I=1,6)
  249. # 122 "linpackd.f"
  250.          WRITE (6, 110) (TIME(4,I), I=1,6)
  251. # 123 "linpackd.f"
  252.   110    format(6(1pe11.3))
  253. c
  254.          call matgen(aa,ldaa,n,b,norma)
  255.          t1 = second()
  256.          call dgefa(aa,ldaa,n,ipvt,info)
  257.          TIME(5,1) = SECOND(  ) - T1
  258. # 129 "linpackd.f"
  259.          t1 = second()
  260.          CALL DGESL (AA,LDAA,N,IPVT,B,0)
  261. # 131 "linpackd.f"
  262.          TIME(5,2) = SECOND(  ) - T1
  263. # 132 "linpackd.f"
  264.          TOTAL = TIME(5,1) + TIME(5,2)
  265. # 133 "linpackd.f"
  266.          TIME(5,3) = TOTAL
  267. # 134 "linpackd.f"
  268.          TIME(5,4) = OPS / (1.0D6 * TOTAL)
  269. # 135 "linpackd.f"
  270.          TIME(5,5) = 2.D0 / TIME(5,4)
  271. # 136 "linpackd.f"
  272.          TIME(5,6) = TOTAL / CRAY
  273. c
  274. # 138 "linpackd.f"
  275.          call matgen(aa,ldaa,n,b,norma)
  276.          t1 = second()
  277.          call dgefa(aa,ldaa,n,ipvt,info)
  278.          TIME(6,1) = SECOND(  ) - T1
  279. # 142 "linpackd.f"
  280.          t1 = second()
  281.          CALL DGESL (AA,LDAA,N,IPVT,B,0)
  282. # 144 "linpackd.f"
  283.          TIME(6,2) = SECOND(  ) - T1
  284. # 145 "linpackd.f"
  285.          TOTAL = TIME(6,1) + TIME(6,2)
  286. # 146 "linpackd.f"
  287.          TIME(6,3) = TOTAL
  288. # 147 "linpackd.f"
  289.          TIME(6,4) = OPS / (1.0D6 * TOTAL)
  290. # 148 "linpackd.f"
  291.          TIME(6,5) = 2.D0 / TIME(6,4)
  292. # 149 "linpackd.f"
  293.          TIME(6,6) = TOTAL / CRAY
  294. c
  295. # 151 "linpackd.f"
  296.          call matgen(aa,ldaa,n,b,norma)
  297.          t1 = second()
  298.          call dgefa(aa,ldaa,n,ipvt,info)
  299.          TIME(7,1) = SECOND(  ) - T1
  300. # 155 "linpackd.f"
  301.          t1 = second()
  302.          CALL DGESL (AA,LDAA,N,IPVT,B,0)
  303. # 157 "linpackd.f"
  304.          TIME(7,2) = SECOND(  ) - T1
  305. # 158 "linpackd.f"
  306.          TOTAL = TIME(7,1) + TIME(7,2)
  307. # 159 "linpackd.f"
  308.          TIME(7,3) = TOTAL
  309. # 160 "linpackd.f"
  310.          TIME(7,4) = OPS / (1.0D6 * TOTAL)
  311. # 161 "linpackd.f"
  312.          TIME(7,5) = 2.D0 / TIME(7,4)
  313. # 162 "linpackd.f"
  314.          TIME(7,6) = TOTAL / CRAY
  315. c
  316. # 165 "linpackd.f"
  317.          TM2 = 0
  318. # 166 "linpackd.f"
  319.          t1 = second()
  320.          DO 120 I=1,10
  321. # 168 "linpackd.f"
  322.             tm = second()
  323.             call matgen(aa,ldaa,n,b,norma)
  324.             tm2 = tm2 + second() - tm
  325.             call dgefa(aa,ldaa,n,ipvt,info)
  326.   120    continue
  327.          TIME(8,1) = (SECOND(  ) - T1 - TM2) / 10
  328. # 174 "linpackd.f"
  329.          t1 = second()
  330.          DO 130 I=1,10
  331. # 176 "linpackd.f"
  332.             CALL DGESL (AA,LDAA,N,IPVT,B,0)
  333. # 177 "linpackd.f"
  334.   130    continue
  335.          TIME(8,2) = (SECOND(  ) - T1) / 10
  336. # 179 "linpackd.f"
  337.          TOTAL = TIME(8,1) + TIME(8,2)
  338. # 180 "linpackd.f"
  339.          TIME(8,3) = TOTAL
  340. # 181 "linpackd.f"
  341.          TIME(8,4) = OPS / (1.0D6 * TOTAL)
  342. # 182 "linpackd.f"
  343.          TIME(8,5) = 2.D0 / TIME(8,4)
  344. # 183 "linpackd.f"
  345.          TIME(8,6) = TOTAL / CRAY
  346. c
  347. # 185 "linpackd.f"
  348.          WRITE (6, 140) LDAA
  349. # 186 "linpackd.f"
  350.   140    format(/' times for array with leading dimension of',i4)
  351.          WRITE (6, 110) (TIME(5,I), I=1,6)
  352. # 188 "linpackd.f"
  353.          WRITE (6, 110) (TIME(6,I), I=1,6)
  354. # 189 "linpackd.f"
  355.          WRITE (6, 110) (TIME(7,I), I=1,6)
  356. # 190 "linpackd.f"
  357.          WRITE (6, 110) (TIME(8,I), I=1,6)
  358. # 191 "linpackd.f"
  359.       stop
  360.       end
  361. # 193 "linpackd.f"
  362.       subroutine matgen(a,lda,n,b,norma)
  363.       DOUBLEPRECISION A(LDA,1), B(1), NORMA
  364. # 195 "linpackd.f"
  365.       INTEGER II1, II2, II3, II4
  366. # 195 "linpackd.f"
  367.       DOUBLE PRECISION DD1
  368. # 195 "linpackd.f"
  369. # 195 "linpackd.f"
  370. c
  371.       INIT = 1325
  372. # 197 "linpackd.f"
  373.       NORMA = 0.
  374. # 198 "linpackd.f"
  375.       DO 4 II1=1,N
  376. # 199 "linpackd.f"
  377.       II2 = MOD (N, 4)
  378. # 199 "linpackd.f"
  379.       DO 3 I=1,II2
  380. # 200 "linpackd.f"
  381.             INIT = MOD (INIT * 3125, 65536)
  382. # 201 "linpackd.f"
  383.             A(I,II1) = (INIT - 32768.) / 16384.
  384. # 202 "linpackd.f"
  385.             NORMA = DMAX1 (DABS (A(I,II1)), NORMA)
  386. # 203 "linpackd.f"
  387.     3       CONTINUE
  388. # 199 "linpackd.f"
  389.          DO 4 I=II2+1,N,4
  390. # 200 "linpackd.f"
  391.             INIT = MOD (INIT * 3125, 65536)
  392. # 201 "linpackd.f"
  393.             A(I,II1) = (INIT - 32768.) / 16384.
  394. # 202 "linpackd.f"
  395.             NORMA = DMAX1 (DABS (A(I,II1)), NORMA)
  396. # 200 "linpackd.f"
  397.             INIT = MOD (INIT * 3125, 65536)
  398. # 201 "linpackd.f"
  399.             A(I+1,II1) = (INIT - 32768.) / 16384.
  400. # 202 "linpackd.f"
  401.             NORMA = DMAX1 (DABS (A(I+1,II1)), NORMA)
  402. # 200 "linpackd.f"
  403.             INIT = MOD (INIT * 3125, 65536)
  404. # 201 "linpackd.f"
  405.             A(I+2,II1) = (INIT - 32768.) / 16384.
  406. # 202 "linpackd.f"
  407.             NORMA = DMAX1 (DABS (A(I+2,II1)), NORMA)
  408. # 200 "linpackd.f"
  409.             INIT = MOD (INIT * 3125, 65536)
  410. # 201 "linpackd.f"
  411.             A(I+3,II1) = (INIT - 32768.) / 16384.
  412. # 202 "linpackd.f"
  413.             NORMA = DMAX1 (DABS (A(I+3,II1)), NORMA)
  414. # 203 "linpackd.f"
  415.     4    CONTINUE
  416. # 198 "linpackd.f"
  417.          II3 = MOD (N, 4)
  418. # 198 "linpackd.f"
  419.          DO 5 II1=1,II3
  420. # 206 "linpackd.f"
  421.           B(II1) = 0.
  422. # 207 "linpackd.f"
  423.     5     CONTINUE
  424. # 198 "linpackd.f"
  425. C$DOACROSS IF(N .GT. 250),SHARE(II3,N,B),LOCAL(II1)
  426. # 198 "linpackd.f"
  427.           DO 6 II1=II3+1,N,4
  428. # 206 "linpackd.f"
  429.           B(II1) = 0.
  430. # 206 "linpackd.f"
  431.           B(II1+1) = 0.
  432. # 206 "linpackd.f"
  433.           B(II1+2) = 0.
  434. # 206 "linpackd.f"
  435.           B(II1+3) = 0.
  436. # 207 "linpackd.f"
  437.     6     CONTINUE
  438. # 209 "linpackd.f"
  439. C$DOACROSS IF(N .GT. 9),SHARE(N,B,A),LOCAL(DD1,II4,I,J)
  440. # 209 "linpackd.f"
  441.          DO 7 I=1,N
  442. # 209 "linpackd.f"
  443.          DD1 = B(I)
  444. # 208 "linpackd.f"
  445.          II4 = MOD (N, 4)
  446. # 208 "linpackd.f"
  447.          DO 2 J=1,II4
  448. # 210 "linpackd.f"
  449.             DD1 = DD1 + A(I,J)
  450. # 212 "linpackd.f"
  451.     2 CONTINUE
  452. # 208 "linpackd.f"
  453.       DO 50 J=II4+1,N,4
  454. # 210 "linpackd.f"
  455.             DD1 = DD1 + A(I,J)
  456. # 210 "linpackd.f"
  457.             DD1 = DD1 + A(I,J+1)
  458. # 210 "linpackd.f"
  459.             DD1 = DD1 + A(I,J+2)
  460. # 210 "linpackd.f"
  461.             DD1 = DD1 + A(I,J+3)
  462. # 212 "linpackd.f"
  463.    50 continue
  464.       B(I) = DD1
  465. # 211 "linpackd.f"
  466.     7    CONTINUE
  467. # 213 "linpackd.f"
  468.       return
  469.       end
  470. # 215 "linpackd.f"
  471.       subroutine dgefa(a,lda,n,ipvt,info)
  472.       INTEGER LDA, N, IPVT(1), INFO
  473. # 217 "linpackd.f"
  474.       DOUBLEPRECISION A(LDA,1)
  475. c
  476. c     dgefa factors a double precision matrix by gaussian elimination.
  477. c
  478. c     dgefa is usually called by dgeco, but it can be called
  479. c     directly with a saving in time if  rcond  is not needed.
  480. c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  481. c
  482. c     on entry
  483. c
  484. c        a       double precision(lda, n)
  485. c                the matrix to be factored.
  486. c
  487. c        lda     integer
  488. c                the leading dimension of the array  a .
  489. c
  490. c        n       integer
  491. c                the order of the matrix  a .
  492. c
  493. c     on return
  494. c
  495. c        a       an upper triangular matrix and the multipliers
  496. c                which were used to obtain it.
  497. c                the factorization can be written  a = l*u  where
  498. c                l  is a product of permutation and unit lower
  499. c                triangular matrices and  u  is upper triangular.
  500. c
  501. c        ipvt    integer(n)
  502. c                an integer vector of pivot indices.
  503. c
  504. c        info    integer
  505. c                = 0  normal value.
  506. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  507. c                     condition for this subroutine, but it does
  508. c                     indicate that dgesl or dgedi will divide by zero
  509. c                     if called.  use  rcond  in dgeco for a reliable
  510. c                     indication of singularity.
  511. c
  512. c     linpack. this version dated 08/14/78 .
  513. c     cleve moler, university of new mexico, argonne national lab.
  514. c
  515. c     subroutines and functions
  516. c
  517. c     blas daxpy,dscal,idamax
  518. c
  519. c     internal variables
  520. c
  521. # 264 "linpackd.f"
  522.       double precision t
  523.       integer idamax,j,k,kp1,l,nm1
  524. # 266 "linpackd.f"
  525.       INTEGER II1, II2, II3, II4
  526. # 266 "linpackd.f"
  527. c
  528. c
  529. c     gaussian elimination with partial pivoting
  530. c
  531.       INFO = 0
  532. # 273 "linpackd.f"
  533.       DO 60 K=1,(N-1)
  534. c
  535. c        find l = pivot index
  536. c
  537. # 278 "linpackd.f"
  538.          L = IDAMAX( N - K + 1, A(K,K), 1 ) + K - 1
  539. # 279 "linpackd.f"
  540.          ipvt(k) = l
  541. c
  542. c        zero pivot implies this column already triangularized
  543. c
  544.          IF (A(L,K) .NE. 0.D0) THEN
  545. c
  546. c           interchange if necessary
  547. c
  548. # 287 "linpackd.f"
  549.             IF (L .NE. K) THEN
  550. # 288 "linpackd.f"
  551.                t = a(l,k)
  552.                a(l,k) = a(k,k)
  553.                a(k,k) = t
  554. # 287 "linpackd.f"
  555.                END IF
  556. c
  557. c           compute multipliers
  558. c
  559. # 295 "linpackd.f"
  560.             T = -1.D0 / A(K,K)
  561. # 296 "linpackd.f"
  562.             CALL DSCAL (N - K,T,A(K+1,K),1)
  563. # 302 "linpackd.f"
  564.                IF (L .NE. K) THEN
  565. # 302 "linpackd.f"
  566.                II1 = K + 1
  567. # 302 "linpackd.f"
  568.                II2 = N - K
  569. c
  570. c           row elimination with column indexing
  571. c
  572. # 300 "linpackd.f"
  573.             DO 30 J=II1,N
  574. # 301 "linpackd.f"
  575.                t = a(l,j)
  576.                   a(l,j) = a(k,j)
  577.                   a(k,j) = t
  578. # 306 "linpackd.f"
  579.                CALL DAXPY (II2,T,A(II1,K),1,A(II1,J),1)
  580. # 307 "linpackd.f"
  581.    30       continue
  582. # 302 "linpackd.f"
  583.             ELSE
  584. # 302 "linpackd.f"
  585.             II3 = K + 1
  586. # 302 "linpackd.f"
  587.             II4 = N - K
  588. # 302 "linpackd.f"
  589.             DO 4 J=II3,N
  590. # 301 "linpackd.f"
  591.                T = A(L,J)
  592. # 306 "linpackd.f"
  593.                CALL DAXPY (II4,T,A(II3,K),1,A(II3,J),1)
  594. # 302 "linpackd.f"
  595.     4          CONTINUE
  596. # 302 "linpackd.f"
  597.                END IF
  598. # 283 "linpackd.f"
  599.                ELSE
  600. # 310 "linpackd.f"
  601.             info = k
  602. # 283 "linpackd.f"
  603.             END IF
  604. # 312 "linpackd.f"
  605.    60 continue
  606. # 314 "linpackd.f"
  607.       ipvt(n) = n
  608.       IF (A(N,N) .EQ. 0.D0) INFO = N
  609. # 316 "linpackd.f"
  610.       return
  611.       end
  612. # 318 "linpackd.f"
  613.       subroutine dgesl(a,lda,n,ipvt,b,job)
  614.       INTEGER LDA, N, IPVT(1), JOB
  615. # 320 "linpackd.f"
  616.       DOUBLEPRECISION A(LDA,1), B(1)
  617. c
  618. c     dgesl solves the double precision system
  619. c     a * x = b  or  trans(a) * x = b
  620. c     using the factors computed by dgeco or dgefa.
  621. c
  622. c     on entry
  623. c
  624. c        a       double precision(lda, n)
  625. c                the output from dgeco or dgefa.
  626. c
  627. c        lda     integer
  628. c                the leading dimension of the array  a .
  629. c
  630. c        n       integer
  631. c                the order of the matrix  a .
  632. c
  633. c        ipvt    integer(n)
  634. c                the pivot vector from dgeco or dgefa.
  635. c
  636. c        b       double precision(n)
  637. c                the right hand side vector.
  638. c
  639. c        job     integer
  640. c                = 0         to solve  a*x = b ,
  641. c                = nonzero   to solve  trans(a)*x = b  where
  642. c                            trans(a)  is the transpose.
  643. c
  644. c     on return
  645. c
  646. c        b       the solution vector  x .
  647. c
  648. c     error condition
  649. c
  650. c        a division by zero will occur if the input factor contains a
  651. c        zero on the diagonal.  technically this indicates singularity
  652. c        but it is often caused by improper arguments or improper
  653. c        setting of lda .  it will not occur if the subroutines are
  654. c        called correctly and if dgeco has set rcond .gt. 0.0
  655. c        or dgefa has set info .eq. 0 .
  656. c
  657. c     to compute  inverse(a) * c  where  c  is a matrix
  658. c     with  p  columns
  659. c           call dgeco(a,lda,n,ipvt,rcond,z)
  660. c           if (rcond is too small) go to ...
  661. c           do 10 j = 1, p
  662. c              call dgesl(a,lda,n,ipvt,c(1,j),0)
  663. c        10 continue
  664. c
  665. c     linpack. this version dated 08/14/78 .
  666. c     cleve moler, university of new mexico, argonne national lab.
  667. c
  668. c     subroutines and functions
  669. c
  670. c     blas daxpy,ddot
  671. c
  672. c     internal variables
  673. c
  674. # 378 "linpackd.f"
  675.       double precision ddot,t
  676.       integer k,kb,l,nm1
  677. c
  678. # 382 "linpackd.f"
  679.       IF (JOB .EQ. 0) THEN
  680. c
  681. c        job = 0 , solve  a * x = b
  682. c        first solve  l*y = b
  683. c
  684. # 388 "linpackd.f"
  685.          DO 20 K=1,(N-1)
  686. # 390 "linpackd.f"
  687.             T = B(IPVT(K))
  688. # 391 "linpackd.f"
  689.             IF (IPVT(K) .NE. K) THEN
  690. # 392 "linpackd.f"
  691.                B(IPVT(K)) = B(K)
  692. # 393 "linpackd.f"
  693.                b(k) = t
  694. # 391 "linpackd.f"
  695.                END IF
  696. # 395 "linpackd.f"
  697.             CALL DAXPY (N - K,T,A(K+1,K),1,B(K+1),1)
  698. # 396 "linpackd.f"
  699.    20    continue
  700. c
  701. c        now solve  u*x = y
  702. c
  703. # 401 "linpackd.f"
  704.          DO 40 KB=1,N
  705. # 402 "linpackd.f"
  706.             K = N - KB + 1
  707. # 403 "linpackd.f"
  708.             b(k) = b(k)/a(k,k)
  709.             T = -B(K)
  710. # 405 "linpackd.f"
  711.             CALL DAXPY (K - 1,T,A(1,K),1,B(1),1)
  712. # 406 "linpackd.f"
  713.    40    continue
  714. # 382 "linpackd.f"
  715.          ELSE
  716. c
  717. c        job = nonzero, solve  trans(a) * x = b
  718. c        first solve  trans(u)*y = b
  719. c
  720. # 413 "linpackd.f"
  721.          DO 60 K=1,N
  722. # 414 "linpackd.f"
  723.             T = DDOT( K - 1, A(1,K), 1, B(1), 1 )
  724. # 415 "linpackd.f"
  725.             b(k) = (b(k) - t)/a(k,k)
  726.    60    continue
  727. c
  728. c        now solve trans(l)*x = y
  729. c
  730. # 421 "linpackd.f"
  731.          DO 80 KB=1,(N-1)
  732. # 423 "linpackd.f"
  733.             B((N-KB)) = B((N-KB)) + DDOT( N - (N - KB), A((N-KB)+1,(N-KB
  734.      X        )), 1, B((N-KB)+1), 1 )
  735. # 425 "linpackd.f"
  736.             IF (IPVT((N-KB)) .NE. (N - KB)) THEN
  737. # 426 "linpackd.f"
  738.                T = B(IPVT((N-KB)))
  739. # 427 "linpackd.f"
  740.                B(IPVT((N-KB))) = B((N-KB))
  741. # 428 "linpackd.f"
  742.                B((N-KB)) = T
  743. # 425 "linpackd.f"
  744.                END IF
  745. # 430 "linpackd.f"
  746.    80    continue
  747. # 382 "linpackd.f"
  748.          END IF
  749. # 433 "linpackd.f"
  750.       return
  751.       end
  752. # 435 "linpackd.f"
  753.       subroutine daxpy(n,da,dx,incx,dy,incy)
  754. c
  755. c     constant times a vector plus a vector.
  756. c     jack dongarra, linpack, 3/11/78.
  757. c
  758.       DOUBLEPRECISION DX(1), DY(1), DA
  759. # 441 "linpackd.f"
  760.       integer i,incx,incy,ix,iy,m,mp1,n
  761. # 442 "linpackd.f"
  762.       INTEGER II1, II2, II3, II4, II5, II6
  763. # 442 "linpackd.f"
  764. c
  765.       IF (N .LE. 0) RETURN
  766. # 444 "linpackd.f"
  767.       IF (DA .EQ. 0.D0) RETURN
  768. # 445 "linpackd.f"
  769.       IF (INCX .NE. 1 .OR. INCY .NE. 1) THEN
  770. c
  771. c        code for unequal increments or equal increments
  772. c          not equal to 1
  773. c
  774. # 450 "linpackd.f"
  775.       IX = 1
  776. # 451 "linpackd.f"
  777.       IY = 1
  778. # 452 "linpackd.f"
  779.       IF (INCX .LT. 0) IX = (1 - N) * INCX + 1
  780. # 453 "linpackd.f"
  781.       IF (INCY .LT. 0) IY = (1 - N) * INCY + 1
  782. # 457 "linpackd.f"
  783.       II1 = IY
  784. # 456 "linpackd.f"
  785.       II2 = IX
  786. # 454 "linpackd.f"
  787.       IF (INCY .EQ. 0) THEN
  788. # 454 "linpackd.f"
  789.       II4 = MOD (N, 4)
  790. # 454 "linpackd.f"
  791.       DO 3 I=1,II4
  792. # 455 "linpackd.f"
  793.         DY(IY) = DY(IY) + DA * DX(IX)
  794. # 456 "linpackd.f"
  795.         IX = IX + INCX
  796. # 457 "linpackd.f"
  797.         IY = IY + INCY
  798. # 458 "linpackd.f"
  799.     3 CONTINUE
  800. # 454 "linpackd.f"
  801.       DO 10 I=II4+1,N,4
  802. # 455 "linpackd.f"
  803.         dy(iy) = dy(iy) + da*dx(ix)
  804.         ix = ix + incx
  805.         iy = iy + incy
  806. # 455 "linpackd.f"
  807.         DY(IY) = DY(IY) + DA * DX(IX)
  808. # 456 "linpackd.f"
  809.         IX = IX + INCX
  810. # 457 "linpackd.f"
  811.         IY = IY + INCY
  812. # 455 "linpackd.f"
  813.         DY(IY) = DY(IY) + DA * DX(IX)
  814. # 456 "linpackd.f"
  815.         IX = IX + INCX
  816. # 457 "linpackd.f"
  817.         IY = IY + INCY
  818. # 455 "linpackd.f"
  819.         DY(IY) = DY(IY) + DA * DX(IX)
  820. # 456 "linpackd.f"
  821.         IX = IX + INCX
  822. # 457 "linpackd.f"
  823.         IY = IY + INCY
  824. # 458 "linpackd.f"
  825.    10 continue
  826. # 454 "linpackd.f"
  827.       ELSE
  828. # 454 "linpackd.f"
  829.       II5 = MOD (N, 4)
  830. # 454 "linpackd.f"
  831.       DO 4 I=1,II5
  832. # 455 "linpackd.f"
  833.         II3 = II1 + (I - 1) * INCY
  834. # 455 "linpackd.f"
  835.         DY(II3) = DY(II3) + DA * DX(II2+(I-1)*INCX)
  836. # 458 "linpackd.f"
  837.     4 CONTINUE
  838. # 454 "linpackd.f"
  839. C$DOACROSS IF(N .GT. 25),SHARE(II5,N,II1,INCY,DY,DA,II2,INCX,DX),LOCAL(
  840. C$&  II3,I)
  841. # 454 "linpackd.f"
  842.       DO 2 I=II5+1,N,4
  843. # 455 "linpackd.f"
  844.         II3 = II1 + (I - 1) * INCY
  845. # 455 "linpackd.f"
  846.         DY(II3) = DY(II3) + DA * DX(II2+(I-1)*INCX)
  847. # 455 "linpackd.f"
  848.         II3 = II1 + I * INCY
  849. # 455 "linpackd.f"
  850.         DY(II3) = DY(II3) + DA * DX(II2+I*INCX)
  851. # 455 "linpackd.f"
  852.         II3 = II1 + (I + 1) * INCY
  853. # 455 "linpackd.f"
  854.         DY(II3) = DY(II3) + DA * DX(II2+(I+1)*INCX)
  855. # 455 "linpackd.f"
  856.         II3 = II1 + (I + 2) * INCY
  857. # 455 "linpackd.f"
  858.         DY(II3) = DY(II3) + DA * DX(II2+(I+2)*INCX)
  859. # 458 "linpackd.f"
  860.     2 CONTINUE
  861. # 454 "linpackd.f"
  862.       END IF
  863. # 459 "linpackd.f"
  864.       return
  865. # 445 "linpackd.f"
  866.       END IF
  867. c
  868. c        code for both increments equal to 1
  869. c
  870. # 464 "linpackd.f"
  871.       II6 = MOD (N, 4)
  872. # 464 "linpackd.f"
  873.       DO 5 I=1,II6
  874. # 465 "linpackd.f"
  875.         DY(I) = DY(I) + DA * DX(I)
  876. # 466 "linpackd.f"
  877.     5   CONTINUE
  878. # 464 "linpackd.f"
  879. C$DOACROSS IF(N .GT. 83),SHARE(II6,N,DY,DA,DX),LOCAL(I)
  880. # 464 "linpackd.f"
  881.       DO 6 I=II6+1,N,4
  882. # 465 "linpackd.f"
  883.         dy(i) = dy(i) + da*dx(i)
  884.         DY(I+1) = DY(I+1) + DA * DX(I+1)
  885. # 465 "linpackd.f"
  886.         DY(I+2) = DY(I+2) + DA * DX(I+2)
  887. # 465 "linpackd.f"
  888.         DY(I+3) = DY(I+3) + DA * DX(I+3)
  889. # 466 "linpackd.f"
  890.     6 CONTINUE
  891. # 467 "linpackd.f"
  892.       return
  893.       end
  894. # 469 "linpackd.f"
  895.       double precision function ddot(n,dx,incx,dy,incy)
  896. c
  897. c     forms the dot product of two vectors.
  898. c     jack dongarra, linpack, 3/11/78.
  899. c
  900.       DOUBLEPRECISION DX(1), DY(1), DTEMP
  901. # 475 "linpackd.f"
  902.       integer i,incx,incy,ix,iy,m,mp1,n
  903. # 476 "linpackd.f"
  904.       INTEGER II3, II4
  905. # 476 "linpackd.f"
  906. c
  907.       DDOT = 0.D0
  908. # 478 "linpackd.f"
  909.       DTEMP = 0.D0
  910. # 479 "linpackd.f"
  911.       IF (N .LE. 0) RETURN
  912. # 480 "linpackd.f"
  913.       IF (INCX .NE. 1 .OR. INCY .NE. 1) THEN
  914. c
  915. c        code for unequal increments or equal increments
  916. c          not equal to 1
  917. c
  918. # 485 "linpackd.f"
  919.       IX = 1
  920. # 486 "linpackd.f"
  921.       IY = 1
  922. # 487 "linpackd.f"
  923.       IF (INCX .LT. 0) IX = (1 - N) * INCX + 1
  924. # 488 "linpackd.f"
  925.       IF (INCY .LT. 0) IY = (1 - N) * INCY + 1
  926. # 489 "linpackd.f"
  927.       II3 = MOD (N, 4)
  928. # 489 "linpackd.f"
  929.       DO 2 I=1,II3
  930. # 490 "linpackd.f"
  931.         DTEMP = DTEMP + DX(IX) * DY(IY)
  932. # 491 "linpackd.f"
  933.         IX = IX + INCX
  934. # 492 "linpackd.f"
  935.         IY = IY + INCY
  936. # 493 "linpackd.f"
  937.     2 CONTINUE
  938. # 489 "linpackd.f"
  939.       DO 10 I=II3+1,N,4
  940. # 490 "linpackd.f"
  941.         dtemp = dtemp + dx(ix)*dy(iy)
  942.         ix = ix + incx
  943.         iy = iy + incy
  944. # 490 "linpackd.f"
  945.         DTEMP = DTEMP + DX(IX) * DY(IY)
  946. # 491 "linpackd.f"
  947.         IX = IX + INCX
  948. # 492 "linpackd.f"
  949.         IY = IY + INCY
  950. # 490 "linpackd.f"
  951.         DTEMP = DTEMP + DX(IX) * DY(IY)
  952. # 491 "linpackd.f"
  953.         IX = IX + INCX
  954. # 492 "linpackd.f"
  955.         IY = IY + INCY
  956. # 490 "linpackd.f"
  957.         DTEMP = DTEMP + DX(IX) * DY(IY)
  958. # 491 "linpackd.f"
  959.         IX = IX + INCX
  960. # 492 "linpackd.f"
  961.         IY = IY + INCY
  962. # 493 "linpackd.f"
  963.    10 continue
  964.       ddot = dtemp
  965.       return
  966. # 480 "linpackd.f"
  967.       END IF
  968. c
  969. c        code for both increments equal to 1
  970. c
  971. # 500 "linpackd.f"
  972.       II4 = MOD (N, 4)
  973. # 500 "linpackd.f"
  974.       DO 3 I=1,II4
  975. # 501 "linpackd.f"
  976.         DTEMP = DTEMP + DX(I) * DY(I)
  977. # 502 "linpackd.f"
  978.     3 CONTINUE
  979. # 500 "linpackd.f"
  980.       DO 30 I=II4+1,N,4
  981. # 501 "linpackd.f"
  982.         dtemp = dtemp + dx(i)*dy(i)
  983.         DTEMP = DTEMP + DX(I+1) * DY(I+1)
  984. # 501 "linpackd.f"
  985.         DTEMP = DTEMP + DX(I+2) * DY(I+2)
  986. # 501 "linpackd.f"
  987.         DTEMP = DTEMP + DX(I+3) * DY(I+3)
  988. # 502 "linpackd.f"
  989.    30 continue
  990.       ddot = dtemp
  991.       return
  992.       end
  993. # 506 "linpackd.f"
  994.       subroutine  dscal(n,da,dx,incx)
  995. c
  996. c     scales a vector by a constant.
  997. c     jack dongarra, linpack, 3/11/78.
  998. c
  999.       DOUBLEPRECISION DA, DX(1)
  1000. # 512 "linpackd.f"
  1001.       integer i,incx,m,mp1,n,nincx
  1002. # 513 "linpackd.f"
  1003.       INTEGER II1, II2
  1004. # 513 "linpackd.f"
  1005. c
  1006.       IF (N .LE. 0) RETURN
  1007. # 515 "linpackd.f"
  1008.       IF (INCX .NE. 1) THEN
  1009. c
  1010. c        code for increment not equal to 1
  1011. c
  1012. # 520 "linpackd.f"
  1013.       II1 = N * INCX
  1014. # 520 "linpackd.f"
  1015. C$DOACROSS IF(II1 / INCX .GT. 125),SHARE(II1,INCX,DX,DA),LOCAL(I)
  1016. # 520 "linpackd.f"
  1017.       DO 2 I=1,II1,INCX
  1018. # 521 "linpackd.f"
  1019.         dx(i) = da*dx(i)
  1020. # 522 "linpackd.f"
  1021.     2 CONTINUE
  1022. # 523 "linpackd.f"
  1023.       return
  1024. # 515 "linpackd.f"
  1025.       END IF
  1026. c
  1027. c        code for increment equal to 1
  1028. c
  1029. # 528 "linpackd.f"
  1030.       II2 = MOD (N, 4)
  1031. # 528 "linpackd.f"
  1032.       DO 3 I=1,II2
  1033. # 529 "linpackd.f"
  1034.         DX(I) = DA * DX(I)
  1035. # 530 "linpackd.f"
  1036.     3   CONTINUE
  1037. # 528 "linpackd.f"
  1038. C$DOACROSS IF(N .GT. 125),SHARE(II2,N,DX,DA),LOCAL(I)
  1039. # 528 "linpackd.f"
  1040.       DO 4 I=II2+1,N,4
  1041. # 529 "linpackd.f"
  1042.         dx(i) = da*dx(i)
  1043.         DX(I+1) = DA * DX(I+1)
  1044. # 529 "linpackd.f"
  1045.         DX(I+2) = DA * DX(I+2)
  1046. # 529 "linpackd.f"
  1047.         DX(I+3) = DA * DX(I+3)
  1048. # 530 "linpackd.f"
  1049.     4 CONTINUE
  1050. # 531 "linpackd.f"
  1051.       return
  1052.       end
  1053. # 533 "linpackd.f"
  1054.       integer function idamax(n,dx,incx)
  1055. c
  1056. c     finds the index of element having max. dabsolute value.
  1057. c     jack dongarra, linpack, 3/11/78.
  1058. c
  1059.       DOUBLEPRECISION DX(1), DMAX
  1060. # 539 "linpackd.f"
  1061.       integer i,incx,ix,n
  1062. # 540 "linpackd.f"
  1063.       INTEGER II1, II2, II3
  1064. # 540 "linpackd.f"
  1065. c
  1066.       IDAMAX = 0
  1067. # 542 "linpackd.f"
  1068.       IF (N .LT. 1) RETURN
  1069. # 543 "linpackd.f"
  1070.       IDAMAX = 1
  1071. # 544 "linpackd.f"
  1072.       IF (N .EQ. 1) RETURN
  1073. # 545 "linpackd.f"
  1074.       IF (INCX .NE. 1) THEN
  1075. c
  1076. c        code for increment not equal to 1
  1077. c
  1078. # 550 "linpackd.f"
  1079.       DMAX = DABS (DX(1))
  1080. # 551 "linpackd.f"
  1081.       IX = INCX + 1
  1082. # 552 "linpackd.f"
  1083.       II2 = MOD (N - 1, 4)
  1084. # 552 "linpackd.f"
  1085.       DO 2 I=2,II2+1
  1086. # 553 "linpackd.f"
  1087.          IF (DABS (DX(IX)) .GT. DMAX) THEN
  1088. # 554 "linpackd.f"
  1089.          IDAMAX = I
  1090. # 555 "linpackd.f"
  1091.          DMAX = DABS (DX(IX))
  1092. # 553 "linpackd.f"
  1093.          END IF
  1094. # 556 "linpackd.f"
  1095.          IX = IX + INCX
  1096. # 557 "linpackd.f"
  1097.     2 CONTINUE
  1098. # 552 "linpackd.f"
  1099.       DO 10 I=II2+2,N,4
  1100. # 553 "linpackd.f"
  1101.          IF (DABS (DX(IX)) .GT. DMAX) THEN
  1102. # 554 "linpackd.f"
  1103.          idamax = i
  1104.          DMAX = DABS (DX(IX))
  1105. # 553 "linpackd.f"
  1106.          END IF
  1107. # 556 "linpackd.f"
  1108.          IX = IX + INCX
  1109. # 553 "linpackd.f"
  1110.          IF (DABS (DX(IX)) .GT. DMAX) THEN
  1111. # 554 "linpackd.f"
  1112.          IDAMAX = I + 1
  1113. # 555 "linpackd.f"
  1114.          DMAX = DABS (DX(IX))
  1115. # 553 "linpackd.f"
  1116.          END IF
  1117. # 556 "linpackd.f"
  1118.          IX = IX + INCX
  1119. # 553 "linpackd.f"
  1120.          IF (DABS (DX(IX)) .GT. DMAX) THEN
  1121. # 554 "linpackd.f"
  1122.          IDAMAX = I + 2
  1123. # 555 "linpackd.f"
  1124.          DMAX = DABS (DX(IX))
  1125. # 553 "linpackd.f"
  1126.          END IF
  1127. # 556 "linpackd.f"
  1128.          IX = IX + INCX
  1129. # 553 "linpackd.f"
  1130.          IF (DABS (DX(IX)) .GT. DMAX) THEN
  1131. # 554 "linpackd.f"
  1132.          IDAMAX = I + 3
  1133. # 555 "linpackd.f"
  1134.          DMAX = DABS (DX(IX))
  1135. # 553 "linpackd.f"
  1136.          END IF
  1137. # 556 "linpackd.f"
  1138.          IX = IX + INCX
  1139. # 557 "linpackd.f"
  1140.    10 continue
  1141.       return
  1142. # 545 "linpackd.f"
  1143.       END IF
  1144. c
  1145. c        code for increment equal to 1
  1146. c
  1147. # 562 "linpackd.f"
  1148.       DMAX = DABS (DX(1))
  1149. # 563 "linpackd.f"
  1150.       II3 = MOD (N - 1, 4)
  1151. # 563 "linpackd.f"
  1152.       DO 3 I=2,II3+1
  1153. # 564 "linpackd.f"
  1154.          IF (DABS (DX(I)) .GT. DMAX) THEN
  1155. # 565 "linpackd.f"
  1156.          IDAMAX = I
  1157. # 566 "linpackd.f"
  1158.          DMAX = DABS (DX(I))
  1159. # 564 "linpackd.f"
  1160.          END IF
  1161. # 567 "linpackd.f"
  1162.     3 CONTINUE
  1163. # 563 "linpackd.f"
  1164.       DO 30 I=II3+2,N,4
  1165. # 564 "linpackd.f"
  1166.          IF (DABS (DX(I)) .GT. DMAX) THEN
  1167. # 565 "linpackd.f"
  1168.          idamax = i
  1169.          dmax = dabs(dx(i))
  1170. # 564 "linpackd.f"
  1171.          END IF
  1172. # 564 "linpackd.f"
  1173.          IF (DABS (DX(I+1)) .GT. DMAX) THEN
  1174. # 565 "linpackd.f"
  1175.          IDAMAX = I + 1
  1176. # 566 "linpackd.f"
  1177.          DMAX = DABS (DX(I+1))
  1178. # 564 "linpackd.f"
  1179.          END IF
  1180. # 564 "linpackd.f"
  1181.          IF (DABS (DX(I+2)) .GT. DMAX) THEN
  1182. # 565 "linpackd.f"
  1183.          IDAMAX = I + 2
  1184. # 566 "linpackd.f"
  1185.          DMAX = DABS (DX(I+2))
  1186. # 564 "linpackd.f"
  1187.          END IF
  1188. # 564 "linpackd.f"
  1189.          IF (DABS (DX(I+3)) .GT. DMAX) THEN
  1190. # 565 "linpackd.f"
  1191.          IDAMAX = I + 3
  1192. # 566 "linpackd.f"
  1193.          DMAX = DABS (DX(I+3))
  1194. # 564 "linpackd.f"
  1195.          END IF
  1196. # 567 "linpackd.f"
  1197.    30 continue
  1198.       return
  1199.       end
  1200. # 570 "linpackd.f"
  1201.       double precision function epslon (x)
  1202.       double precision x
  1203. c
  1204. c     estimate unit roundoff in quantities of size x.
  1205. c
  1206.       double precision a,b,c,eps
  1207. # 576 "linpackd.f"
  1208.       INTEGER II1
  1209. # 576 "linpackd.f"
  1210.       DOUBLE PRECISION DD1
  1211. # 576 "linpackd.f"
  1212. # 576 "linpackd.f"
  1213. c
  1214. c     this program should function properly on all systems
  1215. c     satisfying the following two assumptions,
  1216. c        1.  the base used in representing dfloating point
  1217. c            numbers is not a power of three.
  1218. c        2.  the quantity  a  in statement 10 is represented to
  1219. c            the accuracy used in dfloating point variables
  1220. c            that are stored in memory.
  1221. c     the statement number 10 and the go to 10 are intended to
  1222. c     force optimizing compilers to generate code satisfying
  1223. c     assumption 2.
  1224. c     under these assumptions, it should be true that,
  1225. c            a  is not exactly equal to four-thirds,
  1226. c            b  has a zero for its last bit or digit,
  1227. c            c  is not exactly equal to one,
  1228. c            eps  measures the separation of 1.0 from
  1229. c                 the next larger dfloating point number.
  1230. c     the developers of eispack would appreciate being informed
  1231. c     about any systems where these assumptions do not hold.
  1232. c
  1233. c     *****************************************************************
  1234. c     this routine is one of the auxiliary routines used by eispack iii
  1235. c     to avoid machine dependencies.
  1236. c     *****************************************************************
  1237. c
  1238. c     this version dated 4/6/83.
  1239. c
  1240.       A = 4.D0 / 3.D0
  1241. # 604 "linpackd.f"
  1242.       II1 = 1
  1243. # 604 "linpackd.f"
  1244.       DD1 = A - 1.D0
  1245. # 604 "linpackd.f"
  1246.       DO   WHILE ( II1 .EQ. 1 )
  1247. # 604 "linpackd.f"
  1248.       B = DD1
  1249. # 605 "linpackd.f"
  1250.       c = b + b + b
  1251.       EPS = DABS (C - 1.D0)
  1252. # 607 "linpackd.f"
  1253.       IF (EPS .NE. 0.D0) II1 = 0
  1254. # 607 "linpackd.f"
  1255.       END DO
  1256. # 608 "linpackd.f"
  1257.       epslon = eps*dabs(x)
  1258.       return
  1259.       end
  1260. # 611 "linpackd.f"
  1261.       subroutine dmxpy (n1, y, n2, ldm, x, m)
  1262.       double precision y(*), x(*), m(ldm,*)
  1263. # 613 "linpackd.f"
  1264.       INTEGER II1, II2, II3, II4, II5, II6, II7, II8
  1265. # 613 "linpackd.f"
  1266.       DOUBLE PRECISION DD1, DD2, DD3, DD4, DD5, DD6, DD7, DD8, DD9, DD10
  1267.      X  , DD11, DD12, DD13, DD14, DD15, DD16
  1268. # 613 "linpackd.f"
  1269. # 613 "linpackd.f"
  1270. c
  1271. c   purpose:
  1272. c     multiply matrix m times vector x and add the result to vector y.
  1273. c
  1274. c   parameters:
  1275. c
  1276. c     n1 integer, number of elements in vector y, and number of rows in
  1277. c         matrix m
  1278. c
  1279. c     y double precision(n1), vector of length n1 to which is added
  1280. c         the product m*x
  1281. c
  1282. c     n2 integer, number of elements in vector x, and number of columns
  1283. c         in matrix m
  1284. c
  1285. c     ldm integer, leading dimension of array m
  1286. c
  1287. c     x double precision(n2), vector of length n2
  1288. c
  1289. c     m double precision(ldm,n2), matrix of n1 rows and n2 columns
  1290. c
  1291. c ----------------------------------------------------------------------
  1292. c
  1293. c   cleanup odd vector
  1294. c
  1295.       IF (MOD (N2, 2) .GE. 1) THEN
  1296. # 640 "linpackd.f"
  1297.       II1 = MOD (N2, 2)
  1298. # 640 "linpackd.f"
  1299.       DD1 = X(II1)
  1300. # 640 "linpackd.f"
  1301.       II5 = MOD (N1, 4)
  1302. # 640 "linpackd.f"
  1303.       DO 2 I=1,II5
  1304. # 641 "linpackd.f"
  1305.             Y(I) = Y(I) + DD1 * M(I,II1)
  1306. # 642 "linpackd.f"
  1307.     2       CONTINUE
  1308. # 640 "linpackd.f"
  1309. C$DOACROSS IF(N1 .GT. 41),SHARE(II5,N1,Y,DD1,II1,M),LOCAL(I)
  1310. # 640 "linpackd.f"
  1311.          DO 3 I=II5+1,N1,4
  1312. # 641 "linpackd.f"
  1313.             Y(I) = Y(I) + DD1 * M(I,II1)
  1314. # 641 "linpackd.f"
  1315.             Y(I+1) = Y(I+1) + DD1 * M(I+1,II1)
  1316. # 641 "linpackd.f"
  1317.             Y(I+2) = Y(I+2) + DD1 * M(I+2,II1)
  1318. # 641 "linpackd.f"
  1319.             Y(I+3) = Y(I+3) + DD1 * M(I+3,II1)
  1320. # 642 "linpackd.f"
  1321.     3    CONTINUE
  1322. # 643 "linpackd.f"
  1323.       endif
  1324. c
  1325. c   cleanup odd group of two vectors
  1326. c
  1327. # 648 "linpackd.f"
  1328.       IF (MOD (N2, 4) .GE. 2) THEN
  1329. # 649 "linpackd.f"
  1330.       II2 = MOD (N2, 4)
  1331. # 649 "linpackd.f"
  1332.       DD2 = X(II2)
  1333. # 649 "linpackd.f"
  1334.       DD3 = X(II2-1)
  1335. # 649 "linpackd.f"
  1336.       II6 = MOD (N1, 4)
  1337. # 649 "linpackd.f"
  1338.       DO 4 I=1,II6
  1339. # 650 "linpackd.f"
  1340.             Y(I) = (Y(I) + DD3 * M(I,II2-1)) + DD2 * M(I,II2)
  1341. # 652 "linpackd.f"
  1342.     4       CONTINUE
  1343. # 649 "linpackd.f"
  1344. C$DOACROSS IF(N1 .GT. 21),SHARE(II6,N1,Y,DD3,II2,M,DD2),LOCAL(I)
  1345. # 649 "linpackd.f"
  1346.          DO 5 I=II6+1,N1,4
  1347. # 650 "linpackd.f"
  1348.             Y(I) = (Y(I) + DD3 * M(I,II2-1)) + DD2 * M(I,II2)
  1349. # 650 "linpackd.f"
  1350.             Y(I+1) = (Y(I+1) + DD3 * M(I+1,II2-1)) + DD2 * M(I+1,II2)
  1351. # 650 "linpackd.f"
  1352.             Y(I+2) = (Y(I+2) + DD3 * M(I+2,II2-1)) + DD2 * M(I+2,II2)
  1353. # 650 "linpackd.f"
  1354.             Y(I+3) = (Y(I+3) + DD3 * M(I+3,II2-1)) + DD2 * M(I+3,II2)
  1355. # 652 "linpackd.f"
  1356.     5    CONTINUE
  1357. # 653 "linpackd.f"
  1358.       endif
  1359. c
  1360. c   cleanup odd group of four vectors
  1361. c
  1362. # 658 "linpackd.f"
  1363.       IF (MOD (N2, 8) .GE. 4) THEN
  1364. # 659 "linpackd.f"
  1365.       II3 = MOD (N2, 8)
  1366. # 659 "linpackd.f"
  1367.       DD4 = X(II3)
  1368. # 659 "linpackd.f"
  1369.       DD5 = X(II3-1)
  1370. # 659 "linpackd.f"
  1371.       DD6 = X(II3-2)
  1372. # 659 "linpackd.f"
  1373.       DD7 = X(II3-3)
  1374. # 659 "linpackd.f"
  1375.       II7 = MOD (N1, 4)
  1376. # 659 "linpackd.f"
  1377.       DO 6 I=1,II7
  1378. # 660 "linpackd.f"
  1379.             Y(I) = (((Y(I) + DD7 * M(I,II3-3)) + DD6 * M(I,II3-2)) + DD5
  1380.      X              * M(I,II3-1)) + DD4 * M(I,II3)
  1381. # 663 "linpackd.f"
  1382.     6       CONTINUE
  1383. # 659 "linpackd.f"
  1384. C$DOACROSS IF(N1 .GT. 11),SHARE(II7,N1,Y,DD7,II3,M,DD6,DD5,DD4),LOCAL(I)
  1385. # 659 "linpackd.f"
  1386.          DO 7 I=II7+1,N1,4
  1387. # 660 "linpackd.f"
  1388.             Y(I) = (((Y(I) + DD7 * M(I,II3-3)) + DD6 * M(I,II3-2)) + DD5
  1389.      X              * M(I,II3-1)) + DD4 * M(I,II3)
  1390. # 660 "linpackd.f"
  1391.             Y(I+1) = (((Y(I+1) + DD7 * M(I+1,II3-3)) + DD6 * M(I+1,II3-2
  1392.      X             )) + DD5 * M(I+1,II3-1)) + DD4 * M(I+1,II3)
  1393. # 660 "linpackd.f"
  1394.             Y(I+2) = (((Y(I+2) + DD7 * M(I+2,II3-3)) + DD6 * M(I+2,II3-2
  1395.      X             )) + DD5 * M(I+2,II3-1)) + DD4 * M(I+2,II3)
  1396. # 660 "linpackd.f"
  1397.             Y(I+3) = (((Y(I+3) + DD7 * M(I+3,II3-3)) + DD6 * M(I+3,II3-2
  1398.      X             )) + DD5 * M(I+3,II3-1)) + DD4 * M(I+3,II3)
  1399. # 663 "linpackd.f"
  1400.     7    CONTINUE
  1401. # 664 "linpackd.f"
  1402.       endif
  1403. c
  1404. c   cleanup odd group of eight vectors
  1405. c
  1406. # 669 "linpackd.f"
  1407.       IF (MOD (N2, 16) .GE. 8) THEN
  1408. # 670 "linpackd.f"
  1409.       II4 = MOD (N2, 16)
  1410. # 670 "linpackd.f"
  1411.       DD8 = X(II4)
  1412. # 670 "linpackd.f"
  1413.       DD9 = X(II4-1)
  1414. # 670 "linpackd.f"
  1415.       DD10 = X(II4-2)
  1416. # 670 "linpackd.f"
  1417.       DD11 = X(II4-3)
  1418. # 670 "linpackd.f"
  1419.       DD12 = X(II4-4)
  1420. # 670 "linpackd.f"
  1421.       DD13 = X(II4-5)
  1422. # 670 "linpackd.f"
  1423.       DD14 = X(II4-6)
  1424. # 670 "linpackd.f"
  1425.       DD15 = X(II4-7)
  1426. # 670 "linpackd.f"
  1427.       II8 = MOD (N1, 2)
  1428. # 670 "linpackd.f"
  1429.       DO 8 I=1,II8
  1430. # 671 "linpackd.f"
  1431.             Y(I) = (((((((Y(I) + DD15 * M(I,II4-7)) + DD14 * M(I,II4-6))
  1432.      X              + DD13 * M(I,II4-5)) + DD12 * M(I,II4-4)) + DD11 * M
  1433.      X             (I,II4-3)) + DD10 * M(I,II4-2)) + DD9 * M(I,II4-1)) +
  1434.      X              DD8 * M(I,II4)
  1435. # 676 "linpackd.f"
  1436.     8       CONTINUE
  1437. # 670 "linpackd.f"
  1438. C$DOACROSS SHARE(II8,N1,Y,DD15,II4,M,DD14,DD13,DD12,DD11,DD10,DD9,DD8),
  1439. C$&  LOCAL(I)
  1440. # 670 "linpackd.f"
  1441.          DO 9 I=II8+1,N1,2
  1442. # 671 "linpackd.f"
  1443.             Y(I) = (((((((Y(I) + DD15 * M(I,II4-7)) + DD14 * M(I,II4-6))
  1444.      X              + DD13 * M(I,II4-5)) + DD12 * M(I,II4-4)) + DD11 * M
  1445.      X             (I,II4-3)) + DD10 * M(I,II4-2)) + DD9 * M(I,II4-1)) +
  1446.      X              DD8 * M(I,II4)
  1447. # 671 "linpackd.f"
  1448.             Y(I+1) = (((((((Y(I+1) + DD15 * M(I+1,II4-7)) + DD14 * M(I+1
  1449.      X             ,II4-6)) + DD13 * M(I+1,II4-5)) + DD12 * M(I+1,II4-4)
  1450.      X             ) + DD11 * M(I+1,II4-3)) + DD10 * M(I+1,II4-2)) + DD9
  1451.      X              * M(I+1,II4-1)) + DD8 * M(I+1,II4)
  1452. # 676 "linpackd.f"
  1453.     9    CONTINUE
  1454. # 677 "linpackd.f"
  1455.       endif
  1456. c
  1457. c   main loop - groups of sixteen vectors
  1458. c
  1459.       JMIN = MOD (N2, 16) + 16
  1460. # 683 "linpackd.f"
  1461. C$DOACROSS SHARE(N1,Y,JMIN,N2,X,M),LOCAL(DD16,I,J)
  1462. # 683 "linpackd.f"
  1463.          DO 11 I=1,N1
  1464. # 683 "linpackd.f"
  1465.          DD16 = Y(I)
  1466. # 682 "linpackd.f"
  1467.       DO 60 J=JMIN,N2,16
  1468. # 684 "linpackd.f"
  1469.             DD16 = (((((((((((((((DD16 + X(J-15) * M(I,J-15)) + X(J-14)
  1470.      X             * M(I,J-14)) + X(J-13) * M(I,J-13)) + X(J-12) * M(I,J
  1471.      X             -12)) + X(J-11) * M(I,J-11)) + X(J-10) * M(I,J-10)) +
  1472.      X              X(J-9) * M(I,J-9)) + X(J-8) * M(I,J-8)) + X(J-7) * M
  1473.      X             (I,J-7)) + X(J-6) * M(I,J-6)) + X(J-5) * M(I,J-5)) +
  1474.      X             X(J-4) * M(I,J-4)) + X(J-3) * M(I,J-3)) + X(J-2) * M(
  1475.      X             I,J-2)) + X(J-1) * M(I,J-1)) + X(J) * M(I,J)
  1476. # 694 "linpackd.f"
  1477.    60 continue
  1478.       Y(I) = DD16
  1479. # 693 "linpackd.f"
  1480.    11    CONTINUE
  1481. # 695 "linpackd.f"
  1482.       return
  1483.       end
  1484.  
  1485.  
  1486.