C   
C     This software is copyright (C) 1991,  Regents  of  the
C     Universityof  California.   Anyone may reproduce fmath.f,
C     the software in this distribution, in whole or in part, pro-
C     vided that:
C
C    (1)  Any copy  or  redistribution  of  fmath.f  must  show  the
C         Regents  of  the  University of California, through its
C         Lawrence Berkeley Laboratory, as the source,  and  must
C         include this notice;
C
C    (2)  Any use of this software must reference this  distribu-
C         tion,  state that the software copyright is held by the
C         Regents of the University of California, and  that  the
C         software is used by their permission.
C
C         It is acknowledged that the U.S. Government has  rights
C    in  fmath.f  under  Contract DE-AC03-765F00098 between the U.S.
C    Department of Energy and the University of California.
C
C         fmath.f is provided as a professional  academic  contribu-
C    tion  for  joint exchange.  Thus it is experimental, is pro-
C    vided ``as is'', with no warranties of any kind  whatsoever,
C    no  support,  promise  of updates, or printed documentation.
C    The Regents of the University of California  shall  have  no
C    liability  with respect to the infringement of copyrights by
C    fmath.f, or any part thereof.     
C
C    Author:
C         Wes Bethel   
C         Lawrence Berkeley Laboratory
C         1 Cyclotron Rd.   Mail Stop 50-F
C         Berkeley CA 94720
C         415-486-6626
C         ewbethel@lbl.gov
C
C this fortran code contains many of the numerical, matrix-wise
C routines used by the analytical volume renderer.  fortran was 
C chosen due to the fact that it vectorizes and parallelizes
C easier than C code.
C
C     
C the subroutine m4_pos_mul_3d will do a positionwise multiplication
C of the two 3d matrices p1 & p2.
C
      subroutine m4_pos_mul_3d(p1,p2)
      real p1(1),p2(1)

      integer ii

      do 10 ii=1,64
	p1(ii) = p1(ii)*p2(ii)
 10   continue
      return
      end

C
C     this routine will do a positionwise multiplication of
C     two vectors, i.e. M1[i] = M1[i] * M2[i]
C
      subroutine m4_pos_mul(m1,m2)
      real m1(1),m2(1)
      integer ii

      do ii=1,16
         m1(ii) = m1(ii) * m2(ii)
      end do
      return
      end

C
C this is just a fixed length vector copy routine 
C V1[i] = V2[i]
C
      subroutine matrix_4_copy(p1,p2)
      real p1(1),p2(1)
      integer ii
      do ii=1,16
         p2(ii) = p1(ii)
      end do
      return
      end

C
C
      subroutine pile_mult(m1,m2,m3)
      real m1(1),m2(1),m3(1)
      integer ii

      do ii=1,64
         m3(ii) = m1(ii) * m2(ii)
      end do
      return
      end

C
C
      subroutine f_pile_copy(dest,src)
      real dest(1),src(1)
      integer ii

      do ii=1,64
         dest(ii) = src(ii)
      end do
      return
      end
C
C
      subroutine pile_sum(pile,rval)
      real pile(1),rval
      integer ii

      rval = 0.0
      do ii=1,64
         rval = rval + pile(ii)
      end do
      return
      end
C
C
      subroutine pile_mult_sum(source,weight,work,scalar)
      real source(1),weight(1),work(1),scalar
      integer ii

      scalar = 0.0
      do ii=1,64
         scalar = scalar + source(ii) * weight(ii)
      end do
      return
      end

C
C subroutine to find the min value in a pile
C
      subroutine f_pile_min(p,x)
      real p(1),x

      integer ii
      x = 10000000.0
      do ii=1,64
         if (p(ii) .lt. x) x = p(ii)
      end do
      return
      end
C
C subroutine to find the max value in a pile
C
      subroutine f_pile_max(p,x)
      real p(1),x

      integer ii
      x = -10000000.0
      do ii=1,64
         if (p(ii) .gt. x) x = p(ii)
      end do
      return
      end
         

C
C subroutine to positionwise vector multiply and sum
C
      subroutine f_mult_sum_8(v1,v2,scalar)
      real v1(1),v2(1),scalar

      integer ii
      scalar = 0.0
      do ii=1,8
         scalar = scalar + v1(ii)*v2(ii)
      end do
      return
      end

      subroutine f_8vect_min(v,x)
      real v(1),x
      integer ii
      x = 10000000.0
      do ii=1,8
         if (v(ii) .lt. x) x = v(ii)
      end do
      return
      end

      subroutine f_8vect_max(v,x)
      real v(1),x
      integer ii
      x = -10000000000.0
      do ii=1,8
         if (v(ii) .gt. x) x = v(ii)
      end do
      return
      end

      subroutine load_dvectors(d,v,n)
      double precision d(1),v(1)
      integer n
      integer ii
      do ii=1,n+1
         d(ii) = v(1)
         v(1) = v(1) + v(2)
         v(2) = v(2) + v(3)
         v(3) = v(3) + v(4)
      end do
      return
      end

      subroutine load_lvectors(d,v,n)
      double precision d(1),v(1)
      integer n
      integer ii
      do ii=1,n+1
         d(ii) = v(1)
         v(1) = v(1) + v(2)
      end do
      return
      end