diff --git a/base/newserial/Makefile b/base/newserial/Makefile index 06b297f8..2203c141 100644 --- a/base/newserial/Makefile +++ b/base/newserial/Makefile @@ -1,7 +1,7 @@ include ../../Make.inc -MODULES = psbn_base_mat_mod.o psbn_d_base_mat_mod.o psbn_d_coo_impl.o psbn_mat_mod.o -# psbn_csr_mat.o +MODULES = psbn_base_mat_mod.o psbn_d_base_mat_mod.o psbn_d_coo_impl.o psbn_mat_mod.o\ + psbn_csr_mat.o psbn_d_csr_impl.o LIBMOD= @@ -19,6 +19,7 @@ lib: $(MODULES) $(OBJS) $(LIBMOD) psbn_mat_mod.o: psbn_base_mat_mod.o psbn_coo_mat.o psbn_csr_mat.o: psbn_d_base_mat_mod.o +psbn_d_csr_impl.o: psbn_d_csr_mat_mod.o psbn_d_mat_impl.o: psbn_d_base_mat_mod.o clean: diff --git a/base/newserial/psbn_csr_mat.f03 b/base/newserial/psbn_csr_mat.f03 index 54e87b10..3675e9fc 100644 --- a/base/newserial/psbn_csr_mat.f03 +++ b/base/newserial/psbn_csr_mat.f03 @@ -4,7 +4,6 @@ module psbn_d_csr_sparse_mat_mod type, extends(psbn_d_base_sparse_mat) :: psbn_d_csr_sparse_mat - logical :: sorted integer, allocatable :: irp(:), ja(:) real(psb_dpk_), allocatable :: val(:) @@ -16,9 +15,93 @@ module psbn_d_csr_sparse_mat_mod procedure, pass(a) :: d_base_cssv => d_csr_cssv procedure, pass(a) :: reallocate_nz => d_csr_reallocate_nz procedure, pass(a) :: csins => d_csr_csins + procedure, pass(a) :: allocate_mnnz => d_csr_allocate_mnnz + procedure, pass(a) :: allocate_mn => d_csr_allocate_mn + procedure, pass(a) :: to_coo => d_csr_to_coo + procedure, pass(a) :: from_coo => d_csr_from_coo + procedure, pass(a) :: free => d_csr_free + end type psbn_d_csr_sparse_mat + private :: d_csr_get_nzeros, d_csr_csmm, d_csr_csmv, d_csr_cssm, d_csr_cssv, & + & d_csr_csins, d_csr_reallocate_nz, d_csr_allocate_mnnz, & + & d_csr_allocate_mn, d_csr_to_coo, d_csr_from_coo, & + & d_csr_free + + + interface + subroutine d_csr_to_coo_impl(a,b,info) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine d_csr_to_coo_impl + end interface + + interface + subroutine d_csr_from_coo_impl(a,b,info) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + end subroutine d_csr_from_coo_impl + end interface + + interface + subroutine d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + end subroutine d_csr_csins_impl + end interface + + interface d_csr_cssm_impl + subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_csr_cssv_impl + subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_csr_cssm_impl + end interface + + interface d_csr_csmm_impl + subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_csr_csmv_impl + subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_csr_sparse_mat + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_csr_csmm_impl + end interface - end type psbn_d_csr_sparse_mat contains @@ -108,39 +191,8 @@ contains if (nz == 0) return - - nza = a%get_nzeros() - - if (a%is_bld()) then - ! Build phase should only ever be in COO - info = 1121 - - else if (a%is_upd()) then - if (a%is_sorted()) then - - -!!$#ifdef FIXED_NAG_SEGV -!!$ call d_csr_srch_upd(nz,ia,ja,val,a,& -!!$ & imin,imax,jmin,jmax,info,gtl) -!!$#else - call d_csr_srch_upd(nz,ia,ja,val,& - & a%irp,a%ja,a%val,& - & a%get_dupl(),a%get_nzeros(),a%get_nrows(),& - & info,gtl) -!!$#endif - - else - info = 1121 - end if - - else - ! State is wrong. - info = 1121 - end if - if (info /= 0) then - call psb_errpush(info,name) - goto 9999 - end if + call d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + if (info /= 0) goto 9999 call psb_erractionrestore(err_act) return @@ -153,424 +205,9 @@ contains return end if return - - - contains - - -!!$#ifdef FIXED_NAG_SEGV -!!$ subroutine d_csr_srch_upd(nz,ia,ja,val,a,& -!!$ & imin,imax,jmin,jmax,info,gtl) -!!$ -!!$ use psb_const_mod -!!$ use psb_realloc_mod -!!$ use psb_string_mod -!!$ use psb_serial_mod -!!$ implicit none -!!$ -!!$ class(psbn_d_csr_sparse_mat), intent(inout) :: a -!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax -!!$ integer, intent(in) :: ia(:),ja(:) -!!$ real(psb_dpk_), intent(in) :: val(:) -!!$ integer, intent(out) :: info -!!$ integer, intent(in), optional :: gtl(:) -!!$ integer :: i,ir,ic, ilr, ilc, ip, & -!!$ & i1,i2,nc,nnz,dupl,ng -!!$ integer :: debug_level, debug_unit -!!$ character(len=20) :: name='d_csr_srch_upd' -!!$ -!!$ info = 0 -!!$ debug_unit = psb_get_debug_unit() -!!$ debug_level = psb_get_debug_level() -!!$ -!!$ dupl = a%get_dupl() -!!$ -!!$ if (.not.a%is_sorted()) then -!!$ info = -4 -!!$ return -!!$ end if -!!$ -!!$ ilr = -1 -!!$ ilc = -1 -!!$ nnz = a%get_nzeros() -!!$ -!!$ -!!$ if (present(gtl)) then -!!$ ng = size(gtl) -!!$ -!!$ select case(dupl) -!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_) -!!$ ! Overwrite. -!!$ ! Cannot test for error, should have been caught earlier. -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then -!!$ ir = gtl(ir) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ ic = gtl(ic) -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ a%val(i1+ip-1) = val(i) -!!$ else -!!$ info = i -!!$ return -!!$ end if -!!$ else -!!$ if (debug_level >= psb_debug_serial_) & -!!$ & write(debug_unit,*) trim(name),& -!!$ & ': Discarding row that does not belong to us.' -!!$ endif -!!$ end if -!!$ end do -!!$ case(psbn_dupl_add_) -!!$ ! Add -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then -!!$ ir = gtl(ir) -!!$ ic = gtl(ic) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) -!!$ else -!!$ info = i -!!$ return -!!$ end if -!!$ else -!!$ if (debug_level >= psb_debug_serial_) & -!!$ & write(debug_unit,*) trim(name),& -!!$ & ': Discarding row that does not belong to us.' -!!$ end if -!!$ end if -!!$ end do -!!$ -!!$ case default -!!$ info = -3 -!!$ if (debug_level >= psb_debug_serial_) & -!!$ & write(debug_unit,*) trim(name),& -!!$ & ': Duplicate handling: ',dupl -!!$ end select -!!$ -!!$ else -!!$ -!!$ select case(dupl) -!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_) -!!$ ! Overwrite. -!!$ ! Cannot test for error, should have been caught earlier. -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ a%val(i1+ip-1) = val(i) -!!$ else -!!$ info = i -!!$ return -!!$ end if -!!$ end if -!!$ end do -!!$ -!!$ case(psbn_dupl_add_) -!!$ ! Add -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) -!!$ else -!!$ info = i -!!$ return -!!$ end if -!!$ end if -!!$ end do -!!$ -!!$ case default -!!$ info = -3 -!!$ if (debug_level >= psb_debug_serial_) & -!!$ & write(debug_unit,*) trim(name),& -!!$ & ': Duplicate handling: ',dupl -!!$ end select -!!$ -!!$ end if -!!$ -!!$ end subroutine d_csr_srch_upd -!!$ -!!$#else - subroutine d_csr_srch_upd(nz,ia,ja,val,& - & airp,aja,aval,dupl,nza,nra,& - & info,gtl) - - use psb_error_mod - use psb_sort_mod - implicit none - - integer, intent(inout) :: airp(:),aja(:) - real(psb_dpk_), intent(inout) :: aval(:) - integer, intent(in) :: nz, dupl,nza, nra - integer, intent(in) :: ia(:),ja(:) - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - integer, intent(in), optional :: gtl(:) - - integer :: debug_level, debug_unit - character(len=20) :: name='d_csr_srch_upd' - integer :: i,ir,ic, ilr, ilc, ip, & - & i1,i2,nc,lb,ub,m, ng - - info = 0 - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - - if (present(gtl)) then - ng = size(gtl) - - select case(dupl) - case(psbn_dupl_ovwrt_,psbn_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then - i1 = airp(ir) - i2 = airp(ir+1) - nc=i2-i1 - - ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) - if (ip>0) then - aval(i1+ip-1) = val(i) - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Was searching ',ic,' in: ',i1,i2,& - & ' : ',aja(i1:i2-1) - info = i - return - end if - - else - - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - end if - end do - - case(psbn_dupl_add_) - ! Add - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then - i1 = airp(ir) - i2 = airp(ir+1) - nc = i2-i1 - ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) - if (ip>0) then - aval(i1+ip-1) = aval(i1+ip-1) + val(i) - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Was searching ',ic,' in: ',i1,i2,& - & ' : ',aja(i1:i2-1) - info = i - return - end if - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - - end if - end do - - case default - info = -3 - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Duplicate handling: ',dupl - end select - - else - - select case(dupl) - case(psbn_dupl_ovwrt_,psbn_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - - if ((ir > 0).and.(ir <= a%m)) then - - i1 = airp(ir) - i2 = airp(ir+1) - nc=i2-i1 - - ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) - if (ip>0) then - aval(i1+ip-1) = val(i) - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Was searching ',ic,' in: ',i1,i2,& - & ' : ',aja(i1:i2-1) - info = i - return - end if - - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - - end do - - case(psbn_dupl_add_) - ! Add - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then - i1 = airp(ir) - i2 = airp(ir+1) - nc = i2-i1 - ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) - if (ip>0) then - aval(i1+ip-1) = aval(i1+ip-1) + val(i) - else - info = i - return - end if - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - end do - - case default - info = -3 - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Duplicate handling: ',dupl - end select - - end if - end subroutine d_csr_srch_upd -!!$#endif end subroutine d_csr_csins - - - - subroutine d_csr_csmv(alpha,a,x,beta,y,info,trans) use psb_error_mod class(psbn_d_csr_sparse_mat), intent(in) :: a @@ -578,7 +215,7 @@ contains real(psb_dpk_), intent(inout) :: y(:) integer, intent(out) :: info character, optional, intent(in) :: trans - + character :: trans_ integer :: i,j,k,m,n, nnz, ir, jc real(psb_dpk_) :: acc @@ -588,240 +225,21 @@ contains logical, parameter :: debug=.false. call psb_erractionsave(err_act) - info = 0 - - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if if (.not.a%is_asb()) then - write(0,*) 'Error: csmv called on an unassembled mat' info = 1121 call psb_errpush(info,name) goto 9999 endif - tra = ((trans_=='T').or.(trans_=='t')) - - if (tra) then - m = a%get_ncols() - n = a%get_nrows() - else - n = a%get_ncols() - m = a%get_nrows() - end if - - - if (alpha == dzero) then - if (beta == dzero) then - do i = 1, m - y(i) = dzero - enddo - else - do i = 1, m - y(i) = beta*y(i) - end do - endif - return - end if - - if (.not.tra) then - - if (beta == dzero) then - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = -acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = alpha*acc - end do - - end if - - - else if (beta == done) then - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = y(i) + acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = y(i) -acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = y(i) + alpha*acc - end do - - end if - - else if (beta == -done) then - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = -y(i) + acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = -y(i) -acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = -y(i) + alpha*acc - end do - - end if - - else - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = beta*y(i) + acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = beta*y(i) - acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j)) - enddo - y(i) = beta*y(i) + alpha*acc - end do - - end if - - end if - - else if (tra) then - - if (beta == dzero) then - do i=1, m - y(i) = dzero - end do - else if (beta == done) then - ! Do nothing - else if (beta == -done) then - do i=1, m - y(i) = -y(i) - end do - else - do i=1, m - y(i) = beta*y(i) - end do - end if - - if (alpha.eq.done) then - - do i=1,n - do j=a%irp(i), a%irp(i+1)-1 - ir = a%ja(j) - y(ir) = y(ir) + a%val(j)*x(i) - end do - enddo - - else if (alpha.eq.-done) then - - do i=1,n - do j=a%irp(i), a%irp(i+1)-1 - ir = a%ja(j) - y(ir) = y(ir) - a%val(j)*x(i) - end do - enddo - - else - - do i=1,n - do j=a%irp(i), a%irp(i+1)-1 - ir = a%ja(j) - y(ir) = y(ir) + alpha*a%val(j)*x(i) - end do - enddo - - end if - - endif - - if (a%is_unit()) then - do i=1, min(m,n) - y(i) = y(i) + alpha*x(i) - end do - end if + call d_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) goto 9999 call psb_erractionrestore(err_act) return + 9999 continue call psb_erractionrestore(err_act) @@ -851,243 +269,15 @@ contains call psb_erractionsave(err_act) - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if - - tra = ((trans_=='T').or.(trans_=='t')) - if (.not.a%is_asb()) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - endif - - if (tra) then - m = a%get_ncols() - n = a%get_nrows() - else - n = a%get_ncols() - m = a%get_nrows() - end if - - nc = min(size(x,2) , size(y,2) ) - allocate(acc(nc), stat=info) - if(info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='allocate') - goto 9999 - end if - if (alpha == dzero) then - if (beta == dzero) then - do i = 1, m - y(i,:) = dzero - enddo - else - do i = 1, m - y(i,:) = beta*y(i,:) - end do - endif - return - end if - - if (.not.tra) then - - if (beta == dzero) then - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = -acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = alpha*acc - end do - - end if - - - else if (beta == done) then - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = y(i,:) + acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = y(i,:) -acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = y(i,:) + alpha*acc - end do - - end if - - else if (beta == -done) then - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = -y(i,:) + acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = -y(i,:) -acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = -y(i,:) + alpha*acc - end do - - end if - - else - - if (alpha == done) then - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = beta*y(i,:) + acc - end do - - else if (alpha == -done) then - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = beta*y(i,:) - acc - end do - - else - - do i=1,m - acc = zero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j) * x(a%ja(j),:) - enddo - y(i,:) = beta*y(i,:) + alpha*acc - end do - - end if - - end if - - else if (tra) then - - if (beta == dzero) then - do i=1, m - y(i,:) = dzero - end do - else if (beta == done) then - ! Do nothing - else if (beta == -done) then - do i=1, m - y(i,:) = -y(i,:) - end do - else - do i=1, m - y(i,:) = beta*y(i,:) - end do - end if - - if (alpha.eq.done) then - - do i=1,n - do j=a%irp(i), a%irp(i+1)-1 - ir = a%ja(j) - y(ir,:) = y(ir,:) + a%val(j)*x(i,:) - end do - enddo - - else if (alpha.eq.-done) then - - do i=1,n - do j=a%irp(i), a%irp(i+1)-1 - ir = a%ja(j) - y(ir,:) = y(ir,:) - a%val(j)*x(i,:) - end do - enddo - - else - - do i=1,n - do j=a%irp(i), a%irp(i+1)-1 - ir = a%ja(j) - y(ir,:) = y(ir,:) + alpha*a%val(j)*x(i,:) - end do - enddo - - end if - - endif - - if (a%is_unit()) then - do i=1, min(m,n) - y(i,:) = y(i,:) + alpha*x(i,:) - end do - end if + call d_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + if (info /= 0) goto 9999 + call psb_erractionrestore(err_act) return + 9999 continue call psb_erractionrestore(err_act) @@ -1119,61 +309,26 @@ contains call psb_erractionsave(err_act) - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if if (.not.a%is_asb()) then info = 1121 call psb_errpush(info,name) goto 9999 endif - - tra = ((trans_=='T').or.(trans_=='t')) - m = a%get_nrows() + if (.not. (a%is_triangle())) then + write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 end if - - - if (alpha == dzero) then - if (beta == dzero) then - do i = 1, m - y(i) = dzero - enddo - else - do i = 1, m - y(i) = beta*y(i) - end do - endif - return - end if - - if (beta == dzero) then - call inner_csrsv(tra,a,x,y) - do i = 1, m - y(i) = alpha*y(i) - end do - else - allocate(tmp(m), stat=info) - if (info /= 0) then - write(0,*) 'Memory allocation error in CSRSV ' - return - end if - tmp(1:m) = x(1:m) - call inner_csrsv(tra,a,tmp,y) - do i = 1, m - y(i) = alpha*tmp(i) + beta*y(i) - end do - end if + call d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + call psb_erractionrestore(err_act) return + 9999 continue call psb_erractionrestore(err_act) @@ -1183,108 +338,6 @@ contains end if return - contains - - subroutine inner_csrsv(tra,a,x,y) - logical, intent(in) :: tra - class(psbn_d_csr_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: x(:) - real(psb_dpk_), intent(out) :: y(:) - - integer :: i,j,k,m, ir, jc - real(psb_dpk_) :: acc - - if (.not.tra) then - - if (a%is_lower()) then - if (a%is_unit()) then - do i=1, a%get_nrows() - acc = dzero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j)) - end do - y(i) = x(i) - acc - end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() - acc = dzero - do j=a%irp(i), a%irp(i+1)-2 - acc = acc + a%val(j)*x(a%ja(j)) - end do - y(i) = (x(i) - acc)/a%val(a%irp(i+1)-1) - end do - end if - else if (a%is_upper()) then - - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = dzero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j)) - end do - y(i) = x(i) - acc - end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = dzero - do j=a%irp(i)+1, a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j)) - end do - y(i) = (x(i) - acc)/a%val(a%irp(i)) - end do - end if - - end if - - else if (tra) then - - do i=1, a%get_nrows() - y(i) = x(i) - end do - - if (a%is_lower()) then - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = y(i) - do j=a%irp(i), a%irp(i+1)-1 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - end do - end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 - y(i) = y(i)/a%val(a%irp(i+1)-1) - acc = y(i) - do j=a%irp(i), a%irp(i+1)-2 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - end do - end do - end if - else if (a%is_upper()) then - - if (a%is_unit()) then - do i=1, a%get_nrows() - acc = y(i) - do j=a%irp(i), a%irp(i+1)-1 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - end do - end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() - y(i) = y(i)/a%val(a%irp(i)) - acc = y(i) - do j=a%irp(i)+1, a%irp(i+1)-1 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - end do - end do - end if - - end if - end if - end subroutine inner_csrsv end subroutine d_csr_cssv @@ -1304,26 +357,17 @@ contains real(psb_dpk_), allocatable :: tmp(:,:) logical :: tra Integer :: err_act - character(len=20) :: name='d_base_csmm' + character(len=20) :: name='d_csr_csmm' logical, parameter :: debug=.false. call psb_erractionsave(err_act) - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if if (.not.a%is_asb()) then info = 1121 call psb_errpush(info,name) goto 9999 endif - - tra = ((trans_=='T').or.(trans_=='t')) - m = a%get_nrows() - nc = min(size(x,2) , size(y,2)) if (.not. (a%is_triangle())) then write(0,*) 'Called SM on a non-triangular mat!' @@ -1331,51 +375,144 @@ contains call psb_errpush(info,name) goto 9999 end if + + call d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + call psb_erractionrestore(err_act) + return - if (alpha == dzero) then - if (beta == dzero) then - do i = 1, m - y(i,:) = dzero - enddo - else - do i = 1, m - y(i,:) = beta*y(i,:) - end do - endif +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() return end if + return + + end subroutine d_csr_cssm + + + subroutine d_csr_free(a) + + class(psbn_d_csr_sparse_mat), intent(inout) :: a + + if (allocated(a%irp)) deallocate(a%irp) + if (allocated(a%ja)) deallocate(a%ja) + if (allocated(a%val)) deallocate(a%val) + call a%set_null() + call a%set_nrows(0) + call a%set_ncols(0) + + return + + end subroutine d_csr_free + + + subroutine d_csr_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_csr_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call d_csr_to_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine d_csr_to_coo + + subroutine d_csr_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_csr_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info - if (beta == dzero) then - call inner_csrsm(tra,a,x,y,info) - do i = 1, m - y(i,:) = alpha*y(i,:) - end do - else - allocate(tmp(m,nc), stat=info) - if(info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='allocate') - goto 9999 - end if - - tmp(1:m,:) = x(1:m,:) - call inner_csrsm(tra,a,tmp,y,info) - do i = 1, m - y(i,:) = alpha*tmp(i,:) + beta*y(i,:) - end do + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call d_csr_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() end if + return + + end subroutine d_csr_from_coo + - if(info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='inner_csrsm') + subroutine d_csr_allocate_mnnz(m,n,nz,a) + use psb_error_mod + use psb_realloc_mod + integer, intent(in) :: m,n,nz + class(psbn_d_csr_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='allocate_mnz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + if (m < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) + goto 9999 + endif + if (n < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) goto 9999 + endif + if (nz < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/3,0,0,0,0/)) + goto 9999 + endif + + if (info == 0) call psb_realloc(m+1,a%irp,info) + if (info == 0) call psb_realloc(nz,a%ja,info) + if (info == 0) call psb_realloc(nz,a%val,info) + if (info == 0) then + call a%set_nrows(m) + call a%set_ncols(n) + call a%set_nzeros(0) + call a%set_bld() + call a%set_triangle(.false.) end if call psb_erractionrestore(err_act) return - 9999 continue call psb_erractionrestore(err_act) @@ -1384,119 +521,49 @@ contains return end if return - - contains + end subroutine d_csr_allocate_mnnz - subroutine inner_csrsm(tra,a,x,y,info) - logical, intent(in) :: tra - class(psbn_d_csr_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: x(:,:) - real(psb_dpk_), intent(out) :: y(:,:) - integer, intent(out) :: info - integer :: i,j,k,m, ir, jc - real(psb_dpk_), allocatable :: acc(:) - - allocate(acc(size(x,2)), stat=info) - if(info /= 0) then - info=4010 - return - end if - - - if (.not.tra) then - - if (a%is_lower()) then - if (a%is_unit()) then - do i=1, a%get_nrows() - acc = dzero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j),:) - end do - y(i,:) = x(i,:) - acc - end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() - acc = dzero - do j=a%irp(i), a%irp(i+1)-2 - acc = acc + a%val(j)*x(a%ja(j),:) - end do - y(i,:) = (x(i,:) - acc)/a%val(a%irp(i+1)-1) - end do - end if - else if (a%is_upper()) then - - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = dzero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j),:) - end do - y(i,:) = x(i,:) - acc - end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = dzero - do j=a%irp(i)+1, a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j),:) - end do - y(i,:) = (x(i,:) - acc)/a%val(a%irp(i)) - end do - end if - - end if - - else if (tra) then - - do i=1, a%get_nrows() - y(i,:) = x(i,:) - end do - - if (a%is_lower()) then - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = y(i,:) - do j=a%irp(i), a%irp(i+1)-1 - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - end do - end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 - y(i,:) = y(i,:)/a%val(a%irp(i+1)-1) - acc = y(i,:) - do j=a%irp(i), a%irp(i+1)-2 - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - end do - end do - end if - else if (a%is_upper()) then - - if (a%is_unit()) then - do i=1, a%get_nrows() - acc = y(i,:) - do j=a%irp(i), a%irp(i+1)-1 - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - end do - end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() - y(i,:) = y(i,:)/a%val(a%irp(i)) - acc = y(i,:) - do j=a%irp(i)+1, a%irp(i+1)-1 - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - end do - end do - end if - - end if - end if - end subroutine inner_csrsm - end subroutine d_csr_cssm + subroutine d_csr_allocate_mn(m,n,a) + use psb_error_mod + use psb_realloc_mod + integer, intent(in) :: m,n + class(psbn_d_csr_sparse_mat), intent(inout) :: a + Integer :: err_act, info, nz + character(len=20) :: name='allocate_mn' + logical, parameter :: debug=.false. -end module psbn_d_csr_sparse_mat_mod + call psb_erractionsave(err_act) + info = 0 + if (m < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) + goto 9999 + endif + if (n < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) + goto 9999 + endif + + nz = max(7*m,7*n,1) + call a%allocate(m,n,nz) + + call psb_erractionrestore(err_act) + return +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_csr_allocate_mn + + + +end module psbn_d_csr_sparse_mat_mod diff --git a/base/newserial/psbn_d_base_mat_mod.f03 b/base/newserial/psbn_d_base_mat_mod.f03 index 19df5303..91a22bc4 100644 --- a/base/newserial/psbn_d_base_mat_mod.f03 +++ b/base/newserial/psbn_d_base_mat_mod.f03 @@ -809,17 +809,7 @@ contains endif nz = max(7*m,7*n,1) - - if (info == 0) call psb_realloc(nz,a%ia,info) - if (info == 0) call psb_realloc(nz,a%ja,info) - if (info == 0) call psb_realloc(nz,a%val,info) - if (info == 0) then - call a%set_nrows(m) - call a%set_ncols(n) - call a%set_nzeros(0) - call a%set_bld() - call a%set_triangle(.false.) - end if + call a%allocate(m,n,nz) call psb_erractionrestore(err_act) return diff --git a/base/newserial/psbn_d_csr_impl.f03 b/base/newserial/psbn_d_csr_impl.f03 new file mode 100644 index 00000000..b6439cb0 --- /dev/null +++ b/base/newserial/psbn_d_csr_impl.f03 @@ -0,0 +1,1408 @@ + +subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csmv_impl + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='d_csr_csmv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + if (.not.a%is_asb()) then + write(0,*) 'Error: csmv called on an unassembled mat' + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + tra = ((trans_=='T').or.(trans_=='t')) + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + + + if (alpha == dzero) then + if (beta == dzero) then + do i = 1, m + y(i) = dzero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + end if + + if (.not.tra) then + + if (beta == dzero) then + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = -acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = alpha*acc + end do + + end if + + + else if (beta == done) then + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = y(i) + acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = y(i) -acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = y(i) + alpha*acc + end do + + end if + + else if (beta == -done) then + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = -y(i) + acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = -y(i) -acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = -y(i) + alpha*acc + end do + + end if + + else + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = beta*y(i) + acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = beta*y(i) - acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j)) + enddo + y(i) = beta*y(i) + alpha*acc + end do + + end if + + end if + + else if (tra) then + + if (beta == dzero) then + do i=1, m + y(i) = dzero + end do + else if (beta == done) then + ! Do nothing + else if (beta == -done) then + do i=1, m + y(i) = -y(i) + end do + else + do i=1, m + y(i) = beta*y(i) + end do + end if + + if (alpha.eq.done) then + + do i=1,n + do j=a%irp(i), a%irp(i+1)-1 + ir = a%ja(j) + y(ir) = y(ir) + a%val(j)*x(i) + end do + enddo + + else if (alpha.eq.-done) then + + do i=1,n + do j=a%irp(i), a%irp(i+1)-1 + ir = a%ja(j) + y(ir) = y(ir) - a%val(j)*x(i) + end do + enddo + + else + + do i=1,n + do j=a%irp(i), a%irp(i+1)-1 + ir = a%ja(j) + y(ir) = y(ir) + alpha*a%val(j)*x(i) + end do + enddo + + end if + + endif + + if (a%is_unit()) then + do i=1, min(m,n) + y(i) = y(i) + alpha*x(i) + end do + end if + + call psb_erractionrestore(err_act) + return +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine d_csr_csmv_impl + +subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csmm_impl + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_dpk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_csr_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + tra = ((trans_=='T').or.(trans_=='t')) + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + + nc = min(size(x,2) , size(y,2) ) + + allocate(acc(nc), stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + if (alpha == dzero) then + if (beta == dzero) then + do i = 1, m + y(i,:) = dzero + enddo + else + do i = 1, m + y(i,:) = beta*y(i,:) + end do + endif + return + end if + + if (.not.tra) then + + if (beta == dzero) then + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = -acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = alpha*acc + end do + + end if + + + else if (beta == done) then + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = y(i,:) + acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = y(i,:) -acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = y(i,:) + alpha*acc + end do + + end if + + else if (beta == -done) then + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = -y(i,:) + acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = -y(i,:) -acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = -y(i,:) + alpha*acc + end do + + end if + + else + + if (alpha == done) then + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = beta*y(i,:) + acc + end do + + else if (alpha == -done) then + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = beta*y(i,:) - acc + end do + + else + + do i=1,m + acc = zero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j) * x(a%ja(j),:) + enddo + y(i,:) = beta*y(i,:) + alpha*acc + end do + + end if + + end if + + else if (tra) then + + if (beta == dzero) then + do i=1, m + y(i,:) = dzero + end do + else if (beta == done) then + ! Do nothing + else if (beta == -done) then + do i=1, m + y(i,:) = -y(i,:) + end do + else + do i=1, m + y(i,:) = beta*y(i,:) + end do + end if + + if (alpha.eq.done) then + + do i=1,n + do j=a%irp(i), a%irp(i+1)-1 + ir = a%ja(j) + y(ir,:) = y(ir,:) + a%val(j)*x(i,:) + end do + enddo + + else if (alpha.eq.-done) then + + do i=1,n + do j=a%irp(i), a%irp(i+1)-1 + ir = a%ja(j) + y(ir,:) = y(ir,:) - a%val(j)*x(i,:) + end do + enddo + + else + + do i=1,n + do j=a%irp(i), a%irp(i+1)-1 + ir = a%ja(j) + y(ir,:) = y(ir,:) + alpha*a%val(j)*x(i,:) + end do + enddo + + end if + + endif + + if (a%is_unit()) then + do i=1, min(m,n) + y(i,:) = y(i,:) + alpha*x(i,:) + end do + end if + + call psb_erractionrestore(err_act) + return +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine d_csr_csmm_impl + + +subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_csr_cssv_impl + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_csr_cssv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + tra = ((trans_=='T').or.(trans_=='t')) + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + if (alpha == dzero) then + if (beta == dzero) then + do i = 1, m + y(i) = dzero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + end if + + if (beta == dzero) then + call inner_csrsv(tra,a,x,y) + do i = 1, m + y(i) = alpha*y(i) + end do + else + allocate(tmp(m), stat=info) + if (info /= 0) then + write(0,*) 'Memory allocation error in CSRSV ' + return + end if + tmp(1:m) = x(1:m) + call inner_csrsv(tra,a,tmp,y) + do i = 1, m + y(i) = alpha*tmp(i) + beta*y(i) + end do + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +contains + + subroutine inner_csrsv(tra,a,x,y) + logical, intent(in) :: tra + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(out) :: y(:) + + integer :: i,j,k,m, ir, jc + real(psb_dpk_) :: acc + + if (.not.tra) then + + if (a%is_lower()) then + if (a%is_unit()) then + do i=1, a%get_nrows() + acc = dzero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j)*x(a%ja(j)) + end do + y(i) = x(i) - acc + end do + else if (.not.a%is_unit()) then + do i=1, a%get_nrows() + acc = dzero + do j=a%irp(i), a%irp(i+1)-2 + acc = acc + a%val(j)*x(a%ja(j)) + end do + y(i) = (x(i) - acc)/a%val(a%irp(i+1)-1) + end do + end if + else if (a%is_upper()) then + + if (a%is_unit()) then + do i=a%get_nrows(), 1, -1 + acc = dzero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j)*x(a%ja(j)) + end do + y(i) = x(i) - acc + end do + else if (.not.a%is_unit()) then + do i=a%get_nrows(), 1, -1 + acc = dzero + do j=a%irp(i)+1, a%irp(i+1)-1 + acc = acc + a%val(j)*x(a%ja(j)) + end do + y(i) = (x(i) - acc)/a%val(a%irp(i)) + end do + end if + + end if + + else if (tra) then + + do i=1, a%get_nrows() + y(i) = x(i) + end do + + if (a%is_lower()) then + if (a%is_unit()) then + do i=a%get_nrows(), 1, -1 + acc = y(i) + do j=a%irp(i), a%irp(i+1)-1 + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + end do + end do + else if (.not.a%is_unit()) then + do i=a%get_nrows(), 1, -1 + y(i) = y(i)/a%val(a%irp(i+1)-1) + acc = y(i) + do j=a%irp(i), a%irp(i+1)-2 + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + end do + end do + end if + else if (a%is_upper()) then + + if (a%is_unit()) then + do i=1, a%get_nrows() + acc = y(i) + do j=a%irp(i), a%irp(i+1)-1 + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + end do + end do + else if (.not.a%is_unit()) then + do i=1, a%get_nrows() + y(i) = y(i)/a%val(a%irp(i)) + acc = y(i) + do j=a%irp(i)+1, a%irp(i+1)-1 + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + end do + end do + end if + + end if + end if + end subroutine inner_csrsv + +end subroutine d_csr_cssv_impl + + + +subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_csr_cssm_impl + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_base_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + tra = ((trans_=='T').or.(trans_=='t')) + m = a%get_nrows() + nc = min(size(x,2) , size(y,2)) + + if (.not. (a%is_triangle())) then + write(0,*) 'Called SM on a non-triangular mat!' + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + if (alpha == dzero) then + if (beta == dzero) then + do i = 1, m + y(i,:) = dzero + enddo + else + do i = 1, m + y(i,:) = beta*y(i,:) + end do + endif + return + end if + + if (beta == dzero) then + call inner_csrsm(tra,a,x,y,info) + do i = 1, m + y(i,:) = alpha*y(i,:) + end do + else + allocate(tmp(m,nc), stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + tmp(1:m,:) = x(1:m,:) + call inner_csrsm(tra,a,tmp,y,info) + do i = 1, m + y(i,:) = alpha*tmp(i,:) + beta*y(i,:) + end do + end if + + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='inner_csrsm') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + +contains + + subroutine inner_csrsm(tra,a,x,y,info) + logical, intent(in) :: tra + class(psbn_d_csr_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: x(:,:) + real(psb_dpk_), intent(out) :: y(:,:) + integer, intent(out) :: info + integer :: i,j,k,m, ir, jc + real(psb_dpk_), allocatable :: acc(:) + + allocate(acc(size(x,2)), stat=info) + if(info /= 0) then + info=4010 + return + end if + + + if (.not.tra) then + + if (a%is_lower()) then + if (a%is_unit()) then + do i=1, a%get_nrows() + acc = dzero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j)*x(a%ja(j),:) + end do + y(i,:) = x(i,:) - acc + end do + else if (.not.a%is_unit()) then + do i=1, a%get_nrows() + acc = dzero + do j=a%irp(i), a%irp(i+1)-2 + acc = acc + a%val(j)*x(a%ja(j),:) + end do + y(i,:) = (x(i,:) - acc)/a%val(a%irp(i+1)-1) + end do + end if + else if (a%is_upper()) then + + if (a%is_unit()) then + do i=a%get_nrows(), 1, -1 + acc = dzero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j)*x(a%ja(j),:) + end do + y(i,:) = x(i,:) - acc + end do + else if (.not.a%is_unit()) then + do i=a%get_nrows(), 1, -1 + acc = dzero + do j=a%irp(i)+1, a%irp(i+1)-1 + acc = acc + a%val(j)*x(a%ja(j),:) + end do + y(i,:) = (x(i,:) - acc)/a%val(a%irp(i)) + end do + end if + + end if + + else if (tra) then + + do i=1, a%get_nrows() + y(i,:) = x(i,:) + end do + + if (a%is_lower()) then + if (a%is_unit()) then + do i=a%get_nrows(), 1, -1 + acc = y(i,:) + do j=a%irp(i), a%irp(i+1)-1 + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + end do + end do + else if (.not.a%is_unit()) then + do i=a%get_nrows(), 1, -1 + y(i,:) = y(i,:)/a%val(a%irp(i+1)-1) + acc = y(i,:) + do j=a%irp(i), a%irp(i+1)-2 + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + end do + end do + end if + else if (a%is_upper()) then + + if (a%is_unit()) then + do i=1, a%get_nrows() + acc = y(i,:) + do j=a%irp(i), a%irp(i+1)-1 + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + end do + end do + else if (.not.a%is_unit()) then + do i=1, a%get_nrows() + y(i,:) = y(i,:)/a%val(a%irp(i)) + acc = y(i,:) + do j=a%irp(i)+1, a%irp(i+1)-1 + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + end do + end do + end if + + end if + end if + end subroutine inner_csrsm + +end subroutine d_csr_cssm_impl + + + +subroutine d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psb_error_mod + use psb_realloc_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csins_impl + + class(psbn_d_coo_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='d_csr_csins' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + + nza = a%get_nzeros() + + if (a%is_bld()) then + ! Build phase should only ever be in COO + info = 1121 + + else if (a%is_upd()) then + if (a%is_sorted()) then + + +!!$#ifdef FIXED_NAG_SEGV +!!$ call d_csr_srch_upd(nz,ia,ja,val,a,& +!!$ & imin,imax,jmin,jmax,info,gtl) +!!$#else + call d_csr_srch_upd(nz,ia,ja,val,& + & a%irp,a%ja,a%val,& + & a%get_dupl(),a%get_nzeros(),a%get_nrows(),& + & info,gtl) +!!$#endif + + else + info = 1121 + end if + + else + ! State is wrong. + info = 1121 + end if + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + +contains + + +!!$#ifdef FIXED_NAG_SEGV +!!$ subroutine d_csr_srch_upd(nz,ia,ja,val,a,& +!!$ & imin,imax,jmin,jmax,info,gtl) +!!$ +!!$ use psb_const_mod +!!$ use psb_realloc_mod +!!$ use psb_string_mod +!!$ use psb_serial_mod +!!$ implicit none +!!$ +!!$ class(psbn_d_csr_sparse_mat), intent(inout) :: a +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: gtl(:) +!!$ integer :: i,ir,ic, ilr, ilc, ip, & +!!$ & i1,i2,nc,nnz,dupl,ng +!!$ integer :: debug_level, debug_unit +!!$ character(len=20) :: name='d_csr_srch_upd' +!!$ +!!$ info = 0 +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ dupl = a%get_dupl() +!!$ +!!$ if (.not.a%is_sorted()) then +!!$ info = -4 +!!$ return +!!$ end if +!!$ +!!$ ilr = -1 +!!$ ilc = -1 +!!$ nnz = a%get_nzeros() +!!$ +!!$ +!!$ if (present(gtl)) then +!!$ ng = size(gtl) +!!$ +!!$ select case(dupl) +!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ ic = gtl(ic) +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ a%val(i1+ip-1) = val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ endif +!!$ end if +!!$ end do +!!$ case(psbn_dupl_add_) +!!$ ! Add +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ ic = gtl(ic) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case default +!!$ info = -3 +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Duplicate handling: ',dupl +!!$ end select +!!$ +!!$ else +!!$ +!!$ select case(dupl) +!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ a%val(i1+ip-1) = val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case(psbn_dupl_add_) +!!$ ! Add +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case default +!!$ info = -3 +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Duplicate handling: ',dupl +!!$ end select +!!$ +!!$ end if +!!$ +!!$ end subroutine d_csr_srch_upd +!!$ +!!$#else + subroutine d_csr_srch_upd(nz,ia,ja,val,& + & airp,aja,aval,dupl,nza,nra,& + & info,gtl) + + use psb_error_mod + use psb_sort_mod + implicit none + + integer, intent(inout) :: airp(:),aja(:) + real(psb_dpk_), intent(inout) :: aval(:) + integer, intent(in) :: nz, dupl,nza, nra + integer, intent(in) :: ia(:),ja(:) + real(psb_dpk_), intent(in) :: val(:) + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + integer :: debug_level, debug_unit + character(len=20) :: name='d_csr_srch_upd' + integer :: i,ir,ic, ilr, ilc, ip, & + & i1,i2,nc,lb,ub,m, ng + + info = 0 + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + if (present(gtl)) then + ng = size(gtl) + + select case(dupl) + case(psbn_dupl_ovwrt_,psbn_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= a%m)) then + i1 = airp(ir) + i2 = airp(ir+1) + nc=i2-i1 + + ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) + if (ip>0) then + aval(i1+ip-1) = val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',aja(i1:i2-1) + info = i + return + end if + + else + + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end if + end do + + case(psbn_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= a%m)) then + i1 = airp(ir) + i2 = airp(ir+1) + nc = i2-i1 + ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) + if (ip>0) then + aval(i1+ip-1) = aval(i1+ip-1) + val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',aja(i1:i2-1) + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psbn_dupl_ovwrt_,psbn_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + + if ((ir > 0).and.(ir <= a%m)) then + + i1 = airp(ir) + i2 = airp(ir+1) + nc=i2-i1 + + ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) + if (ip>0) then + aval(i1+ip-1) = val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',aja(i1:i2-1) + info = i + return + end if + + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + + end do + + case(psbn_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= a%m)) then + i1 = airp(ir) + i2 = airp(ir+1) + nc = i2-i1 + ip = psb_ibsrch(ic,nc,aja(i1:i2-1)) + if (ip>0) then + aval(i1+ip-1) = aval(i1+ip-1) + val(i) + else + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + end if + end subroutine d_csr_srch_upd +!!$#endif +end subroutine d_csr_csins_impl