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.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  19.2 KB  |  699 lines
  1. *
  2. *PLEASE NOTE THAT netlib HAS MOVED, THE NEW ADDRESS IS netlib@ornl.gov.
  3. *THE OLD ADDRESS, netlib@mcs.anl.gov, WILL BE TURNED OFF SOON.
  4. *
  5. *** from netlib, Fri Jul 27 14:07:10 EDT 1990 ***
  6.       double precision second
  7.       double precision aa(200,200),a(201,200),b(200),x(200)
  8.       double precision time(8,6),cray,ops,total,norma,normx
  9.       double precision resid,residn,eps,epslon
  10.       integer ipvt(200)
  11.       lda = 201
  12.       ldaa = 200
  13. c
  14.       n = 100
  15.       cray = .056
  16.       write(6,1)
  17.     1 format(' Please send the results of this run to:'//
  18.      $       ' Jack J. Dongarra'/
  19.      $       ' Computer Science Department'/
  20.      $       ' University of Tennessee'/
  21.      $       ' Knoxville, Tennessee 37996-1300'//
  22.      $       ' Fax: 615-974-8296'//
  23.      $       ' Internet: dongarra@cs.utk.edu'/)
  24.       ops = (2.0d0*n**3)/3.0d0 + 2.0d0*n**2
  25. c
  26.          call matgen(a,lda,n,b,norma)
  27.          t1 = second()
  28.          call dgefa(a,lda,n,ipvt,info)
  29.          time(1,1) = second() - t1
  30.          t1 = second()
  31.          call dgesl(a,lda,n,ipvt,b,0)
  32.          time(1,2) = second() - t1
  33.          total = time(1,1) + time(1,2)
  34. c
  35. c     compute a residual to verify results.
  36. c
  37.          do 10 i = 1,n
  38.             x(i) = b(i)
  39.    10    continue
  40.          call matgen(a,lda,n,b,norma)
  41.          do 20 i = 1,n
  42.             b(i) = -b(i)
  43.    20    continue
  44.          call dmxpy(n,b,n,lda,x,a)
  45.          resid = 0.0
  46.          normx = 0.0
  47.          do 30 i = 1,n
  48.             resid = dmax1( resid, dabs(b(i)) )
  49.             normx = dmax1( normx, dabs(x(i)) )
  50.    30    continue
  51.          eps = epslon(1.0d0)
  52.          residn = resid/( n*norma*normx*eps )
  53.          write(6,40)
  54.    40    format('     norm. resid      resid           machep',
  55.      $          '         x(1)          x(n)')
  56.          write(6,50) residn,resid,eps,x(1),x(n)
  57.    50    format(1p5e16.8)
  58. c
  59.          write(6,60) n
  60.    60    format(//'    times are reported for matrices of order ',i5)
  61.          write(6,70)
  62.    70    format(6x,'dgefa',6x,'dgesl',6x,'total',5x,'mflops',7x,'unit',
  63.      $         6x,'ratio')
  64. c
  65.          time(1,3) = total
  66.          time(1,4) = ops/(1.0d6*total)
  67.          time(1,5) = 2.0d0/time(1,4)
  68.          time(1,6) = total/cray
  69.          write(6,80) lda
  70.    80    format(' times for array with leading dimension of',i4)
  71.          write(6,110) (time(1,i),i=1,6)
  72. c
  73.          call matgen(a,lda,n,b,norma)
  74.          t1 = second()
  75.          call dgefa(a,lda,n,ipvt,info)
  76.          time(2,1) = second() - t1
  77.          t1 = second()
  78.          call dgesl(a,lda,n,ipvt,b,0)
  79.          time(2,2) = second() - t1
  80.          total = time(2,1) + time(2,2)
  81.          time(2,3) = total
  82.          time(2,4) = ops/(1.0d6*total)
  83.          time(2,5) = 2.0d0/time(2,4)
  84.          time(2,6) = total/cray
  85. c
  86.          call matgen(a,lda,n,b,norma)
  87.          t1 = second()
  88.          call dgefa(a,lda,n,ipvt,info)
  89.          time(3,1) = second() - t1
  90.          t1 = second()
  91.          call dgesl(a,lda,n,ipvt,b,0)
  92.          time(3,2) = second() - t1
  93.          total = time(3,1) + time(3,2)
  94.          time(3,3) = total
  95.          time(3,4) = ops/(1.0d6*total)
  96.          time(3,5) = 2.0d0/time(3,4)
  97.          time(3,6) = total/cray
  98. c
  99.          ntimes = 10
  100.          tm2 = 0
  101.          t1 = second()
  102.          do 90 i = 1,ntimes
  103.             tm = second()
  104.             call matgen(a,lda,n,b,norma)
  105.             tm2 = tm2 + second() - tm
  106.             call dgefa(a,lda,n,ipvt,info)
  107.    90    continue
  108.          time(4,1) = (second() - t1 - tm2)/ntimes
  109.          t1 = second()
  110.          do 100 i = 1,ntimes
  111.             call dgesl(a,lda,n,ipvt,b,0)
  112.   100    continue
  113.          time(4,2) = (second() - t1)/ntimes
  114.          total = time(4,1) + time(4,2)
  115.          time(4,3) = total
  116.          time(4,4) = ops/(1.0d6*total)
  117.          time(4,5) = 2.0d0/time(4,4)
  118.          time(4,6) = total/cray
  119. c
  120.          write(6,110) (time(2,i),i=1,6)
  121.          write(6,110) (time(3,i),i=1,6)
  122.          write(6,110) (time(4,i),i=1,6)
  123.   110    format(6(1pe11.3))
  124. c
  125.          call matgen(aa,ldaa,n,b,norma)
  126.          t1 = second()
  127.          call dgefa(aa,ldaa,n,ipvt,info)
  128.          time(5,1) = second() - t1
  129.          t1 = second()
  130.          call dgesl(aa,ldaa,n,ipvt,b,0)
  131.          time(5,2) = second() - t1
  132.          total = time(5,1) + time(5,2)
  133.          time(5,3) = total
  134.          time(5,4) = ops/(1.0d6*total)
  135.          time(5,5) = 2.0d0/time(5,4)
  136.          time(5,6) = total/cray
  137. c
  138.          call matgen(aa,ldaa,n,b,norma)
  139.          t1 = second()
  140.          call dgefa(aa,ldaa,n,ipvt,info)
  141.          time(6,1) = second() - t1
  142.          t1 = second()
  143.          call dgesl(aa,ldaa,n,ipvt,b,0)
  144.          time(6,2) = second() - t1
  145.          total = time(6,1) + time(6,2)
  146.          time(6,3) = total
  147.          time(6,4) = ops/(1.0d6*total)
  148.          time(6,5) = 2.0d0/time(6,4)
  149.          time(6,6) = total/cray
  150. c
  151.          call matgen(aa,ldaa,n,b,norma)
  152.          t1 = second()
  153.          call dgefa(aa,ldaa,n,ipvt,info)
  154.          time(7,1) = second() - t1
  155.          t1 = second()
  156.          call dgesl(aa,ldaa,n,ipvt,b,0)
  157.          time(7,2) = second() - t1
  158.          total = time(7,1) + time(7,2)
  159.          time(7,3) = total
  160.          time(7,4) = ops/(1.0d6*total)
  161.          time(7,5) = 2.0d0/time(7,4)
  162.          time(7,6) = total/cray
  163. c
  164.          ntimes = 10
  165.          tm2 = 0
  166.          t1 = second()
  167.          do 120 i = 1,ntimes
  168.             tm = second()
  169.             call matgen(aa,ldaa,n,b,norma)
  170.             tm2 = tm2 + second() - tm
  171.             call dgefa(aa,ldaa,n,ipvt,info)
  172.   120    continue
  173.          time(8,1) = (second() - t1 - tm2)/ntimes
  174.          t1 = second()
  175.          do 130 i = 1,ntimes
  176.             call dgesl(aa,ldaa,n,ipvt,b,0)
  177.   130    continue
  178.          time(8,2) = (second() - t1)/ntimes
  179.          total = time(8,1) + time(8,2)
  180.          time(8,3) = total
  181.          time(8,4) = ops/(1.0d6*total)
  182.          time(8,5) = 2.0d0/time(8,4)
  183.          time(8,6) = total/cray
  184. c
  185.          write(6,140) ldaa
  186.   140    format(/' times for array with leading dimension of',i4)
  187.          write(6,110) (time(5,i),i=1,6)
  188.          write(6,110) (time(6,i),i=1,6)
  189.          write(6,110) (time(7,i),i=1,6)
  190.          write(6,110) (time(8,i),i=1,6)
  191.       stop
  192.       end
  193.       subroutine matgen(a,lda,n,b,norma)
  194.       double precision a(lda,1),b(1),norma
  195. c
  196.       init = 1325
  197.       norma = 0.0
  198.       do 30 j = 1,n
  199.          do 20 i = 1,n
  200.             init = mod(3125*init,65536)
  201.             a(i,j) = (init - 32768.0)/16384.0
  202.             norma = dmax1(dabs(a(i,j)), norma)
  203.    20    continue
  204.    30 continue
  205.       do 35 i = 1,n
  206.           b(i) = 0.0
  207.    35 continue
  208.       do 50 j = 1,n
  209.          do 40 i = 1,n
  210.             b(i) = b(i) + a(i,j)
  211.    40    continue
  212.    50 continue
  213.       return
  214.       end
  215.       subroutine dgefa(a,lda,n,ipvt,info)
  216.       integer lda,n,ipvt(1),info
  217.       double precision a(lda,1)
  218. c
  219. c     dgefa factors a double precision matrix by gaussian elimination.
  220. c
  221. c     dgefa is usually called by dgeco, but it can be called
  222. c     directly with a saving in time if  rcond  is not needed.
  223. c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  224. c
  225. c     on entry
  226. c
  227. c        a       double precision(lda, n)
  228. c                the matrix to be factored.
  229. c
  230. c        lda     integer
  231. c                the leading dimension of the array  a .
  232. c
  233. c        n       integer
  234. c                the order of the matrix  a .
  235. c
  236. c     on return
  237. c
  238. c        a       an upper triangular matrix and the multipliers
  239. c                which were used to obtain it.
  240. c                the factorization can be written  a = l*u  where
  241. c                l  is a product of permutation and unit lower
  242. c                triangular matrices and  u  is upper triangular.
  243. c
  244. c        ipvt    integer(n)
  245. c                an integer vector of pivot indices.
  246. c
  247. c        info    integer
  248. c                = 0  normal value.
  249. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  250. c                     condition for this subroutine, but it does
  251. c                     indicate that dgesl or dgedi will divide by zero
  252. c                     if called.  use  rcond  in dgeco for a reliable
  253. c                     indication of singularity.
  254. c
  255. c     linpack. this version dated 08/14/78 .
  256. c     cleve moler, university of new mexico, argonne national lab.
  257. c
  258. c     subroutines and functions
  259. c
  260. c     blas daxpy,dscal,idamax
  261. c
  262. c     internal variables
  263. c
  264.       double precision t
  265.       integer idamax,j,k,kp1,l,nm1
  266. c
  267. c
  268. c     gaussian elimination with partial pivoting
  269. c
  270.       info = 0
  271.       nm1 = n - 1
  272.       if (nm1 .lt. 1) go to 70
  273.       do 60 k = 1, nm1
  274.          kp1 = k + 1
  275. c
  276. c        find l = pivot index
  277. c
  278.          l = idamax(n-k+1,a(k,k),1) + k - 1
  279.          ipvt(k) = l
  280. c
  281. c        zero pivot implies this column already triangularized
  282. c
  283.          if (a(l,k) .eq. 0.0d0) go to 40
  284. c
  285. c           interchange if necessary
  286. c
  287.             if (l .eq. k) go to 10
  288.                t = a(l,k)
  289.                a(l,k) = a(k,k)
  290.                a(k,k) = t
  291.    10       continue
  292. c
  293. c           compute multipliers
  294. c
  295.             t = -1.0d0/a(k,k)
  296.             call dscal(n-k,t,a(k+1,k),1)
  297. c
  298. c           row elimination with column indexing
  299. c
  300.             do 30 j = kp1, n
  301.                t = a(l,j)
  302.                if (l .eq. k) go to 20
  303.                   a(l,j) = a(k,j)
  304.                   a(k,j) = t
  305.    20          continue
  306.                call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  307.    30       continue
  308.          go to 50
  309.    40    continue
  310.             info = k
  311.    50    continue
  312.    60 continue
  313.    70 continue
  314.       ipvt(n) = n
  315.       if (a(n,n) .eq. 0.0d0) info = n
  316.       return
  317.       end
  318.       subroutine dgesl(a,lda,n,ipvt,b,job)
  319.       integer lda,n,ipvt(1),job
  320.       double precision a(lda,1),b(1)
  321. c
  322. c     dgesl solves the double precision system
  323. c     a * x = b  or  trans(a) * x = b
  324. c     using the factors computed by dgeco or dgefa.
  325. c
  326. c     on entry
  327. c
  328. c        a       double precision(lda, n)
  329. c                the output from dgeco or dgefa.
  330. c
  331. c        lda     integer
  332. c                the leading dimension of the array  a .
  333. c
  334. c        n       integer
  335. c                the order of the matrix  a .
  336. c
  337. c        ipvt    integer(n)
  338. c                the pivot vector from dgeco or dgefa.
  339. c
  340. c        b       double precision(n)
  341. c                the right hand side vector.
  342. c
  343. c        job     integer
  344. c                = 0         to solve  a*x = b ,
  345. c                = nonzero   to solve  trans(a)*x = b  where
  346. c                            trans(a)  is the transpose.
  347. c
  348. c     on return
  349. c
  350. c        b       the solution vector  x .
  351. c
  352. c     error condition
  353. c
  354. c        a division by zero will occur if the input factor contains a
  355. c        zero on the diagonal.  technically this indicates singularity
  356. c        but it is often caused by improper arguments or improper
  357. c        setting of lda .  it will not occur if the subroutines are
  358. c        called correctly and if dgeco has set rcond .gt. 0.0
  359. c        or dgefa has set info .eq. 0 .
  360. c
  361. c     to compute  inverse(a) * c  where  c  is a matrix
  362. c     with  p  columns
  363. c           call dgeco(a,lda,n,ipvt,rcond,z)
  364. c           if (rcond is too small) go to ...
  365. c           do 10 j = 1, p
  366. c              call dgesl(a,lda,n,ipvt,c(1,j),0)
  367. c        10 continue
  368. c
  369. c     linpack. this version dated 08/14/78 .
  370. c     cleve moler, university of new mexico, argonne national lab.
  371. c
  372. c     subroutines and functions
  373. c
  374. c     blas daxpy,ddot
  375. c
  376. c     internal variables
  377. c
  378.       double precision ddot,t
  379.       integer k,kb,l,nm1
  380. c
  381.       nm1 = n - 1
  382.       if (job .ne. 0) go to 50
  383. c
  384. c        job = 0 , solve  a * x = b
  385. c        first solve  l*y = b
  386. c
  387.          if (nm1 .lt. 1) go to 30
  388.          do 20 k = 1, nm1
  389.             l = ipvt(k)
  390.             t = b(l)
  391.             if (l .eq. k) go to 10
  392.                b(l) = b(k)
  393.                b(k) = t
  394.    10       continue
  395.             call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
  396.    20    continue
  397.    30    continue
  398. c
  399. c        now solve  u*x = y
  400. c
  401.          do 40 kb = 1, n
  402.             k = n + 1 - kb
  403.             b(k) = b(k)/a(k,k)
  404.             t = -b(k)
  405.             call daxpy(k-1,t,a(1,k),1,b(1),1)
  406.    40    continue
  407.       go to 100
  408.    50 continue
  409. c
  410. c        job = nonzero, solve  trans(a) * x = b
  411. c        first solve  trans(u)*y = b
  412. c
  413.          do 60 k = 1, n
  414.             t = ddot(k-1,a(1,k),1,b(1),1)
  415.             b(k) = (b(k) - t)/a(k,k)
  416.    60    continue
  417. c
  418. c        now solve trans(l)*x = y
  419. c
  420.          if (nm1 .lt. 1) go to 90
  421.          do 80 kb = 1, nm1
  422.             k = n - kb
  423.             b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
  424.             l = ipvt(k)
  425.             if (l .eq. k) go to 70
  426.                t = b(l)
  427.                b(l) = b(k)
  428.                b(k) = t
  429.    70       continue
  430.    80    continue
  431.    90    continue
  432.   100 continue
  433.       return
  434.       end
  435.       subroutine daxpy(n,da,dx,incx,dy,incy)
  436. c
  437. c     constant times a vector plus a vector.
  438. c     jack dongarra, linpack, 3/11/78.
  439. c
  440.       double precision dx(1),dy(1),da
  441.       integer i,incx,incy,ix,iy,m,mp1,n
  442. c
  443.       if(n.le.0)return
  444.       if (da .eq. 0.0d0) return
  445.       if(incx.eq.1.and.incy.eq.1)go to 20
  446. c
  447. c        code for unequal increments or equal increments
  448. c          not equal to 1
  449. c
  450.       ix = 1
  451.       iy = 1
  452.       if(incx.lt.0)ix = (-n+1)*incx + 1
  453.       if(incy.lt.0)iy = (-n+1)*incy + 1
  454.       do 10 i = 1,n
  455.         dy(iy) = dy(iy) + da*dx(ix)
  456.         ix = ix + incx
  457.         iy = iy + incy
  458.    10 continue
  459.       return
  460. c
  461. c        code for both increments equal to 1
  462. c
  463.    20 continue
  464.       do 30 i = 1,n
  465.         dy(i) = dy(i) + da*dx(i)
  466.    30 continue
  467.       return
  468.       end
  469.       double precision function ddot(n,dx,incx,dy,incy)
  470. c
  471. c     forms the dot product of two vectors.
  472. c     jack dongarra, linpack, 3/11/78.
  473. c
  474.       double precision dx(1),dy(1),dtemp
  475.       integer i,incx,incy,ix,iy,m,mp1,n
  476. c
  477.       ddot = 0.0d0
  478.       dtemp = 0.0d0
  479.       if(n.le.0)return
  480.       if(incx.eq.1.and.incy.eq.1)go to 20
  481. c
  482. c        code for unequal increments or equal increments
  483. c          not equal to 1
  484. c
  485.       ix = 1
  486.       iy = 1
  487.       if(incx.lt.0)ix = (-n+1)*incx + 1
  488.       if(incy.lt.0)iy = (-n+1)*incy + 1
  489.       do 10 i = 1,n
  490.         dtemp = dtemp + dx(ix)*dy(iy)
  491.         ix = ix + incx
  492.         iy = iy + incy
  493.    10 continue
  494.       ddot = dtemp
  495.       return
  496. c
  497. c        code for both increments equal to 1
  498. c
  499.    20 continue
  500.       do 30 i = 1,n
  501.         dtemp = dtemp + dx(i)*dy(i)
  502.    30 continue
  503.       ddot = dtemp
  504.       return
  505.       end
  506.       subroutine  dscal(n,da,dx,incx)
  507. c
  508. c     scales a vector by a constant.
  509. c     jack dongarra, linpack, 3/11/78.
  510. c
  511.       double precision da,dx(1)
  512.       integer i,incx,m,mp1,n,nincx
  513. c
  514.       if(n.le.0)return
  515.       if(incx.eq.1)go to 20
  516. c
  517. c        code for increment not equal to 1
  518. c
  519.       nincx = n*incx
  520.       do 10 i = 1,nincx,incx
  521.         dx(i) = da*dx(i)
  522.    10 continue
  523.       return
  524. c
  525. c        code for increment equal to 1
  526. c
  527.    20 continue
  528.       do 30 i = 1,n
  529.         dx(i) = da*dx(i)
  530.    30 continue
  531.       return
  532.       end
  533.       integer function idamax(n,dx,incx)
  534. c
  535. c     finds the index of element having max. dabsolute value.
  536. c     jack dongarra, linpack, 3/11/78.
  537. c
  538.       double precision dx(1),dmax
  539.       integer i,incx,ix,n
  540. c
  541.       idamax = 0
  542.       if( n .lt. 1 ) return
  543.       idamax = 1
  544.       if(n.eq.1)return
  545.       if(incx.eq.1)go to 20
  546. c
  547. c        code for increment not equal to 1
  548. c
  549.       ix = 1
  550.       dmax = dabs(dx(1))
  551.       ix = ix + incx
  552.       do 10 i = 2,n
  553.          if(dabs(dx(ix)).le.dmax) go to 5
  554.          idamax = i
  555.          dmax = dabs(dx(ix))
  556.     5    ix = ix + incx
  557.    10 continue
  558.       return
  559. c
  560. c        code for increment equal to 1
  561. c
  562.    20 dmax = dabs(dx(1))
  563.       do 30 i = 2,n
  564.          if(dabs(dx(i)).le.dmax) go to 30
  565.          idamax = i
  566.          dmax = dabs(dx(i))
  567.    30 continue
  568.       return
  569.       end
  570.       double precision function epslon (x)
  571.       double precision x
  572. c
  573. c     estimate unit roundoff in quantities of size x.
  574. c
  575.       double precision a,b,c,eps
  576. c
  577. c     this program should function properly on all systems
  578. c     satisfying the following two assumptions,
  579. c        1.  the base used in representing dfloating point
  580. c            numbers is not a power of three.
  581. c        2.  the quantity  a  in statement 10 is represented to 
  582. c            the accuracy used in dfloating point variables
  583. c            that are stored in memory.
  584. c     the statement number 10 and the go to 10 are intended to
  585. c     force optimizing compilers to generate code satisfying 
  586. c     assumption 2.
  587. c     under these assumptions, it should be true that,
  588. c            a  is not exactly equal to four-thirds,
  589. c            b  has a zero for its last bit or digit,
  590. c            c  is not exactly equal to one,
  591. c            eps  measures the separation of 1.0 from
  592. c                 the next larger dfloating point number.
  593. c     the developers of eispack would appreciate being informed
  594. c     about any systems where these assumptions do not hold.
  595. c
  596. c     *****************************************************************
  597. c     this routine is one of the auxiliary routines used by eispack iii
  598. c     to avoid machine dependencies.
  599. c     *****************************************************************
  600. c
  601. c     this version dated 4/6/83.
  602. c
  603.       a = 4.0d0/3.0d0
  604.    10 b = a - 1.0d0
  605.       c = b + b + b
  606.       eps = dabs(c-1.0d0)
  607.       if (eps .eq. 0.0d0) go to 10
  608.       epslon = eps*dabs(x)
  609.       return
  610.       end
  611.       subroutine dmxpy (n1, y, n2, ldm, x, m)
  612.       double precision y(*), x(*), m(ldm,*)
  613. c
  614. c   purpose:
  615. c     multiply matrix m times vector x and add the result to vector y.
  616. c
  617. c   parameters:
  618. c
  619. c     n1 integer, number of elements in vector y, and number of rows in
  620. c         matrix m
  621. c
  622. c     y double precision(n1), vector of length n1 to which is added 
  623. c         the product m*x
  624. c
  625. c     n2 integer, number of elements in vector x, and number of columns
  626. c         in matrix m
  627. c
  628. c     ldm integer, leading dimension of array m
  629. c
  630. c     x double precision(n2), vector of length n2
  631. c
  632. c     m double precision(ldm,n2), matrix of n1 rows and n2 columns
  633. c
  634. c ----------------------------------------------------------------------
  635. c
  636. c   cleanup odd vector
  637. c
  638.       j = mod(n2,2)
  639.       if (j .ge. 1) then
  640.          do 10 i = 1, n1
  641.             y(i) = (y(i)) + x(j)*m(i,j)
  642.    10    continue
  643.       endif
  644. c
  645. c   cleanup odd group of two vectors
  646. c
  647.       j = mod(n2,4)
  648.       if (j .ge. 2) then
  649.          do 20 i = 1, n1
  650.             y(i) = ( (y(i))
  651.      $             + x(j-1)*m(i,j-1)) + x(j)*m(i,j)
  652.    20    continue
  653.       endif
  654. c
  655. c   cleanup odd group of four vectors
  656. c
  657.       j = mod(n2,8)
  658.       if (j .ge. 4) then
  659.          do 30 i = 1, n1
  660.             y(i) = ((( (y(i))
  661.      $             + x(j-3)*m(i,j-3)) + x(j-2)*m(i,j-2))
  662.      $             + x(j-1)*m(i,j-1)) + x(j)  *m(i,j)
  663.    30    continue
  664.       endif
  665. c
  666. c   cleanup odd group of eight vectors
  667. c
  668.       j = mod(n2,16)
  669.       if (j .ge. 8) then
  670.          do 40 i = 1, n1
  671.             y(i) = ((((((( (y(i))
  672.      $             + x(j-7)*m(i,j-7)) + x(j-6)*m(i,j-6))
  673.      $             + x(j-5)*m(i,j-5)) + x(j-4)*m(i,j-4))
  674.      $             + x(j-3)*m(i,j-3)) + x(j-2)*m(i,j-2))
  675.      $             + x(j-1)*m(i,j-1)) + x(j)  *m(i,j)
  676.    40    continue
  677.       endif
  678. c
  679. c   main loop - groups of sixteen vectors
  680. c
  681.       jmin = j+16
  682.       do 60 j = jmin, n2, 16
  683.          do 50 i = 1, n1
  684.             y(i) = ((((((((((((((( (y(i))
  685.      $             + x(j-15)*m(i,j-15)) + x(j-14)*m(i,j-14))
  686.      $             + x(j-13)*m(i,j-13)) + x(j-12)*m(i,j-12))
  687.      $             + x(j-11)*m(i,j-11)) + x(j-10)*m(i,j-10))
  688.      $             + x(j- 9)*m(i,j- 9)) + x(j- 8)*m(i,j- 8))
  689.      $             + x(j- 7)*m(i,j- 7)) + x(j- 6)*m(i,j- 6))
  690.      $             + x(j- 5)*m(i,j- 5)) + x(j- 4)*m(i,j- 4))
  691.      $             + x(j- 3)*m(i,j- 3)) + x(j- 2)*m(i,j- 2))
  692.      $             + x(j- 1)*m(i,j- 1)) + x(j)   *m(i,j)
  693.    50    continue
  694.    60 continue
  695.       return
  696.       end
  697.  
  698.  
  699.