*===== == == == == == == == == == == == == == == == == == == == = ========== *
* cup4.f   for solving Ax=b as [A -b] X = 0 in row forms : Purcell Method
* MAIN Variables/Routines:
*     MAXP -- Max no of processors
*       A  -- The main working vector of size N x N/MAXP
*       U --- The final solution x
*       N --- The overall dimension
*   get_C -- Main routine to supply a row of [A -b] each time as a vector
*  NoProb -- Controls which MATRIX to use inside get_C
* NC_left -- No of columns left to work with on a proc
*===== ==  == == == == == == == == == == == == == == == == = ========== *
       PROGRAM  CUP4_F
       IMPLICIT   INTEGER(I-N)
       IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
       include 'mpif.h'
       double precision  TT1, TT2, mpi_wtime  ! 2^12=4096
       PARAMETER(NMAX=2**12,MAXP=1,NT=100,NA=NMAX/MAXP*(NMAX+2))
       DIMENSION A(NA),U(NMAX),RHS(NMAX+2)  ! A u = b = RHS
       dimension WT2(NMAX+2), IPIVOT(NMAX+2), WT3(NMAX+2), WL(NMAX+2)
       character*256 my_name
       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)

       NLAST = nproc-1 ! This last processor plays the role of Master
       master = 0        ! Some admin (not much)
          RHO = 0.4      ! Only active if NoProb=1 (off for 0, 2)

       do 99699 NoProb = 1, 4
 
c 3 (64) 4 (128) 5 (256) 6 (512) 7 (1024) 8 (2048) 9 (4096) 10 (8192)

       DO 99599 KE = 3,3  ! +/-1 for tests 

       N = 2**(KE+3) !!!!!!!!!!!!!!!!!!!!!   <======= DIM
       IF (N.GT.NMAX .or. N.LE.0) GOTO 99199
       I = N*(N+2)/NPROC
       IF (I.GT.NA) then
        PRINT*,'!!! Check Dimension Setting 4 A = ',I,N,NA
        STOP 'Abort'
       ENDIF

        tt1 = mpi_wtime()

c ==== A=V (not set yet), wt2=Ci=v1, wt3=Ci_p  RHS=CiVi=wk space

        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,NMAX+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 Purcell(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)

*END - END - END - END - END - END - END - END - END - END - END *

      IF (nc_left.ge.1) THEN
*_____________________________ Found Max Cpu _____________________
       WRITE(* ,'(''| cpu = '',G11.4, '' PROC/Prob ='',2I3,
     + '' N '', I4)') TT1A,NPROC,NoProb, N
*-----------------------------------------------------------------
         CALL L2ERROR(U,N, TJ, NoProb)
         TT1=TT1A ! from allreduce !_ Found Max Cpu _____________
         call mpi_get_processor_name(my_name,info,ierr)
      WRITE(*,'(''| Problem number ='',I2, '' On: '',A7)')
     +   NoProb, my_name
      WRITE(*,'(''| --cup4.f--  cpu = '',G11.4, '' PROC ='',I3,
     + '' N ='', I5,'' Er ='',E10.2/)') TT1,NPROC,N,TJ
      ENDIF !!!****************************************** LAST Proc
      
99199 CONTINUE
        call mpi_barrier(mpi_comm_world, ierr)
99599 CONTINUE
99699 CONTINUE !! New prob (NoProb) loop

      call mpi_finalize(ierr)
      IF (nc_left.ge.1) THEN
        print*, '<Exit> from (cup4.f)'
        STOP 
      else
        STOP
      ENDIF
      END

* x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+ x+
* IPIVOT / n+2 are for optimizing vector products (run with n=4,16 for demo)
c       call Purcell(U,A,WT2,WT3,RHS,WL,IPIVOT,
c    +        n,my_id,nproc,k_delta, nc_left)
      subroutine Purcell(u, V,Ci,Ci_p,CiVi, work,
     &           ipivot,n,my_id,nproc,k_delta, nc_left)
      IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
           IMPLICIT   INTEGER(I-N) ! k_delta = n/nproc
      include 'mpif.h'
      character*16 file_test
      dimension u(n), Ci_p(n+1),
     &          V(n+2,k_delta), Ci(n+2), CiVi(n+2),work(n+2) 
      dimension piv_in(2), piv_out(2) ! Pivoting control
      integer ipivot(n+2) ! pivoting history
      integer my_row
* Note: V has an extra column as it may recv a extra col for load balance

      if (n.eq.4 .or. n.eq.16) then
      write(file_test,'(''check_p'',I1,''.txt'')')my_id
c     print*,'Inside Purcell:',file_test
      OPEN(9, FILE=file_test,status='unknown')
      write(9,*)'_______________________ cup4.f _____________________',
     +          '| Non-exact Pivoting -- Narrowing Down Choice of   |',
     +          '| Pivoting for load balancing. No work shifting.   |',
     +          '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
      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

         NLAST = nproc-1
         myused = 0
         my_role = -1
            k_low = 1
            k_up  = k_delta
              TOL = 1.0D-12

