!==============================================================================!
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!==============================================================================!
!==============================================================================!
!                                                                              |
!              WSMP_SETUP1      Oct. 2008                                      |
!              WSMP_SETUP2      Oct. 2008                                      |
!                                                                              |
!==============================================================================!
!==============================================================================!

subroutine wsmp_setup1 (n,n_iproc_st,n_iproc_end,n_iproc,nz_loc,ldb,ndof,vo,osolve,tpl,params,threadinfo,istep,iter)
use definitions
use threads
implicit none

integer n,n_iproc_st,n_iproc_end,n_iproc,nz_loc,ldb
integer ndof
type (void) vo
type (topology) :: tpl(n)
type (octreesolve) osolve
type (parameters) params
type (thread) threadinfo
integer istep,iter


logical, dimension(:), allocatable :: iproc_col
integer icol(100000),kk,nn,n0,jproc,idg_locp,nz
integer ie,k1,k2,j,i2,i1,idof1,idof2,k,err,i3,k3,i
integer iproc,nproc,ierr
character*72  :: shift
integer, dimension(:), allocatable :: irn,jcn

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

!==============================================================================!
!==============================================================================!

! n is the size of the system/matrix
! nz is the total number of non-zero terms in the global matrix
! n_iproc is the number of columns on given thread iproc
! n_iproc_st is the first column on given thread iproc, between 1 and n
! n_iproc_end is the last column on given thread iproc, between 1 and n
! nz_loc is the total number of non-zero terms on given thread iproc
! ia is an integer array of size n_iproc+1 (see wsmp manual)
! ja is an interer array of size nz_loc
! avals is a double precision array of size nz_loc containing the values
! b is a double precision array of size n_iproc containing the rhs

!==============================================================================!
!=====[allocate memory for topology]===========================================!
!==============================================================================!

tpl%nheight=0

do i=1,n
   tpl(i)%icol(:)=0
enddo

!call heap (threadinfo,'tpl(:)%icol(:)','wsmp_setup',wsmp%n*tpl%nheightmax,'int',+1)
 
!==============================================================================!
!=====[find_connectivity_dimension]============================================!
!==============================================================================!
! the list of non-zero terms is computed for every column
! i1 is a column index
! i2 is a row index
! for symmetric systems, we only build and store the lower triangular part 
! of the matrix, i.e. i2.ge.i1

do ie=1,osolve%nleaves
   if (vo%leaf(ie).eq.0) then
      do k1=1,params%mpe
         do idof1=1,ndof
            i1=(vo%ftr(osolve%icon(k1,ie))-1)*ndof+idof1
            do k2=1,params%mpe
               do idof2=1,ndof
                  i2=(vo%ftr(osolve%icon(k2,ie))-1)*ndof+idof2
                  if (i1.gt.i2.and.ndof.gt.1) goto 2
                  do j=1,tpl(i1)%nheight
                     if (i2.eq.tpl(i1)%icol(j)) goto 2
                  enddo
                  if (tpl(i1)%nheight.eq.tpl(i1)%nheightmax) then
                     do k=1,tpl(i1)%nheight
                        icol(k)=tpl(i1)%icol(k)
                     enddo
                     if (tpl(i1)%nheight.gt.0) deallocate (tpl(i1)%icol)
                     tpl(i1)%nheightmax=tpl(i1)%nheightmax+ndof
                     allocate (tpl(i1)%icol(tpl(i1)%nheightmax),stat=err)
                     tpl(i1)%icol(:)=0
                     if (err.ne.0) call stop_run ('Error alloc tpl(i1)%icol$')
                     do k=1,tpl(i1)%nheight
                        tpl(i1)%icol(k)=icol(k)
                     enddo
                  endif
                  tpl(i1)%nheight=tpl(i1)%nheight+1
                  tpl(i1)%icol(tpl(i1)%nheight)=i2
2                 continue
               enddo
            enddo
         enddo
      enddo
   endif
enddo

!print *,iproc,'nface=',osolve%nface

do ie=1,osolve%nface
   if (vo%face(ie).eq.0) then
      do k1=1,9
         do idof1=1,ndof
            i1=(vo%ftr(osolve%iface(k1,ie))-1)*ndof+idof1
            do k2=1,9
               do idof2=1,ndof
                  i2=(vo%ftr(osolve%iface(k2,ie))-1)*ndof+idof2
                  if (i1.gt.i2.and.ndof.gt.1) goto 3             !!!!!!NEW
                  do j=1,tpl(i1)%nheight
                     if (i2.eq.tpl(i1)%icol(j)) goto 3
                  enddo
                  if (tpl(i1)%nheight.eq.tpl(i1)%nheightmax) then
                     do k=1,tpl(i1)%nheight
                        icol(k)=tpl(i1)%icol(k)
                     enddo
                     if (tpl(i1)%nheight.gt.0) deallocate (tpl(i1)%icol)
                     tpl(i1)%nheightmax=tpl(i1)%nheightmax+ndof
                     allocate (tpl(i1)%icol(tpl(i1)%nheightmax),stat=err)
                     if (err.ne.0) call stop_run ('Error alloc tpl(i1)%icol$')
                     do k=1,tpl(i1)%nheight
                        tpl(i1)%icol(k)=icol(k)
                     enddo
                  endif
                  tpl(i1)%nheight=tpl(i1)%nheight+1
                  tpl(i1)%icol(tpl(i1)%nheight)=i2
