c*************
c  Description:  This module contains vector math routines
c
c  Author:       JGS   4/2/90
c  $Revision: 1.1 $
c  $Modtime:   01 Jan 1996 22:44:02  $
c  
c*************
c
c*** 
	function M$VPROMOTE (type1,type2)
	character*1 M$VPROMOTE,type1,type2
	M$VPROMOTE = CHAR( MD$VPROMOTE(ICHAR(type1),ICHAR(type2)) )
	return
	end


	subroutine VECTORINTERFACE ()
	character*1 ctype
	integer*4 type

	entry M$VSTYPE(ctype)
	type = ICHAR(ctype)
	call MD$VSTYPE(type)
	return

	entry M$VTYPE(ctype)
	type = ICHAR(ctype)
	call MD$VTYPE(type)
	return

	entry M$STYPE(ctype)
	type = ICHAR(ctype)
	call MD$STYPE(type)
	return

	entry M$VTYPEQ(ctype)
	call MD$VTYPEQ(type)
	ctype = CHAR(type)
	return

	entry M$STYPEQ(ctype)
	call MD$STYPEQ(type)
	ctype = CHAR(type)
	return

	end	

c***
	function MD$VPROMOTE (type1,type2)
	include 'parameters.inc'
	integer*4 type,type1,type2
	type = MP_LONG
	if (type1.eq.MP_FLOAT.or.type2.eq.MP_FLOAT) type = MP_FLOAT
	if (type1.eq.MP_DOUBLE.or.type2.eq.MP_DOUBLE) type = MP_DOUBLE
	MD$VPROMOTE = type
	return
	end


	subroutine M$VECTORLIB ()
	include 'parameters.inc'
	parameter (L_bbuf=16)
	include 'buffers.inc'
	record /BUFFER/ VA,VB,VC,VD
	include 'values.inc'
	record /VALUE/ SA,SB,SC,S
	integer*4 count,block,flag,type,MV,MS,nbpt,nbps
c
c this is the context for DSES mapping areas of which the first 8 integers
c  are part of the BOC for that DSES.  This keeps the vector mode specific
c  to the DSES as well as to the HOST primitive that may be using this
	common /md$context/ idum(4),MV,MS,nbpt,nbps,idummy(120)
c
	integer*4 nbyt,n,m,i,j,ja,jb,ka,kb,kka,kkb
	integer*4 lr,li,imax,imin
	real*4	fr,fi
	real*8	dr,di

c***
	entry M$VINIT ()
	MV = MP_FLOAT
	MS = MP_FLOAT
	nbpt = 4
	return
