diff --git a/base/modules/psb_c_base_mat_mod.f90 b/base/modules/psb_c_base_mat_mod.f90 index 87e20826..b510621f 100644 --- a/base/modules/psb_c_base_mat_mod.f90 +++ b/base/modules/psb_c_base_mat_mod.f90 @@ -73,6 +73,7 @@ module psb_c_base_mat_mod procedure, pass(a) :: mv_from_fmt => psb_c_base_mv_from_fmt procedure, pass(a) :: mold => psb_c_base_mold procedure, pass(a) :: clone => psb_c_base_clone + procedure, pass(a) :: add_unit_diag => psb_c_base_add_unit_diag ! ! Transpose methods: defined here but not implemented. @@ -430,6 +431,23 @@ module psb_c_base_mat_mod end subroutine psb_c_base_clone end interface + + ! + ! + !> Function add_unit_diag: + !! \memberof psb_c_base_add_unit_diag + !! \brief Given a matrix for which is_unit() is true, explicitly + !! store the unit diagonal and set is_unit() to false. + !! This is needed e.g. when scaling + ! + interface + subroutine psb_c_base_add_unit_diag(a) + import :: psb_c_base_sparse_mat + implicit none + class(psb_c_base_sparse_mat), intent(inout) :: a + end subroutine psb_c_base_add_unit_diag + end interface + ! !> Function cp_to_coo: diff --git a/base/modules/psb_d_base_mat_mod.f90 b/base/modules/psb_d_base_mat_mod.f90 index a54617c2..e190d0f9 100644 --- a/base/modules/psb_d_base_mat_mod.f90 +++ b/base/modules/psb_d_base_mat_mod.f90 @@ -73,6 +73,7 @@ module psb_d_base_mat_mod procedure, pass(a) :: mv_from_fmt => psb_d_base_mv_from_fmt procedure, pass(a) :: mold => psb_d_base_mold procedure, pass(a) :: clone => psb_d_base_clone + procedure, pass(a) :: add_unit_diag => psb_d_base_add_unit_diag ! ! Transpose methods: defined here but not implemented. @@ -430,6 +431,23 @@ module psb_d_base_mat_mod end subroutine psb_d_base_clone end interface + + ! + ! + !> Function add_unit_diag: + !! \memberof psb_d_base_add_unit_diag + !! \brief Given a matrix for which is_unit() is true, explicitly + !! store the unit diagonal and set is_unit() to false. + !! This is needed e.g. when scaling + ! + interface + subroutine psb_d_base_add_unit_diag(a) + import :: psb_d_base_sparse_mat + implicit none + class(psb_d_base_sparse_mat), intent(inout) :: a + end subroutine psb_d_base_add_unit_diag + end interface + ! !> Function cp_to_coo: diff --git a/base/modules/psb_s_base_mat_mod.f90 b/base/modules/psb_s_base_mat_mod.f90 index eac93184..724f0eea 100644 --- a/base/modules/psb_s_base_mat_mod.f90 +++ b/base/modules/psb_s_base_mat_mod.f90 @@ -73,6 +73,7 @@ module psb_s_base_mat_mod procedure, pass(a) :: mv_from_fmt => psb_s_base_mv_from_fmt procedure, pass(a) :: mold => psb_s_base_mold procedure, pass(a) :: clone => psb_s_base_clone + procedure, pass(a) :: add_unit_diag => psb_s_base_add_unit_diag ! ! Transpose methods: defined here but not implemented. @@ -430,6 +431,23 @@ module psb_s_base_mat_mod end subroutine psb_s_base_clone end interface + + ! + ! + !> Function add_unit_diag: + !! \memberof psb_s_base_add_unit_diag + !! \brief Given a matrix for which is_unit() is true, explicitly + !! store the unit diagonal and set is_unit() to false. + !! This is needed e.g. when scaling + ! + interface + subroutine psb_s_base_add_unit_diag(a) + import :: psb_s_base_sparse_mat + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + end subroutine psb_s_base_add_unit_diag + end interface + ! !> Function cp_to_coo: diff --git a/base/modules/psb_z_base_mat_mod.f90 b/base/modules/psb_z_base_mat_mod.f90 index e22ac14d..0b78a0a5 100644 --- a/base/modules/psb_z_base_mat_mod.f90 +++ b/base/modules/psb_z_base_mat_mod.f90 @@ -73,6 +73,7 @@ module psb_z_base_mat_mod procedure, pass(a) :: mv_from_fmt => psb_z_base_mv_from_fmt procedure, pass(a) :: mold => psb_z_base_mold procedure, pass(a) :: clone => psb_z_base_clone + procedure, pass(a) :: add_unit_diag => psb_z_base_add_unit_diag ! ! Transpose methods: defined here but not implemented. @@ -430,6 +431,23 @@ module psb_z_base_mat_mod end subroutine psb_z_base_clone end interface + + ! + ! + !> Function add_unit_diag: + !! \memberof psb_z_base_add_unit_diag + !! \brief Given a matrix for which is_unit() is true, explicitly + !! store the unit diagonal and set is_unit() to false. + !! This is needed e.g. when scaling + ! + interface + subroutine psb_z_base_add_unit_diag(a) + import :: psb_z_base_sparse_mat + implicit none + class(psb_z_base_sparse_mat), intent(inout) :: a + end subroutine psb_z_base_add_unit_diag + end interface + ! !> Function cp_to_coo: diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90 index 81e479fd..683fa266 100644 --- a/base/serial/impl/psb_c_base_mat_impl.F90 +++ b/base/serial/impl/psb_c_base_mat_impl.F90 @@ -580,6 +580,37 @@ subroutine psb_c_base_clone(a,b,info) end subroutine psb_c_base_clone +subroutine psb_c_base_add_unit_diag(a) + use psb_c_base_mat_mod, psb_protect_name => psb_c_base_add_unit_diag + use psb_error_mod + implicit none + class(psb_c_base_sparse_mat), intent(inout) :: a + type(psb_c_coo_sparse_mat) :: tmp + + integer(psb_ipk_) :: i, j, m, n, nz, mnm, info + + if (a%is_unit()) then + call a%mv_to_coo(tmp,info) + if (info /= 0) return + m = tmp%get_nrows() + n = tmp%get_ncols() + mnm = min(m,n) + nz = tmp%get_nzeros() + call tmp%reallocate(nz+mnm) + do i=1, mnm + tmp%val(nz+i) = cone + tmp%ia(nz+i) = i + tmp%ja(nz+i) = i + end do + call tmp%set_nzeros(nz+mnm) + call tmp%set_unit(.false.) + call tmp%fix(info) + if (info /= 0) & + & call a%mv_from_coo(tmp,info) + end if + +end subroutine psb_c_base_add_unit_diag + subroutine psb_c_base_mold(a,b,info) use psb_c_base_mat_mod, psb_protect_name => psb_c_base_mold use psb_error_mod diff --git a/base/serial/impl/psb_c_coo_impl.f90 b/base/serial/impl/psb_c_coo_impl.f90 index b86184b0..0c4a7892 100644 --- a/base/serial/impl/psb_c_coo_impl.f90 +++ b/base/serial/impl/psb_c_coo_impl.f90 @@ -54,11 +54,11 @@ subroutine psb_c_coo_get_diag(a,d,info) call psb_errpush(info,name,i_err=ierr) goto 9999 end if - d(:) = czero - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = cone else + d(1:mnm) = czero do i=1,a%get_nzeros() j=a%ia(i) if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then @@ -101,6 +101,10 @@ subroutine psb_c_coo_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -167,6 +171,9 @@ subroutine psb_c_coo_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d @@ -1313,7 +1320,7 @@ subroutine psb_c_coo_csmv(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == czero) then do i = 1, min(m,n) y(i) = alpha*x(i) @@ -1524,7 +1531,7 @@ subroutine psb_c_coo_csmm(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == czero) then do i = 1, min(m,n) y(i,1:nc) = alpha*x(i,1:nc) @@ -1661,12 +1668,17 @@ function psb_c_coo_maxval(a) result(res) logical, parameter :: debug=.false. - res = szero + if (a%is_unit()) then + res = sone + else + res = szero + end if nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) res = maxval(abs(a%val(1:nnz))) end if + end function psb_c_coo_maxval function psb_c_coo_csnmi(a) result(res) @@ -1696,7 +1708,11 @@ function psb_c_coo_csnmi(a) result(res) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) j = j+1 enddo - acc = szero + if (a%is_unit()) then + acc = sone + else + acc = szero + end if do k=i, j-1 acc = acc + abs(a%val(k)) end do @@ -1707,7 +1723,11 @@ function psb_c_coo_csnmi(a) result(res) m = a%get_nrows() allocate(vt(m),stat=info) if (info /= 0) return - vt(:) = szero + if (a%is_unit()) then + vt = sone + else + vt = szero + end if do j=1, nnz i = a%ia(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1740,10 +1760,14 @@ function psb_c_coo_csnm1(a) result(res) res = szero nnz = a%get_nzeros() - n = a%get_ncols() + n = a%get_ncols() allocate(vt(n),stat=info) if (info /= 0) return - vt(:) = szero + if (a%is_unit()) then + vt = sone + else + vt = szero + end if do j=1, nnz i = a%ja(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1781,18 +1805,17 @@ subroutine psb_c_coo_rowsum(d,a) goto 9999 end if - d = czero + if (a%is_unit()) then + d = cone + else + d = czero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + a%val(j) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + cone - end do - end if return call psb_erractionrestore(err_act) @@ -1835,19 +1858,17 @@ subroutine psb_c_coo_arwsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + sone - end do - end if - return call psb_erractionrestore(err_act) return @@ -1889,18 +1910,17 @@ subroutine psb_c_coo_colsum(d,a) goto 9999 end if - d = czero + if (a%is_unit()) then + d = cone + else + d = czero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + a%val(j) end do - - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + cone - end do - end if return call psb_erractionrestore(err_act) @@ -1943,19 +1963,18 @@ subroutine psb_c_coo_aclsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + sone - end do - end if - return call psb_erractionrestore(err_act) return @@ -1996,15 +2015,15 @@ subroutine psb_c_coo_csgetptn(imin,imax,a,nz,ia,ja,info,& use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_csgetptn implicit none - class(psb_c_coo_sparse_mat), intent(in) :: a - integer(psb_ipk_), intent(in) :: imin,imax - integer(psb_ipk_), intent(out) :: nz - integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) - integer(psb_ipk_),intent(out) :: info - logical, intent(in), optional :: append - integer(psb_ipk_), intent(in), optional :: iren(:) - integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin - logical, intent(in), optional :: rscale,cscale + class(psb_c_coo_sparse_mat), intent(in) :: a + integer(psb_ipk_), intent(in) :: imin,imax + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) + integer(psb_ipk_),intent(out) :: info + logical, intent(in), optional :: append + integer(psb_ipk_), intent(in), optional :: iren(:) + integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i diff --git a/base/serial/impl/psb_c_csc_impl.f90 b/base/serial/impl/psb_c_csc_impl.f90 index 6799fdf2..4fa47574 100644 --- a/base/serial/impl/psb_c_csc_impl.f90 +++ b/base/serial/impl/psb_c_csc_impl.f90 @@ -304,7 +304,7 @@ subroutine psb_c_csc_csmv(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -590,7 +590,7 @@ subroutine psb_c_csc_csmm(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i,:) = y(i,:) + alpha*x(i,:) end do @@ -1081,7 +1081,12 @@ function psb_c_csc_maxval(a) result(res) logical, parameter :: debug=.false. - res = szero + if (a%is_unit()) then + res = sone + else + res = szero + end if + nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) @@ -1112,7 +1117,11 @@ function psb_c_csc_csnmi(a) result(res) if (info /= psb_success_) then return end if - acc(:) = szero + if (a%is_unit()) then + acc = sone + else + acc = szero + end if do i=1, nc do j=a%icp(i),a%icp(i+1)-1 acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) @@ -1149,7 +1158,11 @@ function psb_c_csc_csnm1(a) result(res) m = a%get_nrows() n = a%get_ncols() do j=1, n - acc = szero + if (a%is_unit()) then + acc = sone + else + acc = szero + end if do k=a%icp(j),a%icp(j+1)-1 acc = acc + abs(a%val(k)) end do @@ -1187,19 +1200,17 @@ subroutine psb_c_csc_colsum(d,a) end if do i = 1, a%get_ncols() - d(i) = czero + if (a%is_unit()) then + d(i) = cone + else + d(i) = czero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + (a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_ncols() - d(i) = d(i) + cone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1242,13 +1253,18 @@ subroutine psb_c_csc_aclsum(d,a) do i = 1, a%get_ncols() - d(i) = szero + if (a%is_unit()) then + d(i) = sone + else + d(i) = szero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + abs(a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, a%get_ncols() d(i) = d(i) + sone end do @@ -1295,7 +1311,11 @@ subroutine psb_c_csc_rowsum(d,a) goto 9999 end if - d = czero + if (a%is_unit()) then + d = cone + else + d = czero + end if do i=1, m do j=a%icp(i),a%icp(i+1)-1 @@ -1304,13 +1324,6 @@ subroutine psb_c_csc_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + cone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1352,7 +1365,12 @@ subroutine psb_c_csc_arwsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if + do i=1, m do j=a%icp(i),a%icp(i+1)-1 k = a%ia(j) @@ -1360,13 +1378,6 @@ subroutine psb_c_csc_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + sone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1408,7 +1419,7 @@ subroutine psb_c_csc_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = cone else do i=1, mnm @@ -1445,13 +1456,14 @@ subroutine psb_c_csc_scal(d,a,info,side) use psb_string_mod implicit none class(psb_c_csc_sparse_mat), intent(inout) :: a - complex(psb_spk_), intent(in) :: d(:) - integer(psb_ipk_), intent(out) :: info - character, intent(in), optional :: side + complex(psb_spk_), intent(in) :: d(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side integer(psb_ipk_) :: err_act,mnm, i, j, n + type(psb_c_coo_sparse_mat) :: tmp integer(psb_ipk_) :: ierr(5) - character(len=20) :: name='scal' + character(len=20) :: name='scal' character :: side_ logical :: left logical, parameter :: debug=.false. @@ -1464,6 +1476,10 @@ subroutine psb_c_csc_scal(d,a,info,side) side_ = psb_toupper(side) end if + if (a%is_unit()) then + call a%add_unit_diag() + end if + left = (side_ == 'L') if (left) then @@ -1524,6 +1540,9 @@ subroutine psb_c_csc_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_c_csr_impl.f90 b/base/serial/impl/psb_c_csr_impl.f90 index be6f0441..b541e4f0 100644 --- a/base/serial/impl/psb_c_csr_impl.f90 +++ b/base/serial/impl/psb_c_csr_impl.f90 @@ -380,7 +380,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -733,7 +733,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) end do @@ -1398,7 +1398,7 @@ subroutine psb_c_csr_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + cone end do @@ -1453,7 +1453,7 @@ subroutine psb_c_csr_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + sone end do @@ -1509,7 +1509,7 @@ subroutine psb_c_csr_colsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + cone end do @@ -1566,7 +1566,7 @@ subroutine psb_c_csr_aclsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + sone end do @@ -1613,7 +1613,7 @@ subroutine psb_c_csr_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = cone else do i=1, mnm @@ -1665,6 +1665,10 @@ subroutine psb_c_csr_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -1734,6 +1738,9 @@ subroutine psb_c_csr_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90 index 5f67c547..17192097 100644 --- a/base/serial/impl/psb_d_base_mat_impl.F90 +++ b/base/serial/impl/psb_d_base_mat_impl.F90 @@ -580,6 +580,37 @@ subroutine psb_d_base_clone(a,b,info) end subroutine psb_d_base_clone +subroutine psb_d_base_add_unit_diag(a) + use psb_d_base_mat_mod, psb_protect_name => psb_d_base_add_unit_diag + use psb_error_mod + implicit none + class(psb_d_base_sparse_mat), intent(inout) :: a + type(psb_d_coo_sparse_mat) :: tmp + + integer(psb_ipk_) :: i, j, m, n, nz, mnm, info + + if (a%is_unit()) then + call a%mv_to_coo(tmp,info) + if (info /= 0) return + m = tmp%get_nrows() + n = tmp%get_ncols() + mnm = min(m,n) + nz = tmp%get_nzeros() + call tmp%reallocate(nz+mnm) + do i=1, mnm + tmp%val(nz+i) = done + tmp%ia(nz+i) = i + tmp%ja(nz+i) = i + end do + call tmp%set_nzeros(nz+mnm) + call tmp%set_unit(.false.) + call tmp%fix(info) + if (info /= 0) & + & call a%mv_from_coo(tmp,info) + end if + +end subroutine psb_d_base_add_unit_diag + subroutine psb_d_base_mold(a,b,info) use psb_d_base_mat_mod, psb_protect_name => psb_d_base_mold use psb_error_mod diff --git a/base/serial/impl/psb_d_coo_impl.f90 b/base/serial/impl/psb_d_coo_impl.f90 index 1238ad93..db4f426d 100644 --- a/base/serial/impl/psb_d_coo_impl.f90 +++ b/base/serial/impl/psb_d_coo_impl.f90 @@ -54,11 +54,11 @@ subroutine psb_d_coo_get_diag(a,d,info) call psb_errpush(info,name,i_err=ierr) goto 9999 end if - d(:) = dzero - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = done else + d(1:mnm) = dzero do i=1,a%get_nzeros() j=a%ia(i) if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then @@ -101,6 +101,10 @@ subroutine psb_d_coo_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -167,6 +171,9 @@ subroutine psb_d_coo_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d @@ -1313,7 +1320,7 @@ subroutine psb_d_coo_csmv(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == dzero) then do i = 1, min(m,n) y(i) = alpha*x(i) @@ -1524,7 +1531,7 @@ subroutine psb_d_coo_csmm(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == dzero) then do i = 1, min(m,n) y(i,1:nc) = alpha*x(i,1:nc) @@ -1661,12 +1668,17 @@ function psb_d_coo_maxval(a) result(res) logical, parameter :: debug=.false. - res = dzero + if (a%is_unit()) then + res = done + else + res = dzero + end if nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) res = maxval(abs(a%val(1:nnz))) end if + end function psb_d_coo_maxval function psb_d_coo_csnmi(a) result(res) @@ -1696,7 +1708,11 @@ function psb_d_coo_csnmi(a) result(res) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) j = j+1 enddo - acc = dzero + if (a%is_unit()) then + acc = done + else + acc = dzero + end if do k=i, j-1 acc = acc + abs(a%val(k)) end do @@ -1707,7 +1723,11 @@ function psb_d_coo_csnmi(a) result(res) m = a%get_nrows() allocate(vt(m),stat=info) if (info /= 0) return - vt(:) = dzero + if (a%is_unit()) then + vt = done + else + vt = dzero + end if do j=1, nnz i = a%ia(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1740,10 +1760,14 @@ function psb_d_coo_csnm1(a) result(res) res = dzero nnz = a%get_nzeros() - n = a%get_ncols() + n = a%get_ncols() allocate(vt(n),stat=info) if (info /= 0) return - vt(:) = dzero + if (a%is_unit()) then + vt = done + else + vt = dzero + end if do j=1, nnz i = a%ja(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1781,18 +1805,17 @@ subroutine psb_d_coo_rowsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + a%val(j) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + done - end do - end if return call psb_erractionrestore(err_act) @@ -1835,19 +1858,17 @@ subroutine psb_d_coo_arwsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + done - end do - end if - return call psb_erractionrestore(err_act) return @@ -1889,18 +1910,17 @@ subroutine psb_d_coo_colsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + a%val(j) end do - - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + done - end do - end if return call psb_erractionrestore(err_act) @@ -1943,19 +1963,18 @@ subroutine psb_d_coo_aclsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + done - end do - end if - return call psb_erractionrestore(err_act) return @@ -1996,15 +2015,15 @@ subroutine psb_d_coo_csgetptn(imin,imax,a,nz,ia,ja,info,& use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_csgetptn implicit none - class(psb_d_coo_sparse_mat), intent(in) :: a - integer(psb_ipk_), intent(in) :: imin,imax - integer(psb_ipk_), intent(out) :: nz - integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) - integer(psb_ipk_),intent(out) :: info - logical, intent(in), optional :: append - integer(psb_ipk_), intent(in), optional :: iren(:) - integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin - logical, intent(in), optional :: rscale,cscale + class(psb_d_coo_sparse_mat), intent(in) :: a + integer(psb_ipk_), intent(in) :: imin,imax + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) + integer(psb_ipk_),intent(out) :: info + logical, intent(in), optional :: append + integer(psb_ipk_), intent(in), optional :: iren(:) + integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i diff --git a/base/serial/impl/psb_d_csc_impl.f90 b/base/serial/impl/psb_d_csc_impl.f90 index febab98c..ffdb2afd 100644 --- a/base/serial/impl/psb_d_csc_impl.f90 +++ b/base/serial/impl/psb_d_csc_impl.f90 @@ -304,7 +304,7 @@ subroutine psb_d_csc_csmv(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -590,7 +590,7 @@ subroutine psb_d_csc_csmm(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i,:) = y(i,:) + alpha*x(i,:) end do @@ -1081,7 +1081,12 @@ function psb_d_csc_maxval(a) result(res) logical, parameter :: debug=.false. - res = dzero + if (a%is_unit()) then + res = done + else + res = dzero + end if + nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) @@ -1112,7 +1117,11 @@ function psb_d_csc_csnmi(a) result(res) if (info /= psb_success_) then return end if - acc(:) = dzero + if (a%is_unit()) then + acc = done + else + acc = dzero + end if do i=1, nc do j=a%icp(i),a%icp(i+1)-1 acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) @@ -1149,7 +1158,11 @@ function psb_d_csc_csnm1(a) result(res) m = a%get_nrows() n = a%get_ncols() do j=1, n - acc = dzero + if (a%is_unit()) then + acc = done + else + acc = dzero + end if do k=a%icp(j),a%icp(j+1)-1 acc = acc + abs(a%val(k)) end do @@ -1187,19 +1200,17 @@ subroutine psb_d_csc_colsum(d,a) end if do i = 1, a%get_ncols() - d(i) = dzero + if (a%is_unit()) then + d(i) = done + else + d(i) = dzero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + (a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_ncols() - d(i) = d(i) + done - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1242,13 +1253,18 @@ subroutine psb_d_csc_aclsum(d,a) do i = 1, a%get_ncols() - d(i) = dzero + if (a%is_unit()) then + d(i) = done + else + d(i) = dzero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + abs(a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, a%get_ncols() d(i) = d(i) + done end do @@ -1295,7 +1311,11 @@ subroutine psb_d_csc_rowsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if do i=1, m do j=a%icp(i),a%icp(i+1)-1 @@ -1304,13 +1324,6 @@ subroutine psb_d_csc_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + done - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1352,7 +1365,12 @@ subroutine psb_d_csc_arwsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if + do i=1, m do j=a%icp(i),a%icp(i+1)-1 k = a%ia(j) @@ -1360,13 +1378,6 @@ subroutine psb_d_csc_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + done - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1408,7 +1419,7 @@ subroutine psb_d_csc_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = done else do i=1, mnm @@ -1445,13 +1456,14 @@ subroutine psb_d_csc_scal(d,a,info,side) use psb_string_mod implicit none class(psb_d_csc_sparse_mat), intent(inout) :: a - real(psb_dpk_), intent(in) :: d(:) - integer(psb_ipk_), intent(out) :: info - character, intent(in), optional :: side + real(psb_dpk_), intent(in) :: d(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side integer(psb_ipk_) :: err_act,mnm, i, j, n + type(psb_d_coo_sparse_mat) :: tmp integer(psb_ipk_) :: ierr(5) - character(len=20) :: name='scal' + character(len=20) :: name='scal' character :: side_ logical :: left logical, parameter :: debug=.false. @@ -1464,6 +1476,10 @@ subroutine psb_d_csc_scal(d,a,info,side) side_ = psb_toupper(side) end if + if (a%is_unit()) then + call a%add_unit_diag() + end if + left = (side_ == 'L') if (left) then @@ -1524,6 +1540,9 @@ subroutine psb_d_csc_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_d_csr_impl.f90 b/base/serial/impl/psb_d_csr_impl.f90 index 9ad2db78..007c45cf 100644 --- a/base/serial/impl/psb_d_csr_impl.f90 +++ b/base/serial/impl/psb_d_csr_impl.f90 @@ -380,7 +380,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -733,7 +733,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) end do @@ -1398,7 +1398,7 @@ subroutine psb_d_csr_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + done end do @@ -1453,7 +1453,7 @@ subroutine psb_d_csr_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + done end do @@ -1509,7 +1509,7 @@ subroutine psb_d_csr_colsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + done end do @@ -1566,7 +1566,7 @@ subroutine psb_d_csr_aclsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + done end do @@ -1613,7 +1613,7 @@ subroutine psb_d_csr_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = done else do i=1, mnm @@ -1665,6 +1665,10 @@ subroutine psb_d_csr_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -1734,6 +1738,9 @@ subroutine psb_d_csr_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_s_base_mat_impl.F90 b/base/serial/impl/psb_s_base_mat_impl.F90 index 1483d9b9..b31c1fad 100644 --- a/base/serial/impl/psb_s_base_mat_impl.F90 +++ b/base/serial/impl/psb_s_base_mat_impl.F90 @@ -580,6 +580,37 @@ subroutine psb_s_base_clone(a,b,info) end subroutine psb_s_base_clone +subroutine psb_s_base_add_unit_diag(a) + use psb_s_base_mat_mod, psb_protect_name => psb_s_base_add_unit_diag + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + type(psb_s_coo_sparse_mat) :: tmp + + integer(psb_ipk_) :: i, j, m, n, nz, mnm, info + + if (a%is_unit()) then + call a%mv_to_coo(tmp,info) + if (info /= 0) return + m = tmp%get_nrows() + n = tmp%get_ncols() + mnm = min(m,n) + nz = tmp%get_nzeros() + call tmp%reallocate(nz+mnm) + do i=1, mnm + tmp%val(nz+i) = sone + tmp%ia(nz+i) = i + tmp%ja(nz+i) = i + end do + call tmp%set_nzeros(nz+mnm) + call tmp%set_unit(.false.) + call tmp%fix(info) + if (info /= 0) & + & call a%mv_from_coo(tmp,info) + end if + +end subroutine psb_s_base_add_unit_diag + subroutine psb_s_base_mold(a,b,info) use psb_s_base_mat_mod, psb_protect_name => psb_s_base_mold use psb_error_mod diff --git a/base/serial/impl/psb_s_coo_impl.f90 b/base/serial/impl/psb_s_coo_impl.f90 index 38a8fa76..4a8eaa94 100644 --- a/base/serial/impl/psb_s_coo_impl.f90 +++ b/base/serial/impl/psb_s_coo_impl.f90 @@ -54,11 +54,11 @@ subroutine psb_s_coo_get_diag(a,d,info) call psb_errpush(info,name,i_err=ierr) goto 9999 end if - d(:) = szero - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = sone else + d(1:mnm) = szero do i=1,a%get_nzeros() j=a%ia(i) if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then @@ -101,6 +101,10 @@ subroutine psb_s_coo_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -167,6 +171,9 @@ subroutine psb_s_coo_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d @@ -1313,7 +1320,7 @@ subroutine psb_s_coo_csmv(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == szero) then do i = 1, min(m,n) y(i) = alpha*x(i) @@ -1524,7 +1531,7 @@ subroutine psb_s_coo_csmm(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == szero) then do i = 1, min(m,n) y(i,1:nc) = alpha*x(i,1:nc) @@ -1661,12 +1668,17 @@ function psb_s_coo_maxval(a) result(res) logical, parameter :: debug=.false. - res = szero + if (a%is_unit()) then + res = sone + else + res = szero + end if nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) res = maxval(abs(a%val(1:nnz))) end if + end function psb_s_coo_maxval function psb_s_coo_csnmi(a) result(res) @@ -1696,7 +1708,11 @@ function psb_s_coo_csnmi(a) result(res) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) j = j+1 enddo - acc = szero + if (a%is_unit()) then + acc = sone + else + acc = szero + end if do k=i, j-1 acc = acc + abs(a%val(k)) end do @@ -1707,7 +1723,11 @@ function psb_s_coo_csnmi(a) result(res) m = a%get_nrows() allocate(vt(m),stat=info) if (info /= 0) return - vt(:) = szero + if (a%is_unit()) then + vt = sone + else + vt = szero + end if do j=1, nnz i = a%ia(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1740,10 +1760,14 @@ function psb_s_coo_csnm1(a) result(res) res = szero nnz = a%get_nzeros() - n = a%get_ncols() + n = a%get_ncols() allocate(vt(n),stat=info) if (info /= 0) return - vt(:) = szero + if (a%is_unit()) then + vt = sone + else + vt = szero + end if do j=1, nnz i = a%ja(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1781,18 +1805,17 @@ subroutine psb_s_coo_rowsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + a%val(j) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + sone - end do - end if return call psb_erractionrestore(err_act) @@ -1835,19 +1858,17 @@ subroutine psb_s_coo_arwsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + sone - end do - end if - return call psb_erractionrestore(err_act) return @@ -1889,18 +1910,17 @@ subroutine psb_s_coo_colsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + a%val(j) end do - - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + sone - end do - end if return call psb_erractionrestore(err_act) @@ -1943,19 +1963,18 @@ subroutine psb_s_coo_aclsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + sone - end do - end if - return call psb_erractionrestore(err_act) return @@ -1996,15 +2015,15 @@ subroutine psb_s_coo_csgetptn(imin,imax,a,nz,ia,ja,info,& use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_csgetptn implicit none - class(psb_s_coo_sparse_mat), intent(in) :: a - integer(psb_ipk_), intent(in) :: imin,imax - integer(psb_ipk_), intent(out) :: nz - integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) - integer(psb_ipk_),intent(out) :: info - logical, intent(in), optional :: append - integer(psb_ipk_), intent(in), optional :: iren(:) - integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin - logical, intent(in), optional :: rscale,cscale + class(psb_s_coo_sparse_mat), intent(in) :: a + integer(psb_ipk_), intent(in) :: imin,imax + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) + integer(psb_ipk_),intent(out) :: info + logical, intent(in), optional :: append + integer(psb_ipk_), intent(in), optional :: iren(:) + integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i diff --git a/base/serial/impl/psb_s_csc_impl.f90 b/base/serial/impl/psb_s_csc_impl.f90 index 769670f3..0c5601cb 100644 --- a/base/serial/impl/psb_s_csc_impl.f90 +++ b/base/serial/impl/psb_s_csc_impl.f90 @@ -304,7 +304,7 @@ subroutine psb_s_csc_csmv(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -590,7 +590,7 @@ subroutine psb_s_csc_csmm(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i,:) = y(i,:) + alpha*x(i,:) end do @@ -1081,7 +1081,12 @@ function psb_s_csc_maxval(a) result(res) logical, parameter :: debug=.false. - res = szero + if (a%is_unit()) then + res = sone + else + res = szero + end if + nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) @@ -1112,7 +1117,11 @@ function psb_s_csc_csnmi(a) result(res) if (info /= psb_success_) then return end if - acc(:) = szero + if (a%is_unit()) then + acc = sone + else + acc = szero + end if do i=1, nc do j=a%icp(i),a%icp(i+1)-1 acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) @@ -1149,7 +1158,11 @@ function psb_s_csc_csnm1(a) result(res) m = a%get_nrows() n = a%get_ncols() do j=1, n - acc = szero + if (a%is_unit()) then + acc = sone + else + acc = szero + end if do k=a%icp(j),a%icp(j+1)-1 acc = acc + abs(a%val(k)) end do @@ -1187,19 +1200,17 @@ subroutine psb_s_csc_colsum(d,a) end if do i = 1, a%get_ncols() - d(i) = szero + if (a%is_unit()) then + d(i) = sone + else + d(i) = szero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + (a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_ncols() - d(i) = d(i) + sone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1242,13 +1253,18 @@ subroutine psb_s_csc_aclsum(d,a) do i = 1, a%get_ncols() - d(i) = szero + if (a%is_unit()) then + d(i) = sone + else + d(i) = szero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + abs(a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, a%get_ncols() d(i) = d(i) + sone end do @@ -1295,7 +1311,11 @@ subroutine psb_s_csc_rowsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if do i=1, m do j=a%icp(i),a%icp(i+1)-1 @@ -1304,13 +1324,6 @@ subroutine psb_s_csc_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + sone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1352,7 +1365,12 @@ subroutine psb_s_csc_arwsum(d,a) goto 9999 end if - d = szero + if (a%is_unit()) then + d = sone + else + d = szero + end if + do i=1, m do j=a%icp(i),a%icp(i+1)-1 k = a%ia(j) @@ -1360,13 +1378,6 @@ subroutine psb_s_csc_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + sone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1408,7 +1419,7 @@ subroutine psb_s_csc_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = sone else do i=1, mnm @@ -1445,13 +1456,14 @@ subroutine psb_s_csc_scal(d,a,info,side) use psb_string_mod implicit none class(psb_s_csc_sparse_mat), intent(inout) :: a - real(psb_spk_), intent(in) :: d(:) - integer(psb_ipk_), intent(out) :: info - character, intent(in), optional :: side + real(psb_spk_), intent(in) :: d(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side integer(psb_ipk_) :: err_act,mnm, i, j, n + type(psb_s_coo_sparse_mat) :: tmp integer(psb_ipk_) :: ierr(5) - character(len=20) :: name='scal' + character(len=20) :: name='scal' character :: side_ logical :: left logical, parameter :: debug=.false. @@ -1464,6 +1476,10 @@ subroutine psb_s_csc_scal(d,a,info,side) side_ = psb_toupper(side) end if + if (a%is_unit()) then + call a%add_unit_diag() + end if + left = (side_ == 'L') if (left) then @@ -1524,6 +1540,9 @@ subroutine psb_s_csc_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_s_csr_impl.f90 b/base/serial/impl/psb_s_csr_impl.f90 index 68f0535a..027162bd 100644 --- a/base/serial/impl/psb_s_csr_impl.f90 +++ b/base/serial/impl/psb_s_csr_impl.f90 @@ -380,7 +380,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -733,7 +733,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) end do @@ -1398,7 +1398,7 @@ subroutine psb_s_csr_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + sone end do @@ -1453,7 +1453,7 @@ subroutine psb_s_csr_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + sone end do @@ -1509,7 +1509,7 @@ subroutine psb_s_csr_colsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + sone end do @@ -1566,7 +1566,7 @@ subroutine psb_s_csr_aclsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + sone end do @@ -1613,7 +1613,7 @@ subroutine psb_s_csr_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = sone else do i=1, mnm @@ -1665,6 +1665,10 @@ subroutine psb_s_csr_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -1734,6 +1738,9 @@ subroutine psb_s_csr_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90 index cb615459..58206e22 100644 --- a/base/serial/impl/psb_z_base_mat_impl.F90 +++ b/base/serial/impl/psb_z_base_mat_impl.F90 @@ -580,6 +580,37 @@ subroutine psb_z_base_clone(a,b,info) end subroutine psb_z_base_clone +subroutine psb_z_base_add_unit_diag(a) + use psb_z_base_mat_mod, psb_protect_name => psb_z_base_add_unit_diag + use psb_error_mod + implicit none + class(psb_z_base_sparse_mat), intent(inout) :: a + type(psb_z_coo_sparse_mat) :: tmp + + integer(psb_ipk_) :: i, j, m, n, nz, mnm, info + + if (a%is_unit()) then + call a%mv_to_coo(tmp,info) + if (info /= 0) return + m = tmp%get_nrows() + n = tmp%get_ncols() + mnm = min(m,n) + nz = tmp%get_nzeros() + call tmp%reallocate(nz+mnm) + do i=1, mnm + tmp%val(nz+i) = zone + tmp%ia(nz+i) = i + tmp%ja(nz+i) = i + end do + call tmp%set_nzeros(nz+mnm) + call tmp%set_unit(.false.) + call tmp%fix(info) + if (info /= 0) & + & call a%mv_from_coo(tmp,info) + end if + +end subroutine psb_z_base_add_unit_diag + subroutine psb_z_base_mold(a,b,info) use psb_z_base_mat_mod, psb_protect_name => psb_z_base_mold use psb_error_mod diff --git a/base/serial/impl/psb_z_coo_impl.f90 b/base/serial/impl/psb_z_coo_impl.f90 index 7cbc1739..a1b5c081 100644 --- a/base/serial/impl/psb_z_coo_impl.f90 +++ b/base/serial/impl/psb_z_coo_impl.f90 @@ -54,11 +54,11 @@ subroutine psb_z_coo_get_diag(a,d,info) call psb_errpush(info,name,i_err=ierr) goto 9999 end if - d(:) = zzero - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = zone else + d(1:mnm) = zzero do i=1,a%get_nzeros() j=a%ia(i) if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then @@ -101,6 +101,10 @@ subroutine psb_z_coo_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -167,6 +171,9 @@ subroutine psb_z_coo_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d @@ -1313,7 +1320,7 @@ subroutine psb_z_coo_csmv(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == zzero) then do i = 1, min(m,n) y(i) = alpha*x(i) @@ -1524,7 +1531,7 @@ subroutine psb_z_coo_csmm(alpha,a,x,beta,y,info,trans) endif return else - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then if (beta == zzero) then do i = 1, min(m,n) y(i,1:nc) = alpha*x(i,1:nc) @@ -1661,12 +1668,17 @@ function psb_z_coo_maxval(a) result(res) logical, parameter :: debug=.false. - res = dzero + if (a%is_unit()) then + res = done + else + res = dzero + end if nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) res = maxval(abs(a%val(1:nnz))) end if + end function psb_z_coo_maxval function psb_z_coo_csnmi(a) result(res) @@ -1696,7 +1708,11 @@ function psb_z_coo_csnmi(a) result(res) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) j = j+1 enddo - acc = dzero + if (a%is_unit()) then + acc = done + else + acc = dzero + end if do k=i, j-1 acc = acc + abs(a%val(k)) end do @@ -1707,7 +1723,11 @@ function psb_z_coo_csnmi(a) result(res) m = a%get_nrows() allocate(vt(m),stat=info) if (info /= 0) return - vt(:) = dzero + if (a%is_unit()) then + vt = done + else + vt = dzero + end if do j=1, nnz i = a%ia(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1740,10 +1760,14 @@ function psb_z_coo_csnm1(a) result(res) res = dzero nnz = a%get_nzeros() - n = a%get_ncols() + n = a%get_ncols() allocate(vt(n),stat=info) if (info /= 0) return - vt(:) = dzero + if (a%is_unit()) then + vt = done + else + vt = dzero + end if do j=1, nnz i = a%ja(j) vt(i) = vt(i) + abs(a%val(j)) @@ -1781,18 +1805,17 @@ subroutine psb_z_coo_rowsum(d,a) goto 9999 end if - d = zzero + if (a%is_unit()) then + d = zone + else + d = zzero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + a%val(j) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + zone - end do - end if return call psb_erractionrestore(err_act) @@ -1835,19 +1858,17 @@ subroutine psb_z_coo_arwsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if nnz = a%get_nzeros() do j=1, nnz i = a%ia(j) d(i) = d(i) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, m - d(i) = d(i) + done - end do - end if - return call psb_erractionrestore(err_act) return @@ -1889,18 +1910,17 @@ subroutine psb_z_coo_colsum(d,a) goto 9999 end if - d = zzero + if (a%is_unit()) then + d = zone + else + d = zzero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + a%val(j) end do - - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + zone - end do - end if return call psb_erractionrestore(err_act) @@ -1943,19 +1963,18 @@ subroutine psb_z_coo_aclsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if + nnz = a%get_nzeros() do j=1, nnz k = a%ja(j) d(k) = d(k) + abs(a%val(j)) end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, n - d(i) = d(i) + done - end do - end if - return call psb_erractionrestore(err_act) return @@ -1996,15 +2015,15 @@ subroutine psb_z_coo_csgetptn(imin,imax,a,nz,ia,ja,info,& use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_csgetptn implicit none - class(psb_z_coo_sparse_mat), intent(in) :: a - integer(psb_ipk_), intent(in) :: imin,imax - integer(psb_ipk_), intent(out) :: nz - integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) - integer(psb_ipk_),intent(out) :: info - logical, intent(in), optional :: append - integer(psb_ipk_), intent(in), optional :: iren(:) - integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin - logical, intent(in), optional :: rscale,cscale + class(psb_z_coo_sparse_mat), intent(in) :: a + integer(psb_ipk_), intent(in) :: imin,imax + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) + integer(psb_ipk_),intent(out) :: info + logical, intent(in), optional :: append + integer(psb_ipk_), intent(in), optional :: iren(:) + integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i diff --git a/base/serial/impl/psb_z_csc_impl.f90 b/base/serial/impl/psb_z_csc_impl.f90 index dff25d0c..e5e7e90a 100644 --- a/base/serial/impl/psb_z_csc_impl.f90 +++ b/base/serial/impl/psb_z_csc_impl.f90 @@ -304,7 +304,7 @@ subroutine psb_z_csc_csmv(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -590,7 +590,7 @@ subroutine psb_z_csc_csmm(alpha,a,x,beta,y,info,trans) endif - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, min(m,n) y(i,:) = y(i,:) + alpha*x(i,:) end do @@ -1081,7 +1081,12 @@ function psb_z_csc_maxval(a) result(res) logical, parameter :: debug=.false. - res = dzero + if (a%is_unit()) then + res = done + else + res = dzero + end if + nnz = a%get_nzeros() if (allocated(a%val)) then nnz = min(nnz,size(a%val)) @@ -1112,7 +1117,11 @@ function psb_z_csc_csnmi(a) result(res) if (info /= psb_success_) then return end if - acc(:) = dzero + if (a%is_unit()) then + acc = done + else + acc = dzero + end if do i=1, nc do j=a%icp(i),a%icp(i+1)-1 acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) @@ -1149,7 +1158,11 @@ function psb_z_csc_csnm1(a) result(res) m = a%get_nrows() n = a%get_ncols() do j=1, n - acc = dzero + if (a%is_unit()) then + acc = done + else + acc = dzero + end if do k=a%icp(j),a%icp(j+1)-1 acc = acc + abs(a%val(k)) end do @@ -1187,19 +1200,17 @@ subroutine psb_z_csc_colsum(d,a) end if do i = 1, a%get_ncols() - d(i) = zzero + if (a%is_unit()) then + d(i) = zone + else + d(i) = zzero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + (a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_ncols() - d(i) = d(i) + zone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1242,13 +1253,18 @@ subroutine psb_z_csc_aclsum(d,a) do i = 1, a%get_ncols() - d(i) = dzero + if (a%is_unit()) then + d(i) = done + else + d(i) = dzero + end if + do j=a%icp(i),a%icp(i+1)-1 d(i) = d(i) + abs(a%val(j)) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, a%get_ncols() d(i) = d(i) + done end do @@ -1295,7 +1311,11 @@ subroutine psb_z_csc_rowsum(d,a) goto 9999 end if - d = zzero + if (a%is_unit()) then + d = zone + else + d = zzero + end if do i=1, m do j=a%icp(i),a%icp(i+1)-1 @@ -1304,13 +1324,6 @@ subroutine psb_z_csc_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + zone - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1352,7 +1365,12 @@ subroutine psb_z_csc_arwsum(d,a) goto 9999 end if - d = dzero + if (a%is_unit()) then + d = done + else + d = dzero + end if + do i=1, m do j=a%icp(i),a%icp(i+1)-1 k = a%ia(j) @@ -1360,13 +1378,6 @@ subroutine psb_z_csc_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then - do i=1, a%get_nrows() - d(i) = d(i) + done - end do - end if - - return call psb_erractionrestore(err_act) return @@ -1408,7 +1419,7 @@ subroutine psb_z_csc_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = zone else do i=1, mnm @@ -1445,13 +1456,14 @@ subroutine psb_z_csc_scal(d,a,info,side) use psb_string_mod implicit none class(psb_z_csc_sparse_mat), intent(inout) :: a - complex(psb_dpk_), intent(in) :: d(:) - integer(psb_ipk_), intent(out) :: info - character, intent(in), optional :: side + complex(psb_dpk_), intent(in) :: d(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side integer(psb_ipk_) :: err_act,mnm, i, j, n + type(psb_z_coo_sparse_mat) :: tmp integer(psb_ipk_) :: ierr(5) - character(len=20) :: name='scal' + character(len=20) :: name='scal' character :: side_ logical :: left logical, parameter :: debug=.false. @@ -1464,6 +1476,10 @@ subroutine psb_z_csc_scal(d,a,info,side) side_ = psb_toupper(side) end if + if (a%is_unit()) then + call a%add_unit_diag() + end if + left = (side_ == 'L') if (left) then @@ -1524,6 +1540,9 @@ subroutine psb_z_csc_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d diff --git a/base/serial/impl/psb_z_csr_impl.f90 b/base/serial/impl/psb_z_csr_impl.f90 index 2b7cc374..45781ab6 100644 --- a/base/serial/impl/psb_z_csr_impl.f90 +++ b/base/serial/impl/psb_z_csr_impl.f90 @@ -380,7 +380,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i) = y(i) + alpha*x(i) end do @@ -733,7 +733,7 @@ contains endif - if (is_triangle.and.is_unit) then + if (is_unit) then do i=1, min(m,n) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) end do @@ -1398,7 +1398,7 @@ subroutine psb_z_csr_rowsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + zone end do @@ -1453,7 +1453,7 @@ subroutine psb_z_csr_arwsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, m d(i) = d(i) + done end do @@ -1509,7 +1509,7 @@ subroutine psb_z_csr_colsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + zone end do @@ -1566,7 +1566,7 @@ subroutine psb_z_csr_aclsum(d,a) end do end do - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then do i=1, n d(i) = d(i) + done end do @@ -1613,7 +1613,7 @@ subroutine psb_z_csr_get_diag(a,d,info) end if - if (a%is_triangle().and.a%is_unit()) then + if (a%is_unit()) then d(1:mnm) = zone else do i=1, mnm @@ -1665,6 +1665,10 @@ subroutine psb_z_csr_scal(d,a,info,side) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if + side_ = 'L' if (present(side)) then side_ = psb_toupper(side) @@ -1734,6 +1738,9 @@ subroutine psb_z_csr_scals(d,a,info) info = psb_success_ call psb_erractionsave(err_act) + if (a%is_unit()) then + call a%add_unit_diag() + end if do i=1,a%get_nzeros() a%val(i) = a%val(i) * d