       PROGRAM  DIRICHLET
         IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
           IMPLICIT   INTEGER(I-N)
            include 'mpif.h'
              double precision  TT1, TT2, MPI_WTIME
       PARAMETER(LMAX=2**12,MAXP=1,NA=(LMAX/MAXP+1)*(LMAX+2))
       DIMENSION A(NA),U(LMAX),RHS(LMAX+2)  ! Bare minimum
       dimension WT2(LMAX+2), IPIVOT(LMAX+2), WT3(LMAX+2)
       DIMENSION WL(LMAX+2) ! New work space for communications
       character*256 my_hostname
       COMMON   /WHICHP/RHO,NoProb
 
       call mpi_init(ierr)
       call mpi_comm_rank(mpi_comm_world, my_id, ierr)
       call mpi_comm_size(mpi_comm_world, nproc, ierr)
       call mpi_get_processor_name(my_hostname,info,ierr)
       NLAST = nproc-1 ! This last processor plays the role of Master
            RHO = 0.4D0

      do 99699 NoProb = 1, 9
 
c 3 (64) 4 (128) 5 (256) 6 (512) 7 (1024) 8 (2048) 9 (4096) 10 (8192)
         K_1 = 6
         K_2 = 9
       DO 99599 KE = K_1,K_2 ! 6, 10    !MAX=7
         N = 2**(KE+3) 
         IF (N.GT.LMAX .or. N.LE.0) GOTO 99199
         I = (N+2)*(N/NPROC+1)
       IF (I.GT.NA) then
        PRINT*,'!!! Check Dimension Setting 4 A = ',I,N,NA
        STOP 'Abort'
       ENDIF

          tt1 = mpi_wtime()

        k_delta = n/nproc
        nc_left = 0           ! Only meanful for proc=nlast
        if (my_id .EQ. NLAST) nc_left = N-k_delta*nproc
           k_delta = k_delta + nc_left !Add posssible cols

        do k=1,LMAX+2
           IPIVOT(k)=0  ! Collection of pivots for vector products
           RHS(k)=0.0
           WT2(k)=0.0
           WT3(k)=0.0
            WL(k)=0.0
        ENDDO
        call ZEROUT(A,NA) !! This will be matrix "V" later

        call paralel(U,A,WT2,WT3,RHS,WL,IPIVOT,
     +        n,my_id,nproc,k_delta, nc_left)

         tt2 = mpi_wtime()
         tt1 = tt2 - tt1
 
        call mpi_allreduce(tt1, tt1a, 1, mpi_double_precision,
     +     mpi_max, mpi_comm_world, ierr)

      if (my_id .eq. NLAST) then
       call mpi_get_processor_name(my_hostname,info,ierr)
       WRITE(*,'(''| Problem number ='',I2, '' On: '',A7,
     + '' (NPROC ='',I3,'')'' )') NoPROB, my_hostname, nproc
       CALL L2ERROR(U,N, TJ, NoProb)
       TT1=TT1A ! from allreduce !_ Found Max Cpu _____________
       WRITE(*,'(''| -cup5- MPI cpu = '',G11.4, '' PROC ='',I3,
     + '' N ='', I5,'' Er ='',E10.2/)') TT1,NPROC,N,TJ
      ENDIF !!!****************************************** NLAST
99199 CONTINUE  ! Abort N
        call mpi_barrier(mpi_comm_world, ierr)
99599 CONTINUE  ! New N
99699 CONTINUE  ! New Problem

      call mpi_finalize(ierr)
        IF (my_id .EQ. NLAST) THEN
        STOP '          <Exit> from cup5.f'
        else
        STOP
        ENDIF
      END

* xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx

      subroutine paralel(u, V,Ci,Ci_p,CiVi, work,
     &           ipivot,n,my_id,nproc,k_delta,nc_left)
           IMPLICIT   INTEGER(I-N) ! k_delta = n/nproc
      IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
      include 'mpif.h'
      character*16 file_test
      dimension u(n), Ci_p(n+1), piv_in(2), piv_out(2), ! Pivoting 
     &          V(n+2,k_delta+1), Ci(n+2), CiVi(n+2),work(n+2) 
      integer ipivot(n+2) 
      COMMON    /WHICHP/RHO,NoProb

      if (n.eq.8 .or. n.eq.16) then !----------------------------debug
      write(file_test,'(''check_p'',I1,''.txt'')')my_id
      OPEN(9, FILE=file_test,status='unknown')
      write(9,*)'_______________________ cup5.f ____________________',
     +          '| Exact Pivoting Implemented. Load balancing Ach- |',
     +          '| ieved by work shifting cyclically to Pivot Proc |',
     +          '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
      if (my_id.eq.0) then
        write(9,*)' Matrix C below: ',n,' X', n
        do irow=1,n
        call get_C(Ci,irow,n)
        write(9,'(8F9.3)') (Ci(k),k=1,n+1)
        enddo
      endif
      endif !----------------------------------------------------debug

         NLAST = nproc-1
         mrprof = 0
         mr_count = 0
         my_turn = -1
         my_role = -1
            k_low = 1
            k_up  = k_delta
              TOL = 1.0D-12

C Step irow = 1 ________________________________________________
c ---------  Partial Pivoting Version of Parallel Purcell Method
         irow = 1
              call get_C(Ci,irow,n)
*____ Find the first pivot and set 2:(n+1) with pivot sub by v_1
           my_act=1
           test_m=dabs(Ci(1))
         do k=2,n+1
           if (test_m .LT. dabs(Ci(k))) then
             test_m=dabs(Ci(k))
             my_act=k
           endif
         enddo 
           ipivot(1) = my_act
!----Idea 1 is to swap the pivot col V_k with V_1 (unless V_1 is pivot) 
              k_first = my_id*(k_delta-nc_left) + 1
              pivot_me = 0
              C_first = Ci(my_act) ! pivot - all knows
         do k=1,k_delta
            if (k+k_first .eq. my_act) then
              CiVi(k)=Ci(1)
              V(1,k) = 1.0D0
              V(n+2,k) = 1   !! Store col index of original V_1
              pivot_me = 1
            else
              CiVi(k)=Ci(k+k_first)
              V(k+k_first,k) = 1.0D0
              V(n+2,k) = k+k_first   !! Store col index of original V_j
            endif
              alphak = - CiVi(k)/C_first
              V(my_act,k) = alphak    !! Proper pivot vector V_{my_act}
         enddo
!Also IDea 2 is to advance the counter if k_col => my_act

        if (n.eq.8 .or. N.eq.16) then !---------------------------debug
        if (my_id.eq.0)
     +  write(*,*)'__ Step ',irow,' of ', n,' __cup5._'
        write(9,*)'__ Step ',irow,' of ', n,' __cup5._'
        write(9,*)'Processor =',my_id, '       [ END of Row ',irow,' ]'
        write(9,'(8F10.5)') (CiVi(k),k=1,k_delta)
        write(9,*)'_ New V(:,',1,':',k_delta,') = '
         do j=1,n+2
         write(9,'(I2,2X, 16F9.4)') IPIVOT(j),(V(j,k),k=1,k_delta)
         enddo
         write(9,*)
        endif !---------------------------------------------------debug

        nc_left = N-(k_delta-nc_left)*nproc !nc_left reset for my_assign

C Step irow > 1 ______Full Partial Pivoting Version of Parallel Purcell
c     write(*,*)'Proc ',my_id,' row ',irow,' COLS =',k_up-k_low+1,
c    +            'my_role','my_first','my_count'

         do 999 irow=2, n