C Step irow = 1 ______________________Column Partial Pivoting 
         irow = 1
              call get_C(Ci,irow,n)
*____ Find MAX guy - global pivot MY_ACT i.e. the first pivot (could bcast)
            my_act=1
            test_m=abs(Ci(1))
         do k=2,n+1
            if (test_m .le. abs(Ci(k))) then
             test_m=abs(Ci(k))
             my_act=k
            endif
         enddo 
             my_belong = (my_act - 1) / nproc
         if (my_belong.GT.NLAST) my_belong = NLAST

         k_first = my_id* (k_delta - nc_left)  ! Shift only used once
         C_first = Ci(my_act)      !--> All Know about the pivot

         do k=1,k_delta            !--> k+k_first: Global column index
                k_col = k+k_first
             if (k_col .GT. my_act) THEN
                k_col = k_col + 1  !--> Skip pivoting my_act / reduction
             endif
               V(k_col,k) = 1.0D0  !--> v_k = e_k + alphak * e_{my_act}
                  CiVi(k) = Ci(k_col)
                   alphak = - CiVi(k)/C_first
              V(my_act,k) = alphak !--> Complete the new v_k
                 V(n+2,k) = k_col  !--> Record col index v_k
                ipivot(1) = my_act !--> IPIVOT also nnz positions in v_k
         enddo
      if (n.eq.4 .or. n.eq.16) then!--Only for debugging/checking
       write(9,*)'my_act my_belong k_first k_delta C_first :',
     +            my_act,my_belong,k_first,k_delta,C_first 
       write(9,*)'______________ Step ',irow,' of ', n,' __cup4._'
       write(9,*)'Processor =',my_id, '       [ END of Row ',irow,' ]',
     +  ' Below is CiVi vector ...'
       write(9,'(8F10.5)') (CiVi(k),k=1,k_delta)
       write(9,*)'_____________ New V(:,',1,':',k_delta,') = '
        do j=1,n+2
         write(6,'(I2,2X, 16F9.4)') IPIVOT(j),(V(j,k),k=1,k_delta)
        enddo
         write(9,*)
       if (my_id.eq.0)
     + write(*,*)'______________ Step ',irow,' of ', n,' __cup4._'
      endif

*NOTE: IPIVOT and V_k's last position at n+2
*      determine exactly how many nnz in V_k: useful for product

C Step irow > 1 ________________________________________________
c ---------  Partial Pivoting Version of Parallel Purcell Method
         do 999 irow=2, n

******** my_role for role cycling *******************Control notation****
******** 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 my_assign(my_id,nproc, my_role,my_first,my_count)

* Each time at my_count=1, ALL to fetch a row (so p rows in all)  ! Cool
            if (my_count .eq. 1) then
              my_row1 = irow + my_id
              if (my_row1.gt.n) my_row1 = n
              call get_C(Ci_p,my_row1,n)
            endif 

* Of the p rows fetched, my_role=myused=0 proc is ready to bcast its row
            if (my_role.eq.myused) then
                DO k=1,n+1
                  Ci(k) = Ci_p(k) ! Ci_p: fetched storage / Ci bcast use
                ENDDO 
                  my_row = my_row1
            endif 

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

               k_first = 0

ccc ---- k=1,n+2-irow  !Form Product C^T*v_k for pivoting check (k_up+1=n+1)
         do k = k_low, k_up
	    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) ! local V1 always max *
            piv_in(1) = ABS(CiVi(1))
            piv_in(2) = my_id

         my_dimen = k_up - k_low + 1     !! actual dimension on proc id
         my_speed = n/nproc - (irow-2)/4/nproc**2  ! mean speed (by -1)

* ROUND 1 ______________________________________________________________________
*         With cup4, check speed of dim reduction before making piv_in on/off
      if (k_up.le.0 .OR. (my_dimen.LT.my_speed.AND.k_up.GE.1))
     +      piv_in(1) = -1.0  ! No involvement of IDLE PROC or Reduced Proc

      call mpi_allreduce(piv_in,piv_out,1,mpi_2double_precision, 
     +     mpi_maxloc, mpi_comm_world, ierr) 
      my_act = piv_out(2)     !--> Pivoting Proc ID found, pivot

