!==============================================================================!
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!==============================================================================!
!==============================================================================!
!                                                                              |
!              SOLVE_WITH_WSMP    Oct. 2008                                    |
!                                                                              |
!==============================================================================!
!==============================================================================!

subroutine solve_with_pwssmp (n,nz_loc,n_iproc,n_iproc_st,n_iproc_end,ia,ja,  &
                              ldb,nrhs,avals,b,params,osolve,ov,vo,threadinfo,ndof,istep,iter,iter_nl)

use definitions
use threads

implicit none

INCLUDE 'mpif.h'

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

integer n,nz_loc,n_iproc,n_iproc_st,n_iproc_end
integer ia(n_iproc+1)
integer ja(nz_loc)
integer ldb,nrhs
double precision avals(nz_loc)
double precision b(ldb,nrhs)
type (parameters) params
type (octreesolve) osolve
type (octreev) ov
type (void) vo
type (thread) threadinfo
integer ndof
integer istep,iter,iter_nl

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

integer iproc,nproc,ierr
character*3   :: ciproc
integer i,j,k,iip,ist,iend,iwidth,tag,wsmpdpmem
double precision time1,time2
double precision time_st,time_end
double precision, dimension(:), allocatable :: solution
integer,dimension(MPI_STATUS_SIZE)::statut 
character*72  :: shift

integer iparm(64)
double precision dparm(64)
integer, dimension(:), allocatable :: perm,invp 
double precision aux, diag
integer mrp,naux

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

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

allocate(perm(n),stat=threadinfo%err)                   ; call heap (threadinfo,'perm','solve_with',size(perm),'int',+1)
allocate(invp(n),stat=threadinfo%err)                   ; call heap (threadinfo,'invp','solve_with',size(invp),'int',+1)
allocate (solution(vo%nnode*ndof),stat=threadinfo%err)  ; call heap (threadinfo,'solution','solve_with',size(solution),'dp',+1)

invp=0
perm=0
mrp=0
naux = 0

shift=' '


call int_to_char(ciproc,3,iproc)

!open(unit=770+iproc,file='ia3_'//ciproc//'.dat',status='replace')
!do k=1,n_iproc+1
!   write(770+iproc,*) ia(k)
!end do
!close(770+iproc)
!open(unit=770+iproc,file='ja3_'//ciproc//'.dat',status='replace')
!do k=1,nz_loc
!   write(770+iproc,*) ja(k),avals(k) 
!end do
!close(770+iproc)
!open(unit=770+iproc,file='b3_'//ciproc//'.dat',status='replace')
!do k=1,n_iproc
!   write(770+iproc,*) b(k,1)
!end do
!close(770+iproc)

!open(unit=770+iproc,file='matreconstr3_'//ciproc//'.dat',status='replace')
!do k=1,n_iproc
!   do j=ia(k),ia(k+1)-1
!      write(770+iproc,*) k+n_iproc_st-1,ja(j),avals(j)
!   end do
!end do
!close(770+iproc)

time_st=mpi_wtime()

CALL WSETMAXTHRDS(1)

!CALL PWSMP_INITIALIZE()

if (iproc.eq.0) then
   write(*,'(a,f15.10,a)') shift//'================================='
end if

!==============================================================================!
!=====[WSMP: initialisation]===================================================!
!==============================================================================!

time1=mpi_wtime()

iparm(1)=0
iparm(2)=0
iparm(3)=0

!iparm(20)=2

call pwssmp (n_iproc, ia, ja, avals, diag, perm, invp, b, &
             ldb, nrhs, aux, naux, mrp, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64).ne.0) then
      print *,'WSMP-> initialisation pb',iparm(64)
      call stop_run ('WSMP pb$')
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: initialisation time:   ',time2-time1,' s'
   end if
end if

!==============================================================================!
!=====[WSMP: ordering]=========================================================!
!==============================================================================!

time1=mpi_wtime()

iparm(2)=1
iparm(3)=1

call pwssmp (n_iproc, ia, ja, avals, diag, perm, invp, b, &
             ldb, nrhs, aux, naux, mrp, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64).ne.0) then
      print *,'WSMP-> ordering pb',iparm(64)
      stop 'WSMP pb'
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: ordering time:         ',time2-time1,' s'
   end if
end if

!==============================================================================!
!=====[WSMP: symbolic factorisation]===========================================!
!==============================================================================!

time1=mpi_wtime()

iparm(2)=2
iparm(3)=2

call pwssmp (n_iproc, ia, ja, avals, diag, perm, invp, b, &
             ldb, nrhs, aux, naux, mrp, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64).ne.0) then
      print *,'WSMP-> symbolic factorisation pb',iparm(64)
      stop 
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: symbolic fact. time:   ',time2-time1,' s'
   end if
end if

wsmpdpmem=iparm(23)
call heap (threadinfo,'wsmpdpmem','solve_with',wsmpdpmem,'dp',+1)

!==============================================================================!
!=====[Cholesky Factorization]=================================================!
!==============================================================================!

time1=mpi_wtime()

iparm(2) = 3
iparm(3) = 3

call pwssmp (n_iproc, ia, ja, avals, diag, perm, invp, b, &
             ldb, nrhs, aux, naux, mrp, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64) .ne. 0) then
      print *,'WSMP-> Cholesky factorisation pb',iparm(64)
      stop
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: Cholesky fact. time:   ',time2-time1,' s'
      write(*,'(a,f15.7,a)') shift//'wsmp: Factorisation          ',(dparm(23) * 1.d-6) / (time2-time1),' Mflops'
   end if
end if


!==============================================================================!
!=====[Back substitution]======================================================!
!==============================================================================!