******** my_role for role cycling ***************************************
******** my_first for ID_NO of my_role = 0 processor ********************
******** my_count for next group of "p" rows to be generated ************
******** my_act for ID_NO of pivoting processor (decided by MPI_MAXLOC) *

         call  who_send(my_id,nproc, my_turn,my_tid,my_count)
         call my_assign(my_id,nproc, my_role,my_first,irow,nc_left)

            if (my_count .eq. 1) then
              my_row = irow + my_id
              if (my_row.gt.n) my_row = n
              call get_C(Ci_p,my_row,n)
            endif 

            if (my_turn.eq.mrprof) then
                DO k=1,n+1
                  Ci(k) = Ci_p(k)
                ENDDO 
            endif 

      call mpi_barrier(mpi_comm_world, ierr)
      call mpi_bcast(Ci, n+1, mpi_double_precision,
     +             my_tid, mpi_comm_world, ierr)

           k_first = 0     ! Shift (-0 for no-pivot / =1 for pivot)

      do k = k_low, k_up   ! Product C^T*v_k  (k_up+1=n+1)
       CiVi(k) = scalar_prod(Ci,V(1,k),irow-1,ipivot,n+2)
      enddo
      call my_pivot(CiVi,k_up,V,n+2,k_delta,ip) ! V1 always max *
            piv_in(1) = ABS(CiVi(1))
            piv_in(2) = my_id
            if (k_up.le.0) piv_in(1) = -1.0  ! No involvement of IDLE
      call mpi_allreduce(piv_in,piv_out,1,mpi_2double_precision, 
     +     mpi_maxloc, mpi_comm_world, ierr) 
      my_act = piv_out(2)

      call zerout(WORK,n+2)
      call zerout(  CI,n+2)

      if (my_id.ne.my_act .and. my_role.eq.mrprof) then

         if (n.eq.8 .or. N.eq.16) then !-----------------------debug
         if (ip.GT.1)write(9,*)'Processor =',my_id,' Pivot Proc = ',
     +       my_act
         endif !-----------------------------------------------debug

          call my_pack(V, work, ipivot, irow, n+2, my_row)
          if (my_row.eq.0) print*,'my_row = 0 !!! on P',my_id
          myrow  = my_row
      call mpi_send(myrow,1,mpi_integer, my_act,701,
     +                      mpi_comm_world, ierr) 
      endif

      if (my_id.eq.my_act .and. my_role.ne.mrprof) then
      call mpi_recv(my_row,1,mpi_integer, mpi_any_source,
     +        701, mpi_comm_world, status, ierr)
      endif

      call mpi_barrier(mpi_comm_world, ierr)

***** ******* Sending processor = active role but not pivot ***********
      if (my_id.ne.my_act .and. my_role.eq.mrprof) then
      call mpi_send(CiVi(1),1,mpi_double_precision,my_act,
     +        801, mpi_comm_world, ierr) 
      call mpi_send(work,my_row,mpi_double_precision,my_act,
     +        901, mpi_comm_world, ierr) 
          k_first = 1
          k_up = k_up - 1  ! One V vector less
      endif

***** ******* Recving processor = inactive role but pivot ***********
      if (my_id.eq.my_act .and. my_role.ne.mrprof) then
          k_up = k_up + 1
      call mpi_recv(CiVi(k_up),1,mpi_double_precision, mpi_any_source,
     +        801, mpi_comm_world, status, ierr)
      call mpi_recv(work,my_row, mpi_double_precision, mpi_any_source,
     +        901, mpi_comm_world, status, ierr)
          do k=1,n+2
          V(k,k_up) = 0.0     ! Overall, Number_cols still decreases.
          enddo
       call myn_pack(V(1,k_up), work, ipivot, irow, n+2, my_row)
      endif