c
c***
	entry MD$VSTYPE (type)
	MS = type
	entry MD$VTYPE (type)
	MV = type
	if (MV.eq.MP_DOUBLE) then
	  nbpt = 8
	elseif (MV.eq.MP_LONG .or. MV.eq.MP_FLOAT) then
	  nbpt = 4
	elseif (MV.eq.MP_INT) then
	  nbpt = 2
	elseif (MV.eq.MP_BYTE) then
	  nbpt = 1
	else
	  call M$ERROR('Unsupported vector type in MD$VTYPE: '//CHAR(MV))
	endif
	return
c
c***
	entry MD$STYPE (type)
	MS = type
	return
c
c***
	entry MD$VTYPEQ (type)
	type = MV
	return
c
c***
	entry MD$STYPEQ (type)
	type = MS
	return
c
c***
	entry M$VMOV (VA,ka,VB,kb,count)
	nbyt = nbpt
	goto 101

	entry M$CVMOV (VA,ka,VB,kb,count)
	nbyt = nbpt*2
	goto 101

	entry M$VMOVN (VA,ka,VB,kb,count,block)
	nbyt = block

 101	if (ka.eq.1.and.kb.eq.1) then
	  if (%loc(VA) .ne. %loc(VB) )
     &	  call M$MOVE (VA,VB,count*nbyt)
	  return
	endif
c
	ja = 1
	jb = 1
	if (nbyt.eq.4) then
	  do n = 1,count
	    VB.fbuf(jb) = VA.fbuf(ja)
	    ja = ja+ka
	    jb = jb+kb
	  enddo
	elseif (nbyt.eq.8) then
	  do n = 1,count
	    VB.dbuf(jb) = VA.dbuf(ja)
	    ja = ja+ka
	    jb = jb+kb
	  enddo
	elseif (nbyt.eq.2) then
	  do n = 1,count
	    VB.ibuf(jb) = VA.ibuf(ja)
	    ja = ja+ka
	    jb = jb+kb
	  enddo
	elseif (nbyt.eq.1) then
	  do n = 1,count
	    VB.bbuf(jb) = VA.bbuf(ja)
	    ja = ja+ka
	    jb = jb+kb
	  enddo
	elseif (nbyt.eq.16) then
	  do n = 1,count
	    VB.cdbuf(jb) = VA.cdbuf(ja)
	    ja = ja+ka
	    jb = jb+kb
	  enddo
	elseif (nbyt.gt.0) then
	  ja = 0
	  jb = 0
	  kka = ka*nbyt
	  kkb = kb*nbyt
	  do n = 1,count
	    VB.sbuf(jb+1:jb+nbyt) = VA.sbuf(ja+1:ja+nbyt)
	    ja = ja+kka
	    jb = jb+kkb
	  enddo
	else
	  call M$VERROR('M$xVMOVx')
	endif
	return
c
c***
	entry M$VMOVNB (VA,ka,VB,kb,count,block)
	nbyt = block
	ja = 0
	jb = 0
	do n = 1,count
	  VB.sbuf(jb+1:jb+nbyt) = VA.sbuf(ja+1:ja+nbyt)
	  ja = ja + ka
	  jb = jb + kb
	enddo
	return
c
c***
	entry M$VADD (VA,VB,VC,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    VC.fbuf(n) = VA.fbuf(n) + VB.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    VC.dbuf(n) = VA.dbuf(n) + VB.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,count
	    VC.lbuf(n) = VA.lbuf(n) + VB.lbuf(n)
	  enddo
	else
	  call M$VERROR('M$VADD')
	endif
	return
c
c***
	entry M$VSUB (VA,VB,VC,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    VC.fbuf(n) = VA.fbuf(n) - VB.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    VC.dbuf(n) = VA.dbuf(n) - VB.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,count
	    VC.lbuf(n) = VA.lbuf(n) - VB.lbuf(n)
	  enddo
	else
	  call M$VERROR('M$VSUB')
	endif
	return
c
c***
	entry M$VMUL (VA,VB,VC,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    VC.fbuf(n) = VA.fbuf(n) * VB.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    VC.dbuf(n) = VA.dbuf(n) * VB.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,count
	    VC.lbuf(n) = VA.lbuf(n) * VB.lbuf(n)
	  enddo
	else
	  call M$VERROR('M$VMUL')
	endif
	return
c
c***
	entry M$CVMUL (VA,VB,VC,count,flag)

	if (MV.eq.MP_FLOAT) then
	  if (flag.eq.0) then
	    do n = 1,count
	      VC.cfbuf(n) = VA.cfbuf(n) * VB.cfbuf(n)
	    enddo
	  elseif (flag.eq.1) then
	    do n = 1,count
	      VC.cfbuf(n) = CONJG(VA.cfbuf(n)) * VB.cfbuf(n)
	    enddo
	  elseif (flag.eq.2) then
	    do n = 1,count
	      VC.cfbuf(n) = VA.cfbuf(n) * CONJG(VB.cfbuf(n))
	    enddo
	  else
	    call M$ERROR ('Illegal conjugate flag in M$CVMUL')
	  endif
	elseif (MV.eq.MP_DOUBLE) then
	  if (flag.eq.0) then
	    do n = 1,count
	      VC.cdbuf(n) = VA.cdbuf(n) * VB.cdbuf(n)
	    enddo
	  elseif (flag.eq.1) then
	    do n = 1,count
	      VC.cdbuf(n) = DCONJG(VA.cdbuf(n)) * VB.cdbuf(n)
	    enddo
	  elseif (flag.eq.2) then
	    do n = 1,count
	      VC.cdbuf(n) = VA.cdbuf(n) * DCONJG(VB.cdbuf(n))
	    enddo
	  else
	    call M$ERROR ('Illegal conjugate flag in M$CVMUL')
	  endif
	elseif (MV.eq.MP_LONG) then
	  if (flag.eq.0) then
	    do n = 1,count*2,2
	      m = n+1
	      lr	 = VA.lbuf(n)*VB.lbuf(n) - VA.lbuf(m)*VB.lbuf(m)
	      VC.lbuf(m) = VA.lbuf(n)*VB.lbuf(m) + VA.lbuf(m)*VB.lbuf(n)
	      VC.lbuf(n) = lr ! incase output overwrites either input
	    enddo
	  elseif (flag.eq.1) then
	    do n = 1,count*2,2
	      m = n+1
	      lr	 = VA.lbuf(n)*VB.lbuf(n) + VA.lbuf(m)*VB.lbuf(m)
	      VC.lbuf(m) = VA.lbuf(n)*VB.lbuf(m) - VA.lbuf(m)*VB.lbuf(n)
	      VC.lbuf(n) = lr ! incase output overwrites either input
	    enddo
	  elseif (flag.eq.2) then
	    do n = 1,count*2,2
	      m = n+1
	      lr	 =  VA.lbuf(n)*VB.lbuf(n) + VA.lbuf(m)*VB.lbuf(m)
	      VC.lbuf(m) = -VA.lbuf(n)*VB.lbuf(m) + VA.lbuf(m)*VB.lbuf(n)
	      VC.lbuf(n) = lr ! incase output overwrites either input
	    enddo
	  else
	    call M$ERROR ('Illegal conjugate flag in M$CVMUL')
	  endif
	else
	  call M$VERROR('M$CVMUL')
	endif
	return
c
c***
	entry M$CRVMUL (VA,VB,VC,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    VC.cfbuf(n) = VA.cfbuf(n) * VB.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    VC.cdbuf(n) = VA.cdbuf(n) * VB.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  i = 0
	  do n = 1,count*2,2
	    i = i+1
	    VC.lbuf(n)   = VA.lbuf(n)   * VB.lbuf(i)
	    VC.lbuf(n+1) = VA.lbuf(n+1) * VB.lbuf(i)
	  enddo
	else
	  call M$VERROR('M$CRVMUL')
	endif
	return
c
c***
	entry M$VSADD (VA,SB,VC,count)

	call M$RETYPE(SB,MS,S,MV,1)
	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    VC.fbuf(n) = VA.fbuf(n) + S.f
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    VC.dbuf(n) = VA.dbuf(n) + S.d
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,count
	    VC.lbuf(n) = VA.lbuf(n) + S.l
	  enddo
	else
	  call M$VERROR('M$VSADD')
	endif
	return
c
c***
	entry M$CVSADD (VA,SB,VC,count)

	call M$RETYPE(SB,MS,S,MV,2)
	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    VC.cfbuf(n) = VA.cfbuf(n) + S.cf
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    VC.cdbuf(n) = VA.cdbuf(n) + S.cd
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,count*2,2
	    VC.lbuf(n)   = VA.lbuf(n)   + S.ll(1)
	    VC.lbuf(n+1) = VA.lbuf(n+1) + S.ll(2)
	  enddo
	else
	  call M$VERROR('M$CVSADD')
	endif
	return
c
c***
	entry M$VSMUL (VA,SB,VC,count)

	if (MV.eq.MP_FLOAT) then
	  call M$RETYPE(SB,MS,S,MV,1)
	  do n = 1,count
	    VC.fbuf(n) = VA.fbuf(n) * S.f
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  call M$RETYPE(SB,MS,S,MV,1)
	  do n = 1,count
	    VC.dbuf(n) = VA.dbuf(n) * S.d
	  enddo
	elseif (MV.eq.MP_LONG) then
	  if (MS.eq.MP_FLOAT .and. mod(SB.f,1.0).ne.0.0) then
	    do n = 1,count
	      VC.lbuf(n) = VA.lbuf(n) * SB.f
	    enddo
	  elseif (MS.eq.MP_DOUBLE .and. dmod(SB.d,1.d0).ne.0.d0) then
	    do n = 1,count
	      VC.lbuf(n) = VA.lbuf(n) * SB.d
	    enddo
	  else
	    call M$RETYPE(SB,MS,S,MP_LONG,1)
	    do n = 1,count
	      VC.lbuf(n) = VA.lbuf(n) * S.l
	    enddo
	  endif
	else
	  call M$VERROR('M$VSMUL')
	endif
	return
c
c***
	entry M$CVSMUL (VA,SB,VC,count)

	if (MV.eq.MP_FLOAT) then
	  call M$RETYPE(SB,MS,S,MV,2)
	  do n = 1,count
	    VC.cfbuf(n) = VA.cfbuf(n) * S.cf
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  call M$RETYPE(SB,MS,S,MV,2)
	  do n = 1,count
	    VC.cdbuf(n) = VA.cdbuf(n) * S.cd
	  enddo
	elseif (MV.eq.MP_LONG) then
	  if (MS.eq.MP_FLOAT .and. (mod(SB.ff(1),1.0).ne.0.0 .or. 
     &		mod(SB.ff(2),1.0).ne.0.0)) then
	    do n = 1,count*2,2
	      m = n+1
	      lr	 = VA.lbuf(n)*SB.ff(1) - VA.lbuf(m)*SB.ff(2)
	      VC.lbuf(m) = VA.lbuf(n)*SB.ff(2) + VA.lbuf(m)*SB.ff(1)
	      VC.lbuf(n)   = lr ! incase output overwrites either input
	    enddo
	  elseif (MS.eq.MP_DOUBLE .and. (dmod(SB.dd(1),1.d0).ne.0.d0 .or. 
     &		dmod(SB.dd(2),1.d0).ne.0.d0)) then
	    do n = 1,count*2,2
	      m = n+1
	      lr	 = VA.lbuf(n)*SB.dd(1) - VA.lbuf(m)*SB.dd(2)
	      VC.lbuf(m) = VA.lbuf(n)*SB.dd(2) + VA.lbuf(m)*SB.dd(1)
	      VC.lbuf(n)   = lr ! incase output overwrites either input
	    enddo
	  else
	    call M$RETYPE(SB,MS,S,MP_LONG,2)
	    do n = 1,count*2,2
	      m = n+1
	      lr	 = VA.lbuf(n)*S.ll(1) - VA.lbuf(m)*S.ll(2)
	      VC.lbuf(m) = VA.lbuf(n)*S.ll(2) + VA.lbuf(m)*S.ll(1)
	      VC.lbuf(n)   = lr ! incase output overwrites either input
	    enddo
	  endif
	else
	  call M$VERROR('M$CVSMUL')
	endif
	return
c
c***
	entry M$VSUM (VA,SB,count)

	if (MV.eq.MP_FLOAT) then
	  SB.f = 0.0
	  do n = 1,count
	    SB.f = SB.f + VA.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.d = 0.0
	  do n = 1,count
	    SB.d = SB.d + VA.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  SB.l = 0
	  do n = 1,count
	    SB.l = SB.l + VA.lbuf(n)
	  enddo
	else
	  call M$VERROR('M$VSUM')
	endif
	return
c
c***
	entry M$CVSUM (VA,SB,count)

	if (MV.eq.MP_FLOAT) then
	  SB.cf = 0.0
	  do n = 1,count
	    SB.cf = SB.cf + VA.cfbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.cd = 0.0
	  do n = 1,count
	    SB.cd = SB.cd + VA.cdbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  SB.ll(1) = 0
	  SB.ll(2) = 0
	  do n = 1,count*2,2
	    SB.ll(1) = SB.ll(1) + VA.lbuf(n)
	    SB.ll(2) = SB.ll(2) + VA.lbuf(n+1)
	  enddo
	else
	  call M$VERROR('M$CVSUM')
	endif
	return
c
c***
	entry M$VSUMS (VA,SB,SC,count)

	if (MV.eq.MP_FLOAT) then
	  SB.f = 0.0
	  SC.f = 0.0
	  do n = 1,count
	    SB.f = SB.f + VA.fbuf(n)
	    SC.f = SC.f + VA.fbuf(n)*VA.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.d = 0.0
	  SC.d = 0.0
	  do n = 1,count
	    SB.d = SB.d + VA.dbuf(n)
	    SC.d = SC.d + VA.dbuf(n)*VA.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  SB.l = 0
	  SC.l = 0
	  do n = 1,count
	    SB.l = SB.l + VA.lbuf(n)
	    SC.l = SC.l + VA.lbuf(n)*VA.lbuf(n)
	  enddo
	else
	  call M$VERROR('M$VSUMS')
	endif
	return
c
c***
	entry M$CVSUMS (VA,SB,SC,count)

	if (MV.eq.MP_FLOAT) then
	  SB.cf = 0.0
	  SC.cf = 0.0
	  do n = 1,count
	    SB.cf = SB.cf + VA.cfbuf(n)
	    SC.ff(1) = SC.ff(1) + REAL(VA.cfbuf(n))**2
	    SC.ff(2) = SC.ff(2) + AIMAG(VA.cfbuf(n))**2
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.cd = 0.0
	  SC.cd = 0.0
	  do n = 1,count
	    SB.cd = SB.cd + VA.cdbuf(n)
	    SC.dd(1) = SC.dd(1) + DREAL(VA.cdbuf(n))**2
	    SC.dd(2) = SC.dd(2) + DIMAG(VA.cdbuf(n))**2
	  enddo
	elseif (MV.eq.MP_LONG) then
	  SB.ll(1) = 0
	  SB.ll(2) = 0
	  SC.ll(1) = 0
	  SC.ll(2) = 0
	  do n = 1,count*2,2
	    SB.ll(1) = SB.ll(1) + VA.lbuf(n)
	    SB.ll(2) = SB.ll(2) + VA.lbuf(n+1)
	    SB.ll(1) = SB.ll(1) + VA.lbuf(n)*VA.lbuf(n)
	    SB.ll(2) = SB.ll(2) + VA.lbuf(n+1)*VA.lbuf(n+1)
	  enddo
	else
	  call M$VERROR('M$CVSUMS')
	endif
	return
c
c***
	entry M$VBLK (VA,VB,block,count)

	i = 0
	if (MV.eq.MP_FLOAT) then
	  do n = 1,count/block
	    VB.fbuf(n) = VA.fbuf(i+1)
	    do j = i+2,i+block
	      VB.fbuf(n) = VB.fbuf(n) + VA.fbuf(j)
	    enddo
	    i = i+block
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count/block
	    VB.dbuf(n) = VA.dbuf(i+1)
	    do j = i+2,i+block
	      VB.dbuf(n) = VB.dbuf(n) + VA.dbuf(j)
	    enddo
	    i = i+block
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,count/block
	    VB.lbuf(n) = VA.lbuf(i+1)
	    do j = i+2,i+block
	      VB.lbuf(n) = VB.lbuf(n) + VA.lbuf(j)
	    enddo
	    i = i+block
	  enddo
	else
	  call M$VERROR('M$VBLK')
	endif
	return
c
c***
	entry M$VLOG10 (VA,SA,VB,count)

	if (MV.eq.MP_FLOAT) then
	  call M$RETYPE(SA, MS, fr, MV, 1)
	  do n = 1, count
	    VB.fbuf(n) = log10( amax1( fr, VA.fbuf(n) ) )
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  call M$RETYPE(SA, MS, dr, MV, 1)
	  do n = 1, count
	    VB.dbuf(n) = dlog10( dmax1( dr, VA.dbuf(n) ) )
	  enddo
	else
	  call M$VERROR('M$VLOG10')
	endif
	return
c
c***
	entry M$CVBLK (VA,VB,block,count)

	i = 0
	if (MV.eq.MP_FLOAT) then
	  do n = 1,count/block
	    VB.cfbuf(n) = VA.cfbuf(i+1)
	    do j = i+2,i+block
	      VB.cfbuf(n) = VB.cfbuf(n) + VA.cfbuf(j)
	    enddo
	    i = i+block
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count/block
	    VB.cdbuf(n) = VA.cdbuf(i+1)
	    do j = i+2,i+block
	      VB.cdbuf(n) = VB.cdbuf(n) + VA.cdbuf(j)
	    enddo
	    i = i+block
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 1,2*count/block,2
	    VB.lbuf(n)   = VA.lbuf(i+1)
	    VB.lbuf(n+1) = VA.lbuf(i+2)
	    do j = i+3,i+2*block,2
	      VB.lbuf(n)   = VB.lbuf(n)   + VA.lbuf(j)
	      VB.lbuf(n+1) = VB.lbuf(n+1) + VA.lbuf(j+1)
	    enddo
	    i = i+2*block
	  enddo
	else
	  call M$VERROR('M$CVBLK')
	endif
	return
c
c***
	entry M$CVPHA (VA,VB,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    fr = REAL(VA.cfbuf(n))
	    fi = AIMAG(VA.cfbuf(n))
	    if (fr.eq.0. .and. fi.eq.0.) fr = 1.
	    VB.fbuf(n) = atan2(fi,fr)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    dr = DREAL(VA.cdbuf(n))
	    di = DIMAG(VA.cdbuf(n))
	    if (dr.eq.0. .and. di.eq.0.) dr = 1.
	    VB.dbuf(n) = datan2(di,dr)
	  enddo
	else
	  call M$VERROR('M$CVPHA')
	endif
	return
c
c***
	entry M$CVPHAD (VA,VB,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    fr = REAL(VA.cfbuf(n))
	    fi = AIMAG(VA.cfbuf(n))
	    if (fr.eq.0. .and. fi.eq.0.) fr = 1.
	    VB.fbuf(n) = atan2d(fi,fr)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    dr = DREAL(VA.cdbuf(n))
	    di = DIMAG(VA.cdbuf(n))
	    if (dr.eq.0. .and. di.eq.0.) dr = 1.
	    VB.dbuf(n) = datan2d(di,dr)
	  enddo
	else
	  call M$VERROR('M$CVPHAD')
	endif
	return
c
c***
	entry M$CVMAG2 (VA,VB,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 1,count
	    fr = REAL(VA.cfbuf(n))
	    fi = AIMAG(VA.cfbuf(n))
	    VB.fbuf(n) = fr*fr + fi*fi
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 1,count
	    dr = DREAL(VA.cdbuf(n))
	    di = DIMAG(VA.cdbuf(n))
	    VB.dbuf(n) = dr*dr + di*di
	  enddo
	else
	  call M$VERROR('M$CVMAG2')
	endif
	return
c
c***
	entry M$CVCONJ (VA,VB,count)

	if (MV.eq.MP_FLOAT) then
	  do n = 2,count*2,2
	    VB.fbuf(n) = -VA.fbuf(n)
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  do n = 2,count*2,2
	    VB.dbuf(n) = -VA.dbuf(n)
	  enddo
	elseif (MV.eq.MP_LONG) then
	  do n = 2,count*2,2
	    VB.lbuf(n) = -VA.lbuf(n)
	  enddo
	else
	  call M$VERROR('M$CVCONJ')
	endif
	return
c
c***
	entry M$VMXMN (VA,SB,SC,imax,imin,count)

	imax = 1
	imin = 1
	if (MV.eq.MP_FLOAT) then
	  SB.f = VA.fbuf(1)
	  SC.f = VA.fbuf(1)
	  do n = 2,count
	    if (VA.fbuf(n) .gt. SB.f) then
	      SB.f = VA.fbuf(n)
	      imax = n
	    elseif (VA.fbuf(n) .lt. SC.f) then
	      SC.f = VA.fbuf(n)
	      imin = n
	    endif
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.d = VA.dbuf(1)
	  SC.d = VA.dbuf(1)
	  do n = 2,count
	    if (VA.dbuf(n) .gt. SB.d) then
	      SB.d = VA.dbuf(n)
	      imax = n
	    elseif (VA.dbuf(n) .lt. SC.d) then
	      SC.d = VA.dbuf(n)
	      imin = n
	    endif
	  enddo
	elseif (MV.eq.MP_LONG) then
	  SB.l = VA.lbuf(1)
	  SC.l = VA.lbuf(1)
	  do n = 2,count
	    if (VA.lbuf(n) .gt. SB.l) then
	      SB.l = VA.lbuf(n)
	      imax = n
	    elseif (VA.lbuf(n) .lt. SC.l) then
	      SC.l = VA.lbuf(n)
	      imin = n
	    endif
	  enddo
	else
	  call M$VERROR('M$VMXMN')
	endif
	return

	end

c
c***
	subroutine M$CVMAG (VA, VB, count)
	include 'parameters.inc'
	complex*8 VA(*)
	real*4 VB(*)
	integer*4 count, n
	common /md$context/ idum(4),MV,MS,nbpt,nbps,idummy(120)

	if (MV.eq.MP_FLOAT) then
	  do n = 1, count
	    VB(n) = CABS(VA(n))
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  call CVMAG_D(VA, VB, count)
	else
	  call M$VERROR('M$CVMAG')
	endif
	return
	end

	subroutine CVMAG_D (VA, VB, count)
	implicit none
	complex*16 VA(*)
	real*8 VB(*)
	integer*4 n, count
	do n = 1, count
	  VB(n) = CDABS(VA(n))
	enddo
	return
	end

c
c***
	subroutine M$CVRPHAD (VA, VB, count, telem)
	include 'parameters.inc'
C Arguments
	real*4 VA(*), VB(*)
	integer*4 telem, count
c Locals
	integer*4 mtable
	parameter (mtable=32)
C Saved
	common /md$context/ idum(4),MV,MS,nbpt,nbps,idummy(120)
	save
	integer*4 table(2,mtable), ntable /0/
c Code
	if (telem.le.0) then		! free up phase tables
	  do n = 1, ntable
	    call M$FREE(table(2,n))
	  enddo
	  ntable = 0
	  return
	endif

	do n = 1, ntable		! look for tables already present
	  if (table(1,n).eq.telem) then
	    if (MV.eq.MP_FLOAT) then
	      call CVRPHAD_F(VA, VB, count, %val(table(2,n)), telem)
	    elseif (MV.eq.MP_DOUBLE) then
	      call CVRPHAD_D(VA, VB, count, %val(table(2,n)), telem)
	    else
	      call M$VERROR('M$CVRPHAD')
	    endif
	    return
	  endif
	enddo

	if (ntable.eq.mtable) then		! allocate and init new table
	  call M$FREE(table(2,ntable))
	else
	  ntable = max(1, ntable+1)
	endif
	table(1,ntable) = telem
	call M$MALLOC((telem+1)*nbpt, table(2,ntable))

	if (MV.eq.MP_FLOAT) then		! make call for first buffer
	  call CVRPHAD_INIT_F(%val(table(2,ntable)), telem)
	  call CVRPHAD_F(VA, VB, count, %val(table(2,ntable)), telem)
	elseif (MV.eq.MP_DOUBLE) then
	  call CVRPHAD_INIT_D(%val(table(2,ntable)), telem)
	  call CVRPHAD_D(VA, VB, count, %val(table(2,ntable)), telem)
	else
	  call M$VERROR('M$CVRPHAD')
	endif

	return
	end

	subroutine CVRPHAD_INIT_F (table, telem)
	implicit none
c Arguments
	real*4 table(*)
	integer*4 telem
c Locals
	real*4 scale
	integer*4 ii
c Code
	scale = 1./FLOAT(telem)
	do ii = 0, telem
	  table(ii+1) = ATAND(scale * FLOAT(ii))
	enddo
	return
	end

	subroutine CVRPHAD_F (VA, VB, count, table, telem)
	implicit none
c Arguments
	complex*8 VA(*)
	real*4	  VB(*), table(*)
	integer*4 count, telem
c Locals
	real*4 zr, zi, azr, azi, zz
	integer*4 ind, n
c Code
	do n = 1, count
	  zr = REAL(VA(n))
	  zi = AIMAG(VA(n))

	  azr = abs(zr)
	  azi = abs(zi)
	  if (azr.gt.azi) then
	    ind = float(telem)*azi/azr + 1
	    zz = table(ind)
	  else if (azr.lt.azi) then
	    ind = float(telem)*azr/azi + 1
	    zz = 90. - table(ind)
	  else if (azr.eq.0.) then 
	    zz = 0.
	  else
	    zz = 45.
	  endif

	  if (zr.lt.0.) then
	    if (zi.lt.0.) then
	      VB(n) = zz - 180.
	    else
	      VB(n) = 180. - zz
	    endif
	  else if (zi.lt.0.) then
	    VB(n) = -zz 
	  else
	    VB(n) = zz
	  endif
	enddo

	return
	end

	subroutine CVRPHAD_INIT_D (table, telem)
	implicit none
c Arguments
	real*8 table(*)
	integer*4 telem
c Locals
	real*8 scale
	integer*4 ii
c Code
	scale = 1./DBLE(telem)
	do ii = 0, telem
	  table(ii+1) = DATAND(scale * DBLE(ii))
	enddo
	return
	end

	subroutine CVRPHAD_D (VA, VB, count, table, telem)
	implicit none
c Arguments
	complex*16 VA(*)
	real*8	  VB(*), table(*)
	integer*4 count, telem
c Locals
	real*8 zr, zi, azr, azi, zz
	integer*4 ind, n
c Code
	do n = 1, count
	  zr = DREAL(VA(n))
	  zi = DIMAG(VA(n))

	  azr = dabs(zr)
	  azi = dabs(zi)
	  if (azr.gt.azi) then
	    ind = dble(telem)*azi/azr + 1
	    zz = table(ind)
	  else if (azr.lt.azi) then
	    ind = dble(telem)*azr/azi + 1
	    zz = 90.d0 - table(ind)
	  else if (azr.eq.0.) then 
	    zz = 0.d0
	  else
	    zz = 45.d0
	  endif

	  if (zr.lt.0.) then
	    if (zi.lt.0.) then
	      VB(n) = zz - 180.d0
	    else
	      VB(n) = 180.d0 - zz
	    endif
	  else if (zi.lt.0.d0) then
	    VB(n) = -zz 
	  else
	    VB(n) = zz
	  endif
	enddo

	return
	end

c
c***
	subroutine M$VECTORLIB2 ()
	include 'parameters.inc'
	parameter (L_bbuf=16)
	include 'buffers.inc'
	record /BUFFER/ VA,VB,VC,VD
	include 'values.inc'
	record /VALUE/ SA,SB,SC
	integer*4 count,n,MV,MS,nbpt,nbps
	common /md$context/ idum(4),MV,MS,nbpt,nbps,idummy(120)
c
c***
	entry M$VFILL (VA,SA,count)
	call M$VMOV (SA,0,VA,1,count)
	return
c
c***
	entry M$VNORM (VA,VB,count)
	
	if (count.le.0) return
	call M$VSUMS (VA,SB,SC,count)
	if (MV.eq.MP_FLOAT) then
	  SB.f = SB.f/count
	  SC.f = SC.f/count - SB.f*SB.f
	  if (SC.f .gt. 0.0) SC.f = 1./sqrt(SC.f)
	  do n = 1,count
	    VB.fbuf(n) = (VA.fbuf(n) - SB.f) * SC.f
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.d = SB.d/count
	  SC.d = SC.d/count - SB.d*SB.d
	  if (SC.d .gt. 0.0) SC.d = 1./dsqrt(SC.d)
	  do n = 1,count
	    VB.dbuf(n) = (VA.dbuf(n) - SB.d) * SC.d
	  enddo
	else
	  call M$VERROR('M$VNORM')
	endif
	return
c
c***
	entry M$CVNORM (VA,VB,count)
	
	if (count.le.0) return
	call M$CVSUMS (VA,SB,SC,count)
	if (MV.eq.MP_FLOAT) then
	  SB.cf = SB.cf/count
	  SC.ff(1) = SC.ff(1)/count - SB.ff(1)*SB.ff(1)
	  SC.ff(2) = SC.ff(2)/count - SB.ff(2)*SB.ff(2)
	  SC.f = CABS(SC.cf)
	  if (SC.f .gt. 0.0) SC.f = 1./sqrt(SC.f)
	  do n = 1,count
	    VB.cfbuf(n) = (VA.cfbuf(n) - SB.cf) * SC.f
	  enddo
	elseif (MV.eq.MP_DOUBLE) then
	  SB.cd = SB.cd/count
	  SC.dd(1) = SC.dd(1)/count - SB.dd(1)*SB.dd(1)
	  SC.dd(2) = SC.dd(2)/count - SB.dd(2)*SB.dd(2)
	  SC.d = CDABS(SC.cd)
	  if (SC.d .gt. 0.0) SC.d = 1./dsqrt(SC.d)
	  do n = 1,count
	    VB.cdbuf(n) = (VA.cdbuf(n) - SB.cd) * SC.d
	  enddo
	else
	  call M$VERROR('M$CVNORM')
	endif
	return
	end
c
c***
	subroutine M$VRAND (seed,VA,count)
c
c       The sequence is generated using a linear congruential
c       method based on the recursive formula:
c
c               IS <= a*IS + b, mod 2**26
c
c       a and b are constants.  Pre-scale by 2**26 so arithmetic 
c	need not be done mod 2**26, but rather mod 1 (= take fraction).
c       The starting value, IS, is called the seed; the last value
c       calculated is returned as the new seed.
c
c	Values used for a and b are from the FPSMath library.
c
c
	integer*4 type,count,seed
	real*8 a, b,bi, t26, sis
	parameter (a=67081293.0d0)
	parameter (b=14181771.0d0)
	parameter (t26=67108864.0d0)	! 2**26
	parameter (L_bbuf=16)
	include 'buffers.inc'
	record /BUFFER/ VA
	include 'parameters.inc'
c
	sis = dble(seed) / t26
	bi  = b / t26

	call MD$VTYPEQ(type)
	if (type.eq.MP_FLOAT) then
	  do n = 1,count
	    sis = sis*a + bi
	    sis = sis - INT(sis)
	    VA.fbuf(n) = sis
	  enddo
	elseif (type.eq.MP_DOUBLE) then
	  do n = 1,count
	    sis = sis*a + bi
	    sis = sis - INT(sis)
	    VA.dbuf(n) = sis
	  enddo
	elseif (type.eq.MP_LONG) then
	  do n = 1,count
	    sis = sis*a + bi
	    sis = sis - INT(sis)
	    VA.lbuf(n) = sis*1e6
	  enddo
	else
	  call M$VERROR('M$VRAND')
	endif
	seed = int(sis*t26)
	return
	end


	subroutine M$VERROR (name)
	character*(*) name
	call MD$VTYPEQ(MV)
	call M$ERROR ('Unsupported vector type ['//char(MV)//'] in '//name)
	return
	end