* ROUND 2 ______________________________________________________________________
         if (piv_out(1).LE.TOL) THEN
            piv_in(1) = ABS(CiVi(1))
            piv_in(2) = my_id
            call mpi_allreduce(piv_in,piv_out,1,mpi_2double_precision, 
     +        mpi_maxloc, mpi_comm_world, ierr) 
            my_act = piv_out(2)
         endif

         call zerout(work, n+2)
         call zerout(  ci, n+2)

      if (my_id.eq.my_act) then
        if (n.eq.4 .or. N.eq.16) then !--Only for debugging/checking
        write(9,*)'|||||||||||||| Start',irow,' of ', n,' Pivot proc =',
     +     my_act, ' my work =', k_up+1-k_low, ' xx ID =',my_id
        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)
         my_row = irow  ! Pivoting Proc gets ready to bcast its v_k
        call my_pack(V, work, ipivot, irow, n+2, my_row)
      else
         my_row = irow+1
      endif

        call mpi_bcast(C_first, 1, mpi_double_precision,
     +     my_act, mpi_comm_world, ierr)  ! c^V1
        call mpi_bcast(work, my_row, mpi_double_precision,
     +     my_act, mpi_comm_world, ierr)  ! V1 length=irow (not yet N)

        call myn_pack(Ci, work, ipivot,  n+2, my_row)
             ipivot(irow) = int(Ci(n+2))  ! Source of pivot v_k

        if (n.eq.4 .or. N.eq.16) then!--Only for debugging/checking
         write(9,*)'my_id act role first =',
     +   my_id,my_act,my_role,my_first
         if (k_up.le.0) then
         write(9,*)' ............ Processor ',my_id, ' is idle.'
         else
         write(9,*)'_ _ _ _ _ _ _  _ _ _  This is Mr pivot vec V :'
         write(9,'(16F10.5)') (Ci(k),k=1,n+2)
         write(9,*)'_ _ _ _ _ _ _  _ _ _  This is Product Vector :'
         write(9,'(16F10.5)') (CiVi(k),k=1,k_up+k_first)
         write(9,*)' Pivot position =',ipivot(irow),
     +   ' since C_first=',C_first
         write(9,'(16F10.0)') (V(n+2,k),k=1,k_up+k_first)
         write(9,*)'                      This is C_first =',C_first
         write(9,*)'_____________ PER V(:,',k_low,':',
     +    k_up+k_first,') = ', k_first
          do j=1,n+2
            write(9,'(I2,2X, 16F9.4)') IPIVOT(j),
     +      (V(j,k),k=1,k_up+k_first)
          enddo
         endif 
        endif
 
c Update other V_j below (k_first=1 for pivoting proc and =0 for others)

      do k = k_low, k_up

      if (abs(CiVi(k+k_first)).GT.TOL) then ! Skip small products 
         alphak = - CiVi(k+k_first) / C_first
         Last_C=int(Ci(n+2))          ! The rest nnz's known to IPIVOT
         Last_K=int(V(n+2,k)) 
         if (k_first.eq.1) V(Last_K,k) = 0.0D0  
         Last_V=int(V(n+2,k+k_first)) ! The rest nnz's known to IPIVOT
         Just_C=1
         Just_V=1

         v_col = TOL                  ! Get hold of col max for scaling
      do i=1,irow
         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) ! ----------------- neue

          else  !--------------------------- Skipped small products 
               do j = 1,n+2
               V(j,k) = V(j,k+k_first)
               enddo
          endif ! -------------------------------------------------
      enddo ! ________ col k ends

        if (n.eq.4 .or. N.eq.16) then
        write(9,*)'_ Step ',irow,' of ', n,' __cup4._'
        if (k_up.le.0) then
         write(9,*)' ............ Processor ',my_id, ' is idle.'
        else
        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,'(16F10.5)') (CiVi(k),k=1+k_first,k_up+k_first)
         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_up)
         enddo
         write(9,*)
        endif
        endif
 
  999   continue ! irow

         nc_left = k_up+1-k_low
         v_col = V(n+1,1)
         if (abs(v_col).le.1.e-10) 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
        if ( nc_left.ge.1) then
         write(9,*)'cup4.f ___ Solution from Processor =',my_id
         write(*,*)'cup4.f ___ Solution from Processor =',my_id
         do j=1,n
         write(9,'(I2,2X,2G14.5)') j, u(j), V(j,1)
         write(*,'(I2,2X,2G14.5)') j, u(j), V(j,1)
         enddo
        endif
        close(9)
        endif
      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, 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)
      IF (k_up.LE.0) RETURN
      m_col = 1 ! record pivot col num
c     if (m_col.le.1) goto 9 ! no need to swap - disabled for now
      v_col = CiVi(1)
      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 my_assign(my_id,nproc, my_role,my_first,my_count)
      IMPLICIT   integer(i-n)
* * * my_role toggles between 0 1 2 ... nproc-1 
* * * different from test3b.f in that only my_role=0 is useful * * *
      if (my_role .eq. -1) then
          my_role = my_id
          my_first = 0
          my_count = 1
          return
      endif
          last_p = nproc-1
          my_count = my_count + 1
      if (my_count .ge. nproc+1) my_count=1

          my_role = my_role - 1
          my_first = my_first + 1
      if (my_role .eq. -1) my_role = last_p
      if (my_first .ge. nproc) my_first = my_first-nproc
      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 
* 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 
* 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
*____________________________________ Test Only
c     elseif (n.eq.16) then ! soln = 1:16
c       c(n+1)=0.0d0
c       DO k = 1, n
c         c(n+1)=c(n+1)-c(k)*k
c       ENDDO
      endif
*____________________________________ Test Only
      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 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'
