From 9323d3a7f42d5767f2609349bd39f2877e1075d8 Mon Sep 17 00:00:00 2001 From: Alfredo Buttari Date: Tue, 13 Sep 2005 16:13:33 +0000 Subject: [PATCH] *** empty log message *** --- src/comm/psb_dgather.f90 | 2 + src/comm/psb_dhalo.f90 | 42 +-- src/comm/psb_dovrl.f90 | 59 ++-- src/comm/psb_dscatter.f90 | 2 + src/comm/psb_ihalo.f90 | 34 +-- src/internals/Makefile | 2 +- src/internals/psi_compute_size.f90 | 2 +- src/internals/psi_crea_bnd_elem.f90 | 2 +- src/internals/psi_crea_index.f90 | 2 +- src/internals/psi_crea_ovr_elem.f90 | 2 +- src/internals/psi_desc_index.f90 | 6 +- src/internals/psi_dswapdata.f90 | 76 ++++++ src/internals/psi_dswaptran.f90 | 76 ++++++ src/internals/psi_iswapdata.f90 | 76 ++++++ src/internals/psi_iswaptran.f90 | 76 ++++++ src/internals/srcht.c | 218 +++++++++++++++ src/methd/psb_dcgs.f90 | 2 +- src/modules/Makefile | 2 +- src/modules/parts.fh | 7 + src/modules/psb_check_mod.f90 | 400 ++++++++++++++++++++++++++++ src/modules/psb_comm_mod.f90 | 4 +- src/modules/psb_spmat_type.f90 | 22 +- src/prec/psb_dbldaggrmat.f90 | 14 +- src/prec/psb_dcslu.f90 | 16 ++ src/prec/psb_dprec.f90 | 40 +-- src/prec/psb_dprecbld.f90 | 28 +- src/psblas/Makefile | 1 + src/psblas/psb_damax.f90 | 4 + src/psblas/psb_dasum.f90 | 3 + src/psblas/psb_daxpby.f90 | 6 +- src/psblas/psb_ddot.f90 | 5 +- src/psblas/psb_dnrm2.f90 | 3 + src/psblas/psb_dnrmi.f90 | 1 + src/psblas/psb_dspmm.f90 | 2 + src/psblas/psb_dspsm.f90 | 2 + src/serial/coo/Makefile | 2 +- src/serial/coo/dcoorws.f | 53 ++++ src/serial/coo/dcoosm.f | 6 +- src/serial/csr/Makefile | 3 +- src/serial/csr/dcsrrws.f | 34 +++ src/serial/csr/dcsrsm.f | 6 +- src/serial/dp/dcrjd.f | 4 +- src/serial/f77/Makefile | 3 +- src/serial/f77/dcopy.f | 50 ++++ src/serial/f77/dcsrws.f | 153 +++++++++++ src/serial/f77/dnrm2.f | 60 +++++ src/serial/f77/dscal.f | 43 +++ src/serial/f77/idamax.f | 39 +++ src/serial/f77/smmp.f | 2 +- src/serial/jad/Makefile | 2 +- src/serial/jad/djadrws.f | 48 ++++ src/serial/jad/djadsm.f | 6 +- src/serial/jad/djdrws.f | 30 +++ src/tools/psb_descasb.f90 | 2 +- 54 files changed, 1626 insertions(+), 159 deletions(-) create mode 100644 src/internals/srcht.c create mode 100644 src/modules/parts.fh create mode 100644 src/modules/psb_check_mod.f90 create mode 100644 src/serial/coo/dcoorws.f create mode 100644 src/serial/csr/dcsrrws.f create mode 100644 src/serial/f77/dcopy.f create mode 100644 src/serial/f77/dcsrws.f create mode 100644 src/serial/f77/dnrm2.f create mode 100644 src/serial/f77/dscal.f create mode 100644 src/serial/f77/idamax.f create mode 100644 src/serial/jad/djadrws.f create mode 100644 src/serial/jad/djdrws.f diff --git a/src/comm/psb_dgather.f90 b/src/comm/psb_dgather.f90 index 52000b52..9e18c6fd 100644 --- a/src/comm/psb_dgather.f90 +++ b/src/comm/psb_dgather.f90 @@ -19,6 +19,7 @@ subroutine psb_dgatherm(globx, locx, desc_a, info, iroot,& & iiglobx, ijglobx, iilocx,ijlocx,ik) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -192,6 +193,7 @@ end subroutine psb_dgatherm subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,& & iiglobx, iilocx) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none diff --git a/src/comm/psb_dhalo.f90 b/src/comm/psb_dhalo.f90 index e42ef05c..2f1da7e3 100644 --- a/src/comm/psb_dhalo.f90 +++ b/src/comm/psb_dhalo.f90 @@ -18,6 +18,8 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode) use psb_descriptor_type use psb_const_mod use psi_mod + use psb_check_mod + use psb_realloc_mod use psb_error_mod implicit none @@ -25,7 +27,7 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(kind(1.d0)), intent(in), optional :: alpha - real(kind(1.d0)), intent(inout), target, optional :: work(:) + real(kind(1.d0)), optional, target :: work(:) integer, intent(in), optional :: mode,jx,ik character, intent(in), optional :: tran @@ -116,22 +118,13 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode) end if liwork=ncol - if (present(work)) then - if(size(work).lt.liwork) then - call psrealloc(liwork,work,info) - if(info.ne.0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end if + if (present(work).and.(size(work).ge.liwork)) then iwork => work else - call psrealloc(liwork,iwork,info) + call psb_realloc(liwork,iwork,info) if(info.ne.0) then info=4010 - ch_err='psrealloc' + ch_err='psb_realloc' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -146,7 +139,7 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode) !!$ & size(x,1),desc_a%matrix_data,& !!$ & desc_a%halo_index,iwork,liwork,info) else if((ltran.eq.'T').or.(ltran.eq.'H')) then - call spi_swaptran(imode,k,1.d0,xp,& + call psi_swaptran(imode,k,1.d0,xp,& &desc_a,iwork,info) !!$ call PSI_dSwapTran(imode,k,1.d0,x(1,jjx),& !!$ & size(x,1),desc_a%matrix_data,& @@ -191,6 +184,8 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode) use psb_descriptor_type use psb_const_mod use psi_mod + use psb_check_mod + use psb_realloc_mod use psb_error_mod implicit none @@ -198,7 +193,7 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(kind(1.d0)), intent(in), optional :: alpha - real(kind(1.d0)), intent(inout), target, optional :: work(:) + real(kind(1.d0)), target, optional :: work(:) integer, intent(in), optional :: mode character, intent(in), optional :: tran @@ -271,27 +266,18 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode) end if liwork=ncol - if (present(work)) then - if(size(work).lt.liwork) then - call psrealloc(liwork,work,info) - if(info.ne.0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end if + if (present(work).and.(size(work).ge.liwork)) then iwork => work else - call psrealloc(liwork,iwork,info) + call psb_realloc(liwork,iwork,info) if(info.ne.0) then info=4010 - ch_err='psrealloc' + ch_err='psb_realloc' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if end if - + ! exchange halo elements if(ltran.eq.'N') then call psi_swapdata(imode,0.d0,x(iix:size(x)),& diff --git a/src/comm/psb_dovrl.f90 b/src/comm/psb_dovrl.f90 index 8068cf04..f62bf10b 100644 --- a/src/comm/psb_dovrl.f90 +++ b/src/comm/psb_dovrl.f90 @@ -16,10 +16,13 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,choice,update_type) use psb_descriptor_type use psb_const_mod + use psi_mod + use psb_realloc_mod + use psb_check_mod use psb_error_mod implicit none - real(kind(1.d0)), intent(inout) :: x(:,:) + real(kind(1.d0)), intent(inout), target :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(kind(1.d0)), intent(inout), optional, target :: work(:) @@ -30,7 +33,7 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,choice,update_type) integer :: int_err(5), icontxt, nprow, npcol, myrow, mycol,& & err_act, m, n, iix, jjx, temp(2), ix, ijx, nrow, ncol, k, maxk, iupdate,& & imode, err, liwork, i - real(kind(1.d0)),pointer :: iwork(:) + real(kind(1.d0)),pointer :: iwork(:), xp(:,:) logical :: ichoice character(len=20) :: name, ch_err @@ -109,22 +112,13 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,choice,update_type) ! check for presence/size of a work area liwork=ncol - if (present(work)) then - if(size(work).lt.liwork) then - call psrealloc(liwork,work,info) - if(info.ne.0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end if + if (present(work).and.(size(work).ge.liwork)) then iwork => work else - call psrealloc(liwork,iwork,info) + call psb_realloc(liwork,iwork,info) if(info.ne.0) then info=4010 - ch_err='psrealloc' + ch_err='psb_realloc' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -132,13 +126,13 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,choice,update_type) ! exchange overlap elements if(ichoice) then - call PSI_dSwapData(imode,k,1.d0,x(1,jjx),& - & size(x,1),desc_a%matrix_data,& - & desc_a%halo_index,iwork,liwork,info) + xp => x(iix:size(x,1),jjx:jjx+k-1) + call psi_swapdata(imode,k,1.d0,xp,& + & desc_a,iwork,info) end if if(info.ne.0) then - call psb_errpush(4010,name,a_err='PSI_dSwapData') + call psb_errpush(4010,name,a_err='psi_swapdata') goto 9999 end if @@ -202,11 +196,14 @@ end subroutine psb_dovrlm ! subroutine psb_dovrlv(x,desc_a,info,work,choice,update_type) use psb_descriptor_type + use psi_mod use psb_const_mod + use psb_realloc_mod + use psb_check_mod use psb_error_mod implicit none - real(kind(1.d0)), intent(inout) :: x(:) + real(kind(1.d0)), intent(inout), target :: x(:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(kind(1.d0)), intent(inout), optional, target :: work(:) @@ -264,7 +261,7 @@ subroutine psb_dovrlv(x,desc_a,info,work,choice,update_type) imode = IOR(psb_swap_send_,psb_swap_recv_) ! check vector correctness - call psb_chkvect(m,1,x,1,ix,ijx,desc_a%matrix_data,info,iix,jjx) + call psb_chkvect(m,1,size(x),ix,ijx,desc_a%matrix_data,info,iix,jjx) if(info.ne.0) then info=4010 ch_err='psb_chkvect' @@ -282,22 +279,13 @@ subroutine psb_dovrlv(x,desc_a,info,work,choice,update_type) ! check for presence/size of a work area liwork=ncol - if (present(work)) then - if(size(work).lt.liwork) then - call psrealloc(liwork,work,info) - if(info.ne.0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end if + if (present(work).and.(size(work).ge.liwork)) then iwork => work else - call psrealloc(liwork,iwork,info) + call psb_realloc(liwork,iwork,info) if(info.ne.0) then info=4010 - ch_err='psrealloc' + ch_err='psb_realloc' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -305,13 +293,12 @@ subroutine psb_dovrlv(x,desc_a,info,work,choice,update_type) ! exchange overlap elements if(ichoice) then - call PSI_dSwapData(imode,k,1.d0,x,& - & x,desc_a%matrix_data,& - & desc_a%halo_index,iwork,liwork,info) + call psi_swapdata(imode,1.d0,x(iix:size(x)),& + & desc_a,iwork,info) end if if(info.ne.0) then - call psb_errpush(4010,name,a_err='PSI_dSwapData') + call psb_errpush(4010,name,a_err='PSI_SwapData') goto 9999 end if diff --git a/src/comm/psb_dscatter.f90 b/src/comm/psb_dscatter.f90 index 045b2c2d..f696bae5 100644 --- a/src/comm/psb_dscatter.f90 +++ b/src/comm/psb_dscatter.f90 @@ -21,6 +21,7 @@ subroutine psb_dscatterm(globx, locx, desc_a, info, iroot,& & iiglobx, ijglobx, iilocx,ijlocx,ik) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none include 'mpif.h' @@ -239,6 +240,7 @@ end subroutine psb_dscatterm ! subroutine psb_dscatterv(globx, locx, desc_a, info, iroot) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none include 'mpif.h' diff --git a/src/comm/psb_ihalo.f90 b/src/comm/psb_ihalo.f90 index ff662a5c..483b700a 100644 --- a/src/comm/psb_ihalo.f90 +++ b/src/comm/psb_ihalo.f90 @@ -18,6 +18,8 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode) use psb_descriptor_type use psb_const_mod use psi_mod + use psb_realloc_mod + use psb_check_mod use psb_error_mod implicit none @@ -118,22 +120,13 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode) !!$ end if liwork=ncol - if (present(work)) then - if(size(work).lt.liwork) then - call psrealloc(liwork,work,info) - if(info.ne.0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end if + if (present(work).and.(size(work).ge.liwork)) then iwork => work else - call psrealloc(liwork,iwork,info) + call psb_realloc(liwork,iwork,info) if(info.ne.0) then info=4010 - ch_err='psrealloc' + ch_err='psb_realloc' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -187,6 +180,8 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode) use psb_descriptor_type use psb_const_mod use psi_mod + use psb_realloc_mod + use psb_check_mod use psb_error_mod implicit none @@ -267,22 +262,13 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode) !!$ end if liwork=ncol - if (present(work)) then - if(size(work).lt.liwork) then - call psb_realloc(liwork,work,info) - if(info.ne.0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end if + if (present(work).and.(size(work).ge.liwork)) then iwork => work else - call psrealloc(liwork,iwork,info) + call psb_realloc(liwork,iwork,info) if(info.ne.0) then info=4010 - ch_err='psrealloc' + ch_err='psb_realloc' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if diff --git a/src/internals/Makefile b/src/internals/Makefile index d1bb5e62..441106c8 100644 --- a/src/internals/Makefile +++ b/src/internals/Makefile @@ -4,7 +4,7 @@ FOBJS = psi_compute_size.o psi_crea_bnd_elem.o psi_crea_index.o \ psi_crea_ovr_elem.o psi_dl_check.o \ psi_exist_ovr_elem.o psi_gthsct.o \ psi_list_search.o psi_sort_dl.o srtlist.o -COBJS = avltree.o +COBJS = avltree.o srcht.o MPFOBJS = psi_dswapdata.o psi_dswaptran.o psi_iswapdata.o \ psi_iswaptran.o psi_extrct_dl.o psi_desc_index.o diff --git a/src/internals/psi_compute_size.f90 b/src/internals/psi_compute_size.f90 index 1fe26fb1..eb5f23b5 100644 --- a/src/internals/psi_compute_size.f90 +++ b/src/internals/psi_compute_size.f90 @@ -91,7 +91,7 @@ subroutine psi_compute_size(desc_data,& 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.act_abort) then - call psb_serror(icontxt) + call psb_error(icontxt) return end if return diff --git a/src/internals/psi_crea_bnd_elem.f90 b/src/internals/psi_crea_bnd_elem.f90 index 8b8fffc7..9fab9ed3 100644 --- a/src/internals/psi_crea_bnd_elem.f90 +++ b/src/internals/psi_crea_bnd_elem.f90 @@ -66,7 +66,7 @@ subroutine psi_crea_bnd_elem(desc_a,info) 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.act_abort) then - call psb_serror() + call psb_error() return end if return diff --git a/src/internals/psi_crea_index.f90 b/src/internals/psi_crea_index.f90 index d7daddf4..2fe9a417 100644 --- a/src/internals/psi_crea_index.f90 +++ b/src/internals/psi_crea_index.f90 @@ -81,7 +81,7 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info) 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.act_abort) then - call psb_serror(icontxt) + call psb_error(icontxt) return end if return diff --git a/src/internals/psi_crea_ovr_elem.f90 b/src/internals/psi_crea_ovr_elem.f90 index 395cdfe8..8c1162d6 100644 --- a/src/internals/psi_crea_ovr_elem.f90 +++ b/src/internals/psi_crea_ovr_elem.f90 @@ -45,7 +45,7 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem) if (pnt_new_elem.gt.dim_ovr_elem) then dim=(3*size(ovr_elem))/2+2 write(0,*) 'calling realloc crea_ovr_elem',dim - call psrealloc(dim,ovr_elem,info) + call psb_realloc(dim,ovr_elem,info) endif else ! ....this point already exist in ovr_elem list diff --git a/src/internals/psi_desc_index.f90 b/src/internals/psi_desc_index.f90 index e9b8724b..69dbf3bd 100644 --- a/src/internals/psi_desc_index.f90 +++ b/src/internals/psi_desc_index.f90 @@ -115,10 +115,10 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,& !c$$$ write(0,*) 'potential error on desc_index :', !c$$$ + length_dh, size(desc_index),ntot write(0,*) 'calling irealloc psi_desc_index ',ntot - call psrealloc(ntot,desc_index,info) + call psb_realloc(ntot,desc_index,info) endif if (info /= 0) then - call psb_errpush(4010,name,a_err='psrealloc') + call psb_errpush(4010,name,a_err='psb_realloc') goto 9999 end if @@ -225,7 +225,7 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,& 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.act_abort) then - call psb_serror(icontxt) + call psb_error(icontxt) return end if return diff --git a/src/internals/psi_dswapdata.f90 b/src/internals/psi_dswapdata.f90 index 17b771ae..50c125a8 100644 --- a/src/internals/psi_dswapdata.f90 +++ b/src/internals/psi_dswapdata.f90 @@ -26,6 +26,44 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info) real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_dswap_data' call psb_erractionsave(err_act) @@ -399,6 +437,44 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info) real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_dswap_datav' call psb_erractionsave(err_act) diff --git a/src/internals/psi_dswaptran.f90 b/src/internals/psi_dswaptran.f90 index 8f3d9ec1..b90757ea 100644 --- a/src/internals/psi_dswaptran.f90 +++ b/src/internals/psi_dswaptran.f90 @@ -27,6 +27,44 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info) real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_dswaptranm' call psb_erractionsave(err_act) @@ -397,6 +435,44 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info) real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_dswaptranv' call psb_erractionsave(err_act) diff --git a/src/internals/psi_iswapdata.f90 b/src/internals/psi_iswapdata.f90 index d17efecb..e164a932 100644 --- a/src/internals/psi_iswapdata.f90 +++ b/src/internals/psi_iswapdata.f90 @@ -26,6 +26,44 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info) integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_iswapdata' call psb_erractionsave(err_act) @@ -399,6 +437,44 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info) integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_iswapdatav' call psb_erractionsave(err_act) diff --git a/src/internals/psi_iswaptran.f90 b/src/internals/psi_iswaptran.f90 index 6d733769..0d064bdc 100644 --- a/src/internals/psi_iswaptran.f90 +++ b/src/internals/psi_iswaptran.f90 @@ -27,6 +27,44 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info) integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_dswaptranm' call psb_erractionsave(err_act) @@ -397,6 +435,44 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info) integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name, ch_err + interface psi_gth + subroutine psi_dgthm(n,k,idx,x,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + end subroutine psi_dgthm + subroutine psi_dgthv(n,idx,x,y) + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + end subroutine psi_dgthv + subroutine psi_igthm(n,k,idx,x,y) + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + end subroutine psi_igthm + subroutine psi_igthv(n,idx,x,y) + integer :: n, idx(:) + integer :: x(:), y(:) + end subroutine psi_igthv + end interface + + interface psi_sct + subroutine psi_dsctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + end subroutine psi_dsctm + subroutine psi_dsctv(n,idx,x,beta,y) + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + end subroutine psi_dsctv + subroutine psi_isctm(n,k,idx,x,beta,y) + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + end subroutine psi_isctm + subroutine psi_isctv(n,idx,x,beta,y) + integer :: n, idx(:) + integer :: beta, x(:), y(:) + end subroutine psi_isctv + end interface + info = 0 name='psi_dswaptranv' call psb_erractionsave(err_act) diff --git a/src/internals/srcht.c b/src/internals/srcht.c new file mode 100644 index 00000000..d3ac0e7c --- /dev/null +++ b/src/internals/srcht.c @@ -0,0 +1,218 @@ +/*****************************************************************/ +/* */ +/* srcht.c : specialized insert/search for (key,val) integer */ +/* pairs. Written by: Salvatore Filippone */ +/* */ +/* Last updated: Mar 09 2004 */ +/* */ +/* Uses: avltree */ +/* */ +/* Data types: */ +/* */ +/* KeyType: struct with two integer fields, key and val. */ +/* */ +/* */ +/* User callable functions: */ +/* */ +/* void InitPairSearchTree(int *iret) */ +/* Purpose: initialize a search structure; */ +/* Function value: 0: OK */ +/* -1: failure */ +/* */ +/* */ +/* void SearchInsKeyVal(int *key, int *val, int *res, */ +/* int *iret) */ +/* Purpose: Search for a key, insert it if not present. */ +/* */ +/* Input: 1. key */ +/* Key to be searched for. */ +/* 2. val */ +/* Value to form a (key,val) pair to be */ +/* inserted if key not already present. */ +/* Output: 3. res */ +/* The val part of the pair with key; if the */ +/* key was freshly inserted then res=val */ +/* Function value: 0 Normal termination */ +/* 1 Key was already present. */ +/* -1 Invalid input pointer */ +/* -3 Memory allocation failure */ +/* */ +/* */ +/* void FreePairSearchTree() */ +/* Purpose: free up tree data storage */ +/* */ +/* */ +/*****************************************************************/ + + + +#include +#include +#include +#include "avltree.h" + +#define POOLSIZE 4096 + +#ifdef Add_ +#define InitPairSearchTree initpairsearchtree_ +#define FreePairSearchTree freepairsearchtree_ +#define SearchInsKeyVal searchinskeyval_ +#endif +#ifdef AddDouble_ +#define InitPairSearchTree initpairsearchtree_ +#define FreePairSearchTree freepairsearchtree_ +#define SearchInsKeyVal searchinskeyval_ +#endif +#ifdef NoChange +#define InitPairSearchTree initpairsearchtree +#define FreePairSearchTree freepairsearchtree +#define SearchInsKeyVal searchinskeyval +#endif + + + +typedef struct keypair *KeyPairPtr; +typedef struct keypair { + int key,val; +} KeyPair; + + +typedef struct pairvect *PairVectPtr; +typedef struct pairvect { + KeyPair pool[POOLSIZE]; + int avail; + PairVectPtr previous, next; +} PairVect; + + +static int retval; +static PairVectPtr PairPoolRoot=NULL,PairPoolCrt=NULL; +static AVLTreePtr tree=NULL; + +int CompareKeys(void *key1, void *key2) +{ + if (((KeyPairPtr) key1)->key < ((KeyPairPtr) key2)->key){ + return(-1); + } else if (((KeyPairPtr)key1)->key == ((KeyPairPtr)key2)->key){ + return(0); + } else { + return(1); + } +} + +void InitPairSearchTree(int *iret) +{ + *iret = 0; + + if ((tree = GetAVLTree())==NULL) { + *iret=-1; return; + } + if ((PairPoolRoot=(PairVectPtr)malloc(sizeof(PairVect)))==NULL) { + *iret=-3; + } else { + PairPoolRoot->avail=0; + PairPoolRoot->previous=PairPoolRoot->next=NULL; + PairPoolCrt=PairPoolRoot; + } + return; +} + + +void KeyUpdate( void *key1, void *key2) +{ + retval=((KeyPairPtr) key2)->val; +} + + +void FreePairSearchTree() +{ + PairVectPtr current,next; + + AVLTreeFree(tree,NULL); + + current=PairPoolRoot; + + while (current != NULL) { + next=current->next; + free(current); + current=next; + } + free(tree); + tree = NULL; + return; +} + +int AdvanceKeyPair(PairVectPtr current) +{ + if (current!=NULL) { + current->avail +=1; + return(current->avail); + } + return(-1); +} + + +KeyPairPtr GetKeyPair(PairVectPtr *current) +{ + PairVectPtr new, crt; + KeyPairPtr newnode; + + crt=*current; + if (crt==NULL) { + return(NULL); + } + + if (crt->availpool[crt->avail]); + } else { + if ((new=(PairVectPtr)malloc(sizeof(PairVect)))==NULL) { + fprintf(stderr,"Memory allocation failure\n"); + return(NULL); + } + memset(new,'\0',sizeof(PairVect)); + newnode=&(new->pool[0]); + crt->next=new; + new->previous=crt; + new->next=NULL; + *current=new; + } + return(newnode); +} + + +/* */ +/* void SearchInsKeyVal(int *key, int *val, int *res, */ +/* int *iret) */ +/* Purpose: Search for a key, insert it if not present. */ +/* */ +/* Input: 1. key */ +/* Key to be searched for. */ +/* 2. val */ +/* Value to form a (key,val) pair to be */ +/* inserted if key not already present. */ +/* Output: 3. res */ +/* The val part of the pair with key; if the */ +/* key was freshly inserted then res=val */ +/* Function value: 0 Normal termination */ +/* -1 Invalid input pointer */ +/* -3 Memory allocation failure */ +/* */ + +void SearchInsKeyVal(int *key, int *val, int *res, int *iret) +{ + KeyPairPtr node; int info; + + node = GetKeyPair(&PairPoolCrt); + node->key = *key; + node->val = *val; + + info = AVLTreeInsert(tree,node,CompareKeys,KeyUpdate); + *iret = info; + if (info==0) { + *res = node->val; + AdvanceKeyPair(PairPoolCrt); + } else if (info == 1) { + *res = retval; + } + return; +} diff --git a/src/methd/psb_dcgs.f90 b/src/methd/psb_dcgs.f90 index 6ed6b73f..97764790 100644 --- a/src/methd/psb_dcgs.f90 +++ b/src/methd/psb_dcgs.f90 @@ -307,7 +307,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& End If Deallocate(aux) - Call psb_dsfree(wwrk,desc_a,info) + Call psb_free(wwrk,desc_a,info) ! restore external global coherence behaviour Call blacs_set(icontxt,16,isvch) diff --git a/src/modules/Makefile b/src/modules/Makefile index 692db670..e05d4180 100644 --- a/src/modules/Makefile +++ b/src/modules/Makefile @@ -6,7 +6,7 @@ MODULES = psb_realloc_mod.o psb_string_mod.o psb_spmat_type.o \ psb_prec_type.o psb_error_mod.o psb_prec_mod.o \ psb_methd_mod.o psb_const_mod.o \ psb_comm_mod.o psb_psblas_mod.o psi_mod.o \ - psb_sparse_mod.o + psb_sparse_mod.o psb_check_mod.o OBJS = error.o diff --git a/src/modules/parts.fh b/src/modules/parts.fh new file mode 100644 index 00000000..10e323d0 --- /dev/null +++ b/src/modules/parts.fh @@ -0,0 +1,7 @@ +interface + !.....user passed subroutine..... + subroutine parts(glob_index,nrow,np,pv,nv) + integer, intent (in) :: glob_index,np,nrow + integer, intent (out) :: nv, pv(*) + end subroutine parts +end interface diff --git a/src/modules/psb_check_mod.f90 b/src/modules/psb_check_mod.f90 new file mode 100644 index 00000000..d044c113 --- /dev/null +++ b/src/modules/psb_check_mod.f90 @@ -0,0 +1,400 @@ +! File: psb_check_mod.f90 + +module psb_check_mod + +! interface +! module procedure psb_chkvect +! end interface + +! interface +! module procedure psb_chkglobvect +! end interface + +! interface +! module procedure psb_chkmat +! end interface + +contains + ! Subroutine: psb_chkvect + ! psb_chkvect checks the validity of a descriptor vector desc_dec, the + ! related global indexes ix, jx and the leading dimension lldx. It also + ! eventually computes the starting local indexes (iix,jjx) corresponding + ! to the submatrix starting globally at the entry pointed by (ix,jx). + ! Finally, if an inconsistency is found among its parameters ix, jx, + ! descdec and lldx, the routine returns an error code in info. + ! + ! Parameters: + ! m - integer. The number of rows of the dense matrix X being operated on. + ! n - integer. The number of columns of the dense matrix X being operated on. + ! lldx - integer. The leading dimension of the local dense matrix X. + ! ix - integer. X's global row index, which points to the beginning + ! of the dense submatrix which is to be operated on. + ! jx - integer. X's global column index, which points to the beginning + ! of the dense submatrix which is to be operated on. + ! desc_dec - integer,dimension(:). Is the matrix_data array. + ! info - integer. Eventually returns an error code. + ! iix - integer(optional). The local rows starting index of the submatrix. + ! jjx - integer(optional). The local columns starting index of the submatrix. + subroutine psb_chkvect( m, n, lldx, ix, jx, desc_dec, info, iix, jjx) + + use psb_const_mod + use psb_error_mod + implicit none + + integer, intent(in) :: m,n,ix,jx,lldx + integer, intent(in) :: desc_dec(:) + integer, intent(out) :: info + integer, optional :: iix, jjx + + ! locals + integer :: err_act, int_err(5) + character(len=20) :: name, ch_err + + info=0 + name='psb_chkvect' + call psb_erractionsave(err_act) + + + if (m.lt.0) then + info=10 + int_err(1) = 1 + int_err(2) = m + else if (n.lt.0) then + info=10 + int_err(1) = 3 + int_err(2) = n + else if ((ix.lt.1) .and. (m.ne.0)) then + info=20 + int_err(1) = 4 + int_err(2) = ix + else if ((jx.lt.1) .and. (n.ne.0)) then + info=20 + int_err(1) = 5 + int_err(2) = jx + else if (desc_dec(psb_n_col_).lt.0) then + info=40 + int_err(1) = 6 + int_err(2) = psb_n_col_ + int_err(3) = desc_dec(psb_n_col_) + else if (desc_dec(psb_n_row_).lt.0) then + info=40 + int_err(1) = 6 + int_err(2) = psb_n_row_ + int_err(3) = desc_dec(psb_n_row_) + else if (lldx.lt.desc_dec(psb_n_col_)) then + info=50 + int_err(1) = 3 + int_err(2) = lldx + int_err(3) = 6 + int_err(4) = psb_n_col_ + int_err(5) = desc_dec(psb_n_col_) + else if (desc_dec(psb_n_).lt.m) then + info=60 + int_err(1) = 1 + int_err(2) = m + int_err(3) = 6 + int_err(4) = psb_n_ + int_err(5) = desc_dec(psb_n_) + else if (desc_dec(psb_n_).lt.ix) then + info=60 + int_err(1) = 4 + int_err(2) = ix + int_err(3) = 6 + int_err(4) = psb_n_ + int_err(5) = desc_dec(psb_n_) + else if (desc_dec(psb_m_).lt.jx) then + info=60 + int_err(1) = 5 + int_err(2) = jx + int_err(3) = 6 + int_err(4) = psb_m_ + int_err(5) = desc_dec(psb_m_) + else if (desc_dec(psb_n_).lt.(ix+m-1)) then + info=80 + int_err(1) = 1 + int_err(2) = m + int_err(3) = 4 + int_err(4) = ix + end if + + if (info.ne.0) then + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + ! Compute local indices for submatrix starting + ! at global indices ix and jx + if(present(iix)) iix=ix ! (for our applications iix=ix)) + if(present(jjx)) iix=ix ! (for our applications jjx=jx)) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + + end subroutine psb_chkvect + + ! + ! Subroutine: psb_chkglobvect + ! psb_chkglobvect checks the validity of a descriptor vector desc_dec, the + ! related global indexes ix, jx and the leading dimension lldx. + ! If an inconsistency is found among its parameters ix, jx, + ! descdec and lldx, the routine returns an error code in info. + ! + ! Parameters: + ! m - integer. The number of rows of the dense matrix X being operated on. + ! n - integer. The number of columns of the dense matrix X being operated on. + ! lldx - integer. The leading dimension of the local dense matrix X. + ! ix - integer. X's global row index, which points to the beginning + ! of the dense submatrix which is to be operated on. + ! jx - integer. X's global column index, which points to the beginning + ! of the dense submatrix which is to be operated on. + ! desc_dec - integer,dimension(:). Is the matrix_data array. + ! info - integer. Eventually returns an error code. + ! + subroutine psb_chkglobvect( m, n, lldx, ix, jx, desc_dec, info) + + use psb_const_mod + use psb_error_mod + implicit none + + integer, intent(in) :: m,n,ix,jx,lldx + integer, intent(in) :: desc_dec(:) + integer, intent(out) :: info + + ! locals + integer :: err_act, int_err(5) + character(len=20) :: name, ch_err + + info=0 + name='psb_chkglobvect' + call psb_erractionsave(err_act) + + + if (m.lt.0) then + info=10 + int_err(1) = 1 + int_err(2) = m + else if (n.lt.0) then + info=10 + int_err(1) = 3 + int_err(2) = n + else if ((ix.lt.1) .and. (m.ne.0)) then + info=20 + int_err(1) = 4 + int_err(2) = ix + else if ((jx.lt.1) .and. (n.ne.0)) then + info=20 + int_err(1) = 5 + int_err(2) = jx + else if (desc_dec(psb_n_col_).lt.0) then + info=40 + int_err(1) = 6 + int_err(2) = psb_n_col_ + int_err(3) = desc_dec(psb_n_col_) + else if (desc_dec(psb_n_row_).lt.0) then + info=40 + int_err(1) = 6 + int_err(2) = psb_n_row_ + int_err(3) = desc_dec(psb_n_row_) + else if (lldx.lt.desc_dec(psb_m_)) then + info=50 + int_err(1) = 3 + int_err(2) = lldx + int_err(3) = 6 + int_err(4) = psb_n_col_ + int_err(5) = desc_dec(psb_n_col_) + else if (desc_dec(psb_n_).lt.m) then + info=60 + int_err(1) = 1 + int_err(2) = m + int_err(3) = 6 + int_err(4) = psb_n_ + int_err(5) = desc_dec(psb_n_) + else if (desc_dec(psb_n_).lt.ix) then + info=60 + int_err(1) = 4 + int_err(2) = ix + int_err(3) = 6 + int_err(4) = psb_n_ + int_err(5) = desc_dec(psb_n_) + else if (desc_dec(psb_m_).lt.jx) then + info=60 + int_err(1) = 5 + int_err(2) = jx + int_err(3) = 6 + int_err(4) = psb_m_ + int_err(5) = desc_dec(psb_m_) + else if (desc_dec(psb_n_).lt.(ix+m-1)) then + info=80 + int_err(1) = 1 + int_err(2) = m + int_err(3) = 4 + int_err(4) = ix + end if + + if (info.ne.0) then + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + + end subroutine psb_chkglobvect + + ! + ! Subroutine: psb_chkmat + ! pbmatvect checks the validity of a descriptor vector DESCDEC, the + ! related global indexes IA, JA. It also computes the starting local + ! indexes (IIA,JJA) corresponding to the submatrix starting globally at + ! the entry pointed by (IA,JA). Finally, if an inconsitency is found among + ! its parameters ia, ja and desc_A, the routine returns an error code in + ! info. + ! + ! Parameters: + ! m - integer. The number of rows of the matrix being operated on. + ! n - integer. The number of columns of the matrix being operated on. + ! ia - integer. a's global row index, which points to the beginning + ! of the submatrix which is to be operated on. + ! ja - integer. a's global column index, which points to the beginning + ! of the submatrix which is to be operated on. + ! desc_dec - integer,dimension(:). Is the matrix_data array. + ! info - integer. Eventually returns an error code. + ! iia - integer(optional). The local rows starting index of the submatrix. + ! jja - integer(optional). The local columns starting index of the submatrix. + ! + subroutine psb_chkmat( m, n, ia, ja, desc_dec, info, iia, jja) + + use psb_const_mod + use psb_error_mod + implicit none + + integer, intent(in) :: m,n,ia,ja + integer, intent(in) :: desc_dec(:) + integer, intent(out) :: info + integer, optional :: iia, jja + + ! locals + integer :: err_act, int_err(5) + character(len=20) :: name, ch_err + + info=0 + name='psb_chkmat' + call psb_erractionsave(err_act) + + if (m.lt.0) then + info=10 + int_err(1) = 1 + int_err(2) = m + else if (n.lt.0) then + info=10 + int_err(1) = 3 + int_err(2) = n + else if ((ia.lt.1) .and. (m.ne.0)) then + info=20 + int_err(1) = 4 + int_err(2) = ia + else if ((ja.lt.1) .and. (n.ne.0)) then + info=20 + int_err(1) = 5 + int_err(2) = ja + else if (desc_dec(psb_n_col_).lt.0) then + info=40 + int_err(1) = 6 + int_err(2) = psb_n_col_ + int_err(3) = desc_dec(psb_n_col_) + else if (desc_dec(psb_n_row_).lt.0) then + info=40 + int_err(1) = 6 + int_err(2) = psb_n_row_ + int_err(3) = desc_dec(psb_n_row_) + else if (desc_dec(psb_m_).lt.m) then + info=60 + int_err(1) = 1 + int_err(2) = m + int_err(3) = 5 + int_err(4) = psb_m_ + int_err(5) = desc_dec(psb_m_) + else if (desc_dec(psb_m_).lt.m) then + info=60 + int_err(1) = 2 + int_err(2) = n + int_err(3) = 5 + int_err(4) = psb_m_ + int_err(5) = desc_dec(psb_m_) + else if (desc_dec(psb_m_).lt.ia) then + info=60 + int_err(1) = 3 + int_err(2) = ia + int_err(3) = 5 + int_err(4) = psb_m_ + int_err(5) = desc_dec(psb_m_) + else if (desc_dec(psb_n_).lt.ja) then + info=60 + int_err(1) = 4 + int_err(2) = ja + int_err(3) = 5 + int_err(4) = psb_n_ + int_err(5) = desc_dec(psb_n_) + else if (desc_dec(psb_m_).lt.(ia+m-1)) then + info=80 + int_err(1) = 1 + int_err(2) = m + int_err(3) = 3 + int_err(4) = ia + else if (desc_dec(psb_n_).lt.(ja+n-1)) then + info=80 + int_err(1) = 2 + int_err(2) = n + int_err(3) = 4 + int_err(4) = ja + end if + + if (info.ne.0) then + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + ! Compute local indices for submatrix starting + ! at global indices ix and jx + if(present(iia).and.present(jja)) then + if (desc_dec(psb_n_row_).gt.0) then + iia=1 + jja=1 + else + iia=desc_dec(psb_n_row_)+1 + jja=desc_dec(psb_n_col_)+1 + end if + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + end subroutine psb_chkmat + +end module psb_check_mod diff --git a/src/modules/psb_comm_mod.f90 b/src/modules/psb_comm_mod.f90 index 963e4424..7fe666e8 100644 --- a/src/modules/psb_comm_mod.f90 +++ b/src/modules/psb_comm_mod.f90 @@ -28,7 +28,7 @@ module psb_comm_mod type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(kind(1.d0)), intent(in), optional :: alpha - real(kind(1.d0)), intent(inout), optional :: work(:) + real(kind(1.d0)), target, optional :: work(:) integer, intent(in), optional :: mode,jx,ik character, intent(in), optional :: tran end subroutine psb_dhalom @@ -38,7 +38,7 @@ module psb_comm_mod type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(kind(1.d0)), intent(in), optional :: alpha - real(kind(1.d0)), intent(inout), optional :: work(:) + real(kind(1.d0)), target, optional :: work(:) integer, intent(in), optional :: mode character, intent(in), optional :: tran end subroutine psb_dhalov diff --git a/src/modules/psb_spmat_type.f90 b/src/modules/psb_spmat_type.f90 index bb64931f..eeabde10 100644 --- a/src/modules/psb_spmat_type.f90 +++ b/src/modules/psb_spmat_type.f90 @@ -231,19 +231,19 @@ contains endif if (ifc_ == 1) then - call psrealloc(nnz,a%ia1,a%ia2,a%aspk,info) + call psb_realloc(nnz,a%ia1,a%ia2,a%aspk,info) else - call psrealloc(nnz,a%aspk,info) + call psb_realloc(nnz,a%aspk,info) if (info /= 0) return - call psrealloc(nnz,a%ia2,info) + call psb_realloc(nnz,a%ia2,info) if (info /= 0) return - call psrealloc(ifc*nnz+200,a%ia1,info) + call psb_realloc(ifc*nnz+200,a%ia1,info) if (info /= 0) return end if if (info /= 0) return - call psrealloc(max(1,a%m),a%pl,info) + call psb_realloc(max(1,a%m),a%pl,info) if (info /= 0) return - call psrealloc(max(1,a%k),a%pr,info) + call psb_realloc(max(1,a%k),a%pr,info) if (info /= 0) return Return @@ -261,15 +261,15 @@ contains logical, parameter :: debug=.false. info = 0 - call psrealloc(nd,a%aspk,info) + call psb_realloc(nd,a%aspk,info) if (info /= 0) return - call psrealloc(ni2,a%ia2,info) + call psb_realloc(ni2,a%ia2,info) if (info /= 0) return - call psrealloc(ni1,a%ia1,info) + call psb_realloc(ni1,a%ia1,info) if (info /= 0) return - call psrealloc(max(1,a%m),a%pl,info) + call psb_realloc(max(1,a%m),a%pl,info) if (info /= 0) return - call psrealloc(max(1,a%k),a%pr,info) + call psb_realloc(max(1,a%k),a%pr,info) if (info /= 0) return Return diff --git a/src/prec/psb_dbldaggrmat.f90 b/src/prec/psb_dbldaggrmat.f90 index c6fbc108..e1be9011 100644 --- a/src/prec/psb_dbldaggrmat.f90 +++ b/src/prec/psb_dbldaggrmat.f90 @@ -513,8 +513,8 @@ contains ! Doing it this way means to consider diag(Ai) ! ! - call symbmm90(am3,am4,am1) - call numbmm90(am3,am4,am1) + call psb_symbmm(am3,am4,am1) + call psb_numbmm(am3,am4,am1) call psb_spfree(am4,info) @@ -559,8 +559,8 @@ contains if (test_dump) & & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) - call symbmm90(a,am1,am3) - call numbmm90(a,am1,am3) + call psb_symbmm(a,am1,am3) + call psb_numbmm(a,am1,am3) if (p%iprcparm(smth_kind_) == smth_omg_) then call psb_transp(am1,am2,fmt='COO') @@ -613,8 +613,8 @@ contains end if endif - call symbmm90(am2,am3,b) - call numbmm90(am2,am3,b) + call psb_symbmm(am2,am3,b) + call psb_numbmm(am2,am3,b) !!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.') call psb_spfree(am3,info) @@ -834,7 +834,7 @@ contains end if call psb_dscdec(naggr,icontxt,p%desc_data,info) - call spfree(b,info) + call psb_spfree(b,info) if(info /= 0) then call psb_errpush(4010,name,a_err='spfree') goto 9999 diff --git a/src/prec/psb_dcslu.f90 b/src/prec/psb_dcslu.f90 index d05a62bf..43cfc03e 100644 --- a/src/prec/psb_dcslu.f90 +++ b/src/prec/psb_dcslu.f90 @@ -60,6 +60,22 @@ subroutine psb_dcslu(a,desc_a,p,upd,info) end subroutine psb_dsplu end interface + interface psb_csrsetup + Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) + use psb_serial_mod + Use psb_descriptor_type + Use psb_prec_type + integer, intent(in) :: ptype,novr + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_dspmat_type), Intent(inout) :: blk + Type(psb_desc_type), Intent(inout) :: desc_p + Type(psb_desc_type), Intent(in) :: desc_data + Character, Intent(in) :: upd + integer, intent(out) :: info + character(len=5), optional :: outfmt + end Subroutine psb_dcsrsetup + end interface + info=0 name='psb_dcslu' call psb_erractionsave(err_act) diff --git a/src/prec/psb_dprec.f90 b/src/prec/psb_dprec.f90 index 87430217..8ede592a 100644 --- a/src/prec/psb_dprec.f90 +++ b/src/prec/psb_dprec.f90 @@ -8,12 +8,12 @@ subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) use psb_error_mod implicit none - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(kind(0.d0)),intent(inout), optional, target :: work(:) + type(psb_desc_type),intent(in) :: desc_data + type(psb_dprec_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(kind(0.d0)), optional, target :: work(:) ! Local variables character ::trans_ @@ -97,7 +97,7 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) real(kind(0.d0)),intent(inout) :: x(:), y(:) real(kind(0.d0)),intent(in) :: beta character(len=1) :: trans - real(kind(0.d0)),intent(inout),target :: work(:) + real(kind(0.d0)),target :: work(:) integer, intent(out) :: info ! Local variables @@ -161,7 +161,7 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) case(bja_) - call psb_bjacaply(prec,x,beta,y,desc_data,trans,work,info) + call psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 ch_err='psb_bjacaply' @@ -220,7 +220,7 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) end if endif - call psb_bjacaply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) + call psb_dbjacaply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) if (prec%iprcparm(iren_)>0) then call psb_dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info) @@ -317,7 +317,7 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) real(kind(0.d0)),intent(inout) :: x(:), y(:) real(kind(0.d0)),intent(in) :: beta character(len=1) :: trans - real(kind(0.d0)),intent(inout),target :: work(:) + real(kind(0.d0)),target :: work(:) integer, intent(out) :: info ! Local variables @@ -526,7 +526,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) real(kind(0.d0)),intent(in) :: beta real(kind(0.d0)),intent(inout) :: x(:), y(:) character :: trans - real(kind(0.d0)),intent(inout),target :: work(:) + real(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -567,7 +567,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) t1 = mpi_wtime() n_row=desc_data%matrix_data(psb_n_row_) n_col=desc_data%matrix_data(psb_n_col_) - call psb_baseprcaply(baseprecv(1),x,beta,y,desc_data,trans,work,info) + call psb_dbaseprcaply(baseprecv(1),x,beta,y,desc_data,trans,work,info) if(info /=0) goto 9999 nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) @@ -616,7 +616,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) endif w2l=t2l - call psb_baseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) if (ismth /= no_smth_) then @@ -706,7 +706,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) t6 = mpi_wtime() w2l=t2l - call psb_baseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) if(info /=0) goto 9999 if (ismth /= no_smth_) then @@ -730,7 +730,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) call psb_spmm(-one,baseprecv(2)%aorig,ty,one,tx,desc_data,info,work=work) if(info /=0) goto 9999 - call psb_baseprcaply(baseprecv(1),tx,one,ty,desc_data,trans,work,info) + call psb_dbaseprcaply(baseprecv(1),tx,one,ty,desc_data,trans,work,info) if(info /=0) goto 9999 call psb_axpby(one,ty,beta,y,desc_data,info) @@ -759,7 +759,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) call psb_axpby(one,y,zero,ty,desc_data,info) if(info /=0) goto 9999 - call psb_baseprcaply(baseprecv(1),x,zero,tty,desc_data,trans,work,info) + call psb_dbaseprcaply(baseprecv(1),x,zero,tty,desc_data,trans,work,info) if(info /=0) goto 9999 call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) @@ -796,7 +796,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) t6 = mpi_wtime() w2l=t2l - call psb_baseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) if(info /=0) goto 9999 if (ismth /= no_smth_) then @@ -852,7 +852,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) end subroutine psb_dmlprcaply -subroutine psb_dprec1(prec,x,desc_data,info,trans) +subroutine psb_dprecaply1(prec,x,desc_data,info,trans) use psb_serial_mod use psb_descriptor_type @@ -890,7 +890,7 @@ subroutine psb_dprec1(prec,x,desc_data,info,trans) end if allocate(ww(size(x)),w1(size(x))) - call psb_dprec(prec,x,ww,desc_data,info,trans_,w1) + call psb_dprecaply(prec,x,ww,desc_data,info,trans_,w1) if(info /=0) goto 9999 x(:) = ww(:) deallocate(ww,W1) @@ -906,5 +906,5 @@ subroutine psb_dprec1(prec,x,desc_data,info,trans) return end if return -end subroutine psb_dprec1 +end subroutine psb_dprecaply1 diff --git a/src/prec/psb_dprecbld.f90 b/src/prec/psb_dprecbld.f90 index 93ec8afd..04c046ae 100644 --- a/src/prec/psb_dprecbld.f90 +++ b/src/prec/psb_dprecbld.f90 @@ -448,6 +448,32 @@ subroutine psb_mlprec_bld(a,desc_a,p,info) integer :: i, nrg, nzg, err_act,k character(len=20) :: name, ch_err + + interface psb_genaggrmap + subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) + use psb_spmat_type + use psb_descriptor_type + implicit none + integer, intent(in) :: aggr_type + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer, pointer :: ilaggr(:),nlaggr(:) + integer, intent(out) :: info + end subroutine psb_dgenaggrmap + end interface + + interface psb_bldaggrmat + subroutine psb_dbldaggrmat(a,desc_a,p,info) + use psb_prec_type + use psb_descriptor_type + use psb_spmat_type + type(psb_dspmat_type), intent(in), target :: a + type(psb_dbase_prec), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + end subroutine psb_dbldaggrmat + end interface + name='psb_mlprec_bld' info=0 call psb_erractionsave(err_act) @@ -478,7 +504,7 @@ subroutine psb_mlprec_bld(a,desc_a,p,info) goto 9999 end if - call psb_bld_aggrmat(a,desc_a,p,info) + call psb_bldaggrmat(a,desc_a,p,info) if(info /= 0) then info=4010 ch_err='psb_bld_aggrmat' diff --git a/src/psblas/Makefile b/src/psblas/Makefile index c28234a6..9c724f62 100644 --- a/src/psblas/Makefile +++ b/src/psblas/Makefile @@ -3,6 +3,7 @@ include ../../Make.inc #FCOPT=-O2 OBJS= psb_ddot.o psb_damax.o psb_dasum.o psb_daxpby.o\ psb_dnrm2.o psb_dnrmi.o psb_dspmm.o psb_dspsm.o\ + pdtreecomb.o LIBDIR=../../lib HERE=. diff --git a/src/psblas/psb_damax.f90 b/src/psblas/psb_damax.f90 index 5eda35b9..305a56cc 100644 --- a/src/psblas/psb_damax.f90 +++ b/src/psblas/psb_damax.f90 @@ -17,6 +17,7 @@ function psb_damax (x,desc_a, info, jx) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -119,6 +120,7 @@ end function psb_damax function psb_damaxv (x,desc_a, info) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -218,6 +220,7 @@ end function psb_damaxv subroutine psb_damaxvs (res,x,desc_a, info, jx) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -321,6 +324,7 @@ end subroutine psb_damaxvs subroutine psb_dmamaxs (res,x,desc_a, info) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none diff --git a/src/psblas/psb_dasum.f90 b/src/psblas/psb_dasum.f90 index 1a1acbed..a12b7c29 100644 --- a/src/psblas/psb_dasum.f90 +++ b/src/psblas/psb_dasum.f90 @@ -17,6 +17,7 @@ function psb_dasum (x,desc_a, info, jx) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -137,6 +138,7 @@ function psb_dasumv (x,desc_a, info) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -252,6 +254,7 @@ end function psb_dasumv subroutine psb_dasumvs (res,x,desc_a, info) use psb_serial_mod use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none diff --git a/src/psblas/psb_daxpby.f90 b/src/psblas/psb_daxpby.f90 index 0e6860fa..84df9950 100644 --- a/src/psblas/psb_daxpby.f90 +++ b/src/psblas/psb_daxpby.f90 @@ -21,6 +21,7 @@ ! subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -144,8 +145,9 @@ end subroutine psb_daxpby ! desc_a - type(). The communication descriptor. ! info - integer. Eventually returns an error code. ! -subroutine psb_psdaxpbyv(alpha, x, beta,y,desc_a,info) +subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -217,4 +219,4 @@ subroutine psb_psdaxpbyv(alpha, x, beta,y,desc_a,info) return end if return -end subroutine psb_psdaxpbyv +end subroutine psb_daxpbyv diff --git a/src/psblas/psb_ddot.f90 b/src/psblas/psb_ddot.f90 index 4598c70c..fe516505 100644 --- a/src/psblas/psb_ddot.f90 +++ b/src/psblas/psb_ddot.f90 @@ -19,7 +19,7 @@ ! function psb_ddot(x, y,desc_a, info, jx, jy) use psb_descriptor_type -! use psb_spsb_mod + use psb_check_mod use psb_error_mod implicit none @@ -149,6 +149,7 @@ end function psb_ddot ! function psb_ddotv(x, y,desc_a, info) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -259,6 +260,7 @@ end function psb_ddotv ! subroutine psb_ddotvs(res, x, y,desc_a, info) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -374,6 +376,7 @@ end subroutine psb_ddotvs ! subroutine psb_dmdots(res, x, y, desc_a, info) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none diff --git a/src/psblas/psb_dnrm2.f90 b/src/psblas/psb_dnrm2.f90 index acb7f177..ae70ce7c 100644 --- a/src/psblas/psb_dnrm2.f90 +++ b/src/psblas/psb_dnrm2.f90 @@ -15,6 +15,7 @@ ! function psb_dnrm2(x, desc_a, info, jx) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -127,6 +128,7 @@ end function psb_dnrm2 ! function psb_dnrm2v(x, desc_a, info) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none @@ -236,6 +238,7 @@ end function psb_dnrm2v ! subroutine psb_dnrm2vs(res, x, desc_a, info) use psb_descriptor_type + use psb_check_mod use psb_error_mod implicit none diff --git a/src/psblas/psb_dnrmi.f90 b/src/psblas/psb_dnrmi.f90 index f235c20a..4065ceb0 100644 --- a/src/psblas/psb_dnrmi.f90 +++ b/src/psblas/psb_dnrmi.f90 @@ -13,6 +13,7 @@ function psb_dnrmi(a,desc_a,info) use psb_descriptor_type use psb_serial_mod + use psb_check_mod use psb_error_mod implicit none diff --git a/src/psblas/psb_dspmm.f90 b/src/psblas/psb_dspmm.f90 index 4f7a8faa..9d95fbd7 100644 --- a/src/psblas/psb_dspmm.f90 +++ b/src/psblas/psb_dspmm.f90 @@ -55,6 +55,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,& use psb_descriptor_type use psb_comm_mod use psi_mod + use psb_check_mod use psb_error_mod implicit none @@ -368,6 +369,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,& use psb_descriptor_type use psb_comm_mod use psi_mod + use psb_check_mod use psb_error_mod implicit none diff --git a/src/psblas/psb_dspsm.f90 b/src/psblas/psb_dspsm.f90 index e2f52947..80cc2edc 100644 --- a/src/psblas/psb_dspsm.f90 +++ b/src/psblas/psb_dspsm.f90 @@ -49,6 +49,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& use psb_descriptor_type use psb_comm_mod use psi_mod + use psb_check_mod use psb_error_mod implicit none @@ -325,6 +326,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& use psb_descriptor_type use psb_comm_mod use psi_mod + use psb_check_mod use psb_error_mod real(kind(1.D0)), intent(in) :: alpha, beta diff --git a/src/serial/coo/Makefile b/src/serial/coo/Makefile index 0e4aaeca..859e542f 100644 --- a/src/serial/coo/Makefile +++ b/src/serial/coo/Makefile @@ -3,7 +3,7 @@ include ../../../Make.inc # # The object files # -FOBJS = dcooprt.o dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o +FOBJS = dcooprt.o dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o #zcsrck.o zcrnrmi.o zcsrmm.o zsrmv.o zcsrsm.o zsrsv.o diff --git a/src/serial/coo/dcoorws.f b/src/serial/coo/dcoorws.f new file mode 100644 index 00000000..0c86f881 --- /dev/null +++ b/src/serial/coo/dcoorws.f @@ -0,0 +1,53 @@ + SUBROUTINE DCOORWS(TRANS,M,N,DESCRA,A,IA1,IA2, + & INFOA,ROWSUM,IERROR) + IMPLICIT NONE +C .. Scalar Arguments .. + INTEGER M,N, IERROR + CHARACTER TRANS +C .. Array Arguments .. + INTEGER IA1(*),IA2(*),INFOA(*) + CHARACTER DESCRA*11 + DOUBLE PRECISION A(*), ROWSUM(*) +C .. Local scalars .. + INTEGER I, J, NNZ, K + DOUBLE PRECISION SUM + logical lsame + external lsame + + NNZ = INFOA(1) + IF (lsame(TRANS,'N')) THEN + DO I=1, M + ROWSUM(I) = 0.0D0 + ENDDO + I = 1 + J = I + DO WHILE (I.LE.NNZ) + + DO WHILE ((IA1(J).EQ.IA1(I)).AND. + + (J.LE.NNZ)) + J = J+1 + ENDDO + + SUM = 0.0 + DO K = I, J-1 + SUM = SUM + ABS(A(K)) + ENDDO + ROWSUM(IA1(I)) = ROWSUM(IA1(I)) + SUM + I = J + ENDDO + + ELSE IF (lsame(TRANS,'T').OR.lsame(TRANS,'C')) THEN + DO J = 1, N + ROWSUM(J) = 0.0D0 + ENDDO + DO I = 1, NNZ + ROWSUM(IA2(I)) = ROWSUM(IA2(I)) + ABS(A(I)) + ENDDO + ELSE + ierror = -1 + ENDIF + RETURN + END + + + diff --git a/src/serial/coo/dcoosm.f b/src/serial/coo/dcoosm.f index 927f5682..0f2996ec 100644 --- a/src/serial/coo/dcoosm.f +++ b/src/serial/coo/dcoosm.f @@ -21,7 +21,7 @@ IF((ALPHA.NE.1.D0) .OR. (BETA.NE.0.D0))then IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 ENDIF if (debug) write(*,*) 'DCOOSM ',m @@ -32,14 +32,14 @@ IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') UPLO = 'L' IF (UPLO.EQ.'?') THEN IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 END IF IF (DESCRA(3:3).EQ.'N') DIAG = 'N' IF (DESCRA(3:3).EQ.'U') DIAG = 'U' IF(UNITD.EQ.'B') THEN IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 ENDIF IF (UNITD.EQ.'R') THEN diff --git a/src/serial/csr/Makefile b/src/serial/csr/Makefile index 76984b94..56366126 100644 --- a/src/serial/csr/Makefile +++ b/src/serial/csr/Makefile @@ -5,7 +5,8 @@ include ../../../Make.inc # FOBJS = dcsrck.o dcsrmm.o dcsrsm.o dcsrmv.o dcsrsv.o dcrnrmi.o \ - dcrcrupd.o dcocrupd.o dcsrprt.o dcsrmv4.o dcsrmv2.o dcsrmv3.o\ + dcrcrupd.o dcocrupd.o dcsrprt.o dcsrmv4.o dcsrmv2.o dcsrmv3.o \ + dcsrrws.o OBJS=$(FOBJS) diff --git a/src/serial/csr/dcsrrws.f b/src/serial/csr/dcsrrws.f new file mode 100644 index 00000000..c22ac193 --- /dev/null +++ b/src/serial/csr/dcsrrws.f @@ -0,0 +1,34 @@ + SUBROUTINE DCSRRWS(TRANS,M,N,DESCRA,A,IA1,IA2, + & INFOA,ROWSUM,IERROR) + IMPLICIT NONE +C .. Scalar Arguments .. + INTEGER M,N, IERROR + CHARACTER TRANS +C .. Array Arguments .. + INTEGER IA1(*),IA2(*),INFOA(*) + CHARACTER DESCRA*11 + DOUBLE PRECISION A(*), ROWSUM(*) +C .. Local scalars .. + INTEGER I, J + + IF (TRANS.EQ.'N') THEN + DO I = 1, M + ROWSUM(I) = 0.0D0 + DO J = IA2(I), IA2(I + 1) - 1 + ROWSUM(I) = ROWSUM(I) + ABS(A(J)) + ENDDO + ENDDO + ELSE IF ((TRANS.EQ.'T').OR.(TRANS.EQ.'C')) THEN + DO J = 1, N + ROWSUM(J) = 0.0D0 + ENDDO + DO I = 1, M + DO J = IA2(I), IA2(I + 1) - 1 + ROWSUM(IA1(J)) = ROWSUM(IA1(J)) + ABS(A(J)) + ENDDO + ENDDO + ENDIF + END + + + diff --git a/src/serial/csr/dcsrsm.f b/src/serial/csr/dcsrsm.f index 28a2bf43..1dd14b40 100644 --- a/src/serial/csr/dcsrsm.f +++ b/src/serial/csr/dcsrsm.f @@ -20,7 +20,7 @@ C .. Local Arrays .. IF((ALPHA.NE.1.D0) .OR. (BETA.NE.0.D0))then IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 ENDIF UPLO = '?' @@ -28,14 +28,14 @@ C .. Local Arrays .. IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') UPLO = 'L' IF (UPLO.EQ.'?') THEN IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 END IF IF (DESCRA(3:3).EQ.'N') DIAG = 'N' IF (DESCRA(3:3).EQ.'U') DIAG = 'U' IF(UNITD.EQ.'B') THEN IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 ENDIF IF (UNITD.EQ.'R') THEN diff --git a/src/serial/dp/dcrjd.f b/src/serial/dp/dcrjd.f index b550e1fe..2c26de38 100644 --- a/src/serial/dp/dcrjd.f +++ b/src/serial/dp/dcrjd.f @@ -165,7 +165,7 @@ C INT_VAL(2) = NZ INT_VAL(3) = LAUX LIAN1 = NZ - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 END IF IF (NZ .GT. LARN) THEN @@ -174,7 +174,7 @@ C INT_VAL(2) = NZ INT_VAL(3) = LAUX LIAN1 = NZ - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 END IF diff --git a/src/serial/f77/Makefile b/src/serial/f77/Makefile index d64a2b14..396acbaf 100644 --- a/src/serial/f77/Makefile +++ b/src/serial/f77/Makefile @@ -5,7 +5,8 @@ include ../../../Make.inc # FOBJS = daxpby.o dcsmm.o dcsnmi.o dcsrp.o dcssm.o \ dcsupd.o dgelp.o dlpupd.o dswmm.o dswprt.o \ - dswsm.o smmp.o ddot.o + dswsm.o smmp.o ddot.o dscal.o dcsrws.o idamax.o \ + dnrm2.o dcopy.o #zcsrck.o zcrnrmi.o zcsrmm.o zsrmv.o zcsrsm.o zsrsv.o diff --git a/src/serial/f77/dcopy.f b/src/serial/f77/dcopy.f new file mode 100644 index 00000000..e1689271 --- /dev/null +++ b/src/serial/f77/dcopy.f @@ -0,0 +1,50 @@ + subroutine dcopy(n,dx,incx,dy,incy) +c +c copies a vector, x, to a vector, y. +c uses unrolled loops for increments equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dy(*) + integer i,incx,incy,ix,iy,m,mp1,n +c + if(n.le.0)return + if(incx.eq.1.and.incy.eq.1)go to 20 +c +c code for unequal increments or equal increments +c not equal to 1 +c + ix = 1 + iy = 1 + if(incx.lt.0)ix = (-n+1)*incx + 1 + if(incy.lt.0)iy = (-n+1)*incy + 1 + do 10 i = 1,n + dy(iy) = dx(ix) + ix = ix + incx + iy = iy + incy + 10 continue + return +c +c code for both increments equal to 1 +c +c +c clean-up loop +c + 20 m = mod(n,7) + if( m .eq. 0 ) go to 40 + do 30 i = 1,m + dy(i) = dx(i) + 30 continue + if( n .lt. 7 ) return + 40 mp1 = m + 1 + do 50 i = mp1,n,7 + dy(i) = dx(i) + dy(i + 1) = dx(i + 1) + dy(i + 2) = dx(i + 2) + dy(i + 3) = dx(i + 3) + dy(i + 4) = dx(i + 4) + dy(i + 5) = dx(i + 5) + dy(i + 6) = dx(i + 6) + 50 continue + return + end diff --git a/src/serial/f77/dcsrws.f b/src/serial/f77/dcsrws.f new file mode 100644 index 00000000..650b5eaf --- /dev/null +++ b/src/serial/f77/dcsrws.f @@ -0,0 +1,153 @@ +C SUBROUTINE DCSRS(TRANS,M,N,FIDA,DESCRA,A,IA1,IA2, & +C & INFOA,ROWSUM,IERROR) +C Purpose +C ======= +C +C Computing sum of absolute values for rows of distributed matrix +C ROWSUM(IX) = ASUM(A(IX, 1..N)) +C IX = 1..M +C +C Parameters +C ========== +C +C TRANS - CHARACTER*1 +C On entry TRANS specifies if the routine operates with matrix A +C or with the transpose of A as follows: +C TRANS = 'N' -> use matrix A +C TRANS = 'T' or 'C' -> use A' (transpose of matrix A) +C Unchanged on exit. +C +C M - INTEGER +C On entry: number of rows of matrix A (A') and +C number of rows of matrix C +C Unchanged on exit. +C +C N - INTEGER +C On entry: number of columns of matrix B +C and number of columns of matrix C. +C Unchanged on exit. +C +C FIDA - CHARACTER*5 +C On entry FIDA defines the format of the input sparse matrix. +C Unchanged on exit. +C +C DESCRA - CHARACTER*1 array of DIMENSION (9) +C On entry DESCRA describes the characteristics of the input +C sparse matrix. +C Unchanged on exit. +C +C A - DOUBLE PRECISION array of DIMENSION (*) +C On entry A specifies the values of the input sparse +C matrix. +C Unchanged on exit. +C +C IA1 - INTEGER array of dimension (*) +C On entry IA1 holds integer information on input sparse +C matrix. Actual information will depend on data format used. +C Unchanged on exit. +C +C IA2 - INTEGER array of dimension (*) +C On entry IA2 holds integer information on input sparse +C matrix. Actual information will depend on data format used. +C Unchanged on exit. +C +C INFOA - INTEGER array of length 10. +C On entry can hold auxiliary information on input matrices +C formats or environment of subsequent calls. +C Might be changed on exit. +C +C IERROR - INTEGER +C On exit IERROR contains the value of error flag as follows: +C IERROR = 0 no error +C IERROR > 0 warning +C IERROR < 0 fatal error +C +C ROWSUM - DOUBLE PRECISION array of dimension (*) +C On exit this vector contains the sum of absolute values +C of elements of a row (AMAX of row array). +C It is required that it has dimension: +C ROWSUM(M) if the subroutine in called with the 'N' option +C ROWSUM(N) in other cases ('T' or 'C' options). +C + SUBROUTINE DCSRWS(TRANS,M,N,FIDA,DESCRA,A,IA1,IA2, + & INFOA,ROWSUM,IERROR) + IMPLICIT NONE +C .. Scalar Arguments .. + INTEGER M,N, IERROR + CHARACTER TRANS +C .. Array Arguments .. + INTEGER IA1(*),IA2(*),INFOA(*) + CHARACTER DESCRA*11, FIDA*5 + DOUBLE PRECISION A(*), ROWSUM(*) +C .. Local Array.. + INTEGER INT_VAL(5), ERR_ACT + DOUBLE PRECISION REAL_VAL(5) + CHARACTER*30 NAME,STRINGS(2) +C .. Parameters .. + DOUBLE PRECISION ZERO + INTEGER IONE + PARAMETER (ZERO=0.D0,IONE=1) +C .. Intrinsic Functions .. + INTRINSIC DBLE, IDINT +C .. Executable Statements .. +C +C Check for argument errors +C + IERROR = 0 + NAME = 'DCSRWS\0' + IF (M.LT.0) THEN + IERROR = 10 + INT_VAL(1) = 2 + INT_VAL(2) = M + ELSE IF (N.LT.0) THEN + IERROR = 10 + INT_VAL(1) = 3 + INT_VAL(2) = N + ELSE IF (TRANS.NE.'T' .AND. TRANS.NE.'N' .AND. TRANS.NE.'C') THEN + IERROR = 40 + INT_VAL(1) = 1 + STRINGS(1) = TRANS//'\0' + ENDIF + +C +C Error handling +C + IF(IERROR.NE.0) THEN + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 + ENDIF + + IF(M.LE.0 .OR. N.LE.0) THEN + GOTO 9999 + ENDIF + + IF (FIDA(1:3).EQ.'CSR') THEN + CALL DCSRRWS(TRANS,M,N,DESCRA,A,IA1,IA2, + + INFOA,ROWSUM,IERROR) + ELSE IF (FIDA(1:3).EQ.'COO') THEN + CALL DCOORWS(TRANS,M,N,DESCRA,A,IA1,IA2, + + INFOA,ROWSUM,IERROR) + ELSE IF (FIDA(1:3).EQ.'JAD') THEN + CALL DJDRWS(TRANS,M,N,DESCRA,A,IA1,IA2, + + INFOA,ROWSUM,IERROR) + ELSE +C +C This data structure not yet considered +C + IERROR = 3010 + strings(1) = fida//'\0' + ENDIF + + CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) + RETURN + + 9999 CONTINUE + CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) + + IF ( ERR_ACT .NE. 0 ) THEN + CALL FCPSB_SERROR() + RETURN + ENDIF + + RETURN + END diff --git a/src/serial/f77/dnrm2.f b/src/serial/f77/dnrm2.f new file mode 100644 index 00000000..119d0477 --- /dev/null +++ b/src/serial/f77/dnrm2.f @@ -0,0 +1,60 @@ + DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX ) +* .. Scalar Arguments .. + INTEGER INCX, N +* .. Array Arguments .. + DOUBLE PRECISION X( * ) +* .. +* +* DNRM2 returns the euclidean norm of a vector via the function +* name, so that +* +* DNRM2 := sqrt( x'*x ) +* +* +* +* -- This version written on 25-October-1982. +* Modified on 14-October-1993 to inline the call to DLASSQ. +* Sven Hammarling, Nag Ltd. +* +* +* .. Parameters .. + DOUBLE PRECISION ONE , ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. Local Scalars .. + INTEGER IX + DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ +* .. Intrinsic Functions .. + INTRINSIC ABS, SQRT +* .. +* .. Executable Statements .. + IF( N.LT.1 .OR. INCX.LT.1 )THEN + NORM = ZERO + ELSE IF( N.EQ.1 )THEN + NORM = ABS( X( 1 ) ) + ELSE + SCALE = ZERO + SSQ = ONE +* The following loop is equivalent to this call to the LAPACK +* auxiliary routine: +* CALL DLASSQ( N, X, INCX, SCALE, SSQ ) +* + DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX + IF( X( IX ).NE.ZERO )THEN + ABSXI = ABS( X( IX ) ) + IF( SCALE.LT.ABSXI )THEN + SSQ = ONE + SSQ*( SCALE/ABSXI )**2 + SCALE = ABSXI + ELSE + SSQ = SSQ + ( ABSXI/SCALE )**2 + END IF + END IF + 10 CONTINUE + NORM = SCALE * SQRT( SSQ ) + END IF +* + DNRM2 = NORM + RETURN +* +* End of DNRM2. +* + END diff --git a/src/serial/f77/dscal.f b/src/serial/f77/dscal.f new file mode 100644 index 00000000..e1467faf --- /dev/null +++ b/src/serial/f77/dscal.f @@ -0,0 +1,43 @@ + subroutine dscal(n,da,dx,incx) +c +c scales a vector by a constant. +c uses unrolled loops for increment equal to one. +c jack dongarra, linpack, 3/11/78. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision da,dx(*) + integer i,incx,m,mp1,n,nincx +c + if( n.le.0 .or. incx.le.0 )return + if(incx.eq.1)go to 20 +c +c code for increment not equal to 1 +c + nincx = n*incx + do 10 i = 1,nincx,incx + dx(i) = da*dx(i) + 10 continue + return +c +c code for increment equal to 1 +c +c +c clean-up loop +c + 20 m = mod(n,5) + if( m .eq. 0 ) go to 40 + do 30 i = 1,m + dx(i) = da*dx(i) + 30 continue + if( n .lt. 5 ) return + 40 mp1 = m + 1 + do 50 i = mp1,n,5 + dx(i) = da*dx(i) + dx(i + 1) = da*dx(i + 1) + dx(i + 2) = da*dx(i + 2) + dx(i + 3) = da*dx(i + 3) + dx(i + 4) = da*dx(i + 4) + 50 continue + return + end diff --git a/src/serial/f77/idamax.f b/src/serial/f77/idamax.f new file mode 100644 index 00000000..59d80dc4 --- /dev/null +++ b/src/serial/f77/idamax.f @@ -0,0 +1,39 @@ + integer function idamax(n,dx,incx) +c +c finds the index of element having max. absolute value. +c jack dongarra, linpack, 3/11/78. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) +c + double precision dx(*),dmax + integer i,incx,ix,n +c + idamax = 0 + if( n.lt.1 .or. incx.le.0 ) return + idamax = 1 + if(n.eq.1)return + if(incx.eq.1)go to 20 +c +c code for increment not equal to 1 +c + ix = 1 + dmax = dabs(dx(1)) + ix = ix + incx + do 10 i = 2,n + if(dabs(dx(ix)).le.dmax) go to 5 + idamax = i + dmax = dabs(dx(ix)) + 5 ix = ix + incx + 10 continue + return +c +c code for increment equal to 1 +c + 20 dmax = dabs(dx(1)) + do 30 i = 2,n + if(dabs(dx(i)).le.dmax) go to 30 + idamax = i + dmax = dabs(dx(i)) + 30 continue + return + end diff --git a/src/serial/f77/smmp.f b/src/serial/f77/smmp.f index 24e936d1..3be67e4c 100644 --- a/src/serial/f77/smmp.f +++ b/src/serial/f77/smmp.f @@ -29,7 +29,7 @@ c symbolic matrix multiply c=a*b c c$$$ write(0,*) 'SYMBMM: ',n,m,l,ib(m+1)-1,jb(ib(m+1)-1) if (size(ic) < n+1) then - call psrealloc(n+1,ic,info) + call psb_realloc(n+1,ic,info) endif maxlmn = max(l,m,n) do 10 i=1,maxlmn diff --git a/src/serial/jad/Makefile b/src/serial/jad/Makefile index 2b15d4d8..590719f0 100644 --- a/src/serial/jad/Makefile +++ b/src/serial/jad/Makefile @@ -4,7 +4,7 @@ include ../../../Make.inc # FOBJS = dcojdupd.o djadmm.o djadmv.o djadsm.o djadsv.o djdnrmi.o djadnr.o djadprt.o\ - djadmv2.o djadmv3.o djadmv4.o + djadmv2.o djadmv3.o djadmv4.o djadrws.o djdrws.o OBJS=$(FOBJS) diff --git a/src/serial/jad/djadrws.f b/src/serial/jad/djadrws.f new file mode 100644 index 00000000..fa9138a7 --- /dev/null +++ b/src/serial/jad/djadrws.f @@ -0,0 +1,48 @@ +C ... Compute infinity norma for sparse matrix in CSR Format ... + SUBROUTINE DJADRWS(TRANS,M,N,NG,A,KA,JA,IA, + + INFOA,ROWSUM,IERROR) + IMPLICIT NONE + INCLUDE 'psb_const.fh' +C .. Scalar Arguments .. + INTEGER M,N, IERROR, NG + CHARACTER TRANS +C .. Array Arguments .. + INTEGER KA(*),JA(*),IA(3,*),INFOA(*) + DOUBLE PRECISION A(*), rowsum(*) +C ... Local Scalars .. + DOUBLE PRECISION NRMI + INTEGER I, IR, K, IPG, NPG, IPX + + NRMI = 0.0 + IR = 0 + DO IPG = 1, NG + K = IA(2,IPG) + NPG = JA(K+1)- JA(K) + +C ... ... + DO I = 1, NPG + ROWSUM(IR+I) = 0.0 + ENDDO + + DO K = IA(2,IPG), IA(3,IPG)-1 + IPX = 1 + DO I = JA(K), JA(K+1) - 1 + ROWSUM(IR+IPX) = ROWSUM(IR+IPX) + ABS(A(I)) + IPX = IPX + 1 + ENDDO + ENDDO + +C ... CSR Representation ... + + IPX = 1 + DO K = IA(3,IPG), IA(2,IPG+1)-1 + DO I = JA(K), JA(K+1) - 1 + ROWSUM(IR+IPX) = ROWSUM(IR+IPX) + ABS(A(I)) + ENDDO + IPX = IPX + 1 + ENDDO + + IR = IR + NPG + ENDDO + + END diff --git a/src/serial/jad/djadsm.f b/src/serial/jad/djadsm.f index 6f7a05e5..0628d7e3 100644 --- a/src/serial/jad/djadsm.f +++ b/src/serial/jad/djadsm.f @@ -26,7 +26,7 @@ C IF((ALPHA.NE.1.D0) .OR. (BETA.NE.0.D0))then IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 ENDIF UPLO = '?' @@ -35,13 +35,13 @@ C C IF (UPLO.EQ.'?') THEN IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 END IF IF (DESCRA(3:3).NE.'U') THEN IERROR=5 - CALL PSB_ERRPUSH(IERROR,NAME,INT_VAL) + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 END IF UNITD=DESCRA(3:3) diff --git a/src/serial/jad/djdrws.f b/src/serial/jad/djdrws.f new file mode 100644 index 00000000..da7f6230 --- /dev/null +++ b/src/serial/jad/djdrws.f @@ -0,0 +1,30 @@ +C ... Compute infinity norm for sparse matrix in CSR Format ... + SUBROUTINE DJDRWS(TRANS,M,N,DESCRA,A,JA,IA, + + INFOA,ROWSUM,IERROR) + IMPLICIT NONE +C .. Scalar Arguments .. + INTEGER M,N, IERROR + CHARACTER TRANS +C .. Array Arguments .. + INTEGER JA(*),IA(*),INFOA(*) + CHARACTER DESCRA*11 + DOUBLE PRECISION A(*), ROWSUM(*) +C .. Local scalars .. + INTEGER PNG, PIA, PJA +C .. External routines .. + DOUBLE PRECISION DJADNR + EXTERNAL DJADNR + + IERROR = 0 + PNG = IA(1) + PIA = IA(2) + PJA = IA(3) + + IF (DESCRA(1:1).EQ.'G') THEN + CALL DJADRWS(TRANS,M,N,IA(PNG), + + A,JA,IA(PJA),IA(PIA), + + INFOA,ROWSUM,IERROR) + ELSE + IERROR = 3011 + ENDIF + END diff --git a/src/tools/psb_descasb.f90 b/src/tools/psb_descasb.f90 index 46b06128..348e5c30 100644 --- a/src/tools/psb_descasb.f90 +++ b/src/tools/psb_descasb.f90 @@ -633,7 +633,7 @@ Subroutine psb_descasb(n_ovr,desc_p,desc_a,a,& end if if (.false.) then - call descprt(70+myrow,desc_p,.false.) + call psb_descprt(70+myrow,desc_p,.false.) end if if (debug) write(0,*) myrow,'Done ConvertComm'