diff --git a/amgprec/amg_c_slu_solver.F90 b/amgprec/amg_c_slu_solver.F90 index d312ae0c..902e39f3 100644 --- a/amgprec/amg_c_slu_solver.F90 +++ b/amgprec/amg_c_slu_solver.F90 @@ -63,9 +63,9 @@ module amg_c_slu_solver type(c_ptr) :: lufactors=c_null_ptr integer(c_long_long) :: symbsize=0, numsize=0 contains - procedure, pass(sv) :: build => c_slu_solver_bld - procedure, pass(sv) :: apply_a => c_slu_solver_apply - procedure, pass(sv) :: apply_v => c_slu_solver_apply_vect + procedure, pass(sv) :: build => amg_c_slu_solver_bld + procedure, pass(sv) :: apply_a => amg_c_slu_solver_apply + procedure, pass(sv) :: apply_v => amg_c_slu_solver_apply_vect procedure, pass(sv) :: free => c_slu_solver_free procedure, pass(sv) :: clear_data => c_slu_solver_clear_data procedure, pass(sv) :: descr => c_slu_solver_descr @@ -76,9 +76,8 @@ module amg_c_slu_solver end type amg_c_slu_solver_type - private :: c_slu_solver_bld, c_slu_solver_apply, & - & c_slu_solver_free, c_slu_solver_descr, & - & c_slu_solver_sizeof, c_slu_solver_apply_vect, & + private :: c_slu_solver_free, c_slu_solver_descr, & + & c_slu_solver_sizeof, & & c_slu_solver_get_fmt, c_slu_solver_get_id, & & c_slu_solver_clear_data private :: c_slu_solver_finalize @@ -118,207 +117,264 @@ module amg_c_slu_solver end function amg_cslu_free end interface -contains - - subroutine c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& - & trans,work,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_c_slu_solver_type), intent(inout) :: sv - complex(psb_spk_),intent(inout) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - character, intent(in), optional :: init - complex(psb_spk_),intent(inout), optional :: initu(:) - - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:) - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_slu_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - 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) - 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 - endif - - ww(1:n_row) = x(1:n_row) - select case(trans_) - case('N') - info = amg_cslu_solve(0,n_row,1,ww,n_row,sv%lufactors) - case('T') - info = amg_cslu_solve(1,n_row,1,ww,n_row,sv%lufactors) - case('C') - info = amg_cslu_solve(2,n_row,1,ww,n_row,sv%lufactors) - case default - call psb_errpush(psb_err_internal_error_, & - & name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - if (info == psb_success_) & - & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,& - & name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col > size(work)) then - deallocate(ww) - endif - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine c_slu_solver_apply - - subroutine c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& - & trans,work,wv,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_c_slu_solver_type), intent(inout) :: sv - type(psb_c_vect_type),intent(inout) :: x - type(psb_c_vect_type),intent(inout) :: y - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - type(psb_c_vect_type),intent(inout) :: wv(:) - integer, intent(out) :: info - character, intent(in), optional :: init - type(psb_c_vect_type),intent(inout), optional :: initu - - integer :: err_act - character(len=20) :: name='c_slu_solver_apply_vect' - - call psb_erractionsave(err_act) - - info = psb_success_ - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - call x%v%sync() - call y%v%sync() - call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) - call y%v%set_host() - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine c_slu_solver_apply_vect - - subroutine c_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) - - use psb_base_mod - - Implicit None + interface + subroutine amg_c_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + import amg_c_slu_solver_type + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_c_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_c_slu_solver_bld + end interface - ! Arguments - type(psb_cspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(inout) :: desc_a - class(amg_c_slu_solver_type), intent(inout) :: sv - integer, intent(out) :: info - type(psb_cspmat_type), intent(in), target, optional :: b - class(psb_c_base_sparse_mat), intent(in), optional :: amold - class(psb_c_base_vect_type), intent(in), optional :: vmold - class(psb_i_base_vect_type), intent(in), optional :: imold - ! Local variables - type(psb_cspmat_type) :: atmp - type(psb_c_csc_sparse_mat) :: acsc - type(psb_c_coo_sparse_mat) :: acoo - integer :: n_row,n_col, nrow_a, nztota - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='c_slu_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ctxt = desc_a%get_context() - call psb_info(ctxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - - - call a%cscnv(atmp,info,type='coo') - call psb_rwextd(n_row,atmp,info,b=b) - call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) - nrow_a = atmp%get_nrows() - call atmp%a%csclip(acoo,info,jmax=nrow_a) - call acsc%mv_from_coo(acoo,info) - nztota = acsc%get_nzeros() - ! Fix the entries to call C-base SuperLU - acsc%ia(:) = acsc%ia(:) - 1 - acsc%icp(:) = acsc%icp(:) - 1 - info = amg_cslu_fact(nrow_a,nztota,acsc%val,& - & acsc%icp,acsc%ia,sv%lufactors) - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='amg_cslu_fact' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call acsc%free() - call atmp%free() + interface + subroutine amg_c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + import amg_c_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_slu_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_c_vect_type),intent(inout), optional :: initu + end subroutine amg_c_slu_solver_apply_vect + end interface - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + interface + subroutine amg_c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + import amg_c_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_slu_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_spk_),intent(inout), optional :: initu(:) + end subroutine amg_c_slu_solver_apply + end interface - call psb_erractionrestore(err_act) - return +contains -9999 call psb_error_handler(err_act) - return - end subroutine c_slu_solver_bld +!!$ subroutine c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_c_slu_solver_type), intent(inout) :: sv +!!$ complex(psb_spk_),intent(inout) :: x(:) +!!$ complex(psb_spk_),intent(inout) :: y(:) +!!$ complex(psb_spk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ complex(psb_spk_),target, intent(inout) :: work(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ complex(psb_spk_),intent(inout), optional :: initu(:) +!!$ +!!$ integer :: n_row,n_col +!!$ complex(psb_spk_), pointer :: ww(:) +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act +!!$ character :: trans_ +!!$ character(len=20) :: name='c_slu_solver_apply' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ +!!$ trans_ = psb_toupper(trans) +!!$ select case(trans_) +!!$ case('N') +!!$ case('T','C') +!!$ case default +!!$ call psb_errpush(psb_err_iarg_invalid_i_,name) +!!$ goto 9999 +!!$ end select +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ 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) +!!$ 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 +!!$ endif +!!$ +!!$ ww(1:n_row) = x(1:n_row) +!!$ select case(trans_) +!!$ case('N') +!!$ info = amg_cslu_solve(0,n_row,1,ww,n_row,sv%lufactors) +!!$ case('T') +!!$ info = amg_cslu_solve(1,n_row,1,ww,n_row,sv%lufactors) +!!$ case('C') +!!$ info = amg_cslu_solve(2,n_row,1,ww,n_row,sv%lufactors) +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_, & +!!$ & name,a_err='Invalid TRANS in ILU subsolve') +!!$ goto 9999 +!!$ end select +!!$ +!!$ if (info == psb_success_) & +!!$ & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) +!!$ +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in subsolve') +!!$ goto 9999 +!!$ endif +!!$ +!!$ if (n_col > size(work)) then +!!$ deallocate(ww) +!!$ endif +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine c_slu_solver_apply +!!$ +!!$ subroutine c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,wv,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_c_slu_solver_type), intent(inout) :: sv +!!$ type(psb_c_vect_type),intent(inout) :: x +!!$ type(psb_c_vect_type),intent(inout) :: y +!!$ complex(psb_spk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ complex(psb_spk_),target, intent(inout) :: work(:) +!!$ type(psb_c_vect_type),intent(inout) :: wv(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ type(psb_c_vect_type),intent(inout), optional :: initu +!!$ +!!$ integer :: err_act +!!$ character(len=20) :: name='c_slu_solver_apply_vect' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ call x%v%sync() +!!$ call y%v%sync() +!!$ call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) +!!$ call y%v%set_host() +!!$ if (info /= 0) goto 9999 +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine c_slu_solver_apply_vect +!!$ +!!$ subroutine c_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) +!!$ +!!$ use psb_base_mod +!!$ +!!$ Implicit None +!!$ +!!$ ! Arguments +!!$ type(psb_cspmat_type), intent(in), target :: a +!!$ Type(psb_desc_type), Intent(inout) :: desc_a +!!$ class(amg_c_slu_solver_type), intent(inout) :: sv +!!$ integer, intent(out) :: info +!!$ type(psb_cspmat_type), intent(in), target, optional :: b +!!$ class(psb_c_base_sparse_mat), intent(in), optional :: amold +!!$ class(psb_c_base_vect_type), intent(in), optional :: vmold +!!$ class(psb_i_base_vect_type), intent(in), optional :: imold +!!$ ! Local variables +!!$ type(psb_cspmat_type) :: atmp +!!$ type(psb_c_csc_sparse_mat) :: acsc +!!$ type(psb_c_coo_sparse_mat) :: acoo +!!$ integer :: n_row,n_col, nrow_a, nztota +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act, debug_unit, debug_level +!!$ character(len=20) :: name='c_slu_solver_bld', ch_err +!!$ +!!$ info=psb_success_ +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ ctxt = desc_a%get_context() +!!$ call psb_info(ctxt, me, np) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' start' +!!$ +!!$ +!!$ n_row = desc_a%get_local_rows() +!!$ n_col = desc_a%get_local_cols() +!!$ +!!$ +!!$ call a%cscnv(atmp,info,type='coo') +!!$ call psb_rwextd(n_row,atmp,info,b=b) +!!$ call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ nrow_a = atmp%get_nrows() +!!$ call atmp%a%csclip(acoo,info,jmax=nrow_a) +!!$ call acsc%mv_from_coo(acoo,info) +!!$ nztota = acsc%get_nzeros() +!!$ ! Fix the entries to call C-base SuperLU +!!$ acsc%ia(:) = acsc%ia(:) - 1 +!!$ acsc%icp(:) = acsc%icp(:) - 1 +!!$ info = amg_cslu_fact(nrow_a,nztota,acsc%val,& +!!$ & acsc%icp,acsc%ia,sv%lufactors) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='amg_cslu_fact' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ +!!$ call acsc%free() +!!$ call atmp%free() +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' end' +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ end subroutine c_slu_solver_bld subroutine c_slu_solver_free(sv,info) diff --git a/amgprec/amg_d_slu_solver.F90 b/amgprec/amg_d_slu_solver.F90 index 227d6e18..946dba89 100644 --- a/amgprec/amg_d_slu_solver.F90 +++ b/amgprec/amg_d_slu_solver.F90 @@ -63,9 +63,9 @@ module amg_d_slu_solver type(c_ptr) :: lufactors=c_null_ptr integer(c_long_long) :: symbsize=0, numsize=0 contains - procedure, pass(sv) :: build => d_slu_solver_bld - procedure, pass(sv) :: apply_a => d_slu_solver_apply - procedure, pass(sv) :: apply_v => d_slu_solver_apply_vect + procedure, pass(sv) :: build => amg_d_slu_solver_bld + procedure, pass(sv) :: apply_a => amg_d_slu_solver_apply + procedure, pass(sv) :: apply_v => amg_d_slu_solver_apply_vect procedure, pass(sv) :: free => d_slu_solver_free procedure, pass(sv) :: clear_data => d_slu_solver_clear_data procedure, pass(sv) :: descr => d_slu_solver_descr @@ -76,9 +76,8 @@ module amg_d_slu_solver end type amg_d_slu_solver_type - private :: d_slu_solver_bld, d_slu_solver_apply, & - & d_slu_solver_free, d_slu_solver_descr, & - & d_slu_solver_sizeof, d_slu_solver_apply_vect, & + private :: d_slu_solver_free, d_slu_solver_descr, & + & d_slu_solver_sizeof, & & d_slu_solver_get_fmt, d_slu_solver_get_id, & & d_slu_solver_clear_data private :: d_slu_solver_finalize @@ -118,207 +117,264 @@ module amg_d_slu_solver end function amg_dslu_free end interface -contains - - subroutine d_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& - & trans,work,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_d_slu_solver_type), intent(inout) :: sv - real(psb_dpk_),intent(inout) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - character, intent(in), optional :: init - real(psb_dpk_),intent(inout), optional :: initu(:) - - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:) - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act - character :: trans_ - character(len=20) :: name='d_slu_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - 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) - 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 - endif - - ww(1:n_row) = x(1:n_row) - select case(trans_) - case('N') - info = amg_dslu_solve(0,n_row,1,ww,n_row,sv%lufactors) - case('T') - info = amg_dslu_solve(1,n_row,1,ww,n_row,sv%lufactors) - case('C') - info = amg_dslu_solve(2,n_row,1,ww,n_row,sv%lufactors) - case default - call psb_errpush(psb_err_internal_error_, & - & name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - if (info == psb_success_) & - & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,& - & name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col > size(work)) then - deallocate(ww) - endif - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine d_slu_solver_apply - - subroutine d_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& - & trans,work,wv,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_d_slu_solver_type), intent(inout) :: sv - type(psb_d_vect_type),intent(inout) :: x - type(psb_d_vect_type),intent(inout) :: y - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - type(psb_d_vect_type),intent(inout) :: wv(:) - integer, intent(out) :: info - character, intent(in), optional :: init - type(psb_d_vect_type),intent(inout), optional :: initu - - integer :: err_act - character(len=20) :: name='d_slu_solver_apply_vect' - - call psb_erractionsave(err_act) - - info = psb_success_ - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - call x%v%sync() - call y%v%sync() - call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) - call y%v%set_host() - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine d_slu_solver_apply_vect - - subroutine d_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) - - use psb_base_mod - - Implicit None + interface + subroutine amg_d_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + import amg_d_slu_solver_type + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_d_slu_solver_bld + end interface - ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(inout) :: desc_a - class(amg_d_slu_solver_type), intent(inout) :: sv - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target, optional :: b - class(psb_d_base_sparse_mat), intent(in), optional :: amold - class(psb_d_base_vect_type), intent(in), optional :: vmold - class(psb_i_base_vect_type), intent(in), optional :: imold - ! Local variables - type(psb_dspmat_type) :: atmp - type(psb_d_csc_sparse_mat) :: acsc - type(psb_d_coo_sparse_mat) :: acoo - integer :: n_row,n_col, nrow_a, nztota - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='d_slu_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ctxt = desc_a%get_context() - call psb_info(ctxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - - - call a%cscnv(atmp,info,type='coo') - call psb_rwextd(n_row,atmp,info,b=b) - call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) - nrow_a = atmp%get_nrows() - call atmp%a%csclip(acoo,info,jmax=nrow_a) - call acsc%mv_from_coo(acoo,info) - nztota = acsc%get_nzeros() - ! Fix the entries to call C-base SuperLU - acsc%ia(:) = acsc%ia(:) - 1 - acsc%icp(:) = acsc%icp(:) - 1 - info = amg_dslu_fact(nrow_a,nztota,acsc%val,& - & acsc%icp,acsc%ia,sv%lufactors) - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='amg_dslu_fact' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call acsc%free() - call atmp%free() + interface + subroutine amg_d_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + import amg_d_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_slu_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_d_vect_type),intent(inout), optional :: initu + end subroutine amg_d_slu_solver_apply_vect + end interface - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + interface + subroutine amg_d_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + import amg_d_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_slu_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_dpk_),intent(inout), optional :: initu(:) + end subroutine amg_d_slu_solver_apply + end interface - call psb_erractionrestore(err_act) - return +contains -9999 call psb_error_handler(err_act) - return - end subroutine d_slu_solver_bld +!!$ subroutine d_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_d_slu_solver_type), intent(inout) :: sv +!!$ real(psb_dpk_),intent(inout) :: x(:) +!!$ real(psb_dpk_),intent(inout) :: y(:) +!!$ real(psb_dpk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ real(psb_dpk_),target, intent(inout) :: work(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ real(psb_dpk_),intent(inout), optional :: initu(:) +!!$ +!!$ integer :: n_row,n_col +!!$ real(psb_dpk_), pointer :: ww(:) +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act +!!$ character :: trans_ +!!$ character(len=20) :: name='d_slu_solver_apply' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ +!!$ trans_ = psb_toupper(trans) +!!$ select case(trans_) +!!$ case('N') +!!$ case('T','C') +!!$ case default +!!$ call psb_errpush(psb_err_iarg_invalid_i_,name) +!!$ goto 9999 +!!$ end select +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ 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) +!!$ 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 +!!$ endif +!!$ +!!$ ww(1:n_row) = x(1:n_row) +!!$ select case(trans_) +!!$ case('N') +!!$ info = amg_dslu_solve(0,n_row,1,ww,n_row,sv%lufactors) +!!$ case('T') +!!$ info = amg_dslu_solve(1,n_row,1,ww,n_row,sv%lufactors) +!!$ case('C') +!!$ info = amg_dslu_solve(2,n_row,1,ww,n_row,sv%lufactors) +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_, & +!!$ & name,a_err='Invalid TRANS in ILU subsolve') +!!$ goto 9999 +!!$ end select +!!$ +!!$ if (info == psb_success_) & +!!$ & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) +!!$ +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in subsolve') +!!$ goto 9999 +!!$ endif +!!$ +!!$ if (n_col > size(work)) then +!!$ deallocate(ww) +!!$ endif +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine d_slu_solver_apply +!!$ +!!$ subroutine d_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,wv,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_d_slu_solver_type), intent(inout) :: sv +!!$ type(psb_d_vect_type),intent(inout) :: x +!!$ type(psb_d_vect_type),intent(inout) :: y +!!$ real(psb_dpk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ real(psb_dpk_),target, intent(inout) :: work(:) +!!$ type(psb_d_vect_type),intent(inout) :: wv(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ type(psb_d_vect_type),intent(inout), optional :: initu +!!$ +!!$ integer :: err_act +!!$ character(len=20) :: name='d_slu_solver_apply_vect' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ call x%v%sync() +!!$ call y%v%sync() +!!$ call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) +!!$ call y%v%set_host() +!!$ if (info /= 0) goto 9999 +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine d_slu_solver_apply_vect +!!$ +!!$ subroutine d_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) +!!$ +!!$ use psb_base_mod +!!$ +!!$ Implicit None +!!$ +!!$ ! Arguments +!!$ type(psb_dspmat_type), intent(in), target :: a +!!$ Type(psb_desc_type), Intent(inout) :: desc_a +!!$ class(amg_d_slu_solver_type), intent(inout) :: sv +!!$ integer, intent(out) :: info +!!$ type(psb_dspmat_type), intent(in), target, optional :: b +!!$ class(psb_d_base_sparse_mat), intent(in), optional :: amold +!!$ class(psb_d_base_vect_type), intent(in), optional :: vmold +!!$ class(psb_i_base_vect_type), intent(in), optional :: imold +!!$ ! Local variables +!!$ type(psb_dspmat_type) :: atmp +!!$ type(psb_d_csc_sparse_mat) :: acsc +!!$ type(psb_d_coo_sparse_mat) :: acoo +!!$ integer :: n_row,n_col, nrow_a, nztota +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act, debug_unit, debug_level +!!$ character(len=20) :: name='d_slu_solver_bld', ch_err +!!$ +!!$ info=psb_success_ +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ ctxt = desc_a%get_context() +!!$ call psb_info(ctxt, me, np) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' start' +!!$ +!!$ +!!$ n_row = desc_a%get_local_rows() +!!$ n_col = desc_a%get_local_cols() +!!$ +!!$ +!!$ call a%cscnv(atmp,info,type='coo') +!!$ call psb_rwextd(n_row,atmp,info,b=b) +!!$ call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ nrow_a = atmp%get_nrows() +!!$ call atmp%a%csclip(acoo,info,jmax=nrow_a) +!!$ call acsc%mv_from_coo(acoo,info) +!!$ nztota = acsc%get_nzeros() +!!$ ! Fix the entries to call C-base SuperLU +!!$ acsc%ia(:) = acsc%ia(:) - 1 +!!$ acsc%icp(:) = acsc%icp(:) - 1 +!!$ info = amg_dslu_fact(nrow_a,nztota,acsc%val,& +!!$ & acsc%icp,acsc%ia,sv%lufactors) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='amg_dslu_fact' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ +!!$ call acsc%free() +!!$ call atmp%free() +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' end' +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ end subroutine d_slu_solver_bld subroutine d_slu_solver_free(sv,info) diff --git a/amgprec/amg_d_umf_solver.F90 b/amgprec/amg_d_umf_solver.F90 index a3fc939e..c91cb16f 100644 --- a/amgprec/amg_d_umf_solver.F90 +++ b/amgprec/amg_d_umf_solver.F90 @@ -62,9 +62,9 @@ module amg_d_umf_solver type(c_ptr) :: symbolic=c_null_ptr, numeric=c_null_ptr integer(c_long_long) :: symbsize=0, numsize=0 contains - procedure, pass(sv) :: build => d_umf_solver_bld - procedure, pass(sv) :: apply_a => d_umf_solver_apply - procedure, pass(sv) :: apply_v => d_umf_solver_apply_vect + procedure, pass(sv) :: build => amg_d_umf_solver_bld + procedure, pass(sv) :: apply_a => amg_d_umf_solver_apply + procedure, pass(sv) :: apply_v => amg_d_umf_solver_apply_vect procedure, pass(sv) :: free => d_umf_solver_free procedure, pass(sv) :: clear_data => d_umf_solver_clear_data procedure, pass(sv) :: descr => d_umf_solver_descr @@ -75,9 +75,8 @@ module amg_d_umf_solver end type amg_d_umf_solver_type - private :: d_umf_solver_bld, d_umf_solver_apply, & - & d_umf_solver_free, d_umf_solver_descr, & - & d_umf_solver_sizeof, d_umf_solver_apply_vect, & + private :: d_umf_solver_free, d_umf_solver_descr, & + & d_umf_solver_sizeof, & & d_umf_solver_get_fmt, d_umf_solver_get_id, & & d_umf_solver_clear_data private :: d_umf_solver_finalize @@ -118,208 +117,65 @@ module amg_d_umf_solver end function amg_dumf_free end interface -contains - - subroutine d_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& - & trans,work,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_d_umf_solver_type), intent(inout) :: sv - real(psb_dpk_),intent(inout) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - character, intent(in), optional :: init - real(psb_dpk_),intent(inout), optional :: initu(:) - - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:) - integer(psb_ipk_) :: i, err_act - character :: trans_ - character(len=20) :: name='d_umf_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - 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) - 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 - endif - - select case(trans_) - case('N') - info = amg_dumf_solve(0,n_row,ww,x,n_row,sv%numeric) - case('T') - ! - ! Note: with UMF, 1 meand Ctranspose, 2 means transpose - ! even for complex data. - ! - if (psb_d_is_complex_) then - info = amg_dumf_solve(2,n_row,ww,x,n_row,sv%numeric) - else - info = amg_dumf_solve(1,n_row,ww,x,n_row,sv%numeric) - end if - case('C') - info = amg_dumf_solve(1,n_row,ww,x,n_row,sv%numeric) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col > size(work)) then - deallocate(ww) - endif - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine d_umf_solver_apply - - subroutine d_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& - & trans,work,wv,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_d_umf_solver_type), intent(inout) :: sv - type(psb_d_vect_type),intent(inout) :: x - type(psb_d_vect_type),intent(inout) :: y - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - type(psb_d_vect_type),intent(inout) :: wv(:) - integer, intent(out) :: info - character, intent(in), optional :: init - type(psb_d_vect_type),intent(inout), optional :: initu - - integer :: err_act - character(len=20) :: name='d_umf_solver_apply_vect' - - call psb_erractionsave(err_act) - - info = psb_success_ - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - call x%v%sync() - call y%v%sync() - call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) - call y%v%set_host() - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine d_umf_solver_apply_vect - - subroutine d_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) - - use psb_base_mod - - Implicit None - - ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(inout) :: desc_a - class(amg_d_umf_solver_type), intent(inout) :: sv - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target, optional :: b - class(psb_d_base_sparse_mat), intent(in), optional :: amold - class(psb_d_base_vect_type), intent(in), optional :: vmold - class(psb_i_base_vect_type), intent(in), optional :: imold - ! Local variables - type(psb_dspmat_type) :: atmp - type(psb_d_csc_sparse_mat) :: acsc - integer :: n_row,n_col, nrow_a, nztota - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='d_umf_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ctxt = desc_a%get_context() - call psb_info(ctxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - - call a%cscnv(atmp,info,type='coo') - call psb_rwextd(n_row,atmp,info,b=b) - call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_) - call atmp%mv_to(acsc) - nrow_a = acsc%get_nrows() - nztota = acsc%get_nzeros() - ! Fix the entres to call C-base UMFPACK. - acsc%ia(:) = acsc%ia(:) - 1 - acsc%icp(:) = acsc%icp(:) - 1 - info = amg_dumf_fact(nrow_a,nztota,acsc%val,& - & acsc%ia,acsc%icp,sv%symbolic,sv%numeric,& - & sv%symbsize,sv%numsize) - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='amg_dumf_fact' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call acsc%free() - call atmp%free() + interface + subroutine amg_d_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + import amg_d_umf_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_umf_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_dpk_),intent(inout), optional :: initu(:) + end subroutine amg_d_umf_solver_apply + end interface + + interface + subroutine amg_d_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + import amg_d_umf_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_umf_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_d_vect_type),intent(inout), optional :: initu + end subroutine amg_d_umf_solver_apply_vect + end interface - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + interface + subroutine amg_d_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + import amg_d_umf_solver_type + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_umf_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_d_umf_solver_bld + end interface - call psb_erractionrestore(err_act) - return +contains -9999 call psb_error_handler(err_act) - return - end subroutine d_umf_solver_bld subroutine d_umf_solver_free(sv,info) diff --git a/amgprec/amg_s_slu_solver.F90 b/amgprec/amg_s_slu_solver.F90 index 68092510..6c6d19c6 100644 --- a/amgprec/amg_s_slu_solver.F90 +++ b/amgprec/amg_s_slu_solver.F90 @@ -63,9 +63,9 @@ module amg_s_slu_solver type(c_ptr) :: lufactors=c_null_ptr integer(c_long_long) :: symbsize=0, numsize=0 contains - procedure, pass(sv) :: build => s_slu_solver_bld - procedure, pass(sv) :: apply_a => s_slu_solver_apply - procedure, pass(sv) :: apply_v => s_slu_solver_apply_vect + procedure, pass(sv) :: build => amg_s_slu_solver_bld + procedure, pass(sv) :: apply_a => amg_s_slu_solver_apply + procedure, pass(sv) :: apply_v => amg_s_slu_solver_apply_vect procedure, pass(sv) :: free => s_slu_solver_free procedure, pass(sv) :: clear_data => s_slu_solver_clear_data procedure, pass(sv) :: descr => s_slu_solver_descr @@ -76,9 +76,8 @@ module amg_s_slu_solver end type amg_s_slu_solver_type - private :: s_slu_solver_bld, s_slu_solver_apply, & - & s_slu_solver_free, s_slu_solver_descr, & - & s_slu_solver_sizeof, s_slu_solver_apply_vect, & + private :: s_slu_solver_free, s_slu_solver_descr, & + & s_slu_solver_sizeof, & & s_slu_solver_get_fmt, s_slu_solver_get_id, & & s_slu_solver_clear_data private :: s_slu_solver_finalize @@ -118,207 +117,264 @@ module amg_s_slu_solver end function amg_sslu_free end interface -contains - - subroutine s_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& - & trans,work,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_s_slu_solver_type), intent(inout) :: sv - real(psb_spk_),intent(inout) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - character, intent(in), optional :: init - real(psb_spk_),intent(inout), optional :: initu(:) - - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:) - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act - character :: trans_ - character(len=20) :: name='s_slu_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - 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) - 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 - endif - - ww(1:n_row) = x(1:n_row) - select case(trans_) - case('N') - info = amg_sslu_solve(0,n_row,1,ww,n_row,sv%lufactors) - case('T') - info = amg_sslu_solve(1,n_row,1,ww,n_row,sv%lufactors) - case('C') - info = amg_sslu_solve(2,n_row,1,ww,n_row,sv%lufactors) - case default - call psb_errpush(psb_err_internal_error_, & - & name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - if (info == psb_success_) & - & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,& - & name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col > size(work)) then - deallocate(ww) - endif - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine s_slu_solver_apply - - subroutine s_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& - & trans,work,wv,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_s_slu_solver_type), intent(inout) :: sv - type(psb_s_vect_type),intent(inout) :: x - type(psb_s_vect_type),intent(inout) :: y - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - type(psb_s_vect_type),intent(inout) :: wv(:) - integer, intent(out) :: info - character, intent(in), optional :: init - type(psb_s_vect_type),intent(inout), optional :: initu - - integer :: err_act - character(len=20) :: name='s_slu_solver_apply_vect' - - call psb_erractionsave(err_act) - - info = psb_success_ - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - call x%v%sync() - call y%v%sync() - call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) - call y%v%set_host() - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine s_slu_solver_apply_vect - - subroutine s_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) - - use psb_base_mod - - Implicit None + interface + subroutine amg_s_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + import amg_s_slu_solver_type + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_s_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_s_slu_solver_bld + end interface - ! Arguments - type(psb_sspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(inout) :: desc_a - class(amg_s_slu_solver_type), intent(inout) :: sv - integer, intent(out) :: info - type(psb_sspmat_type), intent(in), target, optional :: b - class(psb_s_base_sparse_mat), intent(in), optional :: amold - class(psb_s_base_vect_type), intent(in), optional :: vmold - class(psb_i_base_vect_type), intent(in), optional :: imold - ! Local variables - type(psb_sspmat_type) :: atmp - type(psb_s_csc_sparse_mat) :: acsc - type(psb_s_coo_sparse_mat) :: acoo - integer :: n_row,n_col, nrow_a, nztota - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='s_slu_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ctxt = desc_a%get_context() - call psb_info(ctxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - - - call a%cscnv(atmp,info,type='coo') - call psb_rwextd(n_row,atmp,info,b=b) - call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) - nrow_a = atmp%get_nrows() - call atmp%a%csclip(acoo,info,jmax=nrow_a) - call acsc%mv_from_coo(acoo,info) - nztota = acsc%get_nzeros() - ! Fix the entries to call C-base SuperLU - acsc%ia(:) = acsc%ia(:) - 1 - acsc%icp(:) = acsc%icp(:) - 1 - info = amg_sslu_fact(nrow_a,nztota,acsc%val,& - & acsc%icp,acsc%ia,sv%lufactors) - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='amg_sslu_fact' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call acsc%free() - call atmp%free() + interface + subroutine amg_s_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + import amg_s_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_slu_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_s_vect_type),intent(inout), optional :: initu + end subroutine amg_s_slu_solver_apply_vect + end interface - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + interface + subroutine amg_s_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + import amg_s_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_slu_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_spk_),intent(inout), optional :: initu(:) + end subroutine amg_s_slu_solver_apply + end interface - call psb_erractionrestore(err_act) - return +contains -9999 call psb_error_handler(err_act) - return - end subroutine s_slu_solver_bld +!!$ subroutine s_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_s_slu_solver_type), intent(inout) :: sv +!!$ real(psb_spk_),intent(inout) :: x(:) +!!$ real(psb_spk_),intent(inout) :: y(:) +!!$ real(psb_spk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ real(psb_spk_),target, intent(inout) :: work(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ real(psb_spk_),intent(inout), optional :: initu(:) +!!$ +!!$ integer :: n_row,n_col +!!$ real(psb_spk_), pointer :: ww(:) +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act +!!$ character :: trans_ +!!$ character(len=20) :: name='s_slu_solver_apply' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ +!!$ trans_ = psb_toupper(trans) +!!$ select case(trans_) +!!$ case('N') +!!$ case('T','C') +!!$ case default +!!$ call psb_errpush(psb_err_iarg_invalid_i_,name) +!!$ goto 9999 +!!$ end select +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ 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) +!!$ 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 +!!$ endif +!!$ +!!$ ww(1:n_row) = x(1:n_row) +!!$ select case(trans_) +!!$ case('N') +!!$ info = amg_sslu_solve(0,n_row,1,ww,n_row,sv%lufactors) +!!$ case('T') +!!$ info = amg_sslu_solve(1,n_row,1,ww,n_row,sv%lufactors) +!!$ case('C') +!!$ info = amg_sslu_solve(2,n_row,1,ww,n_row,sv%lufactors) +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_, & +!!$ & name,a_err='Invalid TRANS in ILU subsolve') +!!$ goto 9999 +!!$ end select +!!$ +!!$ if (info == psb_success_) & +!!$ & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) +!!$ +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in subsolve') +!!$ goto 9999 +!!$ endif +!!$ +!!$ if (n_col > size(work)) then +!!$ deallocate(ww) +!!$ endif +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine s_slu_solver_apply +!!$ +!!$ subroutine s_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,wv,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_s_slu_solver_type), intent(inout) :: sv +!!$ type(psb_s_vect_type),intent(inout) :: x +!!$ type(psb_s_vect_type),intent(inout) :: y +!!$ real(psb_spk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ real(psb_spk_),target, intent(inout) :: work(:) +!!$ type(psb_s_vect_type),intent(inout) :: wv(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ type(psb_s_vect_type),intent(inout), optional :: initu +!!$ +!!$ integer :: err_act +!!$ character(len=20) :: name='s_slu_solver_apply_vect' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ call x%v%sync() +!!$ call y%v%sync() +!!$ call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) +!!$ call y%v%set_host() +!!$ if (info /= 0) goto 9999 +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine s_slu_solver_apply_vect +!!$ +!!$ subroutine s_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) +!!$ +!!$ use psb_base_mod +!!$ +!!$ Implicit None +!!$ +!!$ ! Arguments +!!$ type(psb_sspmat_type), intent(in), target :: a +!!$ Type(psb_desc_type), Intent(inout) :: desc_a +!!$ class(amg_s_slu_solver_type), intent(inout) :: sv +!!$ integer, intent(out) :: info +!!$ type(psb_sspmat_type), intent(in), target, optional :: b +!!$ class(psb_s_base_sparse_mat), intent(in), optional :: amold +!!$ class(psb_s_base_vect_type), intent(in), optional :: vmold +!!$ class(psb_i_base_vect_type), intent(in), optional :: imold +!!$ ! Local variables +!!$ type(psb_sspmat_type) :: atmp +!!$ type(psb_s_csc_sparse_mat) :: acsc +!!$ type(psb_s_coo_sparse_mat) :: acoo +!!$ integer :: n_row,n_col, nrow_a, nztota +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act, debug_unit, debug_level +!!$ character(len=20) :: name='s_slu_solver_bld', ch_err +!!$ +!!$ info=psb_success_ +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ ctxt = desc_a%get_context() +!!$ call psb_info(ctxt, me, np) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' start' +!!$ +!!$ +!!$ n_row = desc_a%get_local_rows() +!!$ n_col = desc_a%get_local_cols() +!!$ +!!$ +!!$ call a%cscnv(atmp,info,type='coo') +!!$ call psb_rwextd(n_row,atmp,info,b=b) +!!$ call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ nrow_a = atmp%get_nrows() +!!$ call atmp%a%csclip(acoo,info,jmax=nrow_a) +!!$ call acsc%mv_from_coo(acoo,info) +!!$ nztota = acsc%get_nzeros() +!!$ ! Fix the entries to call C-base SuperLU +!!$ acsc%ia(:) = acsc%ia(:) - 1 +!!$ acsc%icp(:) = acsc%icp(:) - 1 +!!$ info = amg_sslu_fact(nrow_a,nztota,acsc%val,& +!!$ & acsc%icp,acsc%ia,sv%lufactors) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='amg_sslu_fact' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ +!!$ call acsc%free() +!!$ call atmp%free() +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' end' +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ end subroutine s_slu_solver_bld subroutine s_slu_solver_free(sv,info) diff --git a/amgprec/amg_z_slu_solver.F90 b/amgprec/amg_z_slu_solver.F90 index b5fe5b2a..14f63b11 100644 --- a/amgprec/amg_z_slu_solver.F90 +++ b/amgprec/amg_z_slu_solver.F90 @@ -63,9 +63,9 @@ module amg_z_slu_solver type(c_ptr) :: lufactors=c_null_ptr integer(c_long_long) :: symbsize=0, numsize=0 contains - procedure, pass(sv) :: build => z_slu_solver_bld - procedure, pass(sv) :: apply_a => z_slu_solver_apply - procedure, pass(sv) :: apply_v => z_slu_solver_apply_vect + procedure, pass(sv) :: build => amg_z_slu_solver_bld + procedure, pass(sv) :: apply_a => amg_z_slu_solver_apply + procedure, pass(sv) :: apply_v => amg_z_slu_solver_apply_vect procedure, pass(sv) :: free => z_slu_solver_free procedure, pass(sv) :: clear_data => z_slu_solver_clear_data procedure, pass(sv) :: descr => z_slu_solver_descr @@ -76,9 +76,8 @@ module amg_z_slu_solver end type amg_z_slu_solver_type - private :: z_slu_solver_bld, z_slu_solver_apply, & - & z_slu_solver_free, z_slu_solver_descr, & - & z_slu_solver_sizeof, z_slu_solver_apply_vect, & + private :: z_slu_solver_free, z_slu_solver_descr, & + & z_slu_solver_sizeof, & & z_slu_solver_get_fmt, z_slu_solver_get_id, & & z_slu_solver_clear_data private :: z_slu_solver_finalize @@ -118,207 +117,264 @@ module amg_z_slu_solver end function amg_zslu_free end interface -contains - - subroutine z_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& - & trans,work,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_z_slu_solver_type), intent(inout) :: sv - complex(psb_dpk_),intent(inout) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - character, intent(in), optional :: init - complex(psb_dpk_),intent(inout), optional :: initu(:) - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:) - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_slu_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - 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) - 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 - endif - - ww(1:n_row) = x(1:n_row) - select case(trans_) - case('N') - info = amg_zslu_solve(0,n_row,1,ww,n_row,sv%lufactors) - case('T') - info = amg_zslu_solve(1,n_row,1,ww,n_row,sv%lufactors) - case('C') - info = amg_zslu_solve(2,n_row,1,ww,n_row,sv%lufactors) - case default - call psb_errpush(psb_err_internal_error_, & - & name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - if (info == psb_success_) & - & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,& - & name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col > size(work)) then - deallocate(ww) - endif - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine z_slu_solver_apply - - subroutine z_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& - & trans,work,wv,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_z_slu_solver_type), intent(inout) :: sv - type(psb_z_vect_type),intent(inout) :: x - type(psb_z_vect_type),intent(inout) :: y - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - type(psb_z_vect_type),intent(inout) :: wv(:) - integer, intent(out) :: info - character, intent(in), optional :: init - type(psb_z_vect_type),intent(inout), optional :: initu - - integer :: err_act - character(len=20) :: name='z_slu_solver_apply_vect' - - call psb_erractionsave(err_act) - - info = psb_success_ - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - call x%v%sync() - call y%v%sync() - call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) - call y%v%set_host() - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine z_slu_solver_apply_vect - - subroutine z_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) - - use psb_base_mod - - Implicit None + interface + subroutine amg_z_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + import amg_z_slu_solver_type + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_z_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_z_slu_solver_bld + end interface - ! Arguments - type(psb_zspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(inout) :: desc_a - class(amg_z_slu_solver_type), intent(inout) :: sv - integer, intent(out) :: info - type(psb_zspmat_type), intent(in), target, optional :: b - class(psb_z_base_sparse_mat), intent(in), optional :: amold - class(psb_z_base_vect_type), intent(in), optional :: vmold - class(psb_i_base_vect_type), intent(in), optional :: imold - ! Local variables - type(psb_zspmat_type) :: atmp - type(psb_z_csc_sparse_mat) :: acsc - type(psb_z_coo_sparse_mat) :: acoo - integer :: n_row,n_col, nrow_a, nztota - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='z_slu_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ctxt = desc_a%get_context() - call psb_info(ctxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - - - call a%cscnv(atmp,info,type='coo') - call psb_rwextd(n_row,atmp,info,b=b) - call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) - nrow_a = atmp%get_nrows() - call atmp%a%csclip(acoo,info,jmax=nrow_a) - call acsc%mv_from_coo(acoo,info) - nztota = acsc%get_nzeros() - ! Fix the entries to call C-base SuperLU - acsc%ia(:) = acsc%ia(:) - 1 - acsc%icp(:) = acsc%icp(:) - 1 - info = amg_zslu_fact(nrow_a,nztota,acsc%val,& - & acsc%icp,acsc%ia,sv%lufactors) - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='amg_zslu_fact' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call acsc%free() - call atmp%free() + interface + subroutine amg_z_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + import amg_z_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_slu_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_z_vect_type),intent(inout), optional :: initu + end subroutine amg_z_slu_solver_apply_vect + end interface - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + interface + subroutine amg_z_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + import amg_z_slu_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_slu_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_dpk_),intent(inout), optional :: initu(:) + end subroutine amg_z_slu_solver_apply + end interface - call psb_erractionrestore(err_act) - return +contains -9999 call psb_error_handler(err_act) - return - end subroutine z_slu_solver_bld +!!$ subroutine z_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_z_slu_solver_type), intent(inout) :: sv +!!$ complex(psb_dpk_),intent(inout) :: x(:) +!!$ complex(psb_dpk_),intent(inout) :: y(:) +!!$ complex(psb_dpk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ complex(psb_dpk_),target, intent(inout) :: work(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ complex(psb_dpk_),intent(inout), optional :: initu(:) +!!$ +!!$ integer :: n_row,n_col +!!$ complex(psb_dpk_), pointer :: ww(:) +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act +!!$ character :: trans_ +!!$ character(len=20) :: name='z_slu_solver_apply' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ +!!$ trans_ = psb_toupper(trans) +!!$ select case(trans_) +!!$ case('N') +!!$ case('T','C') +!!$ case default +!!$ call psb_errpush(psb_err_iarg_invalid_i_,name) +!!$ goto 9999 +!!$ end select +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ 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) +!!$ 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 +!!$ endif +!!$ +!!$ ww(1:n_row) = x(1:n_row) +!!$ select case(trans_) +!!$ case('N') +!!$ info = amg_zslu_solve(0,n_row,1,ww,n_row,sv%lufactors) +!!$ case('T') +!!$ info = amg_zslu_solve(1,n_row,1,ww,n_row,sv%lufactors) +!!$ case('C') +!!$ info = amg_zslu_solve(2,n_row,1,ww,n_row,sv%lufactors) +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_, & +!!$ & name,a_err='Invalid TRANS in ILU subsolve') +!!$ goto 9999 +!!$ end select +!!$ +!!$ if (info == psb_success_) & +!!$ & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) +!!$ +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in subsolve') +!!$ goto 9999 +!!$ endif +!!$ +!!$ if (n_col > size(work)) then +!!$ deallocate(ww) +!!$ endif +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine z_slu_solver_apply +!!$ +!!$ subroutine z_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& +!!$ & trans,work,wv,info,init,initu) +!!$ use psb_base_mod +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_z_slu_solver_type), intent(inout) :: sv +!!$ type(psb_z_vect_type),intent(inout) :: x +!!$ type(psb_z_vect_type),intent(inout) :: y +!!$ complex(psb_dpk_),intent(in) :: alpha,beta +!!$ character(len=1),intent(in) :: trans +!!$ complex(psb_dpk_),target, intent(inout) :: work(:) +!!$ type(psb_z_vect_type),intent(inout) :: wv(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: init +!!$ type(psb_z_vect_type),intent(inout), optional :: initu +!!$ +!!$ integer :: err_act +!!$ character(len=20) :: name='z_slu_solver_apply_vect' +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = psb_success_ +!!$ ! +!!$ ! For non-iterative solvers, init and initu are ignored. +!!$ ! +!!$ +!!$ call x%v%sync() +!!$ call y%v%sync() +!!$ call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) +!!$ call y%v%set_host() +!!$ if (info /= 0) goto 9999 +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ +!!$ end subroutine z_slu_solver_apply_vect +!!$ +!!$ subroutine z_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) +!!$ +!!$ use psb_base_mod +!!$ +!!$ Implicit None +!!$ +!!$ ! Arguments +!!$ type(psb_zspmat_type), intent(in), target :: a +!!$ Type(psb_desc_type), Intent(inout) :: desc_a +!!$ class(amg_z_slu_solver_type), intent(inout) :: sv +!!$ integer, intent(out) :: info +!!$ type(psb_zspmat_type), intent(in), target, optional :: b +!!$ class(psb_z_base_sparse_mat), intent(in), optional :: amold +!!$ class(psb_z_base_vect_type), intent(in), optional :: vmold +!!$ class(psb_i_base_vect_type), intent(in), optional :: imold +!!$ ! Local variables +!!$ type(psb_zspmat_type) :: atmp +!!$ type(psb_z_csc_sparse_mat) :: acsc +!!$ type(psb_z_coo_sparse_mat) :: acoo +!!$ integer :: n_row,n_col, nrow_a, nztota +!!$ type(psb_ctxt_type) :: ctxt +!!$ integer :: np,me,i, err_act, debug_unit, debug_level +!!$ character(len=20) :: name='z_slu_solver_bld', ch_err +!!$ +!!$ info=psb_success_ +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ ctxt = desc_a%get_context() +!!$ call psb_info(ctxt, me, np) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' start' +!!$ +!!$ +!!$ n_row = desc_a%get_local_rows() +!!$ n_col = desc_a%get_local_cols() +!!$ +!!$ +!!$ call a%cscnv(atmp,info,type='coo') +!!$ call psb_rwextd(n_row,atmp,info,b=b) +!!$ call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ nrow_a = atmp%get_nrows() +!!$ call atmp%a%csclip(acoo,info,jmax=nrow_a) +!!$ call acsc%mv_from_coo(acoo,info) +!!$ nztota = acsc%get_nzeros() +!!$ ! Fix the entries to call C-base SuperLU +!!$ acsc%ia(:) = acsc%ia(:) - 1 +!!$ acsc%icp(:) = acsc%icp(:) - 1 +!!$ info = amg_zslu_fact(nrow_a,nztota,acsc%val,& +!!$ & acsc%icp,acsc%ia,sv%lufactors) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='amg_zslu_fact' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ +!!$ call acsc%free() +!!$ call atmp%free() +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),' end' +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(err_act) +!!$ return +!!$ end subroutine z_slu_solver_bld subroutine z_slu_solver_free(sv,info) diff --git a/amgprec/amg_z_umf_solver.F90 b/amgprec/amg_z_umf_solver.F90 index 133cc08a..2ca86267 100644 --- a/amgprec/amg_z_umf_solver.F90 +++ b/amgprec/amg_z_umf_solver.F90 @@ -62,9 +62,9 @@ module amg_z_umf_solver type(c_ptr) :: symbolic=c_null_ptr, numeric=c_null_ptr integer(c_long_long) :: symbsize=0, numsize=0 contains - procedure, pass(sv) :: build => z_umf_solver_bld - procedure, pass(sv) :: apply_a => z_umf_solver_apply - procedure, pass(sv) :: apply_v => z_umf_solver_apply_vect + procedure, pass(sv) :: build => amg_z_umf_solver_bld + procedure, pass(sv) :: apply_a => amg_z_umf_solver_apply + procedure, pass(sv) :: apply_v => amg_z_umf_solver_apply_vect procedure, pass(sv) :: free => z_umf_solver_free procedure, pass(sv) :: clear_data => z_umf_solver_clear_data procedure, pass(sv) :: descr => z_umf_solver_descr @@ -75,9 +75,8 @@ module amg_z_umf_solver end type amg_z_umf_solver_type - private :: z_umf_solver_bld, z_umf_solver_apply, & - & z_umf_solver_free, z_umf_solver_descr, & - & z_umf_solver_sizeof, z_umf_solver_apply_vect, & + private :: z_umf_solver_free, z_umf_solver_descr, & + & z_umf_solver_sizeof, & & z_umf_solver_get_fmt, z_umf_solver_get_id, & & z_umf_solver_clear_data private :: z_umf_solver_finalize @@ -118,208 +117,65 @@ module amg_z_umf_solver end function amg_zumf_free end interface -contains - - subroutine z_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& - & trans,work,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_z_umf_solver_type), intent(inout) :: sv - complex(psb_dpk_),intent(inout) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - character, intent(in), optional :: init - complex(psb_dpk_),intent(inout), optional :: initu(:) - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:) - integer(psb_ipk_) :: i, err_act - character :: trans_ - character(len=20) :: name='z_umf_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - 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) - 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 - endif - - select case(trans_) - case('N') - info = amg_zumf_solve(0,n_row,ww,x,n_row,sv%numeric) - case('T') - ! - ! Note: with UMF, 1 meand Ctranspose, 2 means transpose - ! even for complex data. - ! - if (psb_z_is_complex_) then - info = amg_zumf_solve(2,n_row,ww,x,n_row,sv%numeric) - else - info = amg_zumf_solve(1,n_row,ww,x,n_row,sv%numeric) - end if - case('C') - info = amg_zumf_solve(1,n_row,ww,x,n_row,sv%numeric) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col > size(work)) then - deallocate(ww) - endif - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine z_umf_solver_apply - - subroutine z_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& - & trans,work,wv,info,init,initu) - use psb_base_mod - implicit none - type(psb_desc_type), intent(in) :: desc_data - class(amg_z_umf_solver_type), intent(inout) :: sv - type(psb_z_vect_type),intent(inout) :: x - type(psb_z_vect_type),intent(inout) :: y - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - type(psb_z_vect_type),intent(inout) :: wv(:) - integer, intent(out) :: info - character, intent(in), optional :: init - type(psb_z_vect_type),intent(inout), optional :: initu - - integer :: err_act - character(len=20) :: name='z_umf_solver_apply_vect' - - call psb_erractionsave(err_act) - - info = psb_success_ - ! - ! For non-iterative solvers, init and initu are ignored. - ! - - call x%v%sync() - call y%v%sync() - call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) - call y%v%set_host() - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - return - - end subroutine z_umf_solver_apply_vect - - subroutine z_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) - - use psb_base_mod - - Implicit None - - ! Arguments - type(psb_zspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(inout) :: desc_a - class(amg_z_umf_solver_type), intent(inout) :: sv - integer, intent(out) :: info - type(psb_zspmat_type), intent(in), target, optional :: b - class(psb_z_base_sparse_mat), intent(in), optional :: amold - class(psb_z_base_vect_type), intent(in), optional :: vmold - class(psb_i_base_vect_type), intent(in), optional :: imold - ! Local variables - type(psb_zspmat_type) :: atmp - type(psb_z_csc_sparse_mat) :: acsc - integer :: n_row,n_col, nrow_a, nztota - type(psb_ctxt_type) :: ctxt - integer :: np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='z_umf_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ctxt = desc_a%get_context() - call psb_info(ctxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - - call a%cscnv(atmp,info,type='coo') - call psb_rwextd(n_row,atmp,info,b=b) - call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_) - call atmp%mv_to(acsc) - nrow_a = acsc%get_nrows() - nztota = acsc%get_nzeros() - ! Fix the entres to call C-base UMFPACK. - acsc%ia(:) = acsc%ia(:) - 1 - acsc%icp(:) = acsc%icp(:) - 1 - info = amg_zumf_fact(nrow_a,nztota,acsc%val,& - & acsc%ia,acsc%icp,sv%symbolic,sv%numeric,& - & sv%symbsize,sv%numsize) - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='amg_zumf_fact' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call acsc%free() - call atmp%free() + interface + subroutine amg_z_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + import amg_z_umf_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_umf_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_dpk_),intent(inout), optional :: initu(:) + end subroutine amg_z_umf_solver_apply + end interface + + interface + subroutine amg_z_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + import amg_z_umf_solver_type + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_umf_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_z_vect_type),intent(inout), optional :: initu + end subroutine amg_z_umf_solver_apply_vect + end interface - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' + interface + subroutine amg_z_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + import amg_z_umf_solver_type + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_z_umf_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_z_umf_solver_bld + end interface - call psb_erractionrestore(err_act) - return +contains -9999 call psb_error_handler(err_act) - return - end subroutine z_umf_solver_bld subroutine z_umf_solver_free(sv,info) diff --git a/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 index 4a71212e..471e10e9 100644 --- a/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 @@ -80,7 +80,7 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,& & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use amg_c_prec_type - use amg_c_symdec_aggregator_mod, amg_protect_name => amg_c_symdec_aggregator_build_tprol + use amg_c_symdec_aggregator_mod, only : amg_c_symdec_aggregator_type use amg_c_inner_mod implicit none class(amg_c_symdec_aggregator_type), target, intent(inout) :: ag @@ -103,6 +103,16 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,& integer(psb_ipk_) :: debug_level, debug_unit logical :: clean_zeros +!!$ interface +!!$ subroutine psb_caplusat(ain,aout,info) +!!$ use psb_c_mat_mod, only : psb_cspmat_type +!!$ import :: psb_ipk_ +!!$ implicit none +!!$ type(psb_cspmat_type) :: ain, aout +!!$ integer(psb_ipk_) :: info +!!$ end subroutine psb_caplusat +!!$ end interface +!!$ name='amg_c_symdec_aggregator_tprol' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then @@ -123,15 +133,11 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,& call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) nr = a%get_nrows() - call a%csclip(atmp,info,imax=nr,jmax=nr,& + call a%csclip(atrans,info,imax=nr,jmax=nr,& & rscale=.false.,cscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) - if (info == psb_success_) call atmp%transp(atrans) - if (info == psb_success_) call atrans%cscnv(info,type='COO') - if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) + call atrans%set_nrows(nr) + call atrans%set_ncols(nr) + call psb_aplusat(atrans,atmp,info) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') diff --git a/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 index e4b2fd6e..9e4e2179 100644 --- a/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 @@ -80,7 +80,7 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,& & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use amg_d_prec_type - use amg_d_symdec_aggregator_mod, amg_protect_name => amg_d_symdec_aggregator_build_tprol + use amg_d_symdec_aggregator_mod, only : amg_d_symdec_aggregator_type use amg_d_inner_mod implicit none class(amg_d_symdec_aggregator_type), target, intent(inout) :: ag @@ -103,6 +103,16 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,& integer(psb_ipk_) :: debug_level, debug_unit logical :: clean_zeros +!!$ interface +!!$ subroutine psb_daplusat(ain,aout,info) +!!$ use psb_d_mat_mod, only : psb_dspmat_type +!!$ import :: psb_ipk_ +!!$ implicit none +!!$ type(psb_dspmat_type) :: ain, aout +!!$ integer(psb_ipk_) :: info +!!$ end subroutine psb_daplusat +!!$ end interface +!!$ name='amg_d_symdec_aggregator_tprol' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then @@ -123,15 +133,11 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,& call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) nr = a%get_nrows() - call a%csclip(atmp,info,imax=nr,jmax=nr,& + call a%csclip(atrans,info,imax=nr,jmax=nr,& & rscale=.false.,cscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) - if (info == psb_success_) call atmp%transp(atrans) - if (info == psb_success_) call atrans%cscnv(info,type='COO') - if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) + call atrans%set_nrows(nr) + call atrans%set_ncols(nr) + call psb_aplusat(atrans,atmp,info) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') diff --git a/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 index c148b89a..0c1d0011 100644 --- a/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 @@ -80,7 +80,7 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,& & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use amg_s_prec_type - use amg_s_symdec_aggregator_mod, amg_protect_name => amg_s_symdec_aggregator_build_tprol + use amg_s_symdec_aggregator_mod, only : amg_s_symdec_aggregator_type use amg_s_inner_mod implicit none class(amg_s_symdec_aggregator_type), target, intent(inout) :: ag @@ -103,6 +103,16 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,& integer(psb_ipk_) :: debug_level, debug_unit logical :: clean_zeros +!!$ interface +!!$ subroutine psb_saplusat(ain,aout,info) +!!$ use psb_s_mat_mod, only : psb_sspmat_type +!!$ import :: psb_ipk_ +!!$ implicit none +!!$ type(psb_sspmat_type) :: ain, aout +!!$ integer(psb_ipk_) :: info +!!$ end subroutine psb_saplusat +!!$ end interface +!!$ name='amg_s_symdec_aggregator_tprol' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then @@ -123,15 +133,11 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,& call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) nr = a%get_nrows() - call a%csclip(atmp,info,imax=nr,jmax=nr,& + call a%csclip(atrans,info,imax=nr,jmax=nr,& & rscale=.false.,cscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) - if (info == psb_success_) call atmp%transp(atrans) - if (info == psb_success_) call atrans%cscnv(info,type='COO') - if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) + call atrans%set_nrows(nr) + call atrans%set_ncols(nr) + call psb_aplusat(atrans,atmp,info) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') diff --git a/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 index dd4ac4be..02849031 100644 --- a/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 @@ -80,7 +80,7 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,& & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use amg_z_prec_type - use amg_z_symdec_aggregator_mod, amg_protect_name => amg_z_symdec_aggregator_build_tprol + use amg_z_symdec_aggregator_mod, only : amg_z_symdec_aggregator_type use amg_z_inner_mod implicit none class(amg_z_symdec_aggregator_type), target, intent(inout) :: ag @@ -103,6 +103,16 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,& integer(psb_ipk_) :: debug_level, debug_unit logical :: clean_zeros +!!$ interface +!!$ subroutine psb_zaplusat(ain,aout,info) +!!$ use psb_z_mat_mod, only : psb_zspmat_type +!!$ import :: psb_ipk_ +!!$ implicit none +!!$ type(psb_zspmat_type) :: ain, aout +!!$ integer(psb_ipk_) :: info +!!$ end subroutine psb_zaplusat +!!$ end interface +!!$ name='amg_z_symdec_aggregator_tprol' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then @@ -123,15 +133,11 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,& call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) nr = a%get_nrows() - call a%csclip(atmp,info,imax=nr,jmax=nr,& + call a%csclip(atrans,info,imax=nr,jmax=nr,& & rscale=.false.,cscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) - if (info == psb_success_) call atmp%transp(atrans) - if (info == psb_success_) call atrans%cscnv(info,type='COO') - if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - call atmp%set_nrows(nr) - call atmp%set_ncols(nr) + call atrans%set_nrows(nr) + call atrans%set_ncols(nr) + call psb_aplusat(atrans,atmp,info) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') diff --git a/amgprec/impl/amg_c_smoothers_bld.f90 b/amgprec/impl/amg_c_smoothers_bld.f90 index d684e389..34d4dc9b 100644 --- a/amgprec/impl/amg_c_smoothers_bld.f90 +++ b/amgprec/impl/amg_c_smoothers_bld.f90 @@ -81,8 +81,7 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - !use amg_c_inner_mod - use amg_c_prec_mod, amg_protect_name => amg_c_smoothers_bld + use amg_c_prec_type Implicit None diff --git a/amgprec/impl/amg_d_smoothers_bld.f90 b/amgprec/impl/amg_d_smoothers_bld.f90 index c2bf3b99..ed2fd393 100644 --- a/amgprec/impl/amg_d_smoothers_bld.f90 +++ b/amgprec/impl/amg_d_smoothers_bld.f90 @@ -81,8 +81,7 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - !use amg_d_inner_mod - use amg_d_prec_mod, amg_protect_name => amg_d_smoothers_bld + use amg_d_prec_type Implicit None diff --git a/amgprec/impl/amg_s_smoothers_bld.f90 b/amgprec/impl/amg_s_smoothers_bld.f90 index 00385c8a..9cd5f00c 100644 --- a/amgprec/impl/amg_s_smoothers_bld.f90 +++ b/amgprec/impl/amg_s_smoothers_bld.f90 @@ -81,8 +81,7 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - !use amg_s_inner_mod - use amg_s_prec_mod, amg_protect_name => amg_s_smoothers_bld + use amg_s_prec_type Implicit None diff --git a/amgprec/impl/amg_z_smoothers_bld.f90 b/amgprec/impl/amg_z_smoothers_bld.f90 index eeaef05b..f55b3f12 100644 --- a/amgprec/impl/amg_z_smoothers_bld.f90 +++ b/amgprec/impl/amg_z_smoothers_bld.f90 @@ -81,8 +81,7 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - !use amg_z_inner_mod - use amg_z_prec_mod, amg_protect_name => amg_z_smoothers_bld + use amg_z_prec_type Implicit None diff --git a/amgprec/impl/solver/Makefile b/amgprec/impl/solver/Makefile index 56942182..72f22e0a 100644 --- a/amgprec/impl/solver/Makefile +++ b/amgprec/impl/solver/Makefile @@ -54,6 +54,7 @@ amg_c_mumps_solver_apply.o \ amg_c_mumps_solver_apply_vect.o \ amg_c_mumps_solver_bld.o \ amg_c_krm_solver_impl.o \ +amg_c_slu_solver_impl.o \ amg_d_base_solver_apply.o \ amg_d_base_solver_apply_vect.o \ amg_d_base_solver_bld.o \ @@ -101,6 +102,8 @@ amg_d_mumps_solver_apply.o \ amg_d_mumps_solver_apply_vect.o \ amg_d_mumps_solver_bld.o \ amg_d_krm_solver_impl.o \ +amg_d_slu_solver_impl.o \ +amg_d_umf_solver_impl.o \ amg_s_base_solver_apply.o \ amg_s_base_solver_apply_vect.o \ amg_s_base_solver_bld.o \ @@ -148,6 +151,7 @@ amg_s_mumps_solver_apply.o \ amg_s_mumps_solver_apply_vect.o \ amg_s_mumps_solver_bld.o \ amg_s_krm_solver_impl.o \ +amg_s_slu_solver_impl.o \ amg_z_base_solver_apply.o \ amg_z_base_solver_apply_vect.o \ amg_z_base_solver_bld.o \ @@ -303,6 +307,8 @@ amg_z_invk_solver_clone_settings.o \ amg_z_invk_solver_cseti.o \ amg_z_invk_solver_descr.o \ amg_z_krm_solver_impl.o \ +amg_z_slu_solver_impl.o \ +amg_z_umf_solver_impl.o \ amg_c_jac_solver_clone_settings.o \ amg_c_jac_solver_clear_data.o \ amg_c_jac_solver_cnv.o \ diff --git a/amgprec/impl/solver/amg_c_slu_solver_impl.F90 b/amgprec/impl/solver/amg_c_slu_solver_impl.F90 new file mode 100644 index 00000000..1d5c772e --- /dev/null +++ b/amgprec/impl/solver/amg_c_slu_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_c_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + + use psb_base_mod + use amg_c_slu_solver, amg_protect_name => amg_c_slu_solver_bld + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_c_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_cspmat_type) :: atmp + type(psb_c_csc_sparse_mat) :: acsc + type(psb_c_coo_sparse_mat) :: acoo + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_slu_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) + nrow_a = atmp%get_nrows() + call atmp%a%csclip(acoo,info,jmax=nrow_a) + call acsc%mv_from_coo(acoo,info) + nztota = acsc%get_nzeros() + ! Fix the entries to call C-base SuperLU + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_cslu_fact(nrow_a,nztota,acsc%val,& + & acsc%icp,acsc%ia,sv%lufactors) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_cslu_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_c_slu_solver_bld + +subroutine amg_c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_c_slu_solver, amg_protect_name => amg_c_slu_solver_apply_vect + + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_slu_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_c_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='s_slu_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_c_slu_solver_apply_vect + +subroutine amg_c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_c_slu_solver, amg_protect_name => amg_c_slu_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_slu_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_spk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + complex(psb_spk_), pointer :: ww(:) + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_slu_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + ww(1:n_row) = x(1:n_row) + select case(trans_) + case('N') + info = amg_cslu_solve(0,n_row,1,ww,n_row,sv%lufactors) + case('T') + info = amg_cslu_solve(1,n_row,1,ww,n_row,sv%lufactors) + case('C') + info = amg_cslu_solve(2,n_row,1,ww,n_row,sv%lufactors) + case default + call psb_errpush(psb_err_internal_error_, & + & name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) & + & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_c_slu_solver_apply diff --git a/amgprec/impl/solver/amg_c_umf_solver_impl.F90 b/amgprec/impl/solver/amg_c_umf_solver_impl.F90 new file mode 100644 index 00000000..081b23d0 --- /dev/null +++ b/amgprec/impl/solver/amg_c_umf_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_c_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_c_umf_solver, amg_protect_name => amg_c_umf_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_umf_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_spk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + complex(psb_spk_), pointer :: ww(:) + integer(psb_ipk_) :: i, err_act + character :: trans_ + character(len=20) :: name='c_umf_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + select case(trans_) + case('N') + info = amg_cumf_solve(0,n_row,ww,x,n_row,sv%numeric) + case('T') + ! + ! Note: with UMF, 1 meand Ctranspose, 2 means transpose + ! even for complex data. + ! + if (psb_c_is_complex_) then + info = amg_cumf_solve(2,n_row,ww,x,n_row,sv%numeric) + else + info = amg_cumf_solve(1,n_row,ww,x,n_row,sv%numeric) + end if + case('C') + info = amg_cumf_solve(1,n_row,ww,x,n_row,sv%numeric) + case default + call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_c_umf_solver_apply + +subroutine amg_c_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_c_umf_solver, amg_protect_name => amg_c_umf_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_umf_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_c_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='c_umf_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_c_umf_solver_apply_vect + +subroutine amg_c_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + use amg_c_umf_solver, amg_protect_name => amg_c_umf_solver_bld + + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_c_umf_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_cspmat_type) :: atmp + type(psb_c_csc_sparse_mat) :: acsc + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='c_umf_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_) + call atmp%mv_to(acsc) + nrow_a = acsc%get_nrows() + nztota = acsc%get_nzeros() + ! Fix the entres to call C-base UMFPACK. + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_cumf_fact(nrow_a,nztota,acsc%val,& + & acsc%ia,acsc%icp,sv%symbolic,sv%numeric,& + & sv%symbsize,sv%numsize) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_cumf_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_c_umf_solver_bld diff --git a/amgprec/impl/solver/amg_d_slu_solver_impl.F90 b/amgprec/impl/solver/amg_d_slu_solver_impl.F90 new file mode 100644 index 00000000..5f8f18df --- /dev/null +++ b/amgprec/impl/solver/amg_d_slu_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_d_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + + use psb_base_mod + use amg_d_slu_solver, amg_protect_name => amg_d_slu_solver_bld + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_dspmat_type) :: atmp + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_slu_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) + nrow_a = atmp%get_nrows() + call atmp%a%csclip(acoo,info,jmax=nrow_a) + call acsc%mv_from_coo(acoo,info) + nztota = acsc%get_nzeros() + ! Fix the entries to call C-base SuperLU + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_dslu_fact(nrow_a,nztota,acsc%val,& + & acsc%icp,acsc%ia,sv%lufactors) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_dslu_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_d_slu_solver_bld + +subroutine amg_d_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_d_slu_solver, amg_protect_name => amg_d_slu_solver_apply_vect + + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_slu_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_d_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='s_slu_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_d_slu_solver_apply_vect + +subroutine amg_d_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_d_slu_solver, amg_protect_name => amg_d_slu_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_slu_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_dpk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:) + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_slu_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + ww(1:n_row) = x(1:n_row) + select case(trans_) + case('N') + info = amg_dslu_solve(0,n_row,1,ww,n_row,sv%lufactors) + case('T') + info = amg_dslu_solve(1,n_row,1,ww,n_row,sv%lufactors) + case('C') + info = amg_dslu_solve(2,n_row,1,ww,n_row,sv%lufactors) + case default + call psb_errpush(psb_err_internal_error_, & + & name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) & + & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_d_slu_solver_apply diff --git a/amgprec/impl/solver/amg_d_umf_solver_impl.F90 b/amgprec/impl/solver/amg_d_umf_solver_impl.F90 new file mode 100644 index 00000000..615a255c --- /dev/null +++ b/amgprec/impl/solver/amg_d_umf_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_d_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_d_umf_solver, amg_protect_name => amg_d_umf_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_umf_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_dpk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:) + integer(psb_ipk_) :: i, err_act + character :: trans_ + character(len=20) :: name='d_umf_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + select case(trans_) + case('N') + info = amg_dumf_solve(0,n_row,ww,x,n_row,sv%numeric) + case('T') + ! + ! Note: with UMF, 1 meand Ctranspose, 2 means transpose + ! even for complex data. + ! + if (psb_d_is_complex_) then + info = amg_dumf_solve(2,n_row,ww,x,n_row,sv%numeric) + else + info = amg_dumf_solve(1,n_row,ww,x,n_row,sv%numeric) + end if + case('C') + info = amg_dumf_solve(1,n_row,ww,x,n_row,sv%numeric) + case default + call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_d_umf_solver_apply + +subroutine amg_d_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_d_umf_solver, amg_protect_name => amg_d_umf_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_umf_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_d_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='d_umf_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_d_umf_solver_apply_vect + +subroutine amg_d_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + use amg_d_umf_solver, amg_protect_name => amg_d_umf_solver_bld + + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_umf_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_dspmat_type) :: atmp + type(psb_d_csc_sparse_mat) :: acsc + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='d_umf_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_) + call atmp%mv_to(acsc) + nrow_a = acsc%get_nrows() + nztota = acsc%get_nzeros() + ! Fix the entres to call C-base UMFPACK. + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_dumf_fact(nrow_a,nztota,acsc%val,& + & acsc%ia,acsc%icp,sv%symbolic,sv%numeric,& + & sv%symbsize,sv%numsize) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_dumf_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_d_umf_solver_bld diff --git a/amgprec/impl/solver/amg_s_slu_solver_impl.F90 b/amgprec/impl/solver/amg_s_slu_solver_impl.F90 new file mode 100644 index 00000000..c6d0383d --- /dev/null +++ b/amgprec/impl/solver/amg_s_slu_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_s_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + + use psb_base_mod + use amg_s_slu_solver, amg_protect_name => amg_s_slu_solver_bld + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_s_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_sspmat_type) :: atmp + type(psb_s_csc_sparse_mat) :: acsc + type(psb_s_coo_sparse_mat) :: acoo + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_slu_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) + nrow_a = atmp%get_nrows() + call atmp%a%csclip(acoo,info,jmax=nrow_a) + call acsc%mv_from_coo(acoo,info) + nztota = acsc%get_nzeros() + ! Fix the entries to call C-base SuperLU + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_sslu_fact(nrow_a,nztota,acsc%val,& + & acsc%icp,acsc%ia,sv%lufactors) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_sslu_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_s_slu_solver_bld + +subroutine amg_s_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_s_slu_solver, amg_protect_name => amg_s_slu_solver_apply_vect + + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_slu_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_s_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='s_slu_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_s_slu_solver_apply_vect + +subroutine amg_s_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_s_slu_solver, amg_protect_name => amg_s_slu_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_slu_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_spk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + real(psb_spk_), pointer :: ww(:) + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_slu_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + ww(1:n_row) = x(1:n_row) + select case(trans_) + case('N') + info = amg_sslu_solve(0,n_row,1,ww,n_row,sv%lufactors) + case('T') + info = amg_sslu_solve(1,n_row,1,ww,n_row,sv%lufactors) + case('C') + info = amg_sslu_solve(2,n_row,1,ww,n_row,sv%lufactors) + case default + call psb_errpush(psb_err_internal_error_, & + & name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) & + & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_s_slu_solver_apply diff --git a/amgprec/impl/solver/amg_s_umf_solver_impl.F90 b/amgprec/impl/solver/amg_s_umf_solver_impl.F90 new file mode 100644 index 00000000..bb7d904c --- /dev/null +++ b/amgprec/impl/solver/amg_s_umf_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_s_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_s_umf_solver, amg_protect_name => amg_s_umf_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_umf_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_spk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + real(psb_spk_), pointer :: ww(:) + integer(psb_ipk_) :: i, err_act + character :: trans_ + character(len=20) :: name='s_umf_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + select case(trans_) + case('N') + info = amg_sumf_solve(0,n_row,ww,x,n_row,sv%numeric) + case('T') + ! + ! Note: with UMF, 1 meand Ctranspose, 2 means transpose + ! even for complex data. + ! + if (psb_s_is_complex_) then + info = amg_sumf_solve(2,n_row,ww,x,n_row,sv%numeric) + else + info = amg_sumf_solve(1,n_row,ww,x,n_row,sv%numeric) + end if + case('C') + info = amg_sumf_solve(1,n_row,ww,x,n_row,sv%numeric) + case default + call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_s_umf_solver_apply + +subroutine amg_s_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_s_umf_solver, amg_protect_name => amg_s_umf_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_umf_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_s_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='s_umf_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_s_umf_solver_apply_vect + +subroutine amg_s_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + use amg_s_umf_solver, amg_protect_name => amg_s_umf_solver_bld + + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_s_umf_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_sspmat_type) :: atmp + type(psb_s_csc_sparse_mat) :: acsc + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_umf_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_) + call atmp%mv_to(acsc) + nrow_a = acsc%get_nrows() + nztota = acsc%get_nzeros() + ! Fix the entres to call C-base UMFPACK. + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_sumf_fact(nrow_a,nztota,acsc%val,& + & acsc%ia,acsc%icp,sv%symbolic,sv%numeric,& + & sv%symbsize,sv%numsize) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_sumf_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_s_umf_solver_bld diff --git a/amgprec/impl/solver/amg_z_slu_solver_impl.F90 b/amgprec/impl/solver/amg_z_slu_solver_impl.F90 new file mode 100644 index 00000000..fb219935 --- /dev/null +++ b/amgprec/impl/solver/amg_z_slu_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_z_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + + use psb_base_mod + use amg_z_slu_solver, amg_protect_name => amg_z_slu_solver_bld + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_z_slu_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_zspmat_type) :: atmp + type(psb_z_csc_sparse_mat) :: acsc + type(psb_z_coo_sparse_mat) :: acoo + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_slu_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_) + nrow_a = atmp%get_nrows() + call atmp%a%csclip(acoo,info,jmax=nrow_a) + call acsc%mv_from_coo(acoo,info) + nztota = acsc%get_nzeros() + ! Fix the entries to call C-base SuperLU + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_zslu_fact(nrow_a,nztota,acsc%val,& + & acsc%icp,acsc%ia,sv%lufactors) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_zslu_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_z_slu_solver_bld + +subroutine amg_z_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_z_slu_solver, amg_protect_name => amg_z_slu_solver_apply_vect + + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_slu_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_z_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='s_slu_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_z_slu_solver_apply_vect + +subroutine amg_z_slu_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_z_slu_solver, amg_protect_name => amg_z_slu_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_slu_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_dpk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + complex(psb_dpk_), pointer :: ww(:) + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_slu_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + ww(1:n_row) = x(1:n_row) + select case(trans_) + case('N') + info = amg_zslu_solve(0,n_row,1,ww,n_row,sv%lufactors) + case('T') + info = amg_zslu_solve(1,n_row,1,ww,n_row,sv%lufactors) + case('C') + info = amg_zslu_solve(2,n_row,1,ww,n_row,sv%lufactors) + case default + call psb_errpush(psb_err_internal_error_, & + & name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) & + & call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_z_slu_solver_apply diff --git a/amgprec/impl/solver/amg_z_umf_solver_impl.F90 b/amgprec/impl/solver/amg_z_umf_solver_impl.F90 new file mode 100644 index 00000000..3ae4cd75 --- /dev/null +++ b/amgprec/impl/solver/amg_z_umf_solver_impl.F90 @@ -0,0 +1,202 @@ +subroutine amg_z_umf_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_z_umf_solver, amg_protect_name => amg_z_umf_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_umf_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_dpk_),intent(inout), optional :: initu(:) + + integer :: n_row,n_col + complex(psb_dpk_), pointer :: ww(:) + integer(psb_ipk_) :: i, err_act + character :: trans_ + character(len=20) :: name='z_umf_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + 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) + 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 + endif + + select case(trans_) + case('N') + info = amg_zumf_solve(0,n_row,ww,x,n_row,sv%numeric) + case('T') + ! + ! Note: with UMF, 1 meand Ctranspose, 2 means transpose + ! even for complex data. + ! + if (psb_z_is_complex_) then + info = amg_zumf_solve(2,n_row,ww,x,n_row,sv%numeric) + else + info = amg_zumf_solve(1,n_row,ww,x,n_row,sv%numeric) + end if + case('C') + info = amg_zumf_solve(1,n_row,ww,x,n_row,sv%numeric) + case default + call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') + goto 9999 + endif + + if (n_col > size(work)) then + deallocate(ww) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_z_umf_solver_apply + +subroutine amg_z_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + use psb_base_mod + use amg_z_umf_solver, amg_protect_name => amg_z_umf_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_umf_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_z_vect_type),intent(inout), optional :: initu + + integer :: err_act + character(len=20) :: name='z_umf_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + ! + ! For non-iterative solvers, init and initu are ignored. + ! + + call x%v%sync() + call y%v%sync() + call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) + call y%v%set_host() + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine amg_z_umf_solver_apply_vect + +subroutine amg_z_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + use psb_base_mod + use amg_z_umf_solver, amg_protect_name => amg_z_umf_solver_bld + + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_z_umf_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_zspmat_type) :: atmp + type(psb_z_csc_sparse_mat) :: acsc + integer :: n_row,n_col, nrow_a, nztota + type(psb_ctxt_type) :: ctxt + integer :: np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='z_umf_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + + call a%cscnv(atmp,info,type='coo') + call psb_rwextd(n_row,atmp,info,b=b) + call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_) + call atmp%mv_to(acsc) + nrow_a = acsc%get_nrows() + nztota = acsc%get_nzeros() + ! Fix the entres to call C-base UMFPACK. + acsc%ia(:) = acsc%ia(:) - 1 + acsc%icp(:) = acsc%icp(:) - 1 + info = amg_zumf_fact(nrow_a,nztota,acsc%val,& + & acsc%ia,acsc%icp,sv%symbolic,sv%numeric,& + & sv%symbsize,sv%numsize) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_zumf_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call acsc%free() + call atmp%free() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_z_umf_solver_bld diff --git a/cbind/amgprec/amg_dprec_cbind_mod.F90 b/cbind/amgprec/amg_dprec_cbind_mod.F90 index 1c2ca8bd..eddbef41 100644 --- a/cbind/amgprec/amg_dprec_cbind_mod.F90 +++ b/cbind/amgprec/amg_dprec_cbind_mod.F90 @@ -21,17 +21,14 @@ contains !#define MLDC_ERR_FILTER(INFO) min(0,INFO) #define MLDC_ERR_FILTER(INFO) (INFO) #define MLDC_ERR_HANDLE(INFO) if(INFO/=amg_success_)MLDC_ERROR("ERROR!") - function amg_c_dprecinit(cctxt,ph,ptype) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(amg_c_dprec) :: ph type(psb_c_object_type), value :: cctxt character(c_char) :: ptype(*) - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_dprec_type), pointer :: precp character(len=80) :: fptype @@ -41,30 +38,28 @@ contains return end if - allocate(precp,stat=info) - if (info /= 0) return + allocate(precp,stat=iret) + if (iret /= 0) return ph%item = c_loc(precp) call stringc2f(ptype,fptype) - call precp%init(psb_c2f_ctxt(cctxt),fptype,info) + call precp%init(psb_c2f_ctxt(cctxt),fptype,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dprecinit function amg_c_dprecseti(ph,what,val) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph character(c_char) :: what(*) integer(psb_c_ipk_), value :: val - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret character(len=80) :: fwhat type(amg_dprec_type), pointer :: precp @@ -77,24 +72,22 @@ contains call stringc2f(what,fwhat) - call precp%set(fwhat,val,info) + call precp%set(fwhat,val,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dprecseti function amg_c_dprecsetr(ph,what,val) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph character(c_char) :: what(*) real(c_double), value :: val - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret character(len=80) :: fwhat type(amg_dprec_type), pointer :: precp @@ -107,22 +100,20 @@ contains call stringc2f(what,fwhat) - call precp%set(fwhat,val,info) + call precp%set(fwhat,val,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dprecsetr function amg_c_dprecsetc(ph,what,val) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph character(c_char) :: what(*), val(*) - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret character(len=80) :: fwhat,fval type(amg_dprec_type), pointer :: precp @@ -136,25 +127,23 @@ contains call stringc2f(what,fwhat) call stringc2f(val,fval) - call precp%set(fwhat,fval,info) + call precp%set(fwhat,fval,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dprecsetc function amg_c_dprecbld(ah,cdh,ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph,ah,cdh - integer(psb_ipk_) :: info type(amg_dprec_type), pointer :: precp type(psb_dspmat_type), pointer :: ap type(psb_desc_type), pointer :: descp character(len=80) :: fptype + integer(psb_ipk_) :: iret res = -1 @@ -174,22 +163,20 @@ contains return end if - call amg_precbld(ap,descp,precp,info) + call amg_precbld(ap,descp,precp,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dprecbld function amg_c_dhierarchy_build(ah,cdh,ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph,ah,cdh - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_dprec_type), pointer :: precp type(psb_dspmat_type), pointer :: ap type(psb_desc_type), pointer :: descp @@ -213,26 +200,24 @@ contains return end if - call precp%hierarchy_build(ap,descp,info) + call precp%hierarchy_build(ap,descp,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dhierarchy_build function amg_c_dsmoothers_build(ah,cdh,ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph,ah,cdh - integer(psb_ipk_) :: info type(amg_dprec_type), pointer :: precp type(psb_dspmat_type), pointer :: ap type(psb_desc_type), pointer :: descp character(len=80) :: fptype + integer(psb_ipk_) :: iret res = -1 @@ -252,9 +237,9 @@ contains return end if - call precp%smoothers_build(ap,descp,info) + call precp%smoothers_build(ap,descp,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return @@ -302,7 +287,7 @@ contains type(amg_dprec_type), pointer :: precp type(psb_d_vect_type), pointer :: xp, bp - integer(psb_ipk_) :: info,fitmax,fitrace,first,fistop,fiter + integer(psb_ipk_) :: iret,fitmax,fitrace,first,fistop,fiter character(len=20) :: fmethd real(kind(1.d0)) :: feps,ferr @@ -342,23 +327,21 @@ contains fistop = istop call psb_krylov(fmethd, ap, precp, bp, xp, feps, & - & descp, info,& + & descp, iret,& & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,& & irst=first, err=ferr) iter = fiter err = ferr - res = min(info,0) + res = min(iret,0) end function amg_c_dkrylov_opt function amg_c_dprecfree(ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_dprec_type), pointer :: precp character(len=80) :: fptype @@ -370,39 +353,37 @@ contains end if - call precp%free(info) + call precp%free(iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_dprecfree function amg_c_ddescr(ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod - implicit none - - integer(psb_c_ipk_) :: res - type(psb_c_object_type) :: ph - integer(psb_c_ipk_) :: info - type(amg_dprec_type), pointer :: precp - - res = -1 - info = -1 - if (c_associated(ph%item)) then - call c_f_pointer(ph%item,precp) - else - return - end if - - - call precp%descr(info) - call flush(psb_out_unit) - - info = 0 - res = MLDC_ERR_FILTER(info) - MLDC_ERR_HANDLE(res) - return - end function amg_c_ddescr + implicit none + + integer(psb_c_ipk_) :: res + type(psb_c_object_type) :: ph + integer(psb_c_ipk_) :: iret + type(amg_dprec_type), pointer :: precp + + res = -1 + iret = -1 + if (c_associated(ph%item)) then + call c_f_pointer(ph%item,precp) + else + return + end if + + + call precp%descr(iret) + call flush(psb_out_unit) + + iret = 0 + res = MLDC_ERR_FILTER(iret) + MLDC_ERR_HANDLE(res) + return + end function amg_c_ddescr end module amg_dprec_cbind_mod diff --git a/cbind/amgprec/amg_zprec_cbind_mod.F90 b/cbind/amgprec/amg_zprec_cbind_mod.F90 index a58dd4d1..2e8a8b94 100644 --- a/cbind/amgprec/amg_zprec_cbind_mod.F90 +++ b/cbind/amgprec/amg_zprec_cbind_mod.F90 @@ -23,15 +23,13 @@ contains #define MLDC_ERR_HANDLE(INFO) if(INFO/=amg_success_)MLDC_ERROR("ERROR!") function amg_c_zprecinit(cctxt,ph,ptype) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(amg_c_zprec) :: ph type(psb_c_object_type), value :: cctxt character(c_char) :: ptype(*) - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_zprec_type), pointer :: precp character(len=80) :: fptype @@ -41,30 +39,28 @@ contains return end if - allocate(precp,stat=info) - if (info /= 0) return + allocate(precp,stat=iret) + if (iret /= 0) return ph%item = c_loc(precp) call stringc2f(ptype,fptype) - call precp%init(psb_c2f_ctxt(cctxt),fptype,info) + call precp%init(psb_c2f_ctxt(cctxt),fptype,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zprecinit function amg_c_zprecseti(ph,what,val) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph character(c_char) :: what(*) integer(psb_c_ipk_), value :: val - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret character(len=80) :: fwhat type(amg_zprec_type), pointer :: precp @@ -77,24 +73,22 @@ contains call stringc2f(what,fwhat) - call precp%set(fwhat,val,info) + call precp%set(fwhat,val,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zprecseti function amg_c_zprecsetr(ph,what,val) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph character(c_char) :: what(*) real(c_double), value :: val - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret character(len=80) :: fwhat type(amg_zprec_type), pointer :: precp @@ -107,22 +101,20 @@ contains call stringc2f(what,fwhat) - call precp%set(fwhat,val,info) + call precp%set(fwhat,val,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zprecsetr function amg_c_zprecsetc(ph,what,val) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph character(c_char) :: what(*), val(*) - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret character(len=80) :: fwhat,fval type(amg_zprec_type), pointer :: precp @@ -136,21 +128,19 @@ contains call stringc2f(what,fwhat) call stringc2f(val,fval) - call precp%set(fwhat,fval,info) + call precp%set(fwhat,fval,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zprecsetc function amg_c_zprecbld(ah,cdh,ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph,ah,cdh - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_zprec_type), pointer :: precp type(psb_zspmat_type), pointer :: ap type(psb_desc_type), pointer :: descp @@ -174,22 +164,20 @@ contains return end if - call amg_precbld(ap,descp,precp,info) + call amg_precbld(ap,descp,precp,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zprecbld function amg_c_zhierarchy_build(ah,cdh,ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph,ah,cdh - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_zprec_type), pointer :: precp type(psb_zspmat_type), pointer :: ap type(psb_desc_type), pointer :: descp @@ -213,22 +201,20 @@ contains return end if - call precp%hierarchy_build(ap,descp,info) + call precp%hierarchy_build(ap,descp,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zhierarchy_build function amg_c_zsmoothers_build(ah,cdh,ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph,ah,cdh - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_zprec_type), pointer :: precp type(psb_zspmat_type), pointer :: ap type(psb_desc_type), pointer :: descp @@ -252,9 +238,9 @@ contains return end if - call precp%smoothers_build(ap,descp,info) + call precp%smoothers_build(ap,descp,iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return @@ -262,8 +248,6 @@ contains function amg_c_zkrylov(methd,& & ah,ph,bh,xh,cdh,options) bind(c) result(res) - use psb_base_mod - use psb_prec_mod use psb_linsolve_mod use psb_prec_cbind_mod use psb_zkrylov_cbind_mod @@ -283,11 +267,8 @@ contains function amg_c_zkrylov_opt(methd,& & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res) - use psb_base_mod - use psb_prec_mod use psb_linsolve_mod use psb_objhandle_mod - use psb_prec_cbind_mod use psb_base_string_cbind_mod implicit none integer(psb_c_ipk_) :: res @@ -302,7 +283,7 @@ contains type(amg_zprec_type), pointer :: precp type(psb_z_vect_type), pointer :: xp, bp - integer(psb_ipk_) :: info,fitmax,fitrace,first,fistop,fiter + integer(psb_ipk_) :: iret,fitmax,fitrace,first,fistop,fiter character(len=20) :: fmethd real(kind(1.d0)) :: feps,ferr @@ -342,23 +323,21 @@ contains fistop = istop call psb_krylov(fmethd, ap, precp, bp, xp, feps, & - & descp, info,& + & descp, iret,& & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,& & irst=first, err=ferr) iter = fiter err = ferr - res = min(info,0) + res = min(iret,0) end function amg_c_zkrylov_opt function amg_c_zprecfree(ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph - integer(psb_ipk_) :: info + integer(psb_ipk_) :: iret type(amg_zprec_type), pointer :: precp character(len=80) :: fptype @@ -370,25 +349,23 @@ contains end if - call precp%free(info) + call precp%free(iret) - res = MLDC_ERR_FILTER(info) + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zprecfree function amg_c_zdescr(ph) bind(c) result(res) - use psb_base_mod - use amg_prec_mod implicit none integer(psb_c_ipk_) :: res type(psb_c_object_type) :: ph - integer(psb_c_ipk_) :: info + integer(psb_c_ipk_) :: iret type(amg_zprec_type), pointer :: precp res = -1 - info = -1 + iret = -1 if (c_associated(ph%item)) then call c_f_pointer(ph%item,precp) else @@ -396,11 +373,11 @@ contains end if - call precp%descr(info) + call precp%descr(iret) call flush(psb_out_unit) - info = 0 - res = MLDC_ERR_FILTER(info) + iret = 0 + res = MLDC_ERR_FILTER(iret) MLDC_ERR_HANDLE(res) return end function amg_c_zdescr