***** ******* ***** ******* ***** ******* ***** ******* ******* *******

        call mpi_barrier(mpi_comm_world, ierr)

      if (my_id.eq.my_act) then
        if (n.eq.8 .or. N.eq.16) then !---------------------------------
        if (ip.GT.1)write(9,*)'Processor =',my_id,' ___ Pivot = ',my_act
        endif !---------------------------------------------------------
             k_first = 1
               k_up  = k_up-1
             C_first = CiVi(1)
             call my_pack(V, work, ipivot, irow, n+2, my_row)
             piv_in(1)=C_first
             piv_in(2)=my_row ! save 1 bcast
      endif

        call mpi_bcast(piv_in, 2, mpi_double_precision,
     +     my_act, mpi_comm_world, ierr)  ! c^V1 product + my_row
             C_first = piv_in(1)
              my_row = INT(piv_in(2))

        call mpi_bcast(work, my_row, mpi_double_precision,
     +     my_act, mpi_comm_world, ierr)  !V1 pivot with length my_row

        call myn_pack(Ci, work, ipivot, irow, n+2, my_row)
             ipivot(irow) = int(Ci(n+2))

        if (n.eq.8 .or. N.eq.16) then !------------------------debug
       write(9,*)'my_id act role first =',my_id,my_act,my_role,my_first
        write(9,*)'   xxxx ',irow,' of ', n,' Per/Send_cv'
        write(9,*)'__ Role =',my_role,' Root_ID = ',
     +             my_first, ' Column_shift =',k_first,' (prod below)'
        write(9,'(8F10.5)') (CiVi(k),k=1+k_first,k_up+k_first)
        write(9,*)'__ ______ This is Mr pivot vec V :'
        write(9,'(8F10.5)') (Ci(k),k=1,n+2)
        if (k_up.le.0) then
         write(9,*)'... Processor ',my_id, ' is idle.'
        else
         write(9,*)'_ PER V(:,',k_low,':',k_up,') = '
         do j=1,n+2
         write(9,'(I2,2X, 16F9.4)') IPIVOT(j),
     +   (V(j,k),k=1+k_first,k_up+k_first)
         enddo
         write(9,*)'                      This is C_first =',C_first
        endif 
        endif !------------------------------------------------debug
 
      do k = k_low, k_up  ! Local Columns of V_k (excluding pivot v)

         alphak = - CiVi(k+k_first) / C_first
      if (dabs(alphak).GT.TOL) then 
         Last_K=int(V(n+2,k))
         if (k_first.eq.1) V(Last_K,k) = 0.0d0
         Last_C=int(Ci(n+2))          ! Other nnz positions in IPIVOT
         Last_V=int(V(n+2,k+k_first)) ! Other nnz positions in IPIVOT
         Just_C=1
         Just_V=1
         v_col = TOL
      do i=1,irow    !! Update column v_k at its nnz positions only
         j = ipivot(i)
         V(j,k) = alphak*Ci(j) + V(j,k+k_first)
         Just_C=min(Just_C,abs(j-Last_C))
         Just_V=min(Just_V,abs(j-Last_V))
          v_col=max(v_col,abs(V(j,k)))
      enddo
       j=Last_C
         if (Just_C.NE.0) V(j,k) = alphak*Ci(j) + V(j,k+k_first)
           v_col=max(v_col,abs(V(j,k)))
       j=Last_V
         if (Just_V.NE.0) V(j,k) = alphak*Ci(j) + V(j,k+k_first)
            v_col=max(v_col,abs(V(j,k)))

         if (v_col .GE. 2.0d0) then !__Found a col max in abs_
            do i=1,irow
               j = ipivot(i)
               V(j,k) = V(j,k) / v_col
            enddo
              j=Last_C
            if (Just_C.NE.0) V(j,k) = V(j,k) / v_col
              j=Last_V
            if (Just_V.NE.0) V(j,k) = V(j,k) / v_col
         endif !_Scaling Case_
               V(n+2,k) = V(n+2,k+k_first) 
          else !--------------------------- Skip small \alpha_k
               do j = 1,n+2
               V(j,k) = V(j,k+k_first)
               enddo
          endif ! ---- if TOL
      enddo ! ________ do k ends

        if (n.eq.8 .or. N.eq.16) then !--------------------------debug
        write(9,*)'__ Step ',irow,' of ', n,' __cup5._'
        write(9,*)'Processor =',my_id, ' [END of Row',irow,']',
     +  ' Scaling has occurred for vectors V'
        write(9,*)'__ Role =',my_role,' Root_ID = ',
     +             my_first, ' Column_shift =',k_first,' (prod below)'
        write(9,'(8F10.5)') (CiVi(k),k=1+k_first,k_up+k_first)
        if (k_up.le.0) then
         write(9,*)'... Processor ',my_id, ' is idle.'
        else
         write(9,*)'__ New V(:,',k_low,':',k_up,') = '
         do j=1,n+2
         write(9,'(I2,2X, 16F9.4)') IPIVOT(j),
     +   (V(j,k),k=1+k_first,k_up+k_first)
         enddo
         write(9,*)
        endif
        endif !-------------------------------------------------debug
 
