home *** CD-ROM | disk | FTP | other *** search
/ Power Programming / powerprogramming1994.iso / progtool / progjrn / pj_6_1.arc / WHETLIN.ARC / SLINPAKR.FOR < prev    next >
Text File  |  1987-09-20  |  22KB  |  785 lines
  1. *     a TIME function for Ryan/McFarland Fortran and Microsoft Version 4.0
  2.  
  3. *    Author:    M. Steven Baker
  4. *    Date:    September 20, 1986
  5. *
  6.        real function second()
  7.        integer*2 hh,mm,ss,hd
  8.        call gettim(hh,mm,ss,hd)
  9.        second = float(hh)*3600 + float(mm*60+ss) + float(hd)/100
  10.        end
  11.  
  12. *      real function second()
  13.       
  14. *      external msec
  15. *      second = msec()*0.001
  16. *      end
  17. *$system
  18.  
  19.       real aa(200,200),a(201,200),b(200),x(200)
  20.       real time(8,6),cray,ops,total,norma,normx
  21.       real resid,residn,eps,epslon
  22.       integer ipvt(200)
  23.       lda = 201
  24.       ldaa = 200
  25. c
  26.       n = 100
  27.       cray = .056
  28.       write(*,*)
  29.     1 format(' Please send the results of this run to:'//
  30.      $       ' Jack J. Dongarra'/
  31.      $       ' Mathematics and Computer Science Division'/
  32.      $       ' Argonne National Laboratory'/
  33.      $       ' Argonne, Illinois 60439'//
  34.      $       ' Telephone: 312-972-7246'//
  35.      $       ' ARPAnet: DONGARRA@ANL-MCS'/)
  36.       ops = (2.0e0*n**3)/3.0e0 + 2.0e0*n**2
  37. c
  38.          call matgen(a,lda,n,b,norma)
  39.          t1 = second()
  40.          call sgefa(a,lda,n,ipvt,info)
  41.          time(1,1) = second() - t1
  42.          t1 = second()
  43.          call sgesl(a,lda,n,ipvt,b,0)
  44.          time(1,2) = second() - t1
  45.          total = time(1,1) + time(1,2)
  46. C
  47. C     COMPUTE A RESIDUAL TO VERIFY RESULTS.
  48. C
  49.          do 10 i = 1,n
  50.             x(i) = b(i)
  51.    10    continue
  52.          call matgen(a,lda,n,b,norma)
  53.          do 20 i = 1,n
  54.             b(i) = -b(i)
  55.    20    continue
  56.          CALL SMXPY(n,b,n,lda,x,a)
  57.          RESID = 0.0
  58.          NORMX = 0.0
  59.          DO 30 I = 1,N
  60.             RESID = amax1( RESID, ABS(b(i)) )
  61.             NORMX = amax1( NORMX, ABS(X(I)) )
  62.    30    CONTINUE
  63.          eps = epslon(1.0)
  64.          RESIDn = RESID/( N*NORMA*NORMX*EPS )
  65.          write(*,40)
  66.    40    format('     norm. resid      resid           machep',
  67.      $          '         x(1)          x(n)')
  68.          write(*,50) residn,resid,eps,x(1),x(n)
  69.    50    format(1p5e16.8)
  70. c
  71.          write(*,60) n
  72.    60    format(//'    times are reported for matrices of order ',i5)
  73.          write(*,70)
  74.    70    format(6x,'sgefa',6x,'sgesl',6x,'total',5x,'mflops',7x,'unit',
  75.      $         6x,'ratio')
  76. c
  77.          time(1,3) = total
  78.          time(1,4) = ops/(1.0e6*total)
  79.          time(1,5) = 2.0e0/time(1,4)
  80.          time(1,6) = total/cray
  81.          write(*,80) lda
  82.    80    format(' times for array with leading dimension of',i4)
  83.          write(*,110) (time(1,i),i=1,6)
  84. c
  85.          call matgen(a,lda,n,b,norma)
  86.          t1 = second()
  87.          call sgefa(a,lda,n,ipvt,info)
  88.          time(2,1) = second() - t1
  89.          t1 = second()
  90.          call sgesl(a,lda,n,ipvt,b,0)
  91.          time(2,2) = second() - t1
  92.          total = time(2,1) + time(2,2)
  93.          time(2,3) = total
  94.          time(2,4) = ops/(1.0e6*total)
  95.          time(2,5) = 2.0e0/time(2,4)
  96.          time(2,6) = total/cray
  97. c
  98.          call matgen(a,lda,n,b,norma)
  99.          t1 = second()
  100.          call sgefa(a,lda,n,ipvt,info)
  101.          time(3,1) = second() - t1
  102.          t1 = second()
  103.          call sgesl(a,lda,n,ipvt,b,0)
  104.          time(3,2) = second() - t1
  105.          total = time(3,1) + time(3,2)
  106.          time(3,3) = total
  107.          time(3,4) = ops/(1.0e6*total)
  108.          time(3,5) = 2.0e0/time(3,4)
  109.          time(3,6) = total/cray
  110. c
  111.          ntimes = 10
  112.          tm2 = 0
  113.          t1 = second()
  114.          do 90 i = 1,ntimes
  115.             tm = second()
  116.             call matgen(a,lda,n,b,norma)
  117.             tm2 = tm2 + second() - tm
  118.             call sgefa(a,lda,n,ipvt,info)
  119.    90    continue
  120.          time(4,1) = (second() - t1 - tm2)/ntimes
  121.          t1 = second()
  122.          do 100 i = 1,ntimes
  123.             call sgesl(a,lda,n,ipvt,b,0)
  124.   100    continue
  125.          time(4,2) = (second() - t1)/ntimes
  126.          total = time(4,1) + time(4,2)
  127.          time(4,3) = total
  128.          time(4,4) = ops/(1.0e6*total)
  129.          time(4,5) = 2.0e0/time(4,4)
  130.          time(4,6) = total/cray
  131. c
  132.          write(*,110) (time(2,i),i=1,6)
  133.          write(*,110) (time(3,i),i=1,6)
  134.          write(*,110) (time(4,i),i=1,6)
  135.   110    format(6(1pe11.3))
  136. c
  137.          call matgen(aa,ldaa,n,b,norma)
  138.          t1 = second()
  139.          call sgefa(aa,ldaa,n,ipvt,info)
  140.          time(5,1) = second() - t1
  141.          t1 = second()
  142.          call sgesl(aa,ldaa,n,ipvt,b,0)
  143.          time(5,2) = second() - t1
  144.          total = time(5,1) + time(5,2)
  145.          time(5,3) = total
  146.          time(5,4) = ops/(1.0e6*total)
  147.          time(5,5) = 2.0e0/time(5,4)
  148.          time(5,6) = total/cray
  149. c
  150.          call matgen(aa,ldaa,n,b,norma)
  151.          t1 = second()
  152.          call sgefa(aa,ldaa,n,ipvt,info)
  153.          time(6,1) = second() - t1
  154.          t1 = second()
  155.          call sgesl(aa,ldaa,n,ipvt,b,0)
  156.          time(6,2) = second() - t1
  157.          total = time(6,1) + time(6,2)
  158.          time(6,3) = total
  159.          time(6,4) = ops/(1.0e6*total)
  160.          time(6,5) = 2.0e0/time(6,4)
  161.          time(6,6) = total/cray
  162. c
  163.          call matgen(aa,ldaa,n,b,norma)
  164.          t1 = second()
  165.          call sgefa(aa,ldaa,n,ipvt,info)
  166.          time(7,1) = second() - t1
  167.          t1 = second()
  168.          call sgesl(aa,ldaa,n,ipvt,b,0)
  169.          time(7,2) = second() - t1
  170.          total = time(7,1) + time(7,2)
  171.          time(7,3) = total
  172.          time(7,4) = ops/(1.0e6*total)
  173.          time(7,5) = 2.0e0/time(7,4)
  174.          time(7,6) = total/cray
  175. c
  176.          ntimes = 10
  177.          tm2 = 0
  178.          t1 = second()
  179.          do 120 i = 1,ntimes
  180.             tm = second()
  181.             call matgen(aa,ldaa,n,b,norma)
  182.             tm2 = tm2 + second() - tm
  183.             call sgefa(aa,ldaa,n,ipvt,info)
  184.   120    continue
  185.          time(8,1) = (second() - t1 - tm2)/ntimes
  186.          t1 = second()
  187.          do 130 i = 1,ntimes
  188.             call sgesl(aa,ldaa,n,ipvt,b,0)
  189.   130    continue
  190.          time(8,2) = (second() - t1)/ntimes
  191.          total = time(8,1) + time(8,2)
  192.          time(8,3) = total
  193.          time(8,4) = ops/(1.0e6*total)
  194.          time(8,5) = 2.0e0/time(8,4)
  195.          time(8,6) = total/cray
  196. c
  197.          write(*,140) ldaa
  198.   140    format(/' times for array with leading dimension of',i4)
  199.          write(*,110) (time(5,i),i=1,6)
  200.          write(*,110) (time(6,i),i=1,6)
  201.          write(*,110) (time(7,i),i=1,6)
  202.          write(*,110) (time(8,i),i=1,6)
  203.       stop
  204.       end
  205.       subroutine matgen(a,lda,n,b,norma)
  206.       real a(lda,1),b(1),norma
  207. c
  208.       init = 1325
  209.       norma = 0.0
  210.       do 30 j = 1,n
  211.          do 20 i = 1,n
  212.             init = mod(3125*init,65536)
  213.             a(i,j) = (init - 32768.0)/16384.0
  214.             norma = amax1(a(i,j), norma)
  215.    20    continue
  216.    30 continue
  217.       do 35 i = 1,n
  218.           b(i) = 0.0
  219.    35 continue
  220.       do 50 j = 1,n
  221.          do 40 i = 1,n
  222.             b(i) = b(i) + a(i,j)
  223.    40    continue
  224.    50 continue
  225.       return
  226.       end
  227.       subroutine sgefa(a,lda,n,ipvt,info)
  228.       integer lda,n,ipvt(1),info
  229.       real a(lda,1)
  230. c
  231. c     sgefa factors a real matrix by gaussian elimination.
  232. c
  233. c     sgefa is usually called by dgeco, but it can be called
  234. c     directly with a saving in time if  rcond  is not needed.
  235. c     (time for dgeco) = (1 + 9/n)*(time for sgefa) .
  236. c
  237. c     on entry
  238. c
  239. c        a       real(lda, n)
  240. c                the matrix to be factored.
  241. c
  242. c        lda     integer
  243. c                the leading dimension of the array  a .
  244. c
  245. c        n       integer
  246. c                the order of the matrix  a .
  247. c
  248. c     on return
  249. c
  250. c        a       an upper triangular matrix and the multipliers
  251. c                which were used to obtain it.
  252. c                the factorization can be written  a = l*u  where
  253. c                l  is a product of permutation and unit lower
  254. c                triangular matrices and  u  is upper triangular.
  255. c
  256. c        ipvt    integer(n)
  257. c                an integer vector of pivot indices.
  258. c
  259. c        info    integer
  260. c                = 0  normal value.
  261. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  262. c                     condition for this subroutine, but it does
  263. c                     indicate that sgesl or dgedi will divide by zero
  264. c                     if called.  use  rcond  in dgeco for a reliable
  265. c                     indication of singularity.
  266. c
  267. c     linpack. this version dated 08/14/78 .
  268. c     cleve moler, university of new mexico, argonne national lab.
  269. c
  270. c     subroutines and functions
  271. c
  272. c     blas saxpy,sscal,isamax
  273. c
  274. c     internal variables
  275. c
  276.       real t
  277.       integer isamax,j,k,kp1,l,nm1
  278. c
  279. c
  280. c     gaussian elimination with partial pivoting
  281. c
  282.       info = 0
  283.       nm1 = n - 1
  284.       if (nm1 .lt. 1) go to 70
  285.       do 60 k = 1, nm1
  286.          kp1 = k + 1
  287. c
  288. c        find l = pivot index
  289. c
  290.          l = isamax(n-k+1,a(k,k),1) + k - 1
  291.          ipvt(k) = l
  292. c
  293. c        zero pivot implies this column already triangularized
  294. c
  295.          if (a(l,k) .eq. 0.0e0) go to 40
  296. c
  297. c           interchange if necessary
  298. c
  299.             if (l .eq. k) go to 10
  300.                t = a(l,k)
  301.                a(l,k) = a(k,k)
  302.                a(k,k) = t
  303.    10       continue
  304. c
  305. c           compute multipliers
  306. c
  307.             t = -1.0e0/a(k,k)
  308.             call sscal(n-k,t,a(k+1,k),1)
  309. c
  310. c           row elimination with column indexing
  311. c
  312.             do 30 j = kp1, n
  313.                t = a(l,j)
  314.                if (l .eq. k) go to 20
  315.                   a(l,j) = a(k,j)
  316.                   a(k,j) = t
  317.    20          continue
  318.                call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  319.    30       continue
  320.          go to 50
  321.    40    continue
  322.             info = k
  323.    50    continue
  324.    60 continue
  325.    70 continue
  326.       ipvt(n) = n
  327.       if (a(n,n) .eq. 0.0e0) info = n
  328.       return
  329.       end
  330.       subroutine sgesl(a,lda,n,ipvt,b,job)
  331.       integer lda,n,ipvt(1),job
  332.       real a(lda,1),b(1)
  333. c
  334. c     sgesl solves the real system
  335. c     a * x = b  or  trans(a) * x = b
  336. c     using the factors computed by dgeco or sgefa.
  337. c
  338. c     on entry
  339. c
  340. c        a       real(lda, n)
  341. c                the output from dgeco or sgefa.
  342. c
  343. c        lda     integer
  344. c                the leading dimension of the array  a .
  345. c
  346. c        n       integer
  347. c                the order of the matrix  a .
  348. c
  349. c        ipvt    integer(n)
  350. c                the pivot vector from dgeco or sgefa.
  351. c
  352. c        b       real(n)
  353. c                the right hand side vector.
  354. c
  355. c        job     integer
  356. c                = 0         to solve  a*x = b ,
  357. c                = nonzero   to solve  trans(a)*x = b  where
  358. c                            trans(a)  is the transpose.
  359. c
  360. c     on return
  361. c
  362. c        b       the solution vector  x .
  363. c
  364. c     error condition
  365. c
  366. c        a division by zero will occur if the input factor contains a
  367. c        zero on the diagonal.  technically this indicates singularity
  368. c        but it is often caused by improper arguments or improper
  369. c        setting of lda .  it will not occur if the subroutines are
  370. c        called correctly and if dgeco has set rcond .gt. 0.0
  371. c        or sgefa has set info .eq. 0 .
  372. c
  373. c     to compute  inverse(a) * c  where  c  is a matrix
  374. c     with  p  columns
  375. c           call dgeco(a,lda,n,ipvt,rcond,z)
  376. c           if (rcond is too small) go to ...
  377. c           do 10 j = 1, p
  378. c              call sgesl(a,lda,n,ipvt,c(1,j),0)
  379. c        10 continue
  380. c
  381. c     linpack. this version dated 08/14/78 .
  382. c     cleve moler, university of new mexico, argonne national lab.
  383. c
  384. c     subroutines and functions
  385. c
  386. c     blas saxpy,sdot
  387. c
  388. c     internal variables
  389. c
  390.       real sdot,t
  391.       integer k,kb,l,nm1
  392. c
  393.       nm1 = n - 1
  394.       if (job .ne. 0) go to 50
  395. c
  396. c        job = 0 , solve  a * x = b
  397. c        first solve  l*y = b
  398. c
  399.          if (nm1 .lt. 1) go to 30
  400.          do 20 k = 1, nm1
  401.             l = ipvt(k)
  402.             t = b(l)
  403.             if (l .eq. k) go to 10
  404.                b(l) = b(k)
  405.                b(k) = t
  406.    10       continue
  407.             call saxpy(n-k,t,a(k+1,k),1,b(k+1),1)
  408.    20    continue
  409.    30    continue
  410. c
  411. c        now solve  u*x = y
  412. c
  413.          do 40 kb = 1, n
  414.             k = n + 1 - kb
  415.             b(k) = b(k)/a(k,k)
  416.             t = -b(k)
  417.             call saxpy(k-1,t,a(1,k),1,b(1),1)
  418.    40    continue
  419.       go to 100
  420.    50 continue
  421. c
  422. c        job = nonzero, solve  trans(a) * x = b
  423. c        first solve  trans(u)*y = b
  424. c
  425.          do 60 k = 1, n
  426.             t = sdot(k-1,a(1,k),1,b(1),1)
  427.             b(k) = (b(k) - t)/a(k,k)
  428.    60    continue
  429. c
  430. c        now solve trans(l)*x = y
  431. c
  432.          if (nm1 .lt. 1) go to 90
  433.          do 80 kb = 1, nm1
  434.             k = n - kb
  435.             b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1)
  436.             l = ipvt(k)
  437.             if (l .eq. k) go to 70
  438.                t = b(l)
  439.                b(l) = b(k)
  440.                b(k) = t
  441.    70       continue
  442.    80    continue
  443.    90    continue
  444.   100 continue
  445.       return
  446.       end
  447.       subroutine saxpy(n,da,dx,incx,dy,incy)
  448. c
  449. c     constant times a vector plus a vector.
  450. c     uses unrolled loops for increments equal to one.
  451. c     jack dongarra, linpack, 3/11/78.
  452. c
  453.       real dx(1),dy(1),da
  454.       integer i,incx,incy,ix,iy,m,mp1,n
  455. c
  456.       if(n.le.0)return
  457.       if (da .eq. 0.0e0) return
  458.       if(incx.eq.1.and.incy.eq.1)go to 20
  459. c
  460. c        code for unequal increments or equal increments
  461. c          not equal to 1
  462. c
  463.       ix = 1
  464.       iy = 1
  465.       if(incx.lt.0)ix = (-n+1)*incx + 1
  466.       if(incy.lt.0)iy = (-n+1)*incy + 1
  467.       do 10 i = 1,n
  468.         dy(iy) = dy(iy) + da*dx(ix)
  469.         ix = ix + incx
  470.         iy = iy + incy
  471.    10 continue
  472.       return
  473. c
  474. c        code for both increments equal to 1
  475. c
  476. c
  477. c        clean-up loop
  478. c
  479.    20 m = mod(n,4)
  480.       if( m .eq. 0 ) go to 40
  481.       do 30 i = 1,m
  482.         dy(i) = dy(i) + da*dx(i)
  483.    30 continue
  484.       if( n .lt. 4 ) return
  485.    40 mp1 = m + 1
  486.       do 50 i = mp1,n,4
  487.         dy(i) = dy(i) + da*dx(i)
  488.         dy(i + 1) = dy(i + 1) + da*dx(i + 1)
  489.         dy(i + 2) = dy(i + 2) + da*dx(i + 2)
  490.         dy(i + 3) = dy(i + 3) + da*dx(i + 3)
  491.    50 continue
  492.       return
  493.       end
  494.       real function sdot(n,dx,incx,dy,incy)
  495. c
  496. c     forms the dot product of two vectors.
  497. c     uses unrolled loops for increments equal to one.
  498. c     jack dongarra, linpack, 3/11/78.
  499. c
  500.       real dx(1),dy(1),dtemp
  501.       integer i,incx,incy,ix,iy,m,mp1,n
  502. c
  503.       sdot = 0.0e0
  504.       dtemp = 0.0e0
  505.       if(n.le.0)return
  506.       if(incx.eq.1.and.incy.eq.1)go to 20
  507. c
  508. c        code for unequal increments or equal increments
  509. c          not equal to 1
  510. c
  511.       ix = 1
  512.       iy = 1
  513.       if(incx.lt.0)ix = (-n+1)*incx + 1
  514.       if(incy.lt.0)iy = (-n+1)*incy + 1
  515.       do 10 i = 1,n
  516.         dtemp = dtemp + dx(ix)*dy(iy)
  517.         ix = ix + incx
  518.         iy = iy + incy
  519.    10 continue
  520.       sdot = dtemp
  521.       return
  522. c
  523. c        code for both increments equal to 1
  524. c
  525. c
  526. c        clean-up loop
  527. c
  528.    20 m = mod(n,5)
  529.       if( m .eq. 0 ) go to 40
  530.       do 30 i = 1,m
  531.         dtemp = dtemp + dx(i)*dy(i)
  532.    30 continue
  533.       if( n .lt. 5 ) go to 60
  534.    40 mp1 = m + 1
  535.       do 50 i = mp1,n,5
  536.         dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
  537.      *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
  538.    50 continue
  539.    60 sdot = dtemp
  540.       return
  541.       end
  542.       subroutine  sscal(n,da,dx,incx)
  543. c
  544. c     scales a vector by a constant.
  545. c     uses unrolled loops for increment equal to one.
  546. c     jack dongarra, linpack, 3/11/78.
  547. c
  548.       real da,dx(1)
  549.       integer i,incx,m,mp1,n,nincx
  550. c
  551.       if(n.le.0)return
  552.       if(incx.eq.1)go to 20
  553. c
  554. c        code for increment not equal to 1
  555. c
  556.       nincx = n*incx
  557.       do 10 i = 1,nincx,incx
  558.         dx(i) = da*dx(i)
  559.    10 continue
  560.       return
  561. c
  562. c        code for increment equal to 1
  563. c
  564. c
  565. c        clean-up loop
  566. c
  567.    20 m = mod(n,5)
  568.       if( m .eq. 0 ) go to 40
  569.       do 30 i = 1,m
  570.         dx(i) = da*dx(i)
  571.    30 continue
  572.       if( n .lt. 5 ) return
  573.    40 mp1 = m + 1
  574.       do 50 i = mp1,n,5
  575.         dx(i) = da*dx(i)
  576.         dx(i + 1) = da*dx(i + 1)
  577.         dx(i + 2) = da*dx(i + 2)
  578.         dx(i + 3) = da*dx(i + 3)
  579.         dx(i + 4) = da*dx(i + 4)
  580.    50 continue
  581.       return
  582.       end
  583.       integer function isamax(n,dx,incx)
  584. c
  585. c     finds the index of element having max. absolute value.
  586. c     jack dongarra, linpack, 3/11/78.
  587. c
  588.       real dx(1),dmax
  589.       integer i,incx,ix,n
  590. c
  591.       isamax = 0
  592.       if( n .lt. 1 ) return
  593.       isamax = 1
  594.       if(n.eq.1)return
  595.       if(incx.eq.1)go to 20
  596. c
  597. c        code for increment not equal to 1
  598. c
  599.       ix = 1
  600.       dmax = abs(dx(1))
  601.       ix = ix + incx
  602.       do 10 i = 2,n
  603.          if(abs(dx(ix)).le.dmax) go to 5
  604.          isamax = i
  605.          dmax = abs(dx(ix))
  606.     5    ix = ix + incx
  607.    10 continue
  608.       return
  609. c
  610. c        code for increment equal to 1
  611. c
  612.    20 dmax = abs(dx(1))
  613.       do 30 i = 2,n
  614.          if(abs(dx(i)).le.dmax) go to 30
  615.          isamax = i
  616.          dmax = abs(dx(i))
  617.    30 continue
  618.       return
  619.       end
  620.       REAL FUNCTION EPSLON (X)
  621.       REAL X
  622. C
  623. C     ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
  624. C
  625.       REAL A,B,C,EPS
  626. C
  627. C     THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS
  628. C     SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
  629. C        1.  THE BASE USED IN REPRESENTING FLOATING POINT
  630. C            NUMBERS IS NOT A POWER OF THREE.
  631. C        2.  THE QUANTITY  A  IN STATEMENT 10 IS REPRESENTED TO 
  632. C            THE ACCURACY USED IN FLOATING POINT VARIABLES
  633. C            THAT ARE STORED IN MEMORY.
  634. C     THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
  635. C     FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING 
  636. C     ASSUMPTION 2.
  637. C     UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
  638. C            A  IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
  639. C            B  HAS A ZERO FOR ITS LAST BIT OR DIGIT,
  640. C            C  IS NOT EXACTLY EQUAL TO ONE,
  641. C            EPS  MEASURES THE SEPARATION OF 1.0 FROM
  642. C                 THE NEXT LARGER FLOATING POINT NUMBER.
  643. C     THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
  644. C     ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
  645. C
  646. C     *****************************************************************
  647. C     THIS ROUTINE IS ONE OF THE AUXILIARY ROUTINES USED BY EISPACK III
  648. C     TO AVOID MACHINE DEPENDENCIES.
  649. C     *****************************************************************
  650. C
  651. C     THIS VERSION DATED 4/6/83.
  652. C
  653.       A = 4.0E0/3.0E0
  654.    10 B = A - 1.0E0
  655.       C = B + B + B
  656.       EPS = ABS(C-1.0E0)
  657.       IF (EPS .EQ. 0.0E0) GO TO 10
  658.       EPSLON = EPS*ABS(X)
  659.       RETURN
  660.       END
  661.       SUBROUTINE MM (A, LDA, N1, N3, B, LDB, N2, C, LDC)
  662.       REAL A(LDA,*), B(LDB,*), C(LDC,*)
  663. C
  664. C   PURPOSE:
  665. C     MULTIPLY MATRIX B TIMES MATRIX C AND STORE THE RESULT IN MATRIX A.
  666. C
  667. C   PARAMETERS:
  668. C
  669. C     A REAL(LDA,N3), MATRIX OF N1 ROWS AND N3 COLUMNS
  670. C
  671. C     LDA INTEGER, LEADING DIMENSION OF ARRAY A
  672. C
  673. C     N1 INTEGER, NUMBER OF ROWS IN MATRICES A AND B
  674. C
  675. C     N3 INTEGER, NUMBER OF COLUMNS IN MATRICES A AND C
  676. C
  677. C     B REAL(LDB,N2), MATRIX OF N1 ROWS AND N2 COLUMNS
  678. C
  679. C     LDB INTEGER, LEADING DIMENSION OF ARRAY B
  680. C
  681. C     N2 INTEGER, NUMBER OF COLUMNS IN MATRIX B, AND NUMBER OF ROWS IN
  682. C         MATRIX C
  683. C
  684. C     C REAL(LDC,N3), MATRIX OF N2 ROWS AND N3 COLUMNS
  685. C
  686. C     LDC INTEGER, LEADING DIMENSION OF ARRAY C
  687. C
  688. C ----------------------------------------------------------------------
  689. C
  690.       DO 20 J = 1, N3
  691.          DO 10 I = 1, N1
  692.             A(I,J) = 0.0
  693.    10    CONTINUE
  694.          CALL SMXPY (N2,A(1,J),N1,LDB,C(1,J),B)
  695.    20 CONTINUE
  696. C
  697.       RETURN
  698.       END
  699.       SUBROUTINE SMXPY (N1, Y, N2, LDM, X, M)
  700.       REAL Y(*), X(*), M(LDM,*)
  701. C
  702. C   PURPOSE:
  703. C     MULTIPLY MATRIX M TIMES VECTOR X AND ADD THE RESULT TO VECTOR Y.
  704. C
  705. C   PARAMETERS:
  706. C
  707. C     N1 INTEGER, NUMBER OF ELEMENTS IN VECTOR Y, AND NUMBER OF ROWS IN
  708. C         MATRIX M
  709. C
  710. C     Y REAL(N1), VECTOR OF LENGTH N1 TO WHICH IS ADDED THE PRODUCT M*X
  711. C
  712. C     N2 INTEGER, NUMBER OF ELEMENTS IN VECTOR X, AND NUMBER OF COLUMNS
  713. C         IN MATRIX M
  714. C
  715. C     LDM INTEGER, LEADING DIMENSION OF ARRAY M
  716. C
  717. C     X REAL(N2), VECTOR OF LENGTH N2
  718. C
  719. C     M REAL(LDM,N2), MATRIX OF N1 ROWS AND N2 COLUMNS
  720. C
  721. C ----------------------------------------------------------------------
  722. C
  723. C   CLEANUP ODD VECTOR
  724. C
  725.       J = MOD(N2,2)
  726.       IF (J .GE. 1) THEN
  727.          DO 10 I = 1, N1
  728.             Y(I) = (Y(I)) + X(J)*M(I,J)
  729.    10    CONTINUE
  730.       ENDIF
  731. C
  732. C   CLEANUP ODD GROUP OF TWO VECTORS
  733. C
  734.       J = MOD(N2,4)
  735.       IF (J .GE. 2) THEN
  736.          DO 20 I = 1, N1
  737.             Y(I) = ( (Y(I))
  738.      $             + X(J-1)*M(I,J-1)) + X(J)*M(I,J)
  739.    20    CONTINUE
  740.       ENDIF
  741. C
  742. C   CLEANUP ODD GROUP OF FOUR VECTORS
  743. C
  744.       J = MOD(N2,8)
  745.       IF (J .GE. 4) THEN
  746.          DO 30 I = 1, N1
  747.             Y(I) = ((( (Y(I))
  748.      $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
  749.      $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
  750.    30    CONTINUE
  751.       ENDIF
  752. C
  753. C   CLEANUP ODD GROUP OF EIGHT VECTORS
  754. C
  755.       J = MOD(N2,16)
  756.       IF (J .GE. 8) THEN
  757.          DO 40 I = 1, N1
  758.             Y(I) = ((((((( (Y(I))
  759.      $             + X(J-7)*M(I,J-7)) + X(J-6)*M(I,J-6))
  760.      $             + X(J-5)*M(I,J-5)) + X(J-4)*M(I,J-4))
  761.      $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
  762.      $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
  763.    40    CONTINUE
  764.       ENDIF
  765. C
  766. C   MAIN LOOP - GROUPS OF SIXTEEN VECTORS
  767. C
  768.       JMIN = J+16
  769.       DO 60 J = JMIN, N2, 16
  770.          DO 50 I = 1, N1
  771.             Y(I) = ((((((((((((((( (Y(I))
  772.      $             + X(J-15)*M(I,J-15)) + X(J-14)*M(I,J-14))
  773.      $             + X(J-13)*M(I,J-13)) + X(J-12)*M(I,J-12))
  774.      $             + X(J-11)*M(I,J-11)) + X(J-10)*M(I,J-10))
  775.      $             + X(J- 9)*M(I,J- 9)) + X(J- 8)*M(I,J- 8))
  776.      $             + X(J- 7)*M(I,J- 7)) + X(J- 6)*M(I,J- 6))
  777.      $             + X(J- 5)*M(I,J- 5)) + X(J- 4)*M(I,J- 4))
  778.      $             + X(J- 3)*M(I,J- 3)) + X(J- 2)*M(I,J- 2))
  779.      $             + X(J- 1)*M(I,J- 1)) + X(J)   *M(I,J)
  780.    50    CONTINUE
  781.    60 CONTINUE
  782.       RETURN
  783.       END
  785.