mirror of
https://github.com/sfilippone/amg4psblas.git
synced 2026-10-06 22:55:12 +00:00
Merge branch 'development' into MixedI8
This commit is contained in:
@@ -82,25 +82,37 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
! Running in local mode?
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
gx = x
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='complex(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='complex(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
else
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='complex(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
sv%id%icntl(9) = 1
|
||||
@@ -120,10 +132,13 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
sv%id%icntl(4)=-1
|
||||
sv%id%job = 3
|
||||
call cmumps(sv%id)
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_geaxpby(alpha,gx,beta,y,desc_data,info)
|
||||
else
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
end if
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
@@ -132,9 +147,7 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (nglob > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
if (allocated(ww)) deallocate(ww)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -58,135 +58,145 @@
|
||||
! Local variables
|
||||
type(psb_cspmat_type) :: atmp
|
||||
type(psb_c_coo_sparse_mat), target :: acoo
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='c_mumps_solver_bld', ch_err
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, iam, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='c_mumps_solver_bld', ch_err
|
||||
|
||||
#if defined(HAVE_MUMPS_) && !defined(LPK8)
|
||||
|
||||
info=psb_success_
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) < 0 ) then
|
||||
call psb_info(ictxt, me, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/me/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt1
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, me, np)
|
||||
npr = np
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_info(ictxt, iam, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/iam/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,sv%local_ictxt
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, iam, np)
|
||||
me = iam
|
||||
npr = np
|
||||
else
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_cmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call cmumps(sv%id)
|
||||
!WARNING: CALLING cMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) < 0) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) >= 0 ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)'calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call cmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_cmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_dmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
end if
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call cmumps(sv%id)
|
||||
!WARNING: CALLING dMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) then
|
||||
!write(0,*) 'MUMPS_BLD: setting entry ',i,' to ', sv%icntl(i)%item
|
||||
sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. '
|
||||
write(*,*)'A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)iam, ' calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call dmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_dmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) iam,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
#else
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
#endif
|
||||
end subroutine c_mumps_solver_bld
|
||||
end subroutine c_mumps_solver_bld
|
||||
|
||||
|
||||
@@ -82,25 +82,37 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
! Running in local mode?
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
gx = x
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
else
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
sv%id%icntl(9) = 1
|
||||
@@ -120,10 +132,13 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
sv%id%icntl(4)=-1
|
||||
sv%id%job = 3
|
||||
call dmumps(sv%id)
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_geaxpby(alpha,gx,beta,y,desc_data,info)
|
||||
else
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
end if
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
@@ -132,9 +147,7 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (nglob > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
if (allocated(ww)) deallocate(ww)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -58,135 +58,145 @@
|
||||
! Local variables
|
||||
type(psb_dspmat_type) :: atmp
|
||||
type(psb_d_coo_sparse_mat), target :: acoo
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='d_mumps_solver_bld', ch_err
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, iam, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='d_mumps_solver_bld', ch_err
|
||||
|
||||
#if defined(HAVE_MUMPS_) && !defined(LPK8)
|
||||
|
||||
info=psb_success_
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) < 0 ) then
|
||||
call psb_info(ictxt, me, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/me/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt1
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, me, np)
|
||||
npr = np
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_info(ictxt, iam, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/iam/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,sv%local_ictxt
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, iam, np)
|
||||
me = iam
|
||||
npr = np
|
||||
else
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_dmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call dmumps(sv%id)
|
||||
!WARNING: CALLING dMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) < 0) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) >= 0 ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)'calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call dmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_dmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_dmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
end if
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call dmumps(sv%id)
|
||||
!WARNING: CALLING dMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) then
|
||||
!write(0,*) 'MUMPS_BLD: setting entry ',i,' to ', sv%icntl(i)%item
|
||||
sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. '
|
||||
write(*,*)'A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)iam, ' calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call dmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_dmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) iam,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
#else
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
#endif
|
||||
end subroutine d_mumps_solver_bld
|
||||
end subroutine d_mumps_solver_bld
|
||||
|
||||
|
||||
@@ -82,25 +82,37 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
! Running in local mode?
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
gx = x
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='real(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
else
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='real(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
sv%id%icntl(9) = 1
|
||||
@@ -120,10 +132,13 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
sv%id%icntl(4)=-1
|
||||
sv%id%job = 3
|
||||
call smumps(sv%id)
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_geaxpby(alpha,gx,beta,y,desc_data,info)
|
||||
else
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
end if
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
@@ -132,9 +147,7 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (nglob > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
if (allocated(ww)) deallocate(ww)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -58,135 +58,145 @@
|
||||
! Local variables
|
||||
type(psb_sspmat_type) :: atmp
|
||||
type(psb_s_coo_sparse_mat), target :: acoo
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_mumps_solver_bld', ch_err
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, iam, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_mumps_solver_bld', ch_err
|
||||
|
||||
#if defined(HAVE_MUMPS_) && !defined(LPK8)
|
||||
|
||||
info=psb_success_
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) < 0 ) then
|
||||
call psb_info(ictxt, me, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/me/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt1
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, me, np)
|
||||
npr = np
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_info(ictxt, iam, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/iam/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,sv%local_ictxt
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, iam, np)
|
||||
me = iam
|
||||
npr = np
|
||||
else
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_smumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call smumps(sv%id)
|
||||
!WARNING: CALLING sMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) < 0) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) >= 0 ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)'calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call smumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_smumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_dmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
end if
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call smumps(sv%id)
|
||||
!WARNING: CALLING dMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) then
|
||||
!write(0,*) 'MUMPS_BLD: setting entry ',i,' to ', sv%icntl(i)%item
|
||||
sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. '
|
||||
write(*,*)'A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)iam, ' calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call dmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_dmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) iam,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
#else
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
#endif
|
||||
end subroutine s_mumps_solver_bld
|
||||
end subroutine s_mumps_solver_bld
|
||||
|
||||
|
||||
@@ -82,25 +82,37 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
! Running in local mode?
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
gx = x
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww = work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='complex(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='complex(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
else
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
allocate(gx(nglob),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/nglob,0,0,0,0/),&
|
||||
& a_err='complex(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_gather(gx, x, desc_data, info, root=0)
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
sv%id%icntl(9) = 1
|
||||
@@ -120,10 +132,13 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
sv%id%icntl(4)=-1
|
||||
sv%id%job = 3
|
||||
call zmumps(sv%id)
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_geaxpby(alpha,gx,beta,y,desc_data,info)
|
||||
else
|
||||
call psb_scatter(gx, ww, desc_data, info, root=0)
|
||||
if (info == psb_success_) then
|
||||
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
end if
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
@@ -132,9 +147,7 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (nglob > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
if (allocated(ww)) deallocate(ww)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -58,135 +58,145 @@
|
||||
! Local variables
|
||||
type(psb_zspmat_type) :: atmp
|
||||
type(psb_z_coo_sparse_mat), target :: acoo
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='z_mumps_solver_bld', ch_err
|
||||
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nglobrec, nzt, npr, npc
|
||||
integer(psb_ipk_) :: ifrst, ibcheck
|
||||
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, iam, me, i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='z_mumps_solver_bld', ch_err
|
||||
|
||||
#if defined(HAVE_MUMPS_) && !defined(LPK8)
|
||||
|
||||
info=psb_success_
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) < 0 ) then
|
||||
call psb_info(ictxt, me, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/me/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt1
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, me, np)
|
||||
npr = np
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = desc_a%get_context()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
call psb_info(ictxt, iam, np)
|
||||
call psb_init(ictxt1,np=1,basectxt=ictxt,ids=(/iam/))
|
||||
call psb_get_mpicomm(ictxt1, icomm)
|
||||
allocate(sv%local_ictxt,stat=info)
|
||||
sv%local_ictxt = ictxt1
|
||||
write(*,*)'mumps_bld: +++++>',icomm,sv%local_ictxt
|
||||
call psb_info(ictxt1, me, np)
|
||||
npr = np
|
||||
else if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_get_mpicomm(ictxt,icomm)
|
||||
write(*,*)'mumps_bld: +++++>',icomm,ictxt
|
||||
call psb_info(ictxt, iam, np)
|
||||
me = iam
|
||||
npr = np
|
||||
else
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='Invalid local/global solver in MUMPS')
|
||||
goto 9999
|
||||
end if
|
||||
npc = 1
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
! if (allocated(sv%id)) then
|
||||
! call sv%free(info)
|
||||
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_zmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call zmumps(sv%id)
|
||||
!WARNING: CALLING zMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) < 0) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) >= 0 ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)'calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call zmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
! deallocate(sv%id)
|
||||
! end if
|
||||
if(.not.allocated(sv%id)) then
|
||||
allocate(sv%id,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_zmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name,a_err='mld_dmumps_default')
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
end if
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
sv%id%comm = icomm
|
||||
sv%id%job = -1
|
||||
sv%id%par = 1
|
||||
call zmumps(sv%id)
|
||||
!WARNING: CALLING dMUMPS WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
|
||||
if (allocated(sv%icntl)) then
|
||||
do i=1,mld_mumps_icntl_size
|
||||
if (allocated(sv%icntl(i)%item)) then
|
||||
!write(0,*) 'MUMPS_BLD: setting entry ',i,' to ', sv%icntl(i)%item
|
||||
sv%id%icntl(i) = sv%icntl(i)%item
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if (allocated(sv%rcntl)) then
|
||||
do i=1,mld_mumps_rcntl_size
|
||||
if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
|
||||
end do
|
||||
end if
|
||||
sv%id%icntl(3)=sv%ipar(2)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
nglob = desc_a%get_global_rows()
|
||||
if (sv%ipar(1) == mld_local_solver_ ) then
|
||||
nglobrec=desc_a%get_local_rows()
|
||||
call a%csclip(c,info,jmax=a%get_nrows())
|
||||
call c%cp_to(acoo)
|
||||
nglob = c%get_nrows()
|
||||
if (nglobrec /= nglob) then
|
||||
write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. '
|
||||
write(*,*)'A zero-overlap is used instead'
|
||||
end if
|
||||
else
|
||||
call a%cp_to(acoo)
|
||||
end if
|
||||
nztota = acoo%get_nzeros()
|
||||
|
||||
! switch to global numbering
|
||||
if (sv%ipar(1) == mld_global_solver_ ) then
|
||||
call psb_loc_to_glob(acoo%ja(1:nztota), desc_a, info, iact='I')
|
||||
call psb_loc_to_glob(acoo%ia(1:nztota), desc_a, info, iact='I')
|
||||
end if
|
||||
sv%id%irn_loc => acoo%ia
|
||||
sv%id%jcn_loc => acoo%ja
|
||||
sv%id%a_loc => acoo%val
|
||||
sv%id%icntl(18) = 3
|
||||
if(acoo%is_upper() .or. acoo%is_lower()) then
|
||||
sv%id%sym = 2
|
||||
else
|
||||
sv%id%sym = 0
|
||||
end if
|
||||
sv%id%n = nglob
|
||||
! there should be a better way for this
|
||||
sv%id%nz_loc = acoo%get_nzeros()
|
||||
sv%id%nz = acoo%get_nzeros()
|
||||
sv%id%job = 4
|
||||
!call psb_barrier(ictxt)
|
||||
write(*,*)iam, ' calling mumps N,nz,nz_loc',sv%id%n,sv%id%nz,sv%id%nz_loc
|
||||
call dmumps(sv%id)
|
||||
!call psb_barrier(ictxt)
|
||||
info = sv%id%infog(1)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_dmumps_fact '
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
nullify(sv%id%irn)
|
||||
nullify(sv%id%jcn)
|
||||
nullify(sv%id%a)
|
||||
|
||||
call acoo%free()
|
||||
sv%built=.true.
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) iam,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
#else
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
write(psb_err_unit,*) "MUMPS Not Configured, fix make.inc and recompile "
|
||||
#endif
|
||||
end subroutine z_mumps_solver_bld
|
||||
end subroutine z_mumps_solver_bld
|
||||
|
||||
|
||||
@@ -286,15 +286,15 @@ module mld_base_prec_type
|
||||
integer(psb_ipk_), parameter :: mld_coarse_iluthrs_ = 4
|
||||
integer(psb_ipk_), parameter :: mld_solver_eps_ = 6
|
||||
integer(psb_ipk_), parameter :: mld_rfpsz_ = 8
|
||||
|
||||
!
|
||||
! Is the current solver local or global
|
||||
!
|
||||
integer(psb_ipk_), parameter :: mld_local_solver_ = 0
|
||||
integer(psb_ipk_), parameter :: mld_global_solver_ = 1
|
||||
|
||||
!
|
||||
! Entries for mumps
|
||||
!
|
||||
!parameter controling the sequential/parallel building of MUMPS
|
||||
integer(psb_ipk_), parameter :: mld_as_sequential_ = 40
|
||||
!parameter regulating the error printing of MUMPS
|
||||
integer(psb_ipk_), parameter :: mld_mumps_print_err_ = 41
|
||||
! Size of the control vectors
|
||||
integer, parameter :: mld_mumps_icntl_size=40
|
||||
integer, parameter :: mld_mumps_rcntl_size=15
|
||||
@@ -483,6 +483,10 @@ contains
|
||||
val = mld_no_filter_mat_
|
||||
case('OUTER_SWEEPS')
|
||||
val = mld_outer_sweeps_
|
||||
case('LOCAL_SOLVER')
|
||||
val = mld_local_solver_
|
||||
case('GLOBAL_SOLVER')
|
||||
val = mld_global_solver_
|
||||
case default
|
||||
val = -1
|
||||
end select
|
||||
|
||||
@@ -291,7 +291,7 @@ contains
|
||||
|
||||
select case(psb_toupper(what))
|
||||
#if defined(HAVE_MUMPS_)
|
||||
case('MUMPS_AS_SEQUENTIAL')
|
||||
case('MUMPS_LOC_GLOB')
|
||||
sv%ipar(1)=val
|
||||
case('MUMPS_PRINT_ERR')
|
||||
sv%ipar(2)=val
|
||||
@@ -400,8 +400,9 @@ contains
|
||||
! INSTANTIATION OF sv%id needed to set parmater but mpi communicator needed
|
||||
! sv%id%job = -1
|
||||
! sv%id%par=1
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar(1)=2
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar = 0
|
||||
sv%ipar(1) = mld_global_solver_
|
||||
!sv%ipar(10)=6
|
||||
!sv%ipar(11)=0
|
||||
!sv%ipar(12)=6
|
||||
|
||||
@@ -291,7 +291,7 @@ contains
|
||||
|
||||
select case(psb_toupper(what))
|
||||
#if defined(HAVE_MUMPS_)
|
||||
case('MUMPS_AS_SEQUENTIAL')
|
||||
case('MUMPS_LOC_GLOB')
|
||||
sv%ipar(1)=val
|
||||
case('MUMPS_PRINT_ERR')
|
||||
sv%ipar(2)=val
|
||||
@@ -400,8 +400,9 @@ contains
|
||||
! INSTANTIATION OF sv%id needed to set parmater but mpi communicator needed
|
||||
! sv%id%job = -1
|
||||
! sv%id%par=1
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar(1)=2
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar = 0
|
||||
sv%ipar(1) = mld_global_solver_
|
||||
!sv%ipar(10)=6
|
||||
!sv%ipar(11)=0
|
||||
!sv%ipar(12)=6
|
||||
|
||||
@@ -291,7 +291,7 @@ contains
|
||||
|
||||
select case(psb_toupper(what))
|
||||
#if defined(HAVE_MUMPS_)
|
||||
case('MUMPS_AS_SEQUENTIAL')
|
||||
case('MUMPS_LOC_GLOB')
|
||||
sv%ipar(1)=val
|
||||
case('MUMPS_PRINT_ERR')
|
||||
sv%ipar(2)=val
|
||||
@@ -400,8 +400,9 @@ contains
|
||||
! INSTANTIATION OF sv%id needed to set parmater but mpi communicator needed
|
||||
! sv%id%job = -1
|
||||
! sv%id%par=1
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar(1)=2
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar = 0
|
||||
sv%ipar(1) = mld_global_solver_
|
||||
!sv%ipar(10)=6
|
||||
!sv%ipar(11)=0
|
||||
!sv%ipar(12)=6
|
||||
|
||||
@@ -291,7 +291,7 @@ contains
|
||||
|
||||
select case(psb_toupper(what))
|
||||
#if defined(HAVE_MUMPS_)
|
||||
case('MUMPS_AS_SEQUENTIAL')
|
||||
case('MUMPS_LOC_GLOB')
|
||||
sv%ipar(1)=val
|
||||
case('MUMPS_PRINT_ERR')
|
||||
sv%ipar(2)=val
|
||||
@@ -400,8 +400,9 @@ contains
|
||||
! INSTANTIATION OF sv%id needed to set parmater but mpi communicator needed
|
||||
! sv%id%job = -1
|
||||
! sv%id%par=1
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar(1)=2
|
||||
! call dmumps(sv%id)
|
||||
sv%ipar = 0
|
||||
sv%ipar(1) = mld_global_solver_
|
||||
!sv%ipar(10)=6
|
||||
!sv%ipar(11)=0
|
||||
!sv%ipar(12)=6
|
||||
|
||||
Reference in New Issue
Block a user