c c   write(9,*)'Proc ',my_id,' row ',irow,' cols =',k_up-k_low+1,
c c  +            my_role,my_first,my_count

  999   continue ! irow

         v_col = abs(V(n+1,1))
         if (v_col.le.TOL) v_col = 1.0D0
         do j = 1, n
            u(j) = V(j,1)/v_col
         enddo
        if (n.eq.4 .or. N.eq.16) then !----------------------debug
        if (my_id.eq.NLAST) then
         write(9,*)'cup5.f ___ Solution from Processor =',my_id
         do j=1,n
         write(9,'(I2,2X,1G14.5)') j, u(j)
         enddo
        endif
        close(9)
        endif !----------------------------------------------debug
      return
      end

C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine my_pack(Ci, pack, ipivot, irow, n2, ip)
      IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT   integer(i-n)
      DIMENSION Ci(n2),ipivot(n2),pack(n2)
         Last=Ci(n2)
      do k=1,irow-1
         j=ipivot(k)
         pack(k)=Ci(j)
      enddo
          j=ipivot(irow)
            ip=irow
         if (j.eq.0) then
            pack(ip)=Ci(Last)
         else
            pack(ip)=Ci(j)
         endif
            ip=irow+1
            pack(ip)=Last
      return
      end
C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine myn_pack(Ci, pack, ipivot, irow, n2, ip)
      IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT   integer(i-n)
      DIMENSION Ci(n2),ipivot(n2),pack(n2)
         Ci(n2)=pack(ip)
         Last=Ci(n2)
      do k=1,ip-2
         j=ipivot(k)
         Ci(j)=pack(k)
      enddo
         j=ipivot(ip-1)
         if (j.eq.0) then
           Ci(Last)=pack(ip-1)
         else
           Ci(j)=pack(ip-1)
         endif
      return
      end

C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine my_pivot(CiVi,k_up,V,n1,k_delta,m_col)
      IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT   integer(i-n)
      DIMENSION CiVi(k_up),V(n1,k_delta)
      m_col = 1 ! record pivot col num
      v_col = 1.0E-10
      do k = 1, k_up  ! not including col V(:,n+1)
      if (abs(CiVi(k)) .gt. v_col) then
          m_col = k
          v_col = abs(CiVi(k))
      endif
      enddo
      if (m_col.eq.0) stop 'Pivoting fails'
      if (m_col.le.1) goto 9 ! no need to swap
      v_col = CiVi(1)
      CiVi(1)=CiVi(m_col)
      CiVi(m_col)=v_col
      DO k = 1, n1
      v_col = V(k,1)
      V(k,1) = V(k,m_col)
      V(k,m_col) = v_col
      ENDDO ! swap cols 1 and m_col in Matrix V
    9 return
      end

C234567891123456789212345678931234567894123456789512345678961234567897
      double precision function scalar_prod(C,vk,irow,ipivot,n)
      IMPLICIT  DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT  integer (I-N)
      DIMENSION ipivot(n), C(n),vk(n) !n=n+2 usual
         Last=vk(n)
         Just=1
         scalar_prod = 0.0
      do k=1,irow
         j=ipivot(k)
	 scalar_prod = scalar_prod + C(j)*vk(j)
         Just=min(Just,abs(j-Last))
      enddo
         j=Last
         if (Just.NE.0) scalar_prod = scalar_prod + C(j)*vk(j)
      return
      end 

