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