      subroutine bsmbk(n,lda,nb,a,y,x,work,info,me,nprocs,
     &            proclb,nbks,colsthisproc,blksthisproc)
      implicit none

*
*     -- NPAC, Syracuse University
*
*     written by David Koester
*     November 1993
*     
*     block forward substitution for dense matrices
*
      integer n,lda,nb,info,j,jb,nprocs,proclb,me
      integer proc,mydone,mycols,mypanel
      integer i1,i2,i3,j1,j2,j3,k
      integer previnring
      integer nbks,colsthisproc,blksthisproc
      integer PanelMsgType
      integer t3,t4
      real*8  a(*),y(*),x(*),work(*),one
      parameter (one=1.0e0)

      COMMON /LOCN/ PROC, MYDONE, MYCOLS, MYPANEL

      data PanelMsgType /3000/

      mypanel = 0
      mydone  = 0
      proc    = proclb
      mycols  = colsthisproc
      jb = n - ((n/nb)*nb)
      if (jb .eq. 0) jb = nb
      i2 = (n-jb)
*     ** set the value of j2 for all procs but the last proc
      j2 = (colsthisproc-nb)

*     ** for one processor,set nb = n
*     ** for Choleski factorization only
      if (nb .eq. n) then
        do i1 = n,1,-1
          x(i1) = y(i1) / a( (((i1-1)*n)+i1) )
          do j1 = 1,(i1-1)
            y(j1) = y(j1) - (x(i1) * a( (((j1-1)*n)+i1) ) )
          enddo
        enddo
      else
*        ** use blocked code **
         if (me .eq. proc) then
*         ** set the value of j2 for the last proc
          j2 = (colsthisproc-jb)
          do i1 = jb,1,-1
            x(i2+i1) = y(i2+i1) / a( (((j2+i1-1)*n)+(i2+i1)) )
            do j1 = 1,(i1-1)
              y(i2+j1) = y(i2+j1) - (x(i2+i1) * 
     &                     a((((j2+j1-1)*n)+(i2+i1))))
            enddo
          enddo
*         ** column block finished - so update j2 **
          j2 = (j2-nb)
*        ** broadcast the x-vector  **
	  call bcastdoublesd(jb,1,x(i2+1),'P',(PanelMsgType+1))
         endif

      DO j = 1,nbks
            if (j .gt. 1) jb = nb

	    if (proc .ne. me) then
*              ** receive the x-vector **
               call bcastdoublerc(t3,t4,x(i2+1),'P',
     &                       (PanelMsgType+j))
            endif

            proc = previnring(proc)

            if (j2 .ge. 0) then
               if (proc .eq. me) then
*                 ** i'm the next processor to calculate values of x **
*                 ** so only update the y-vector for the first panel **
                  i3 = (i2-nb)
                  do i1 = jb,1,-1
                    do j1 = 1,nb
                      y(i3+j1) = y(i3+j1) - (x(i2+i1) * 
     &                            a((((j2+j1-1)*n)+(i2+i1))))
                    enddo
                  enddo
*                 ** then substitute the values of the x-vector and send it **
                  do i1 = nb,1,-1
                    x(i3+i1) = y(i3+i1) / a( (((j2+i1-1)*n)+(i3+i1)) )
                    do j1 = 1,(i1-1)
                      y(j2+j1) = y(j2+j1) - (x(i3+i1) * 
     &                            a((((j2+j1-1)*n)+(i3+i1))))
                    enddo
                  enddo
*                 ** column block finished - so update j2 **
                  j2 = (j2-nb)
                  call bcastdoublesd(nb,1,x(i2-nb+1),'P',
     &                                       (PanelMsgType+j+1))
               endif

                j3 = j2
                do while (j3 .ge. 0)
                  k = (j3*nprocs) + (me*nb)
                  do i1 = jb,1,-1
                    do j1 = 1,nb
                      y(k+j1) = y(k+j1) - (x(i2+i1) * 
     &                            a((((j3+j1-1)*n)+(i2+i1))))
                    enddo
                  enddo
                  j3 = (j3-nb)
                enddo
*               ** row block finished - so update i2 **
                i2 = (i2-nb)
            else
             i2 = (i2-nb)
            endif
          enddo

      ENDIF

      RETURN
*
*     End of FRMBK .
*
      END