C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine who_send(my_id,nproc, my_turn,my_tid,my_count)
      IMPLICIT   integer(i-n)
* * * my_turn toggles between 0 1 2 ... p-1 / only my_turn=0 is useful
      if (my_turn .eq. -1) then
          my_turn = my_id
          my_tid = 0
          my_count = 1
          return
      endif
          NLAST = nproc-1
          my_count = my_count + 1
      if (my_count .ge. nproc+1) my_count=1

          my_turn = my_turn - 1
          my_tid = my_tid + 1
      if (my_turn .eq. -1) my_turn = NLAST
      if (my_tid .ge. nproc) my_tid = my_tid-nproc
      return
      end

C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine my_assign(my_id,nproc, my_role,my_first,irow,nc_left)
      IMPLICIT   integer(i-n)
* * * my_role toggles between 0 1 2 ... p-1 / only my_role=0 is useful
          NLAST = nproc-1
      if (my_role .eq. -1 .and. nc_left.eq.0) then
          my_role = my_id
          my_first = 0
          return
      endif
      if (irow .LT. nc_left+2) then ! Need to reduce load in NLAST
          my_role = my_id+1
          if (my_role .gt. nlast) my_role = 0
          my_first = NLAST
          return
      endif
      if (irow .EQ. nc_left+2) then
          my_role = my_id
          my_first = 0
          return
      endif
          my_role = my_role - 1
          my_first = my_first + 1
      if (my_role .eq. -1) my_role = NLAST
      if (my_first .ge. nproc) my_first = my_first-nproc
      return
      end 
C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine L2ERROR(U,N,TJ,NoProb)
      IMPLICIT  NONE
      integer n,j,NoProb
      double precision U(n),TJ
        TJ = 0.0D0
        do j = 1, N
         TJ = TJ + (U(j)-1.0d0)**2
        enddo
        TJ = DSQRT(TJ / DBLE(N))
      return
      end 

C234567891123456789212345678931234567894123456789512345678961234567897
      subroutine get_C(C,irow,n)
      IMPLICIT  NONE
      integer n,irow,jcol,NoProb
      double precision C(n+1),RHO
      COMMON    /WHICHP/RHO,NoProb