time1=mpi_wtime()

iparm(2) = 4
iparm(3) = 4

call pwssmp (n_iproc, ia, ja, avals, diag, perm, invp, b, &
             ldb, nrhs, aux, naux, mrp, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64) .ne. 0) then
      print *,'WSMP-> back substitution pb',iparm(64)
      return
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: back substitution time:',time2-time1,' s'
   end if
end if

call heap (threadinfo,'perm','solve_with',size(perm),'int',-1)            ; deallocate(perm)
call heap (threadinfo,'invp','solve_with',size(invp),'int',-1)            ; deallocate(invp)
call heap (threadinfo,'wsmpdpmem','solve_with',wsmpdpmem,'dp',-1)

!==============================================================================!
!=====[transfer b's in solution vector]========================================!
!==============================================================================!

time1=mpi_wtime()

tag=111
solution=1.d31

if (iproc.eq.0) solution(1:n_iproc)=b(:,1)

do iip=1,nproc-1
  call mpi_sendrecv(n_iproc_st,1,mpi_integer,0,tag,ist,1,mpi_integer,iip,tag,&
       mpi_comm_world,statut,ierr)
  call mpi_sendrecv(n_iproc_end,1,mpi_integer,0,tag,iend,1,mpi_integer,iip,tag,&
       mpi_comm_world,statut,ierr)
  call mpi_sendrecv(n_iproc,1,mpi_integer,0,tag,iwidth,1,mpi_integer,iip,tag,&
       mpi_comm_world,statut,ierr)
  call mpi_sendrecv(b(:,1),n_iproc,mpi_double_precision,0,tag,&
       solution(ist:iend),iwidth,mpi_double_precision,iip,tag,mpi_comm_world,&
       statut,ierr)
!   if (iproc==iip) then
!      call mpi_send (n_iproc_st,1,mpi_integer,0,tag,mpi_comm_world,ierr)  
!      call mpi_send (n_iproc_end,1,mpi_integer,0,tag,mpi_comm_world,ierr)
!      call mpi_send (n_iproc,1,mpi_integer,0,tag,mpi_comm_world,ierr)
!      call mpi_send (b(:,1),n_iproc,mpi_double_precision,0,tag,mpi_comm_world,ierr)
!   elseif (iproc==0) then
!      call mpi_recv (ist,1,mpi_integer,iip,tag,mpi_comm_world,statut,ierr)
!      call mpi_recv (iend,1,mpi_integer,iip,tag,mpi_comm_world,statut,ierr)
!      call mpi_recv (iwidth,1,mpi_integer,iip,tag,mpi_comm_world,statut,ierr)
!      call mpi_recv (solution(ist:iend),iwidth,mpi_double_precision,iip,tag,mpi_comm_world,statut,ierr)
!   end if
end do
call mpi_bcast (solution,n,mpi_double_precision,0,mpi_comm_world,ierr)

time2=mpi_wtime()

if (iproc.eq.0) then
   write(*,'(a,f15.7,a)') shift//'wsmp: from b to rhs vector:  ',time2-time1,' s'
end if

!write(*,*) iproc,'solution',minval(solution),maxval(solution)

time_end=mpi_wtime()

if (iproc.eq.0) then
   write(*,'(a,f15.7,a)') shift//'---------------------------------'
   write(*,'(a,f15.7,a)') shift//'wsmp: total solve time:      ',time_end-time_st,' s' 
end if

if (iproc.eq.0) then
   write(*,'(a,f15.7,a)') shift//'================================='
end if


!==============================================================================!
!=====[transfer the solution onto ov for next step or iteration]===============!
!==============================================================================!

ov%unode(1:osolve%nnode)=0.d0
ov%vnode(1:osolve%nnode)=0.d0
ov%wnode(1:osolve%nnode)=0.d0
do i=1,vo%nnode
   j=vo%rtf(i)
   ov%unode(j)=solution((i-1)*ndof+1)
   ov%vnode(j)=solution((i-1)*ndof+2)
   ov%wnode(j)=solution((i-1)*ndof+3)
enddo

if (iproc.eq.0) then
   if (params%debug>=1) then
   write(*,'(a,2e13.5)') shift//'u :',minval(ov%unode),maxval(ov%unode)
   write(*,'(a,2e13.5)') shift//'v :',minval(ov%vnode),maxval(ov%vnode)
   write(*,'(a,2e13.5)') shift//'w :',minval(ov%wnode),maxval(ov%wnode)
   write(*,'(a,f15.7,a)') shift//'================================='
   end if
   write (8,*) 'raw solution',minval(solution),maxval(solution)
   write (8,'(a,f30.20,a,f30.20)') 'min ov%unode ',minval(ov%unode),'maxval ov%unode',maxval(ov%unode)
   write (8,'(a,f30.20,a,f30.20)') 'min ov%vnode ',minval(ov%vnode),'maxval ov%vnode',maxval(ov%vnode)
   write (8,'(a,f30.20,a,f30.20)') 'min ov%wnode ',minval(ov%wnode),'maxval ov%wnode',maxval(ov%wnode)
end if

call heap (threadinfo,'solution','solve_with',size(solution),'dp',-1)                   ; deallocate (solution)

call pwsmp_clear()

!==============================================================================!
!=====[DoRuRe: output octree and solver statistics]============================!
!==============================================================================!

call DoRuRe_octree_stats (osolve,params,istep,iter,iter_nl)
!call DoRuRe_solver_stats (params%doDoRuRe,istep,iter,iter_nl,id%n,id%nz,id%rinfog(3),id%infog(3),time1,time2,time3,osolve%octree(2))

end subroutine solve_with_pwssmp


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