Newer
Older
!==============================================================================!
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!==============================================================================!
!==============================================================================!
! |
! 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, &
!use mpi
include 'mpif.h'
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
!==============================================================================!
!==============================================================================!
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
integer badnode,cnt,curbadnode,idof,rbadnode
integer,dimension(:),allocatable :: bad_elems
!==============================================================================!
!==============================================================================!
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
allocate(perm(n),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'perm','solve_with',size(perm),'int',+1)
allocate(invp(n),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'invp','solve_with',size(invp),'int',+1)
allocate (solution(vo%nnode*ndof),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'solution','solve_with',size(solution),'dp',+1)
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
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)
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)
if (params%debug.gt.1) 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
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
call write_streepjes(6,1)
write (*,*) 'Cholesky factorization problem at location ',iparm(64),' in WSMP'
call write_streepjes(6,1)
if (params%debug.gt.0) then
write (*,*) 'Debug output:'
write (*,*) 'Bad main diagonal value (avals): ',avals(ia(iparm(64)))
idof=mod(iparm(64),ndof)
if (idof.eq.0) idof=3
rbadnode = ((iparm(64)-idof+1)/(ndof))
badnode = vo%rtf(rbadnode)
allocate (bad_elems(8))
cnt=0
do i=1,osolve%nleaves
do j=1,params%mpe
if (osolve%icon(j,i).eq.badnode) then
cnt=cnt+1
bad_elems(cnt)=i
endif
enddo
enddo
call write_streepjes(6,1)
do i=1,cnt
write (*,*) '-- Elemental output --'
write (*,*) 'Values used in element ',bad_elems(cnt),':'
write (*,*) 'Pressure: ',osolve%pressure(bad_elems(cnt))
write (*,*) 'Smoothed pressure: ',osolve%spressure(bad_elems(cnt))
write (*,*) 'Effective viscosity: ',osolve%eviscosity(bad_elems(cnt))
write (*,*) 'Effective strain rate: ',osolve%e2d(bad_elems(cnt))
do j=1,params%mpe
curbadnode=osolve%icon(j,bad_elems(cnt))
write (*,*) '-- Nodal output --'
write (*,*) 'Values for node ',curbadnode
write (*,*) 'x-coorindate: ',osolve%x(curbadnode)
write (*,*) 'y-coorindate: ',osolve%y(curbadnode)
write (*,*) 'z-coorindate: ',osolve%z(curbadnode)
write (*,*) 'u-velocity: ',osolve%u(curbadnode)
write (*,*) 'v-velocity: ',osolve%v(curbadnode)
write (*,*) 'w-velocity: ',osolve%w(curbadnode)
enddo
enddo
endif
call write_streepjes(6,1)
!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
if (params%debug.gt.1) call heap (threadinfo,'perm','solve_with',size(perm),'int',-1)
deallocate(perm)
if (params%debug.gt.1) call heap (threadinfo,'invp','solve_with',size(invp),'int',-1)
deallocate(invp)
if (params%debug.gt.1) 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
Dave Whipp
committed
! 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_ssend (n_iproc_st,1,mpi_integer,0,tag,mpi_comm_world,ierr)
call mpi_ssend (n_iproc_end,1,mpi_integer,0,tag,mpi_comm_world,ierr)
call mpi_ssend (n_iproc,1,mpi_integer,0,tag,mpi_comm_world,ierr)
call mpi_ssend (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
! 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
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
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
if (params%debug.gt.1) call heap (threadinfo,'solution','solve_with',size(solution),'dp',-1)
deallocate (solution)
!==============================================================================!
!=====[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
!==============================================================================!
!==============================================================================!