diff --git a/mlprec/impl/solver/mld_c_mumps_solver_apply.F90 b/mlprec/impl/solver/mld_c_mumps_solver_apply.F90 index 436c3bf9..a0f9da5e 100644 --- a/mlprec/impl/solver/mld_c_mumps_solver_apply.F90 +++ b/mlprec/impl/solver/mld_c_mumps_solver_apply.F90 @@ -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 diff --git a/mlprec/impl/solver/mld_c_mumps_solver_bld.F90 b/mlprec/impl/solver/mld_c_mumps_solver_bld.F90 index 31e9d6b2..d5f62827 100644 --- a/mlprec/impl/solver/mld_c_mumps_solver_bld.F90 +++ b/mlprec/impl/solver/mld_c_mumps_solver_bld.F90 @@ -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_) 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 + ! 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 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' + 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 - else - call a%cp_to(acoo) - end if - nztota = acoo%get_nzeros() + 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) - ! 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') + 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 - 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) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mld_cmumps_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) + else + call a%cp_to(acoo) + end if + nztota = acoo%get_nzeros() - call acoo%free() - sv%built=.true. + ! 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) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + call acoo%free() + sv%built=.true. - call psb_erractionrestore(err_act) - return + 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 diff --git a/mlprec/impl/solver/mld_d_mumps_solver_apply.F90 b/mlprec/impl/solver/mld_d_mumps_solver_apply.F90 index 3900d211..0facc887 100644 --- a/mlprec/impl/solver/mld_d_mumps_solver_apply.F90 +++ b/mlprec/impl/solver/mld_d_mumps_solver_apply.F90 @@ -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 diff --git a/mlprec/impl/solver/mld_d_mumps_solver_bld.F90 b/mlprec/impl/solver/mld_d_mumps_solver_bld.F90 index 8ec378e0..44da6b05 100644 --- a/mlprec/impl/solver/mld_d_mumps_solver_bld.F90 +++ b/mlprec/impl/solver/mld_d_mumps_solver_bld.F90 @@ -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_) 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 + ! 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' + 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 - else - call a%cp_to(acoo) - end if - nztota = acoo%get_nzeros() + 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) - ! 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') + 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 - 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) - 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) + else + call a%cp_to(acoo) + end if + nztota = acoo%get_nzeros() - call acoo%free() - sv%built=.true. + ! 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) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + call acoo%free() + sv%built=.true. - call psb_erractionrestore(err_act) - return + 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 diff --git a/mlprec/impl/solver/mld_s_mumps_solver_apply.F90 b/mlprec/impl/solver/mld_s_mumps_solver_apply.F90 index f0fae77d..8c5aebdb 100644 --- a/mlprec/impl/solver/mld_s_mumps_solver_apply.F90 +++ b/mlprec/impl/solver/mld_s_mumps_solver_apply.F90 @@ -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 diff --git a/mlprec/impl/solver/mld_s_mumps_solver_bld.F90 b/mlprec/impl/solver/mld_s_mumps_solver_bld.F90 index e7f46fcc..caeea837 100644 --- a/mlprec/impl/solver/mld_s_mumps_solver_bld.F90 +++ b/mlprec/impl/solver/mld_s_mumps_solver_bld.F90 @@ -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_) 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 + ! 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 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' + 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 - else - call a%cp_to(acoo) - end if - nztota = acoo%get_nzeros() + 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) - ! 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') + 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 - 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) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mld_smumps_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) + else + call a%cp_to(acoo) + end if + nztota = acoo%get_nzeros() - call acoo%free() - sv%built=.true. + ! 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) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + call acoo%free() + sv%built=.true. - call psb_erractionrestore(err_act) - return + 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 diff --git a/mlprec/impl/solver/mld_z_mumps_solver_apply.F90 b/mlprec/impl/solver/mld_z_mumps_solver_apply.F90 index 252e4c3a..71742830 100644 --- a/mlprec/impl/solver/mld_z_mumps_solver_apply.F90 +++ b/mlprec/impl/solver/mld_z_mumps_solver_apply.F90 @@ -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 diff --git a/mlprec/impl/solver/mld_z_mumps_solver_bld.F90 b/mlprec/impl/solver/mld_z_mumps_solver_bld.F90 index 467530e8..bbf1fe52 100644 --- a/mlprec/impl/solver/mld_z_mumps_solver_bld.F90 +++ b/mlprec/impl/solver/mld_z_mumps_solver_bld.F90 @@ -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_) 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 + ! 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 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' + 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 - else - call a%cp_to(acoo) - end if - nztota = acoo%get_nzeros() + 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) - ! 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') + 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 - 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) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mld_zmumps_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) + else + call a%cp_to(acoo) + end if + nztota = acoo%get_nzeros() - call acoo%free() - sv%built=.true. + ! 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) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + call acoo%free() + sv%built=.true. - call psb_erractionrestore(err_act) - return + 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 diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index f07d30a8..6bd3b515 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -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 diff --git a/mlprec/mld_c_mumps_solver.F90 b/mlprec/mld_c_mumps_solver.F90 index 3728eec6..47e1082c 100644 --- a/mlprec/mld_c_mumps_solver.F90 +++ b/mlprec/mld_c_mumps_solver.F90 @@ -291,7 +291,7 @@ contains select case(psb_toupper(what)) #if defined(HAVE_MUMPS_) - case('MUMPS_AS_SEQUENTIAL') + case('MUMPS_SOLVER_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 diff --git a/mlprec/mld_d_mumps_solver.F90 b/mlprec/mld_d_mumps_solver.F90 index a8e225f2..39360791 100644 --- a/mlprec/mld_d_mumps_solver.F90 +++ b/mlprec/mld_d_mumps_solver.F90 @@ -291,7 +291,7 @@ contains select case(psb_toupper(what)) #if defined(HAVE_MUMPS_) - case('MUMPS_AS_SEQUENTIAL') + case('MUMPS_SOLVER_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 diff --git a/mlprec/mld_s_mumps_solver.F90 b/mlprec/mld_s_mumps_solver.F90 index 8798d90e..d93ca4dc 100644 --- a/mlprec/mld_s_mumps_solver.F90 +++ b/mlprec/mld_s_mumps_solver.F90 @@ -291,7 +291,7 @@ contains select case(psb_toupper(what)) #if defined(HAVE_MUMPS_) - case('MUMPS_AS_SEQUENTIAL') + case('MUMPS_SOLVER_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 diff --git a/mlprec/mld_z_mumps_solver.F90 b/mlprec/mld_z_mumps_solver.F90 index 6732bf9f..1c9f2995 100644 --- a/mlprec/mld_z_mumps_solver.F90 +++ b/mlprec/mld_z_mumps_solver.F90 @@ -291,7 +291,7 @@ contains select case(psb_toupper(what)) #if defined(HAVE_MUMPS_) - case('MUMPS_AS_SEQUENTIAL') + case('MUMPS_SOLVER_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