cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c Gregor von Laszewski * North East Parallel Architectures Center
c                      * Center for Computational Science, Syracuse University
c                      * 111 College Place * Syracuse, NY 13244-4100
c                      * (315) 443-4883 * gregor@npac.syr.edu 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      PROGRAM NODE
      
#include "lu_const.f"

      real LU(mLDA,mLDA)
      real a(mLDA,mLDA)
      real b(mLDA)	
#ifdef SAXWPY
      real L(mLDA,MAX_BlockSize)
#endif
      real h,one,sum

      double precision t_start,t_stop,dt,clock0

      integer ipiv (mLDA)
      integer sdouble	
      integer ione
      integer N, M, processors, PanelSize

      integer msg(6), MsgLength,MsgType
      
      integer iseed(4) 

      integer me, host, nprocs, iarch, iord, idir,lda

      real MsgBlock(mLDA*MAX_BlockSize + 5)

      external open0, setarc0,bcast0

      lda = mLDA
      call open0 ( nprocs, me, host )
      call setarc0 (iprocs,iarch,iord,idir)



      print *, ' Node ', me, ' of ', nprocs, ' started'
      call SizeDouble (sdouble)
      call SizeInteger (MsgLength)
      Msglength = MsgLength * 5
      MsgType = HostInitMsg
      call bcast0 (msg,MsgLength,MsgType, host)

      m = msg(1)
      n = msg(2) 
      m = msg (3)
      processors = msg(4) 
      PanelSize = msg(5) 


      print *, me, ': Host Init msg recieved ', n,m,processors,PanelSize

      iseed(1) = 617
      iseed(2) = 1553
      iseed(3) = 2191
      iseed(4) = 3771
#ifdef GAXPY
      print *, me,' Strategy: GAXPY'
c      call MatGen(a,lda,m,n,iseed)
      call StdProblem(a,lda,m,n,b)
      call CopyMatrix (A,LU,n,lda)
#ifdef DEBUG
      call PrintMatrix (LU,m,n,lda,'A    ',me)	
      call PrintMatrix (B,m,1,lda,'B    ',me)	
#endif
#else
      print *, me,' Strategy: SAXPY'
      call MatrixGenerateSAXPY(A,lda,m,n,iseed,Processors,
     $     PanelSize, me)
#endif
      print *, me,' Matrix generated ', m,n

#ifdef DEBUG
      call tracenode (150000,0,1)
      call tracelevel (6,6,0)
#endif

      t_start = clock0()

#ifdef SAXPY
	
      call G_saxpy (A,L,lda,m,n,ipiv, PanelSize, Processors,
     $            me,host, MsgBlock)  
#else
      call G_gaxpy (LU,lda,m,n,ipiv, PanelSize, Processors,
     $            me,host, MsgBlock)  
#endif

      t_stop = clock0()

#ifdef DEBUG
      call tracelevel (0,0,0)
      call traceexit()
#endif
#if DEBUG
#ifdef SAXPY
      call PrintMatrix(a,n,m/Processors,m,'A  ',me)
#else
      call PrintMatrix(LU,n,m,lda,'A  ',me)
#endif
#endif

      call SendResult (m,n,PanelSize,lda,ipiv,LU,
     $     me,processors,MsgBlock)


      if (me .eq. 0) then
        call PrintMatrix(LU,n,m,lda,'A           ',me)
        call PrintMatrix(LU,n,m,lda,'LU          ',me)	
        call PrintTIVector (ipiv,n,'Pivt         ',0)

        do i=1,m
           if (i .ne. ipiv(i)) then
              h = b(i)
              b(i) = b(ipiv(i)) 
              b(ipiv(i)) = h
           end if
        end do
        call PrintMatrix (B,m,1,lda,'B           ',0)
	call strsv ('l','n','u',m,lu,lda,b,1)
        call strsv ('u','n','n',m,lu,lda,b,1)


        call PrintMatrix (B,m,1,lda,'X           ',0)
        sum = 0.0
        do i=1,m
           sum = sum + abs(b(i))
        end do
	write (*, 101) 'Relative error:', sum / (float(m)*one)
 101    format ( A,  E )

      end if

      dt = t_stop - t_start
      print *, 'Processor ',me, ' finished in ', dt

      ione = 1
      call send0 (dt, sdouble, 3, host)
      call sync0()

      call close0(ione)
      stop

      END


      subroutine GenMsgBlock (a,LDA,x0,y0,x1,y1,ipiv,PivotLength,
     $                        MsgLength,MsgBlock)

      integer  LDA, x0,y0,x1,y1
      real A(LDA, *)
      real MsgBlock (*)
      integer ipiv (*)
      integer MsgLength, PivotLength
      external sizereal
      integer index, i,j
      integer rsize

      MsgBlock(1) = float(x0)
      MsgBlock(2) = float(y0)
      MsgBlock(3) = float(x1)
      MsgBlock(4) = float(y1)
      index = 5

      do j=x0,x1
         do i=y0,y1
            MsgBlock(index) = A(i,j)
            index = index + 1
         end do
      end do
      do i=y0,y0+PivotLength - 1
         MsgBlock(index) = ipiv(i)
         index = index +1 
      end do

      call SizeReal   (rsize)

      MsgLength = rsize * (index - 1)


      return
      end

      subroutine GetMsgBlock (a,LDA, x0,y0,x1,y1,MsgBlock,
     $                        ipiv,PivotLength)

      real MsgBlock(*)
      integer LDA
      real A(LDA, *)
      integer ipiv (*)
      integer  x0,y0,x1,y1,index, PivotLength, i,j

      x0 = int(MsgBlock(1))
      y0 = int(MsgBlock(2))
      x1 = int(MsgBlock(3))
      y1 = int(MsgBlock(4))
      index = 5

      do j=x0,x1
         do i=y0,y1
            A(i,j) = MsgBlock(index) 
            index = index + 1
         end do
      end do
      do i=y0,y0 + PivotLength- 1
         ipiv(i) = int(MsgBlock(index))
         index = index +1 
      end do

      return
      end



