diff --git a/test/kernel/vecoperation.f90 b/test/kernel/vecoperation.f90 index e13fc37a..5c87768b 100644 --- a/test/kernel/vecoperation.f90 +++ b/test/kernel/vecoperation.f90 @@ -33,102 +33,815 @@ ! module unittestvector_mod - use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& - & psb_dspmat_type, psb_d_vect_type, dzero, psb_ctxt_type,& - & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type,& - & psb_d_base_multivect_type + use psb_base_mod interface psb_gen_const - module procedure psb_d_gen_const, psb_d_gen_const_multi + module procedure psb_d_gen_const, psb_s_gen_const, & + & psb_c_gen_const, psb_z_gen_const, & + & psb_d_gen_const_multi, psb_s_gen_const_multi end interface psb_gen_const interface psb_check_ans - module procedure psb_check_ans_v, psb_check_ans_mv + module procedure psb_d_check_ans_v, psb_c_check_ans_v, & + & psb_z_check_ans_v, psb_s_check_ans_v, & + & psb_d_check_ans_mv, psb_s_check_ans_mv, & + & psb_c_check_ans_mv, psb_z_check_ans_mv end interface psb_check_ans -contains +contains + + function psb_d_check_ans_v(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_d_vect_type) :: v + real(psb_dpk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_dpk_) :: check + real(psb_dpk_), allocatable :: va(:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(va); + + call psb_sum(ctxt,check) + + if(check == 0.d0) then + ans = .true. + else + ans = .false. + end if + + end function psb_d_check_ans_v + + function psb_s_check_ans_v(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_s_vect_type) :: v + real(psb_spk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_spk_) :: check + real(psb_spk_), allocatable :: va(:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(va); + + call psb_sum(ctxt,check) + + if(check == 0.e0) then + ans = .true. + else + ans = .false. + end if + + end function psb_s_check_ans_v + + function psb_c_check_ans_v(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_c_vect_type) :: v + complex(psb_spk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_spk_) :: check + complex(psb_spk_), allocatable :: va(:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(abs(va)); + + call psb_sum(ctxt,check) + + if(check == 0.e0) then + ans = .true. + else + ans = .false. + end if + + end function psb_c_check_ans_v + + function psb_z_check_ans_v(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_z_vect_type) :: v + complex(psb_dpk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_dpk_) :: check + complex(psb_dpk_), allocatable :: va(:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(abs(va)); + + call psb_sum(ctxt,check) + + if(check == 0.d0) then + ans = .true. + else + ans = .false. + end if + + end function psb_z_check_ans_v + + function psb_d_check_ans_mv(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_d_multivect_type) :: v + real(psb_dpk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_dpk_) :: check + real(psb_dpk_), allocatable :: va(:,:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(va); + + call psb_sum(ctxt,check) + + if(check == 0.d0) then + ans = .true. + else + ans = .false. + end if + + end function psb_d_check_ans_mv + + function psb_s_check_ans_mv(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_s_multivect_type) :: v + real(psb_spk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_spk_) :: check + real(psb_spk_), allocatable :: va(:,:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(va); + + call psb_sum(ctxt,check) + + if(check == 0.e0) then + ans = .true. + else + ans = .false. + end if + + end function psb_s_check_ans_mv + + function psb_c_check_ans_mv(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_c_multivect_type) :: v + complex(psb_spk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_spk_) :: check + complex(psb_spk_), allocatable :: va(:,:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(abs(va)); + + call psb_sum(ctxt,check) + + if(check == 0.e0) then + ans = .true. + else + ans = .false. + end if + + end function psb_c_check_ans_mv + + function psb_z_check_ans_mv(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_z_multivect_type) :: v + complex(psb_dpk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_dpk_) :: check + complex(psb_dpk_), allocatable :: va(:,:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(abs(va)); + + call psb_sum(ctxt,check) + + if(check == 0.d0) then + ans = .true. + else + ans = .false. + end if + + end function psb_z_check_ans_mv + ! + ! subroutine to fill a vector with constant entries + ! + subroutine psb_z_gen_const(v,val,idim,ctxt,desc_a,info) + use psb_base_mod + implicit none + + type(psb_z_vect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: info + complex(psb_dpk_) :: val + + ! Local variables + integer(psb_ipk_), parameter :: nb=20 + complex(psb_dpk_) :: zt(nb) + character(len=20) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) + + + info = psb_success_ + name = 'create_constant_vector' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + n = idim*np ! The global dimension is the number of process times + ! the input size + + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr + + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),v,desc_a,info) + if(info /= psb_success_) exit + end do + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + end subroutine psb_z_gen_const + ! + ! subroutine to fill a vector with constant entries + ! + subroutine psb_c_gen_const(v,val,idim,ctxt,desc_a,info) + use psb_base_mod + implicit none + + type(psb_c_vect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: info + complex(psb_spk_) :: val + + ! Local variables + integer(psb_ipk_), parameter :: nb=20 + complex(psb_spk_) :: zt(nb) + character(len=20) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) + + + info = psb_success_ + name = 'create_constant_vector' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + n = idim*np ! The global dimension is the number of process times + ! the input size + + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr + + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),v,desc_a,info) + if(info /= psb_success_) exit + end do + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + end subroutine psb_c_gen_const + ! + ! subroutine to fill a vector with constant entries + ! + subroutine psb_s_gen_const(v,val,idim,ctxt,desc_a,info) + use psb_base_mod + implicit none + + type(psb_s_vect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: info + real(psb_spk_) :: val + + ! Local variables + integer(psb_ipk_), parameter :: nb=20 + real(psb_spk_) :: zt(nb) + character(len=20) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) + + + info = psb_success_ + name = 'create_constant_vector' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + n = idim*np ! The global dimension is the number of process times + ! the input size + + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr + + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),v,desc_a,info) + if(info /= psb_success_) exit + end do + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + end subroutine psb_s_gen_const + ! + ! subroutine to fill a vector with constant entries + ! + subroutine psb_d_gen_const(v,val,idim,ctxt,desc_a,info) + use psb_base_mod + implicit none + + type(psb_d_vect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: info + real(psb_dpk_) :: val + + ! Local variables + integer(psb_ipk_), parameter :: nb=20 + real(psb_dpk_) :: zt(nb) + character(len=20) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) + + + info = psb_success_ + name = 'create_constant_vector' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + n = idim*np ! The global dimension is the number of process times + ! the input size + + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr + + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),v,desc_a,info) + if(info /= psb_success_) exit + end do + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + end subroutine psb_d_gen_const + ! + ! subroutine to fill a multivectorvector with constant entries + ! + subroutine psb_d_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info) + use psb_base_mod + implicit none + + type(psb_d_multivect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + integer(psb_ipk_) :: jdim !! Number of columns of the multivector + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: info + real(psb_dpk_) :: val + + ! Local variables + integer(psb_ipk_), parameter :: nb=20 + real(psb_dpk_) :: zt(nb,jdim) ! Temporary array to fill the vector + character(len=40) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) + + + info = psb_success_ + name = 'create_constant_multivector' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + n = idim*np ! The global dimension is the number of process times + ! the input size + + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr + + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info,n=jdim) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:,:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib,1:jdim),v,desc_a,info) + if(info /= psb_success_) exit + end do + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) + + call psb_erractionrestore(err_act) + return - function psb_check_ans_v(v,val,ctxt) result(ans) - use psb_base_mod +9999 call psb_error_handler(ctxt,err_act) + return + end subroutine psb_d_gen_const_multi + ! + ! subroutine to fill a multivectorvector with constant entries + ! + subroutine psb_s_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info) + use psb_base_mod implicit none - type(psb_d_vect_type) :: v - real(psb_dpk_) :: val + type(psb_s_multivect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + integer(psb_ipk_) :: jdim !! Number of columns of the multivector type(psb_ctxt_type) :: ctxt - logical :: ans + integer(psb_ipk_) :: info + real(psb_spk_) :: val ! Local variables - integer(psb_ipk_) :: np, iam, info - real(psb_dpk_) :: check - real(psb_dpk_), allocatable :: va(:) - - call psb_info(ctxt,iam,np) - - va = v%get_vect() - va = va - val; - - check = maxval(va); + integer(psb_ipk_), parameter :: nb=20 + real(psb_spk_) :: zt(nb,jdim) ! Temporary array to fill the vector + character(len=40) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) - call psb_sum(ctxt,check) - if(check == 0.d0) then - ans = .true. - else - ans = .false. - end if + info = psb_success_ + name = 'create_constant_multivector' + call psb_erractionsave(err_act) - end function psb_check_ans_v + call psb_info(ctxt, iam, np) - function psb_check_ans_mv(v,val,ctxt) result(ans) - use psb_base_mod + n = idim*np ! The global dimension is the number of process times + ! the input size - implicit none + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr - type(psb_d_multivect_type) :: v - real(psb_dpk_) :: val - type(psb_ctxt_type) :: ctxt - logical :: ans + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info,n=jdim) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - ! Local variables - integer(psb_ipk_) :: np, iam, info - real(psb_dpk_) :: check - real(psb_dpk_), allocatable :: va(:,:) + myidx = desc_a%get_global_indices() + nlr = size(myidx) - call psb_info(ctxt,iam,np) + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:,:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib,1:jdim),v,desc_a,info) + if(info /= psb_success_) exit + end do - va = v%get_vect() - va = va - val; + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - check = maxval(va); + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) - call psb_sum(ctxt,check) + call psb_erractionrestore(err_act) + return - if(check == 0.d0) then - ans = .true. - else - ans = .false. - end if +9999 call psb_error_handler(ctxt,err_act) - end function psb_check_ans_mv - ! - ! subroutine to fill a vector with constant entries + return + end subroutine psb_s_gen_const_multi + ! + ! subroutine to fill a multivectorvector with constant entries ! - subroutine psb_d_gen_const(v,val,idim,ctxt,desc_a,info) + subroutine psb_c_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info) use psb_base_mod implicit none - type(psb_d_vect_type) :: v + type(psb_c_multivect_type) :: v type(psb_desc_type) :: desc_a integer(psb_lpk_) :: idim + integer(psb_ipk_) :: jdim !! Number of columns of the multivector type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: info - real(psb_dpk_) :: val + complex(psb_spk_) :: val ! Local variables integer(psb_ipk_), parameter :: nb=20 - real(psb_dpk_) :: zt(nb) - character(len=20) :: name, ch_err + complex(psb_spk_) :: zt(nb,jdim) ! Temporary array to fill the vector + character(len=40) :: name, ch_err integer(psb_ipk_) :: np, iam, nr, nt integer(psb_ipk_) :: n,nlr,ib,ii integer(psb_ipk_) :: err_act @@ -136,7 +849,7 @@ contains info = psb_success_ - name = 'create_constant_vector' + name = 'create_constant_multivector' call psb_erractionsave(err_act) call psb_info(ctxt, iam, np) @@ -160,7 +873,7 @@ contains ! Allocate the descriptor with simple minded data distribution call psb_cdall(ctxt,desc_a,info,nl=nr) ! Allocate the vector on the recently build descriptor - if (info == psb_success_) call psb_geall(v,desc_a,info) + if (info == psb_success_) call psb_geall(v,desc_a,info,n=jdim) ! Check that allocation has gone good if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -174,8 +887,8 @@ contains do ii=1,nlr,nb ib = min(nb,nlr-ii+1) - zt(:) = val - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),v,desc_a,info) + zt(:,:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib,1:jdim),v,desc_a,info) if(info /= psb_success_) exit end do @@ -196,25 +909,25 @@ contains 9999 call psb_error_handler(ctxt,err_act) return - end subroutine psb_d_gen_const + end subroutine psb_c_gen_const_multi ! ! subroutine to fill a multivectorvector with constant entries ! - subroutine psb_d_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info) + subroutine psb_z_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info) use psb_base_mod implicit none - type(psb_d_multivect_type) :: v + type(psb_z_multivect_type) :: v type(psb_desc_type) :: desc_a integer(psb_lpk_) :: idim integer(psb_ipk_) :: jdim !! Number of columns of the multivector type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: info - real(psb_dpk_) :: val + complex(psb_dpk_) :: val ! Local variables integer(psb_ipk_), parameter :: nb=20 - real(psb_dpk_) :: zt(nb,jdim) ! Temporary array to fill the vector + complex(psb_dpk_) :: zt(nb,jdim) ! Temporary array to fill the vector character(len=40) :: name, ch_err integer(psb_ipk_) :: np, iam, nr, nt integer(psb_ipk_) :: n,nlr,ib,ii @@ -283,7 +996,7 @@ contains 9999 call psb_error_handler(ctxt,err_act) return - end subroutine psb_d_gen_const_multi + end subroutine psb_z_gen_const_multi end module unittestvector_mod @@ -305,12 +1018,36 @@ program vecoperation real(psb_dpk_), parameter :: negativeone = -1.d0 real(psb_dpk_), parameter :: negativetwo = -2.d0 real(psb_dpk_), parameter :: negativeonehalf = -0.5_psb_dpk_ + + real(psb_spk_), parameter :: stwo = 2.e0 + real(psb_spk_), parameter :: sonehalf = 0.5_psb_spk_ + real(psb_spk_), parameter :: snegativeone = -1.e0 + real(psb_spk_), parameter :: snegativetwo = -2.e0 + real(psb_spk_), parameter :: snegativeonehalf = -0.5_psb_spk_ + + complex(psb_spk_), parameter :: ctwo = (2.e0, 0.e0) + complex(psb_spk_), parameter :: cnegativeone = (-1.e0, 0.e0) + complex(psb_spk_), parameter :: cnegativetwo = (-2.e0, 0.e0) + complex(psb_spk_), parameter :: conehalf = (0.5e0, 0.e0) + + complex(psb_dpk_), parameter :: ztwo = (2.d0, 0.d0) + complex(psb_dpk_), parameter :: znegativeone = (-1.d0, 0.d0) + complex(psb_dpk_), parameter :: znegativetwo = (-2.d0, 0.d0) + complex(psb_dpk_), parameter :: zonehalf = (0.5d0, 0.d0) + + ! descriptor type(psb_desc_type) :: desc_a ! vector - type(psb_d_vect_type) :: x,y,z + type(psb_d_vect_type) :: x, y, z + type(psb_s_vect_type) :: sx, sy, sz + type(psb_c_vect_type) :: cx, cy, cz + type(psb_z_vect_type) :: zx, zy, zz ! multivector type(psb_d_multivect_type) :: mv1, mv2 + type(psb_s_multivect_type) :: smv1, smv2 + type(psb_c_multivect_type) :: cmv1, cmv2 + type(psb_z_multivect_type) :: zmv1, zmv2 ! blacs parameters type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np @@ -354,8 +1091,10 @@ program vecoperation ! ! Test of standard vector operation ! + if (iam == psb_root_) write(psb_out_unit,'("---- Standard Vector Operations ----")') if (iam == psb_root_) write(psb_out_unit,'(" ")') - if (iam == psb_root_) write(psb_out_unit,'("Standard Vector Operations")') + if (iam == psb_root_) write(psb_out_unit,'(" ")') + if (iam == psb_root_) write(psb_out_unit,'(" Real Double Precision")') if (iam == psb_root_) write(psb_out_unit,'(" ")') ! X = 1 call psb_d_gen_const(x,one,idim,ctxt,desc_a,info) @@ -411,16 +1150,6 @@ program vecoperation if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = X - Y")') if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = X - Y ")') end if - ! X = 2 , Y = 1, Z = 0, Z = -X + Y = -2 + 1 = -1 - call psb_d_gen_const(x,two,idim,ctxt,desc_a,info) - call psb_d_gen_const(y,one,idim,ctxt,desc_a,info) - call psb_d_gen_const(z,dzero,idim,ctxt,desc_a,info) - call psb_geaxpby(negativeone,x,one,y,z,desc_a,info) - hasitnotfailed = psb_check_ans(z,negativeone,ctxt) - if (iam == psb_root_) then - if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = -X + Y")') - if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = -X + Y ")') - end if ! X = 2 , Y = -0.5, Z = 0, Z = X*Y = 2*(-0.5) = -1 call psb_d_gen_const(x,two,idim,ctxt,desc_a,info) call psb_d_gen_const(y,negativeonehalf,idim,ctxt,desc_a,info) @@ -431,42 +1160,229 @@ program vecoperation if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> mlt Z = X*Y")') if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- mlt Z = X*Y ")') end if - ! X = 1 , Y = 2, Z = 0, Z = X/Y = 1/2 = 0.5 - call psb_d_gen_const(x,one,idim,ctxt,desc_a,info) - call psb_d_gen_const(y,two,idim,ctxt,desc_a,info) - call psb_d_gen_const(z,dzero,idim,ctxt,desc_a,info) - call psb_gediv(x,y,z,desc_a,info) - hasitnotfailed = psb_check_ans(z,onehalf,ctxt) + + ! Single Precision Real + if (iam == psb_root_) write(psb_out_unit,'(" ")') + if (iam == psb_root_) write(psb_out_unit,'(" Real Single Precision")') + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! X = 1 + call psb_s_gen_const(sx,sone,idim,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(sx,sone,ctxt) if (iam == psb_root_) then - if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> div Z = X/Y")') - if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- div Z = X/Y ")') + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant vector (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant vector (Single Precision)")') end if - ! X = -1 , Z = 0, Z = |X| = |-1| = 1 - call psb_d_gen_const(x,negativeone,idim,ctxt,desc_a,info) - call psb_d_gen_const(z,dzero,idim,ctxt,desc_a,info) - call psb_geabs(x,z,desc_a,info) - hasitnotfailed = psb_check_ans(z,one,ctxt) + ! X = 2 , Y = -2, Y = 0.5*X + Y = 1 - 2 = -1 + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,snegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(sonehalf,sx,sone,sy,desc_a,info) + hasitnotfailed = psb_check_ans(sy,snegativeone,ctxt) if (iam == psb_root_) then - if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> abs Z = |X|")') - if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- abs Z = |X| ")') + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = 0.5 X + Y (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = 0.5 X + Y (Single Precision)")') end if - ! X = 2 , Z = 0, Z = 1/X = 1/2 = 0.5 - call psb_d_gen_const(x,two,idim,ctxt,desc_a,info) - call psb_d_gen_const(z,dzero,idim,ctxt,desc_a,info) - call psb_geinv(x,z,desc_a,info) - hasitnotfailed = psb_check_ans(z,onehalf,ctxt) + ! X = 2 , Y = -2, Y = X + Y = 0 + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,snegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(sone,sx,sone,sy,desc_a,info) + hasitnotfailed = psb_check_ans(sy,(0.0_psb_spk_),ctxt) if (iam == psb_root_) then - if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> inv Z = 1/X")') - if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- inv Z = 1/X ")') + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = X + Y (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = X + Y (Single Precision)")') end if - ! X = 1, Z = 0, c = -2, Z = X + c = -1 - call psb_d_gen_const(x,one,idim,ctxt,desc_a,info) - call psb_d_gen_const(z,dzero,idim,ctxt,desc_a,info) - call psb_geaddconst(x,negativetwo,z,desc_a,info) - hasitnotfailed = psb_check_ans(z,negativeone,ctxt) + ! X = 2 , Y = -2, Y = 0.5*X + Y = 1 - 2 = -1 + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,snegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(sonehalf,sx,sone,sy,desc_a,info) + hasitnotfailed = psb_check_ans(sy,snegativeone,ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = 0.5 X + Y (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = 0.5 X + Y (Single Precision)")') + end if + ! X = 2 , Y = -2, Y = 0.5*X + Y = 1 - 2 = -1 + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,snegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(sonehalf,sx,sone,sy,desc_a,info) + hasitnotfailed = psb_check_ans(sy,snegativeone,ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = 0.5 X + Y (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = 0.5 X + Y (Single Precision)")') + end if + ! X = 2 , Y = 1, Z = 0, Z = X - Y = 2 - 1 = 1 + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,sone,idim,ctxt,desc_a,info) + call psb_s_gen_const(sz,szero,idim,ctxt,desc_a,info) + call psb_geaxpby(sone,sx,snegativeone,sy,sz,desc_a,info) + hasitnotfailed = psb_check_ans(sz,sone,ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = X - Y (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = X - Y (Single Precision)")') + end if + ! X = 2 , Y = -0.5, Z = 0, Z = X*Y = -1 + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,snegativeonehalf,idim,ctxt,desc_a,info) + call psb_s_gen_const(sz,szero,idim,ctxt,desc_a,info) + call psb_gemlt(sone,sx,sy,szero,sz,desc_a,info) + hasitnotfailed = psb_check_ans(sz,snegativeone,ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> mlt Z = X*Y (Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- mlt Z = X*Y (Single Precision)")') + end if + + ! Single Precision Complex + if (iam == psb_root_) write(psb_out_unit,'(" ")') + if (iam == psb_root_) write(psb_out_unit,'(" Complex Single Precision")') + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! X = 1 + 0i + call psb_c_gen_const(cx,(1.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(cx,(1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant vector (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant vector (Complex Single Precision)")') + end if + ! X = 2+0i , Y = -2+0i, Y = X + Y = 0+0i + call psb_c_gen_const(cx,ctwo,idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,cnegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(cone,cx,cone,cy,desc_a,info) + hasitnotfailed = psb_check_ans(cy,(0.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = X + Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = X + Y (Complex Single Precision)")') + end if + ! X = 1 , Y = -2, Y = X + Y = 1 -2 = -1 + call psb_c_gen_const(cx,(1.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,cnegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(cone,cx,cone,cy,desc_a,info) + hasitnotfailed = psb_check_ans(cy,(-1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = X + Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = X + Y (Complex Single Precision)")') + end if + ! X = 1 , Y = 2, Y = -X + Y = -1 +2 = 1 + call psb_c_gen_const(cx,(1.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,(2.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_geaxpby(cnegativeone,cx,cone,cy,desc_a,info) + hasitnotfailed = psb_check_ans(cy,(1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = -X + Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = -X + Y (Complex Single Precision)")') + end if + ! X = 2 , Y = -2, Y = 0.5*X + Y = 1 - 2 = -1 + call psb_c_gen_const(cx,ctwo,idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,cnegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(conehalf,cx,cone,cy,desc_a,info) + hasitnotfailed = psb_check_ans(cy,(-1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = 0.5 X + Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = 0.5 X + Y (Complex Single Precision)")') + end if + ! X = -2 , Y = 1, Z = 0, Z = X + Y = -2 + 1 = -1 + call psb_c_gen_const(cx,cnegativetwo,idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,(1.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_c_gen_const(cz,(0.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_geaxpby(cone,cx,cone,cy,cz,desc_a,info) + hasitnotfailed = psb_check_ans(cz,(-1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = X + Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = X + Y (Complex Single Precision)")') + end if + ! X = 2 , Y = 1, Z = 0, Z = X - Y = 2 - 1 = 1 + call psb_c_gen_const(cx,ctwo,idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,(1.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_c_gen_const(cz,(0.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_geaxpby(cone,cx,cnegativeone,cy,cz,desc_a,info) + hasitnotfailed = psb_check_ans(cz,(1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = X - Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = X - Y (Complex Single Precision)")') + end if + ! X = 2 , Y = -0.5, Z = 0, Z = X*Y = 2*(-0.5) = -1 + call psb_c_gen_const(cx,ctwo,idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,(snegativeonehalf,0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_c_gen_const(cz,(0.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_gemlt(cone,cx,cy,(0.0_psb_spk_, 0.0_psb_spk_),cz,desc_a,info) + hasitnotfailed = psb_check_ans(cz,(-1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> mlt Z = X*Y (Complex Single Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- mlt Z = X*Y (Complex Single Precision)")') + end if + + ! Double Precision Complex + if (iam == psb_root_) write(psb_out_unit,'(" ")') + if (iam == psb_root_) write(psb_out_unit,'(" Complex Double Precision")') + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! X = 1 + 0i + call psb_z_gen_const(zx,(1.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(zx,(1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant vector (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant vector (Complex Double Precision)")') + end if + ! X = 2+0i , Y = -2+0i, Y = X + Y = 0+0i + call psb_z_gen_const(zx,ztwo,idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,znegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(zone,zx,zone,zy,desc_a,info) + hasitnotfailed = psb_check_ans(zy,(0.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = X + Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = X + Y (Complex Double Precision)")') + end if + ! X = 1 , Y = -2, Y = X + Y = 1 -2 = -1 + call psb_z_gen_const(zx,(1.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,znegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(zone,zx,zone,zy,desc_a,info) + hasitnotfailed = psb_check_ans(zy,(-1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = X + Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = X + Y (Complex Double Precision)")') + end if + ! X = 1 , Y = 2, Y = -X + Y = -1 +2 = 1 + call psb_z_gen_const(zx,(1.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,(2.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_geaxpby(znegativeone,zx,zone,zy,desc_a,info) + hasitnotfailed = psb_check_ans(zy,(1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) if (iam == psb_root_) then - if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Add constant Z = X + c")') - if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Add constant Z = X + c")') + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = -X + Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = -X + Y (Complex Double Precision)")') + end if + ! X = 2 , Y = -2, Y = 0.5*X + Y = 1 - 2 = -1 + call psb_z_gen_const(zx,ztwo,idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,znegativetwo,idim,ctxt,desc_a,info) + call psb_geaxpby(zonehalf,zx,zone,zy,desc_a,info) + hasitnotfailed = psb_check_ans(zy,(-1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Y = 0.5 X + Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Y = 0.5 X + Y (Complex Double Precision)")') + end if + ! X = -2 , Y = 1, Z = 0, Z = X + Y = -2 + 1 = -1 + call psb_z_gen_const(zx,znegativetwo,idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,(1.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_z_gen_const(zz,(0.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_geaxpby(zone,zx,zone,zy,zz,desc_a,info) + hasitnotfailed = psb_check_ans(zz,(-1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = X + Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = X + Y (Complex Double Precision)")') + end if + ! X = 2 , Y = 1, Z = 0, Z = X - Y = 2 - 1 = 1 + call psb_z_gen_const(zx,ztwo,idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,(1.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_z_gen_const(zz,(0.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_geaxpby(zone,zx,znegativeone,zy,zz,desc_a,info) + hasitnotfailed = psb_check_ans(zz,(1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> axpby Z = X - Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- axpby Z = X - Y (Complex Double Precision)")') + end if + ! X = 2 , Y = -0.5, Z = 0, Z = X*Y = 2*(-0.5) = -1 + call psb_z_gen_const(zx,ztwo,idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,(snegativeonehalf,0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_z_gen_const(zz,(0.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_gemlt(zone,zx,zy,(0.0_psb_dpk_, 0.0_psb_dpk_),zz,desc_a,info) + hasitnotfailed = psb_check_ans(zz,(-1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> mlt Z = X*Y (Complex Double Precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- mlt Z = X*Y (Complex Double Precision)")') end if ! @@ -476,28 +1392,74 @@ program vecoperation if (iam == psb_root_) write(psb_out_unit,'("Vector to Field Operations")') if (iam == psb_root_) write(psb_out_unit,'(" ")') - ! Dot product + ! Dot product (double real) call psb_d_gen_const(x,two,idim,ctxt,desc_a,info) call psb_d_gen_const(y,onehalf,idim,ctxt,desc_a,info) ans = psb_gedot(x,y,desc_a,info) if (iam == psb_root_) then - if(ans == np*idim) write(psb_out_unit,'("TEST PASSED >>> Dot product")') - if(ans /= np*idim) write(psb_out_unit,'("TEST FAILED --- Dot product")') + if(ans == np*idim) write(psb_out_unit,'("TEST PASSED >>> Dot product (double real)")') + if(ans /= np*idim) write(psb_out_unit,'("TEST FAILED --- Dot product (double real)")') + end if + ! Dot product (single precision real) + call psb_s_gen_const(sx,stwo,idim,ctxt,desc_a,info) + call psb_s_gen_const(sy,sonehalf,idim,ctxt,desc_a,info) + ans = psb_gedot(sx,sy,desc_a,info) + if (iam == psb_root_) then + if(ans == np*idim) write(psb_out_unit,'("TEST PASSED >>> Dot product (single precision real)")') + if(ans /= np*idim) write(psb_out_unit,'("TEST FAILED --- Dot product (single precision real)")') + end if + ! Dot product (single precision complex) + call psb_c_gen_const(cx,(2.0_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + call psb_c_gen_const(cy,(0.5_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + ans = psb_gedot(cx,cy,desc_a,info) + if (iam == psb_root_) then + if(ans == np*idim) write(psb_out_unit,'("TEST PASSED >>> Dot product (single precision complex)")') + if(ans /= np*idim) write(psb_out_unit,'("TEST FAILED --- Dot product (single precision complex)")') + end if + ! Dot product (double precision complex) + call psb_z_gen_const(zx,(2.0_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + call psb_z_gen_const(zy,(0.5_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + ans = psb_gedot(zx,zy,desc_a,info) + if (iam == psb_root_) then + if(ans == np*idim) write(psb_out_unit,'("TEST PASSED >>> Dot product (double precision complex)")') + if(ans /= np*idim) write(psb_out_unit,'("TEST FAILED --- Dot product (double precision complex)")') end if - ! MaxNorm + ! MaxNorm (double real) call psb_d_gen_const(x,negativeonehalf,idim,ctxt,desc_a,info) ans = psb_geamax(x,desc_a,info) if (iam == psb_root_) then - if(ans == onehalf) write(psb_out_unit,'("TEST PASSED >>> MaxNorm")') - if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm")') + if(ans == onehalf) write(psb_out_unit,'("TEST PASSED >>> MaxNorm (double real)")') + if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm (double real)")') + end if + ! MaxNorm (single real) + call psb_s_gen_const(sx,snegativeonehalf,idim,ctxt,desc_a,info) + ans = psb_geamax(sx,desc_a,info) + if (iam == psb_root_) then + if(ans == onehalf) write(psb_out_unit,'("TEST PASSED >>> MaxNorm (single precision real)")') + if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm (single precision real)")') + end if + ! MaxNorm (single complex) + call psb_c_gen_const(cx,(snegativeonehalf,0.0_psb_spk_),idim,ctxt,desc_a,info) + ans = psb_geamax(cx,desc_a,info) + if (iam == psb_root_) then + if(ans == onehalf) write(psb_out_unit,'("TEST PASSED >>> MaxNorm (single precision complex)")') + if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm (single precision complex)")') + end if + ! MaxNorm (double complex) + call psb_z_gen_const(zx,(snegativeonehalf,0.0_psb_dpk_),idim,ctxt,desc_a,info) + ans = psb_geamax(zx,desc_a,info) + if (iam == psb_root_) then + if(ans == onehalf) write(psb_out_unit,'("TEST PASSED >>> MaxNorm (double precision complex)")') + if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm (double precision complex)")') end if + ! ! Test of multivector operation ! if (iam == psb_root_) write(psb_out_unit,'(" ")') - if (iam == psb_root_) write(psb_out_unit,'("Multivector Operations")') + if (iam == psb_root_) write(psb_out_unit,'("---- Multivector Operations ----")') if (iam == psb_root_) write(psb_out_unit,'(" ")') ! X = 1 call psb_d_gen_const_multi(mv1,one,idim,nmv,ctxt,desc_a,info) @@ -506,6 +1468,27 @@ program vecoperation if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant multivector ")') if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant multivector ")') end if + ! X = 1 (single precision) + call psb_s_gen_const_multi(smv1,sone,idim,nmv,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(smv1,sone,ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant multivector (single precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant multivector (single precision)")') + end if + ! X = 1 (complex single precision) + call psb_c_gen_const_multi(cmv1,(1.0_psb_spk_, 0.0_psb_spk_),idim,nmv,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(cmv1,(1.0_psb_spk_, 0.0_psb_spk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant multivector (complex single precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant multivector (complex single precision)")') + end if + ! X = 1 (complex double precision) + call psb_z_gen_const_multi(zmv1,(1.0_psb_dpk_, 0.0_psb_dpk_),idim,nmv,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(zmv1,(1.0_psb_dpk_, 0.0_psb_dpk_),ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant multivector (complex double precision)")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant multivector (complex double precision)")') + end if ! ! Multivector to field operation @@ -514,23 +1497,77 @@ program vecoperation if (iam == psb_root_) write(psb_out_unit,'("Multivector to Field Operations")') if (iam == psb_root_) write(psb_out_unit,'(" ")') - ! Dot product: multivector vs multivector + ! Dot product: multivector vs multivector (double real) call psb_d_gen_const_multi(mv1,two,idim,nmv,ctxt,desc_a,info) call psb_d_gen_const_multi(mv2,onehalf,idim,nmv,ctxt,desc_a,info) ansmv = psb_gedot(mv1,mv2,desc_a,info) if (iam == psb_root_) then ! write ansmv to check - if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product")') - if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product")') + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs mv) (double real)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs mv) (double real)")') + end if + ! Dot product: multivector vs multivector (single real) + call psb_s_gen_const_multi(smv1,stwo,idim,nmv,ctxt,desc_a,info) + call psb_s_gen_const_multi(smv2,sonehalf,idim,nmv,ctxt,desc_a,info) + ansmv = psb_gedot(smv1,smv2,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs mv) (single real)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs mv) (single real)")') + end if + ! Dot product: multivector vs multivector (single complex) + call psb_c_gen_const_multi(cmv1,(2.0_psb_spk_, 0.0_psb_spk_),idim,nmv,ctxt,desc_a,info) + call psb_c_gen_const_multi(cmv2,(0.5_psb_spk_, 0.0_psb_spk_),idim,nmv,ctxt,desc_a,info) + ansmv = psb_gedot(cmv1,cmv2,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs mv) (single complex)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs mv) (single complex)")') end if - ! Dot product: multivector vs vector + ! Dot product: multivector vs multivector (double complex) + call psb_z_gen_const_multi(zmv1,(2.0_psb_dpk_, 0.0_psb_dpk_),idim,nmv,ctxt,desc_a,info) + call psb_z_gen_const_multi(zmv2,(0.5_psb_dpk_, 0.0_psb_dpk_),idim,nmv,ctxt,desc_a,info) + ansmv = psb_gedot(zmv1,zmv2,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs mv) (double complex)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs mv) (double complex)")') + end if + ! Dot product: multivector vs vector (double real) call psb_d_gen_const_multi(mv1,two,idim,nmv,ctxt,desc_a,info) call psb_d_gen_const(x,onehalf,idim,ctxt,desc_a,info) ansmv = psb_gedot(mv1,x,desc_a,info) if (iam == psb_root_) then ! write ansmv to check - if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product")') - if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product")') + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs vector) (double real)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs vector) (double real)")') + end if + ! Dot product: multivector vs vector (single real) + call psb_s_gen_const_multi(smv1,stwo,idim,nmv,ctxt,desc_a,info) + call psb_s_gen_const(sx,sonehalf,idim,ctxt,desc_a,info) + ansmv = psb_gedot(smv1,sx,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs vector) (single real)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs vector) (single real)")') + end if + ! Dot product: multivector vs vector (single complex) + call psb_c_gen_const_multi(cmv1,(2.0_psb_spk_, 0.0_psb_spk_),idim,nmv,ctxt,desc_a,info) + call psb_c_gen_const(cx,(0.5_psb_spk_, 0.0_psb_spk_),idim,ctxt,desc_a,info) + ansmv = psb_gedot(cmv1,cx,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs vector) (single complex)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs vector) (single complex)")') + end if + ! Dot product: multivector vs vector (double complex) + call psb_z_gen_const_multi(zmv1,(2.0_psb_dpk_, 0.0_psb_dpk_),idim,nmv,ctxt,desc_a,info) + call psb_z_gen_const(zx,(0.5_psb_dpk_, 0.0_psb_dpk_),idim,ctxt,desc_a,info) + ansmv = psb_gedot(zmv1,zx,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product (mv vs vector) (double complex)")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product (mv vs vector) (double complex)")') end if @@ -538,8 +1575,23 @@ program vecoperation call psb_gefree(x,desc_a,info) call psb_gefree(y,desc_a,info) call psb_gefree(z,desc_a,info) + call psb_gefree(sx,desc_a,info) + call psb_gefree(sy,desc_a,info) + call psb_gefree(sz,desc_a,info) + call psb_gefree(cx,desc_a,info) + call psb_gefree(cy,desc_a,info) + call psb_gefree(cz,desc_a,info) + call psb_gefree(zx,desc_a,info) + call psb_gefree(zy,desc_a,info) + call psb_gefree(zz,desc_a,info) call psb_gefree(mv1,desc_a,info) call psb_gefree(mv2,desc_a,info) + call psb_gefree(smv1,desc_a,info) + call psb_gefree(smv2,desc_a,info) + call psb_gefree(cmv1,desc_a,info) + call psb_gefree(cmv2,desc_a,info) + call psb_gefree(zmv1,desc_a,info) + call psb_gefree(zmv2,desc_a,info) call psb_cdfree(desc_a,info) if(info /= psb_success_) then info=psb_err_from_subroutine_