C
C Calculate IBS sharing for all pair of relatives
C
      subroutine clcibs(nobs,ncol,ped,nrel,npairs,relist,
     &                  relate,mibs,vibs,typloc)
      double precision NA  
      parameter (NA=0.0d0)
      integer ncol, npairs, nobs, nrel
      integer relist(nrel)
      double precision ped(nobs,ncol)
      integer relate(npairs), typloc(npairs)
      double precision mibs(npairs), vibs(npairs)


      integer i,j,k,l,rel1,rel2,typed
      logical samefa, samemo
      double precision ibs, m, v


      k=0
      do 10 i=2,nrel
      do 10 j=1,(i-1)
        k=k+1
        rel1=relist(i)
        rel2=relist(j)
        relate(k)=0
        mibs(k)=0.0d0
        vibs(k)=0.0d0
        if (ped(rel1,1).eq.ped(rel2,1)) then
          samefa=(ped(rel1,3).ne.NA .and. ped(rel1,3).eq.ped(rel2,3))
          samemo=(ped(rel1,4).ne.NA .and. ped(rel1,4).eq.ped(rel2,4))
          if (samefa .and. samemo) then
            relate(k)=2
          else if (samefa .or. samemo) then
            relate(k)=3
          else if (ped(rel1,2).eq.ped(rel2,3) .or. 
     2             ped(rel1,2).eq.ped(rel2,4) .or.
     3             ped(rel2,2).eq.ped(rel1,3) .or. 
     4             ped(rel2,2).eq.ped(rel1,4)) then
            relate(k)=1
          end if
        end if
        typed=0
        m=0.0d0
        v=0.0d0
        do 20 l=6,ncol-1,2
        if (ped(rel1,l).ne.NA .and. ped(rel2,l).ne.NA) then
          typed=typed+1
          l2=l+1
          ibs=0.5d0
          if (ped(rel1,l).eq.ped(rel2,l) .and. 
     &        ped(rel1,l2).eq.ped(rel2,l2)) then
            ibs=1.0d0
          else if (ped(rel1,l).eq.ped(rel2,l2) .and. 
     &             ped(rel1,l2).eq.ped(rel2,l)) then
            ibs=2.0d0
          elseif (ped(rel1,l).ne.ped(rel2,l) .and. 
     2        ped(rel1,l).ne.ped(rel2,l2) .and.
     3        ped(rel1,l2).ne.ped(rel2,l) .and.
     4        ped(rel1,l2).ne.ped(rel2,l2)) then
            ibs=0.0d0
          end if
          call moment(typed,ibs,m,v)
        end if
   20   continue
        if (typed.gt.0) then
          mibs(k)=m
          vibs(k)=v/max(typed-1,1)
        end if
        typloc(k)=typed
   10 continue
      return
      end
C end-of-clcibs
C
C accumulate mean and sum-of-squares following AS41
C
      subroutine moment(n,x,mean,ss)
      integer n
      double precision mean, ss, x
C local variables
      double precision dev

      dev=x-mean
      mean=mean+dev/dfloat(n)
      ss=ss+dev*dev*dfloat(n-1)/dfloat(n)
      return
      end
C end-of-moment
C