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