* NoProb is used for 1 -- new rho=0.4  (1.4)^n-1 growth / Wilkinson/ GV89
*                    2 -- new rho=0.4  (1.4)^n-1  A^t
*                    3 -- new rho=1        2^n-1 growth worst (Higham's 89)
*                    4 -- new rho=1               A^t
* --------------- / -------------------------------------\
c                |  1     0     0     0     0     0     1 |
c                | -1     1     0     0     0     0     1 |
c                | -1    -1     1     0     0     0     1 |
c A=gfpp(7)=A =  | -1    -1    -1     1     0     0     1 |
c                | -1    -1    -1    -1     1     0     1 |
c                | -1    -1    -1    -1    -1     1     1 |
c                | -1    -1    -1    -1    -1    -1     1 |
* --------------- \ -------------------------------------/
c A=round( gfpp(12,.4)*10 )
c    10     0     0     0     0     0     0     0     0     0     0    10
c    -4    10     0     0     0     0     0     0     0     0     0    10
c    -4    -4    10     0     0     0     0     0     0     0     0    10
c    -4    -4    -4    10     0     0     0     0     0     0     0    10
c    -4    -4    -4    -4    10     0     0     0     0     0     0    10
c    -4    -4    -4    -4    -4    10     0     0     0     0     0    10
c    -4    -4    -4    -4    -4    -4    10     0     0     0     0    10
c    -4    -4    -4    -4    -4    -4    -4    10     0     0     0    10
c    -4    -4    -4    -4    -4    -4    -4    -4    10     0     0    10
c    -4    -4    -4    -4    -4    -4    -4    -4    -4    10     0    10
c    -4    -4    -4    -4    -4    -4    -4    -4    -4    -4    10    10
c    -4    -4    -4    -4    -4    -4    -4    -4    -4    -4    -4    10
C---------------------------------------------------------------------------

          IF (NoProb.Eq.1) THEN !________Rho case 1 of 0<Rho<1 for (1+Rho)^n-1
* Assume that x_j = 1 to work out C(n+1) or RHS easily ...
        DO jcol = 1, N-1
          C(jcol) = -RHO
          IF (jcol.GT.irow) C(jcol) = 0.0D0
          IF (jcol.EQ.irow) C(jcol) = 1.0D0
        ENDDO
          C(N) = 1.0D0
          C(N+1) = 2.0D0 - RHO*DBLE(irow-1)
          IF (irow.GE.N) C(N+1) = 1.0D0 - RHO*DBLE(irow-1)
          C(N+1) = -C(N+1)

      ELSEIF (NoProb.Eq.2) THEN !________Rho   case 2  A^T=Transpose 1
* Assume that x_j = 1 to work out C(n+1) or RHS easily ...
        DO jcol = 1, N
          C(jcol) = -RHO
          IF (jcol.LT.irow) C(jcol) = 0.0D0
          IF (jcol.EQ.irow) C(jcol) = 1.0D0
          IF (irow.EQ.N   ) C(jcol) = 1.0D0
        ENDDO
          C(N+1) = 1.0D0 - RHO*DBLE(N-irow)
          IF (irow.GE.N) C(N+1) = DBLE(N)
          C(N+1) = -C(N+1)
      ELSEIF (NoProb.Eq.3) THEN !________
* Assume that x_j = 1 to work out C(n+1) or RHS easily ...
        DO jcol = 1, N-1
          C(jcol) = -1.0D0
          IF (jcol.GT.irow) C(jcol) = 0.0D0
          IF (jcol.EQ.irow) C(jcol) = 1.0D0
        ENDDO
          C(N) = 1.0D0
          C(N+1) = 2.0D0 - DBLE(irow-1)
          IF (irow.GE.N) C(N+1) = 1.0D0 - DBLE(irow-1)
          C(N+1) = -C(N+1) 
      ELSEIF (NoProb.Eq.4) THEN !________Rho=1 case 4  A^T=Transpose 3 
* Assume that x_j = 1 to work out C(n+1) or RHS easily ...
        DO jcol = 1, N
          C(jcol) = -1.0D0
          IF (jcol.LT.irow) C(jcol) = 0.0D0
          IF (jcol.EQ.irow) C(jcol) = 1.0D0
          IF (irow.EQ.N   ) C(jcol) = 1.0D0
        ENDDO
          C(N+1) = 1.0D0 - DBLE(N-irow)
          IF (irow.GE.N) C(N+1) = DBLE(N)
          C(N+1) = -C(N+1) 
      ELSE !_____________________
          stop 'Check NoProb setting'
      ENDIF
*____________________________________ Test Only
      IF (n.eq.4) then
      IF (irow.eq.1) then
      C(1) = 5
      C(2) = 1
      C(3) = 2
      C(4) = 1
      C(5) =     -17.0
      else IF (irow.eq.2) then
      C(1) = 2
      C(2) = 10
      C(3) = 3
      C(4) = 1
      C(5) =     -35.0
      else IF (irow.eq.3) then
      C(1) = 1
      C(2) = 4
      C(3) = 8
      C(4) = 2
      C(5) =     -41.0
      else IF (irow.eq.4) then
      C(1) = 6
      C(2) = 2
      C(3) = 4
      C(4) = 20
      C(5) =     -102.0
      endif
      endif
      return
      end

C234567891123456789212345678931234567894123456789512345678961234567897 

      subroutine zerout(X,N)
      integer N,I
      real*8 X(N)
      do i=1,n
         X(i)=0.0d0
      enddo
      RETURN
      END

c     include 'mpif.h'
