SUBROUTINE DMRRRR (NRA, NCA, A, LDA, NRB, NCB, B, LDB, NRC, NCC,
     &                   C, LDC)
      include 'mpif.h'
      parameter (MASTER = 0)
      parameter (FROM_MASTER=1)
      parameter (FROM_WORKER=2)
      integer nra,nca,lda,nrb,ncb,ldb,nrc,ncc,ldc
      double precision a(nra,nca),b(nrb,ncb),c(nrc,ncc)
c     double precision a(lda,*),b(ldb,*),c(ldc,*)
      integer ierr, myrank, nprocs
      integer   numworkers,source,dest,nbytes,mtype,
     &          cols,avecol,extra, offset

      integer status(MPI_STATUS_SIZE)

      call MPI_INIT(ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)
      
      if (myrank==0) then 
      if (nca.ne.nrb) print*, '#Columns of A must equal #rows of B'
      if (ncb.ne.ncc) print*, '#Columns of B must equal #columns of C'
      if (nra.ne.nrc) print*, '#Rows of A must equal #rows of C'

      endif

      numworkers = nprocs-1

C * * * * * * Master Task * * * *
      if (myrank .eq. MASTER) then
      avecol = NCB/numworkers
      extra = mod(NCB, numworkers)
      offset = 1
      mtype = FROM_MASTER
      do 50 dest=1, numworkers
      if (dest .le. extra) then
        cols = avecol + 1
      else
        cols = avecol
      endif

      write(*,*)'   sending',cols,' cols to task',dest
      call MPI_SEND(offset,1, MPI_INTEGER, dest, mtype,
     .    MPI_COMM_WORLD, ierr)
      call MPI_SEND(cols, 1, MPI_INTEGER, dest, mtype,
     .    MPI_COMM_WORLD, ierr)
      call MPI_SEND(a, NRA*NCA, MPI_DOUBLE_PRECISION, dest,
     .    mtype, MPI_COMM_WORLD, ierr)
      call MPI_SEND(b(1,offset), cols*NCA, MPI_DOUBLE_PRECISION,
     .    dest, mtype, MPI_COMM_WORLD, ierr)

      offset = offset + cols
 50   continue

      mtype = FROM_WORKER
      do 60 i=1, numworkers
      source = i
      call MPI_RECV(offset, 1, MPI_INTEGER, source, mtype,
     .    MPI_COMM_WORLD, status, ierr)
      call MPI_RECV(cols, 1, MPI_INTEGER, source, mtype,
     .    MPI_COMM_WORLD, status, ierr)
      call MPI_RECV(c(1,offset), cols*NRA, MPI_DOUBLE_PRECISION,
     .    source, mtype, MPI_COMM_WORLD, status, ierr)
 60   continue

      endif



C *************************** worker task *************************************
      if (myrank > MASTER) then
C     Receive matrix data from master task
      mtype = FROM_MASTER
      call MPI_RECV(offset, 1, MPI_INTEGER, MASTER,
     .     mtype, MPI_COMM_WORLD, status, ierr)
      call MPI_RECV(cols, 1, MPI_INTEGER, MASTER,
     .     mtype, MPI_COMM_WORLD, status, ierr)
      call MPI_RECV(a, NRA*NCA, MPI_DOUBLE_PRECISION,
     .     MASTER, mtype, MPI_COMM_WORLD, status, ierr)
      call MPI_RECV(b, cols*NCA, MPI_DOUBLE_PRECISION,
     .     MASTER, mtype, MPI_COMM_WORLD, status, ierr)

C     Do matrix multiply
       do 100 k=1, cols
         do 100 i=1, NRA
         c(i,k) = 0.0
           do 100 j=1, NCA
           c(i,k) = c(i,k) + a(i,j) * b(j,k)
  100       continue

C     Send results back to master task

      mtype = FROM_WORKER
       call MPI_SEND(offset, 1, MPI_INTEGER, MASTER,
     .     mtype, MPI_COMM_WORLD, ierr)
       call MPI_SEND(cols, 1, MPI_INTEGER, MASTER,
     .     mtype, MPI_COMM_WORLD, ierr)
       call MPI_SEND(c, cols*NRA, MPI_DOUBLE_PRECISION,
     .     MASTER, mtype, MPI_COMM_WORLD, ierr)
      endif

      if (myrank==0) then
      print*, 'nprocs=',nprocs
      endif

      call MPI_FINALIZE(ierr)


c      do 10 i=1,nrc
c      do 10 j=1,ncc
c         c(i,j)=0.0
c         do 15 k=1,nca
c            c(i,j)=c(i,j)+a(i,k)*b(k,j)
c15       continue
c10    continue

      return
      end


      SUBROUTINE DEVCSF (N, A, LDA, EVAL, EVEC, LDEVEC)

      integer n,lda,ldevec
      double precision a(lda,lda),eval(lda),evec(lda,lda)
      double precision diag(lda),sdiag(lda)

      external tred2, tql2

      call tred2(lda,lda,a,diag,sdiag,evec)
      call tql2(lda,lda,diag,sdiag,evec,ierr)

      do 20 l=1,lda
         eval(l)=diag(l)
20    continue

      return
      end