3                 continue
               enddo
            enddo
         enddo
      enddo
   endif
enddo

!print *,iproc,'min/max tpl%nheight',minval(tpl(:)%nheight),maxval(tpl(:)%nheight)
!print *,iproc,'avrg tpl%nheight',sum(tpl(:)%nheight)/dble(wsmp%n)

!==============================================================================!
!=====[computing total number of non-zero terms]===============================!
!==============================================================================!
! the total number of non-zero terms is the sum of the height of each column

nz=sum(tpl(1:vo%nnode*ndof)%nheight)

if (iproc.eq.0) write(*,'(a,i9)') shift//'n =',n
if (iproc.eq.0) write(*,'(a,i9)') shift//'nz=',nz

!==============================================================================!
!=====[sorting icol's]=========================================================!
!==============================================================================!
! wsmp requires the indices of the row contained in each column to be sorted

do i=1,n
   call iqsort_s (tpl(i)%icol,tpl(i)%nheight)
end do

!==============================================================================!
!=====[find_processors]========================================================!
!==============================================================================!

nn  =(n-1)/nproc+1
n_iproc_st     = nn*iproc+1
n_iproc_end    = min(n_iproc_st+nn-1,n)
n_iproc        = n_iproc_end - n_iproc_st +1

ldb  = n_iproc

allocate(iproc_col(n))
iproc_col=.false.
iproc_col(n_iproc_st:n_iproc_end)=.true.
nz_loc=sum(tpl(:)%nheight,iproc_col)
deallocate(iproc_col) 

!==============================================================================!
!=====[outputs visual reprenstation of matrix]=================================!
!==============================================================================!

if (params%visualise_matrix) then
   allocate(irn(nz))      ; call heap (threadinfo,'irn','wsmp_setup',size(irn),'int',+1)
   allocate(jcn(nz))      ; call heap (threadinfo,'jcn','wsmp_setup',size(jcn),'int',+1)
   irn=0
   jcn=0
   ! the coordinates (irn,jcn) of the non zero terms in the global matrix are 
   ! computed, as well as the indices of the diagonal terms in the list of dof's.
   kk=0
   do i=1,n
      do j=1,tpl(i)%nheight
         kk=kk+1
         if (kk.gt.nz) call stop_run ('error in find_connectivity$')
         irn(kk)=tpl(i)%icol(j)
         jcn(kk)=i
      enddo
   enddo
!   call visualise_matrix (nz,irn,jcn,n,istep,iter,ndof)
   call heap (threadinfo,'irn','wsmp_setup',size(irn),'int',-1)    ; deallocate (irn)
   call heap (threadinfo,'jcn','wsmp_setup',size(jcn),'int',-1)    ; deallocate (jcn)
end if

end subroutine wsmp_setup1




!==============================================================================!
!==============================================================================!
!
!              WSMP_SETUP2      Oct. 2008                                      |
!
!==============================================================================!
!==============================================================================!

subroutine wsmp_setup2 (n,n_iproc,n_iproc_st,n_iproc_end,nz_loc,iproc_col,     &
                        ia,ja,ndof,vo,osolve,tpl,params,threadinfo,istep,iter)

use definitions
use threads
implicit none

integer n,nz_loc,n_iproc,n_iproc_st,n_iproc_end
logical iproc_col(n)
integer ia(n_iproc+1)
integer ja(nz_loc)
integer ndof
type (void) vo
type (topology) :: tpl(n)
type (octreesolve) osolve
type (parameters) params
type (thread) threadinfo
integer istep,iter

integer iloc,iatemp,nheightp,i,k
integer iproc,nproc,ierr

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

iproc_col=.false.
iproc_col(n_iproc_st:n_iproc_end)=.true.

!==============================================================================!
!=====[find_connectivity_local]================================================!
!==============================================================================!

iloc=1

do i=1,n
   if (iproc_col(i)) then
      do k=1,tpl(i)%nheight
         ja(iloc+k-1)=tpl(i)%icol(k)
      enddo
      iloc=iloc+tpl(i)%nheight
   endif
enddo

!==============================================================================!
!=====[build ia,ja]============================================================!
!==============================================================================!

iatemp=0
k=0
do i=n_iproc_st,n_iproc_end
   k=k+1
   if (iatemp==0) then
      ia(k)=1
   else
      ia(k)=iatemp+nheightp
   end if
   iatemp=ia(k)
   nheightp=tpl(i)%nheight
end do 
k=k+1
ia(k)=iatemp+nheightp

end subroutine wsmp_setup2

!==============================================================================!
!==============================================================================!



