Enable MUMPS at LPK8, by checking if actual size fits in IPK4

stopcriterion
Salvatore Filippone 5 years ago
parent f59c816914
commit 6ef5ea9748

@ -55,7 +55,8 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
character, intent(in), optional :: init
complex(psb_spk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row, n_col, nglob
integer(psb_ipk_) :: n_row, n_col
integer(psb_lpk_) :: nglob
complex(psb_spk_), allocatable :: ww(:)
complex(psb_spk_), allocatable, target :: gx(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@ -64,7 +65,7 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
call psb_erractionsave(err_act)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
@ -93,7 +94,7 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,i_err=(/n_col/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
@ -101,7 +102,7 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,e_err=(/nglob/),&
& a_err='complex(psb_spk_)')
goto 9999
end if

@ -59,7 +59,7 @@ subroutine c_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
integer(psb_ipk_) :: err_act
character(len=20) :: name='c_mumps_solver_apply_vect'
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
call psb_erractionsave(err_act)

@ -39,7 +39,7 @@
! Ambra Abdullahi Hassan
!
!
subroutine c_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
subroutine c_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_c_mumps_solver
@ -58,12 +58,16 @@
! 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
#if defined(IPK4) && defined(LPK8)
integer(psb_lpk_), allocatable :: gia(:), gja(:)
#endif
integer(psb_ipk_) :: n_row,n_col, nrow_a, nza, npr, npc
integer(psb_lpk_) :: nglob, nglobrec, nzt
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)
#if defined(HAVE_MUMPS_)
info=psb_success_
@ -113,6 +117,12 @@
sv%id%comm = icomm
sv%id%job = -1
sv%id%par = 1
if (sv%ipar(3) == 2) then
sv%id%sym = 2
else
sv%id%sym = 0
end if
call cmumps(sv%id)
!WARNING: CALLING cmumps WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
if (allocated(sv%icntl)) then
@ -133,7 +143,13 @@
nglob = desc_a%get_global_rows()
if (sv%ipar(1) == mld_local_solver_ ) then
nglobrec=desc_a%get_local_rows()
if (sv%ipar(3) == 2) then
! Always pass the upper triangle to MUMPS
call a%triu(c,info,jmax=a%get_nrows())
call c%set_symmetric()
else
call a%csclip(c,info,jmax=a%get_nrows())
end if
call c%cp_to(acoo)
nglob = c%get_nrows()
if (nglobrec /= nglob) then
@ -143,22 +159,63 @@
else
call a%cp_to(acoo)
end if
nztota = acoo%get_nzeros()
nza = 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')
#if defined(IPK4) && defined(LPK8)
!
! Strategy here is as follows: because a call to MUMPS
! as a gobal solver is mostly done at the coarsest level,
! even if we start from a problem requiring 8 bytes, chances
! are that the global size will be suitable for 4 bytes
! anyway, so we hope for the best, and throw an error
! if something goes wrong.
!
if (nglob > huge(1_psb_ipk_)) then
write(0,*) iam,' ',trim(name),': Error: overflow of local indices '
info=psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
gia = acoo%ia(1:nza)
gja = acoo%ja(1:nza)
call psb_loc_to_glob(gia(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(gja(1:nza), desc_a, info, iact='I')
acoo%ia(1:nza) = gia(1:nza)
acoo%ja(1:nza) = gja(1:nza)
#else
!
! Here global and local numbers have the same size, so this must work.
!
call psb_loc_to_glob(acoo%ja(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(acoo%ia(1:nza), desc_a, info, iact='I')
#endif
if (sv%ipar(3) == 2 ) then
! Always pass the upper triangle to MUMPS
block
integer(psb_ipk_) :: j,nz
nz = 0
do j=1,nza
if (acoo%ja(j) >= acoo%ia(j)) then
nz = nz + 1
acoo%ia(nz) = acoo%ia(j)
acoo%ja(nz) = acoo%ja(j)
acoo%val(nz) = acoo%val(j)
end if
end do
call acoo%set_nzeros(nz)
call acoo%set_triangle()
call acoo%set_upper()
call acoo%set_symmetric()
end block
end if
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%nnz_loc = acoo%get_nzeros()
@ -174,7 +231,7 @@
info = sv%id%infog(1)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='mld_dmumps_fact '
ch_err='mld_cmumps_fact '
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if

@ -55,7 +55,8 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
character, intent(in), optional :: init
real(psb_dpk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row, n_col, nglob
integer(psb_ipk_) :: n_row, n_col
integer(psb_lpk_) :: nglob
real(psb_dpk_), allocatable :: ww(:)
real(psb_dpk_), allocatable, target :: gx(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@ -64,7 +65,7 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
call psb_erractionsave(err_act)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
@ -93,7 +94,7 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,i_err=(/n_col/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
@ -101,7 +102,7 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,e_err=(/nglob/),&
& a_err='real(psb_dpk_)')
goto 9999
end if

@ -59,7 +59,7 @@ subroutine d_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
integer(psb_ipk_) :: err_act
character(len=20) :: name='d_mumps_solver_apply_vect'
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
call psb_erractionsave(err_act)

@ -39,7 +39,7 @@
! Ambra Abdullahi Hassan
!
!
subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_d_mumps_solver
@ -58,12 +58,16 @@
! 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
#if defined(IPK4) && defined(LPK8)
integer(psb_lpk_), allocatable :: gia(:), gja(:)
#endif
integer(psb_ipk_) :: n_row,n_col, nrow_a, nza, npr, npc
integer(psb_lpk_) :: nglob, nglobrec, nzt
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)
#if defined(HAVE_MUMPS_)
info=psb_success_
@ -113,6 +117,12 @@
sv%id%comm = icomm
sv%id%job = -1
sv%id%par = 1
if (sv%ipar(3) == 2) then
sv%id%sym = 2
else
sv%id%sym = 0
end if
call dmumps(sv%id)
!WARNING: CALLING dmumps WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
if (allocated(sv%icntl)) then
@ -133,7 +143,13 @@
nglob = desc_a%get_global_rows()
if (sv%ipar(1) == mld_local_solver_ ) then
nglobrec=desc_a%get_local_rows()
if (sv%ipar(3) == 2) then
! Always pass the upper triangle to MUMPS
call a%triu(c,info,jmax=a%get_nrows())
call c%set_symmetric()
else
call a%csclip(c,info,jmax=a%get_nrows())
end if
call c%cp_to(acoo)
nglob = c%get_nrows()
if (nglobrec /= nglob) then
@ -143,22 +159,63 @@
else
call a%cp_to(acoo)
end if
nztota = acoo%get_nzeros()
nza = 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')
#if defined(IPK4) && defined(LPK8)
!
! Strategy here is as follows: because a call to MUMPS
! as a gobal solver is mostly done at the coarsest level,
! even if we start from a problem requiring 8 bytes, chances
! are that the global size will be suitable for 4 bytes
! anyway, so we hope for the best, and throw an error
! if something goes wrong.
!
if (nglob > huge(1_psb_ipk_)) then
write(0,*) iam,' ',trim(name),': Error: overflow of local indices '
info=psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
gia = acoo%ia(1:nza)
gja = acoo%ja(1:nza)
call psb_loc_to_glob(gia(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(gja(1:nza), desc_a, info, iact='I')
acoo%ia(1:nza) = gia(1:nza)
acoo%ja(1:nza) = gja(1:nza)
#else
!
! Here global and local numbers have the same size, so this must work.
!
call psb_loc_to_glob(acoo%ja(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(acoo%ia(1:nza), desc_a, info, iact='I')
#endif
if (sv%ipar(3) == 2 ) then
! Always pass the upper triangle to MUMPS
block
integer(psb_ipk_) :: j,nz
nz = 0
do j=1,nza
if (acoo%ja(j) >= acoo%ia(j)) then
nz = nz + 1
acoo%ia(nz) = acoo%ia(j)
acoo%ja(nz) = acoo%ja(j)
acoo%val(nz) = acoo%val(j)
end if
end do
call acoo%set_nzeros(nz)
call acoo%set_triangle()
call acoo%set_upper()
call acoo%set_symmetric()
end block
end if
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%nnz_loc = acoo%get_nzeros()

@ -55,7 +55,8 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
character, intent(in), optional :: init
real(psb_spk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row, n_col, nglob
integer(psb_ipk_) :: n_row, n_col
integer(psb_lpk_) :: nglob
real(psb_spk_), allocatable :: ww(:)
real(psb_spk_), allocatable, target :: gx(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@ -64,7 +65,7 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
call psb_erractionsave(err_act)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
@ -93,7 +94,7 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,i_err=(/n_col/),&
& a_err='real(psb_spk_)')
goto 9999
end if
@ -101,7 +102,7 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,e_err=(/nglob/),&
& a_err='real(psb_spk_)')
goto 9999
end if

@ -59,7 +59,7 @@ subroutine s_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
integer(psb_ipk_) :: err_act
character(len=20) :: name='s_mumps_solver_apply_vect'
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
call psb_erractionsave(err_act)

@ -39,7 +39,7 @@
! Ambra Abdullahi Hassan
!
!
subroutine s_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
subroutine s_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_s_mumps_solver
@ -58,12 +58,16 @@
! 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
#if defined(IPK4) && defined(LPK8)
integer(psb_lpk_), allocatable :: gia(:), gja(:)
#endif
integer(psb_ipk_) :: n_row,n_col, nrow_a, nza, npr, npc
integer(psb_lpk_) :: nglob, nglobrec, nzt
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)
#if defined(HAVE_MUMPS_)
info=psb_success_
@ -113,6 +117,12 @@
sv%id%comm = icomm
sv%id%job = -1
sv%id%par = 1
if (sv%ipar(3) == 2) then
sv%id%sym = 2
else
sv%id%sym = 0
end if
call smumps(sv%id)
!WARNING: CALLING smumps WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
if (allocated(sv%icntl)) then
@ -133,7 +143,13 @@
nglob = desc_a%get_global_rows()
if (sv%ipar(1) == mld_local_solver_ ) then
nglobrec=desc_a%get_local_rows()
if (sv%ipar(3) == 2) then
! Always pass the upper triangle to MUMPS
call a%triu(c,info,jmax=a%get_nrows())
call c%set_symmetric()
else
call a%csclip(c,info,jmax=a%get_nrows())
end if
call c%cp_to(acoo)
nglob = c%get_nrows()
if (nglobrec /= nglob) then
@ -143,22 +159,63 @@
else
call a%cp_to(acoo)
end if
nztota = acoo%get_nzeros()
nza = 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')
#if defined(IPK4) && defined(LPK8)
!
! Strategy here is as follows: because a call to MUMPS
! as a gobal solver is mostly done at the coarsest level,
! even if we start from a problem requiring 8 bytes, chances
! are that the global size will be suitable for 4 bytes
! anyway, so we hope for the best, and throw an error
! if something goes wrong.
!
if (nglob > huge(1_psb_ipk_)) then
write(0,*) iam,' ',trim(name),': Error: overflow of local indices '
info=psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
gia = acoo%ia(1:nza)
gja = acoo%ja(1:nza)
call psb_loc_to_glob(gia(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(gja(1:nza), desc_a, info, iact='I')
acoo%ia(1:nza) = gia(1:nza)
acoo%ja(1:nza) = gja(1:nza)
#else
!
! Here global and local numbers have the same size, so this must work.
!
call psb_loc_to_glob(acoo%ja(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(acoo%ia(1:nza), desc_a, info, iact='I')
#endif
if (sv%ipar(3) == 2 ) then
! Always pass the upper triangle to MUMPS
block
integer(psb_ipk_) :: j,nz
nz = 0
do j=1,nza
if (acoo%ja(j) >= acoo%ia(j)) then
nz = nz + 1
acoo%ia(nz) = acoo%ia(j)
acoo%ja(nz) = acoo%ja(j)
acoo%val(nz) = acoo%val(j)
end if
end do
call acoo%set_nzeros(nz)
call acoo%set_triangle()
call acoo%set_upper()
call acoo%set_symmetric()
end block
end if
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%nnz_loc = acoo%get_nzeros()
@ -174,7 +231,7 @@
info = sv%id%infog(1)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='mld_dmumps_fact '
ch_err='mld_smumps_fact '
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if

@ -55,7 +55,8 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
character, intent(in), optional :: init
complex(psb_dpk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row, n_col, nglob
integer(psb_ipk_) :: n_row, n_col
integer(psb_lpk_) :: nglob
complex(psb_dpk_), allocatable :: ww(:)
complex(psb_dpk_), allocatable, target :: gx(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@ -64,7 +65,7 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
call psb_erractionsave(err_act)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
@ -93,7 +94,7 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,i_err=(/n_col/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
@ -101,7 +102,7 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
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/),&
call psb_errpush(info,name,e_err=(/nglob/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if

@ -59,7 +59,7 @@ subroutine z_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
integer(psb_ipk_) :: err_act
character(len=20) :: name='z_mumps_solver_apply_vect'
#if defined(HAVE_MUMPS_) && !defined(LPK8)
#if defined(HAVE_MUMPS_)
call psb_erractionsave(err_act)

@ -39,7 +39,7 @@
! Ambra Abdullahi Hassan
!
!
subroutine z_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
subroutine z_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_z_mumps_solver
@ -58,12 +58,16 @@
! 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
#if defined(IPK4) && defined(LPK8)
integer(psb_lpk_), allocatable :: gia(:), gja(:)
#endif
integer(psb_ipk_) :: n_row,n_col, nrow_a, nza, npr, npc
integer(psb_lpk_) :: nglob, nglobrec, nzt
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)
#if defined(HAVE_MUMPS_)
info=psb_success_
@ -113,6 +117,12 @@
sv%id%comm = icomm
sv%id%job = -1
sv%id%par = 1
if (sv%ipar(3) == 2) then
sv%id%sym = 2
else
sv%id%sym = 0
end if
call zmumps(sv%id)
!WARNING: CALLING zmumps WITH JOB=-1 DESTROY THE SETTING OF DEFAULT:TO FIX
if (allocated(sv%icntl)) then
@ -133,7 +143,13 @@
nglob = desc_a%get_global_rows()
if (sv%ipar(1) == mld_local_solver_ ) then
nglobrec=desc_a%get_local_rows()
if (sv%ipar(3) == 2) then
! Always pass the upper triangle to MUMPS
call a%triu(c,info,jmax=a%get_nrows())
call c%set_symmetric()
else
call a%csclip(c,info,jmax=a%get_nrows())
end if
call c%cp_to(acoo)
nglob = c%get_nrows()
if (nglobrec /= nglob) then
@ -143,22 +159,63 @@
else
call a%cp_to(acoo)
end if
nztota = acoo%get_nzeros()
nza = 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')
#if defined(IPK4) && defined(LPK8)
!
! Strategy here is as follows: because a call to MUMPS
! as a gobal solver is mostly done at the coarsest level,
! even if we start from a problem requiring 8 bytes, chances
! are that the global size will be suitable for 4 bytes
! anyway, so we hope for the best, and throw an error
! if something goes wrong.
!
if (nglob > huge(1_psb_ipk_)) then
write(0,*) iam,' ',trim(name),': Error: overflow of local indices '
info=psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
gia = acoo%ia(1:nza)
gja = acoo%ja(1:nza)
call psb_loc_to_glob(gia(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(gja(1:nza), desc_a, info, iact='I')
acoo%ia(1:nza) = gia(1:nza)
acoo%ja(1:nza) = gja(1:nza)
#else
!
! Here global and local numbers have the same size, so this must work.
!
call psb_loc_to_glob(acoo%ja(1:nza), desc_a, info, iact='I')
call psb_loc_to_glob(acoo%ia(1:nza), desc_a, info, iact='I')
#endif
if (sv%ipar(3) == 2 ) then
! Always pass the upper triangle to MUMPS
block
integer(psb_ipk_) :: j,nz
nz = 0
do j=1,nza
if (acoo%ja(j) >= acoo%ia(j)) then
nz = nz + 1
acoo%ia(nz) = acoo%ia(j)
acoo%ja(nz) = acoo%ja(j)
acoo%val(nz) = acoo%val(j)
end if
end do
call acoo%set_nzeros(nz)
call acoo%set_triangle()
call acoo%set_upper()
call acoo%set_symmetric()
end block
end if
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%nnz_loc = acoo%get_nzeros()
@ -174,7 +231,7 @@
info = sv%id%infog(1)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='mld_dmumps_fact '
ch_err='mld_zmumps_fact '
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if

@ -109,12 +109,13 @@ module mld_c_base_solver_mod
procedure, nopass :: get_fmt => c_base_solver_get_fmt
procedure, nopass :: get_id => c_base_solver_get_id
procedure, nopass :: is_iterative => c_base_solver_is_iterative
procedure, pass(sv) :: is_global => c_base_solver_is_global
end type mld_c_base_solver_type
private :: c_base_solver_sizeof, c_base_solver_default,&
& c_base_solver_get_nzeros, c_base_solver_get_fmt, &
& c_base_solver_is_iterative, c_base_solver_get_id, &
& c_base_solver_get_wrksize
& c_base_solver_get_wrksize, c_base_solver_is_global
interface
@ -365,6 +366,17 @@ contains
val = .false.
end function c_base_solver_is_iterative
!
! Is the solver acting globally? In most cases
! not, SuperLU_Dist does, MUMPS can do either.
!
function c_base_solver_is_global(sv) result(val)
implicit none
class(mld_c_base_solver_type), intent(in) :: sv
logical :: val
val = .false.
end function c_base_solver_is_global
function c_base_solver_get_id() result(val)
implicit none

@ -1,3 +1,4 @@
!
!
! MLD2P4 version 2.2
@ -58,12 +59,6 @@ module mld_c_mumps_solver
include 'cmumps_struc.h'
#endif
#if defined(LPK8)
type, extends(mld_c_base_solver_type) :: mld_c_mumps_solver_type
end type mld_c_mumps_solver_type
#else
type :: mld_c_mumps_icntl_item
integer(psb_ipk_), allocatable :: item
@ -80,7 +75,13 @@ module mld_c_mumps_solver
#endif
type(mld_c_mumps_icntl_item), allocatable :: icntl(:)
type(mld_c_mumps_rcntl_item), allocatable :: rcntl(:)
integer(psb_ipk_), dimension(2) :: ipar
!
! Controls to be set before MUMPS instantiation:
!
! IPAR(1) : MUMPS_LOC_GLOB 0==mld_local_solver_: LOCAL 1==mld_global_solver_: GLOBAL
! IPAR(2) : MUMPS_PRINT_ERR print verbosity (see MUMPS)
! IPAR(3) : MUMPS_SYM 0: non-symmetric 2: symmetric
integer(psb_ipk_), dimension(3) :: ipar
integer(psb_ipk_), allocatable :: local_ictxt
logical :: built = .false.
contains
@ -95,8 +96,8 @@ module mld_c_mumps_solver
procedure, pass(sv) :: default => c_mumps_solver_default
procedure, nopass :: get_fmt => c_mumps_solver_get_fmt
procedure, nopass :: get_id => c_mumps_solver_get_id
procedure, pass(sv) :: is_global => c_mumps_solver_is_global
#if defined(HAVE_FINAL)
final :: c_mumps_solver_finalize
#endif
end type mld_c_mumps_solver_type
@ -107,7 +108,7 @@ module mld_c_mumps_solver
& c_mumps_solver_sizeof, c_mumps_solver_apply_vect,&
& c_mumps_solver_cseti, c_mumps_solver_csetr, &
& c_mumps_solver_default, c_mumps_solver_get_fmt, &
& c_mumps_solver_get_id
& c_mumps_solver_get_id, c_mumps_solver_is_global
#if defined(HAVE_FINAL)
private :: c_mumps_solver_finalize
#endif
@ -295,9 +296,11 @@ contains
select case(psb_toupper(what))
#if defined(HAVE_MUMPS_)
case('MUMPS_LOC_GLOB')
sv%ipar(1)=val
sv%ipar(1) = val
case('MUMPS_PRINT_ERR')
sv%ipar(2)=val
sv%ipar(2) = val
case('MUMPS_SYM')
sv%ipar(3) = val
case('MUMPS_IPAR_ENTRY')
if(present(idx)) then
! Note: this will allocate %item
@ -455,7 +458,15 @@ contains
val = mld_mumps_
end function c_mumps_solver_get_id
#endif
function c_mumps_solver_is_global(sv) result(val)
implicit none
class(mld_c_mumps_solver_type), intent(in) :: sv
logical :: val
val = (sv%ipar(1) == mld_global_solver_ )
end function c_mumps_solver_is_global
end module mld_c_mumps_solver

@ -109,12 +109,13 @@ module mld_d_base_solver_mod
procedure, nopass :: get_fmt => d_base_solver_get_fmt
procedure, nopass :: get_id => d_base_solver_get_id
procedure, nopass :: is_iterative => d_base_solver_is_iterative
procedure, pass(sv) :: is_global => d_base_solver_is_global
end type mld_d_base_solver_type
private :: d_base_solver_sizeof, d_base_solver_default,&
& d_base_solver_get_nzeros, d_base_solver_get_fmt, &
& d_base_solver_is_iterative, d_base_solver_get_id, &
& d_base_solver_get_wrksize
& d_base_solver_get_wrksize, d_base_solver_is_global
interface
@ -365,6 +366,17 @@ contains
val = .false.
end function d_base_solver_is_iterative
!
! Is the solver acting globally? In most cases
! not, SuperLU_Dist does, MUMPS can do either.
!
function d_base_solver_is_global(sv) result(val)
implicit none
class(mld_d_base_solver_type), intent(in) :: sv
logical :: val
val = .false.
end function d_base_solver_is_global
function d_base_solver_get_id() result(val)
implicit none

@ -1,3 +1,4 @@
!
!
! MLD2P4 version 2.2
@ -58,12 +59,6 @@ module mld_d_mumps_solver
include 'dmumps_struc.h'
#endif
#if defined(LPK8)
type, extends(mld_d_base_solver_type) :: mld_d_mumps_solver_type
end type mld_d_mumps_solver_type
#else
type :: mld_d_mumps_icntl_item
integer(psb_ipk_), allocatable :: item
@ -80,7 +75,13 @@ module mld_d_mumps_solver
#endif
type(mld_d_mumps_icntl_item), allocatable :: icntl(:)
type(mld_d_mumps_rcntl_item), allocatable :: rcntl(:)
integer(psb_ipk_), dimension(2) :: ipar
!
! Controls to be set before MUMPS instantiation:
!
! IPAR(1) : MUMPS_LOC_GLOB 0==mld_local_solver_: LOCAL 1==mld_global_solver_: GLOBAL
! IPAR(2) : MUMPS_PRINT_ERR print verbosity (see MUMPS)
! IPAR(3) : MUMPS_SYM 0: non-symmetric 2: symmetric
integer(psb_ipk_), dimension(3) :: ipar
integer(psb_ipk_), allocatable :: local_ictxt
logical :: built = .false.
contains
@ -95,8 +96,8 @@ module mld_d_mumps_solver
procedure, pass(sv) :: default => d_mumps_solver_default
procedure, nopass :: get_fmt => d_mumps_solver_get_fmt
procedure, nopass :: get_id => d_mumps_solver_get_id
procedure, pass(sv) :: is_global => d_mumps_solver_is_global
#if defined(HAVE_FINAL)
final :: d_mumps_solver_finalize
#endif
end type mld_d_mumps_solver_type
@ -107,7 +108,7 @@ module mld_d_mumps_solver
& d_mumps_solver_sizeof, d_mumps_solver_apply_vect,&
& d_mumps_solver_cseti, d_mumps_solver_csetr, &
& d_mumps_solver_default, d_mumps_solver_get_fmt, &
& d_mumps_solver_get_id
& d_mumps_solver_get_id, d_mumps_solver_is_global
#if defined(HAVE_FINAL)
private :: d_mumps_solver_finalize
#endif
@ -295,9 +296,11 @@ contains
select case(psb_toupper(what))
#if defined(HAVE_MUMPS_)
case('MUMPS_LOC_GLOB')
sv%ipar(1)=val
sv%ipar(1) = val
case('MUMPS_PRINT_ERR')
sv%ipar(2)=val
sv%ipar(2) = val
case('MUMPS_SYM')
sv%ipar(3) = val
case('MUMPS_IPAR_ENTRY')
if(present(idx)) then
! Note: this will allocate %item
@ -455,7 +458,15 @@ contains
val = mld_mumps_
end function d_mumps_solver_get_id
#endif
function d_mumps_solver_is_global(sv) result(val)
implicit none
class(mld_d_mumps_solver_type), intent(in) :: sv
logical :: val
val = (sv%ipar(1) == mld_global_solver_ )
end function d_mumps_solver_is_global
end module mld_d_mumps_solver

@ -70,6 +70,7 @@ module mld_d_sludist_solver
procedure, pass(sv) :: sizeof => d_sludist_solver_sizeof
procedure, nopass :: get_fmt => d_sludist_solver_get_fmt
procedure, nopass :: get_id => d_sludist_solver_get_id
procedure, pass(sv) :: is_global => d_sludist_solver_is_global
#if defined(HAVE_FINAL)
final :: d_sludist_solver_finalize
#endif
@ -79,7 +80,8 @@ module mld_d_sludist_solver
private :: d_sludist_solver_bld, d_sludist_solver_apply, &
& d_sludist_solver_free, d_sludist_solver_descr, &
& d_sludist_solver_sizeof, d_sludist_solver_apply_vect, &
& d_sludist_solver_get_fmt, d_sludist_solver_get_id
& d_sludist_solver_get_fmt, d_sludist_solver_get_id, &
& d_sludist_solver_is_global
#if defined(HAVE_FINAL)
private :: d_sludist_solver_finalize
#endif
@ -355,6 +357,15 @@ contains
return
end subroutine d_sludist_solver_free
!
function d_sludist_solver_is_global(sv) result(val)
implicit none
class(mld_d_sludist_solver_type), intent(in) :: sv
logical :: val
val = .true.
end function d_sludist_solver_is_global
#if defined(HAVE_FINAL)
subroutine d_sludist_solver_finalize(sv)

@ -109,12 +109,13 @@ module mld_s_base_solver_mod
procedure, nopass :: get_fmt => s_base_solver_get_fmt
procedure, nopass :: get_id => s_base_solver_get_id
procedure, nopass :: is_iterative => s_base_solver_is_iterative
procedure, pass(sv) :: is_global => s_base_solver_is_global
end type mld_s_base_solver_type
private :: s_base_solver_sizeof, s_base_solver_default,&
& s_base_solver_get_nzeros, s_base_solver_get_fmt, &
& s_base_solver_is_iterative, s_base_solver_get_id, &
& s_base_solver_get_wrksize
& s_base_solver_get_wrksize, s_base_solver_is_global
interface
@ -365,6 +366,17 @@ contains
val = .false.
end function s_base_solver_is_iterative
!
! Is the solver acting globally? In most cases
! not, SuperLU_Dist does, MUMPS can do either.
!
function s_base_solver_is_global(sv) result(val)
implicit none
class(mld_s_base_solver_type), intent(in) :: sv
logical :: val
val = .false.
end function s_base_solver_is_global
function s_base_solver_get_id() result(val)
implicit none

@ -1,3 +1,4 @@
!
!
! MLD2P4 version 2.2
@ -58,12 +59,6 @@ module mld_s_mumps_solver
include 'smumps_struc.h'
#endif
#if defined(LPK8)
type, extends(mld_s_base_solver_type) :: mld_s_mumps_solver_type
end type mld_s_mumps_solver_type
#else
type :: mld_s_mumps_icntl_item
integer(psb_ipk_), allocatable :: item
@ -80,7 +75,13 @@ module mld_s_mumps_solver
#endif
type(mld_s_mumps_icntl_item), allocatable :: icntl(:)
type(mld_s_mumps_rcntl_item), allocatable :: rcntl(:)
integer(psb_ipk_), dimension(2) :: ipar
!
! Controls to be set before MUMPS instantiation:
!
! IPAR(1) : MUMPS_LOC_GLOB 0==mld_local_solver_: LOCAL 1==mld_global_solver_: GLOBAL
! IPAR(2) : MUMPS_PRINT_ERR print verbosity (see MUMPS)
! IPAR(3) : MUMPS_SYM 0: non-symmetric 2: symmetric
integer(psb_ipk_), dimension(3) :: ipar
integer(psb_ipk_), allocatable :: local_ictxt
logical :: built = .false.
contains
@ -95,8 +96,8 @@ module mld_s_mumps_solver
procedure, pass(sv) :: default => s_mumps_solver_default
procedure, nopass :: get_fmt => s_mumps_solver_get_fmt
procedure, nopass :: get_id => s_mumps_solver_get_id
procedure, pass(sv) :: is_global => s_mumps_solver_is_global
#if defined(HAVE_FINAL)
final :: s_mumps_solver_finalize
#endif
end type mld_s_mumps_solver_type
@ -107,7 +108,7 @@ module mld_s_mumps_solver
& s_mumps_solver_sizeof, s_mumps_solver_apply_vect,&
& s_mumps_solver_cseti, s_mumps_solver_csetr, &
& s_mumps_solver_default, s_mumps_solver_get_fmt, &
& s_mumps_solver_get_id
& s_mumps_solver_get_id, s_mumps_solver_is_global
#if defined(HAVE_FINAL)
private :: s_mumps_solver_finalize
#endif
@ -295,9 +296,11 @@ contains
select case(psb_toupper(what))
#if defined(HAVE_MUMPS_)
case('MUMPS_LOC_GLOB')
sv%ipar(1)=val
sv%ipar(1) = val
case('MUMPS_PRINT_ERR')
sv%ipar(2)=val
sv%ipar(2) = val
case('MUMPS_SYM')
sv%ipar(3) = val
case('MUMPS_IPAR_ENTRY')
if(present(idx)) then
! Note: this will allocate %item
@ -455,7 +458,15 @@ contains
val = mld_mumps_
end function s_mumps_solver_get_id
#endif
function s_mumps_solver_is_global(sv) result(val)
implicit none
class(mld_s_mumps_solver_type), intent(in) :: sv
logical :: val
val = (sv%ipar(1) == mld_global_solver_ )
end function s_mumps_solver_is_global
end module mld_s_mumps_solver

@ -109,12 +109,13 @@ module mld_z_base_solver_mod
procedure, nopass :: get_fmt => z_base_solver_get_fmt
procedure, nopass :: get_id => z_base_solver_get_id
procedure, nopass :: is_iterative => z_base_solver_is_iterative
procedure, pass(sv) :: is_global => z_base_solver_is_global
end type mld_z_base_solver_type
private :: z_base_solver_sizeof, z_base_solver_default,&
& z_base_solver_get_nzeros, z_base_solver_get_fmt, &
& z_base_solver_is_iterative, z_base_solver_get_id, &
& z_base_solver_get_wrksize
& z_base_solver_get_wrksize, z_base_solver_is_global
interface
@ -365,6 +366,17 @@ contains
val = .false.
end function z_base_solver_is_iterative
!
! Is the solver acting globally? In most cases
! not, SuperLU_Dist does, MUMPS can do either.
!
function z_base_solver_is_global(sv) result(val)
implicit none
class(mld_z_base_solver_type), intent(in) :: sv
logical :: val
val = .false.
end function z_base_solver_is_global
function z_base_solver_get_id() result(val)
implicit none

@ -1,3 +1,4 @@
!
!
! MLD2P4 version 2.2
@ -58,12 +59,6 @@ module mld_z_mumps_solver
include 'zmumps_struc.h'
#endif
#if defined(LPK8)
type, extends(mld_z_base_solver_type) :: mld_z_mumps_solver_type
end type mld_z_mumps_solver_type
#else
type :: mld_z_mumps_icntl_item
integer(psb_ipk_), allocatable :: item
@ -80,7 +75,13 @@ module mld_z_mumps_solver
#endif
type(mld_z_mumps_icntl_item), allocatable :: icntl(:)
type(mld_z_mumps_rcntl_item), allocatable :: rcntl(:)
integer(psb_ipk_), dimension(2) :: ipar
!
! Controls to be set before MUMPS instantiation:
!
! IPAR(1) : MUMPS_LOC_GLOB 0==mld_local_solver_: LOCAL 1==mld_global_solver_: GLOBAL
! IPAR(2) : MUMPS_PRINT_ERR print verbosity (see MUMPS)
! IPAR(3) : MUMPS_SYM 0: non-symmetric 2: symmetric
integer(psb_ipk_), dimension(3) :: ipar
integer(psb_ipk_), allocatable :: local_ictxt
logical :: built = .false.
contains
@ -95,8 +96,8 @@ module mld_z_mumps_solver
procedure, pass(sv) :: default => z_mumps_solver_default
procedure, nopass :: get_fmt => z_mumps_solver_get_fmt
procedure, nopass :: get_id => z_mumps_solver_get_id
procedure, pass(sv) :: is_global => z_mumps_solver_is_global
#if defined(HAVE_FINAL)
final :: z_mumps_solver_finalize
#endif
end type mld_z_mumps_solver_type
@ -107,7 +108,7 @@ module mld_z_mumps_solver
& z_mumps_solver_sizeof, z_mumps_solver_apply_vect,&
& z_mumps_solver_cseti, z_mumps_solver_csetr, &
& z_mumps_solver_default, z_mumps_solver_get_fmt, &
& z_mumps_solver_get_id
& z_mumps_solver_get_id, z_mumps_solver_is_global
#if defined(HAVE_FINAL)
private :: z_mumps_solver_finalize
#endif
@ -295,9 +296,11 @@ contains
select case(psb_toupper(what))
#if defined(HAVE_MUMPS_)
case('MUMPS_LOC_GLOB')
sv%ipar(1)=val
sv%ipar(1) = val
case('MUMPS_PRINT_ERR')
sv%ipar(2)=val
sv%ipar(2) = val
case('MUMPS_SYM')
sv%ipar(3) = val
case('MUMPS_IPAR_ENTRY')
if(present(idx)) then
! Note: this will allocate %item
@ -455,7 +458,15 @@ contains
val = mld_mumps_
end function z_mumps_solver_get_id
#endif
function z_mumps_solver_is_global(sv) result(val)
implicit none
class(mld_z_mumps_solver_type), intent(in) :: sv
logical :: val
val = (sv%ipar(1) == mld_global_solver_ )
end function z_mumps_solver_is_global
end module mld_z_mumps_solver

@ -70,6 +70,7 @@ module mld_z_sludist_solver
procedure, pass(sv) :: sizeof => z_sludist_solver_sizeof
procedure, nopass :: get_fmt => z_sludist_solver_get_fmt
procedure, nopass :: get_id => z_sludist_solver_get_id
procedure, pass(sv) :: is_global => z_sludist_solver_is_global
#if defined(HAVE_FINAL)
final :: z_sludist_solver_finalize
#endif
@ -79,7 +80,8 @@ module mld_z_sludist_solver
private :: z_sludist_solver_bld, z_sludist_solver_apply, &
& z_sludist_solver_free, z_sludist_solver_descr, &
& z_sludist_solver_sizeof, z_sludist_solver_apply_vect, &
& z_sludist_solver_get_fmt, z_sludist_solver_get_id
& z_sludist_solver_get_fmt, z_sludist_solver_get_id, &
& z_sludist_solver_is_global
#if defined(HAVE_FINAL)
private :: z_sludist_solver_finalize
#endif
@ -355,6 +357,15 @@ contains
return
end subroutine z_sludist_solver_free
!
function z_sludist_solver_is_global(sv) result(val)
implicit none
class(mld_z_sludist_solver_type), intent(in) :: sv
logical :: val
val = .true.
end function z_sludist_solver_is_global
#if defined(HAVE_FINAL)
subroutine z_sludist_solver_finalize(sv)

Loading…
Cancel
Save