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

subroutine solve_with_pwgsmp (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'

!==============================================================================!
!==============================================================================!
! This subroutine solves the system obtained in the temperature calculations
! It relies on the pwgsmp routine, which is the parallel wsmp solver for 
! asymmetric matrices. 

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

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
double precision time1,time2
double precision time_st,time_end
double precision, dimension(:), allocatable :: solution
integer,dimension(MPI_STATUS_SIZE)::statut 
character*72  :: shift
character*35  :: equal

integer iparm(64)
double precision dparm(64)
double precision, dimension(:,:), allocatable :: rmisc

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

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

shift=' ' ; equal='==================================='

allocate (solution(vo%nnode*ndof),stat=threadinfo%err)   ; call heap (threadinfo,'solution','solve_with',size(solution),'dp',+1)
allocate (rmisc(ldb,nrhs),stat=threadinfo%err)            ; call heap (threadinfo,'rmisc','solve_with',size(rmisc),'dp',+1)

!call int_to_char(ciproc,3,iproc)
!open(unit=770+iproc,file='ia1_'//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='ja1_'//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='b1_'//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='matreconstr1_'//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) write(*,'(a)') shift//equal

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

time1=mpi_wtime()

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

iparm(4)=1 !CSC format

call pwgsmp (n_iproc, ia, ja, avals, b, &
             ldb, nrhs, rmisc, iparm, dparm)

time2=mpi_wtime()

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


!==============================================================================!
!=====[WSMP: analyis]==========================================================!
!==============================================================================!

time1=mpi_wtime()

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

call pwgsmp (n_iproc, ia, ja, avals, b, &
             ldb, nrhs, rmisc, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64).ne.0) then
      print *,'WSMP-> analysis pb',iparm(64)
      stop 'WSMP pb'
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: analysis time:         ',time2-time1,' s'
      write(*,'(a,i12,a)')   shift//'wsmp: nb of exp. non-z in factors',iparm(24)
      write(*,'(a,f15.7)')   shift//'wsmp: nb of exp. Mflops in factorisation',dparm(24)*1.d-6
   end if
end if


!==============================================================================!
!=====[WSMP: LU factorisation]=================================================!
!==============================================================================!

time1=mpi_wtime()

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

call pwgsmp (n_iproc, ia, ja, avals, b, &
             ldb, nrhs, rmisc, iparm, dparm)

time2=mpi_wtime()

if (iproc.eq.0) then
   if (iparm(64).ne.0) then
      print *,'WSMP-> LU factorisation pb',iparm(64)
      stop 'WSMP pb' 
   else
      write(*,'(a,f15.7,a)') shift//'wsmp: LU time:               ',time2-time1,' s'
      write(*,'(a,i12,a)')   shift//'wsmp: actual nb of non-z in factors',iparm(23)
      write(*,'(a,f15.7)')   shift//'wsmp: actual nb of Mflops in factorisation',(dparm(23)*1.d-6)/(time2-time1)
   end if
end if


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

time1=mpi_wtime()

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

call pwgsmp (n_iproc, ia, ja, avals, b, &
             ldb, nrhs, rmisc, iparm, dparm)

time2=mpi_wtime()


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


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

!==============================================================================!
!=====[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%temp(1:osolve%nnode)=0.d0
do i=1,vo%nnode
   j=vo%rtf(i)
   ov%temp(j)=solution(i)
   if (.not.vo%influid(j)) ov%temp(j)=0.d0
enddo


if (iproc.eq.0) then
   if (params%debug>=1) then
   write(*,'(a,2e13.5)') shift//'temp :',minval(ov%temp),maxval(ov%temp)
   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%temp ',minval(ov%temp),'maxval ov%temp',maxval(ov%temp)
end if

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


call pwsmp_clear ()

end subroutine solve_with_pwgsmp

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