diff --git a/base/modules/serial/psb_d_base_mat_mod.F90 b/base/modules/serial/psb_d_base_mat_mod.F90 index 24c2097f7..92c0d1dd0 100644 --- a/base/modules/serial/psb_d_base_mat_mod.F90 +++ b/base/modules/serial/psb_d_base_mat_mod.F90 @@ -375,6 +375,8 @@ module psb_d_base_mat_mod procedure, pass(a) :: mv_from_fmt => psb_ld_mv_coo_from_fmt procedure, pass(a) :: cp_to_icoo => psb_ld_cp_coo_to_icoo procedure, pass(a) :: cp_from_icoo => psb_ld_cp_coo_from_icoo + procedure, pass(a) :: tril => psb_ld_coo_tril + procedure, pass(a) :: triu => psb_ld_coo_triu procedure, pass(a) :: csput_a => psb_ld_coo_csput_a procedure, pass(a) :: get_diag => psb_ld_coo_get_diag @@ -3548,6 +3550,92 @@ module psb_d_base_mat_mod integer(psb_ipk_), intent(out) :: info end subroutine psb_ld_coo_mold end interface + ! + !> Function tril: + !! \memberof psb_d_coo_sparse_mat + !! \brief Copy the lower triangle, i.e. all entries + !! A(I,J) such that J-I <= DIAG + !! default value is DIAG=0, i.e. lower triangle up to + !! the main diagonal. + !! DIAG=-1 means copy the strictly lower triangle + !! DIAG= 1 means copy the lower triangle plus the first diagonal + !! of the upper triangle. + !! Moreover, apply a clipping by copying entries A(I,J) only if + !! IMIN<=I<=IMAX + !! JMIN<=J<=JMAX + !! + !! \param l the output (sub)matrix + !! \param info return code + !! \param diag [0] the last diagonal (J-I) to be considered. + !! \param imin [1] the minimum row index we are interested in + !! \param imax [a\%get_nrows()] the minimum row index we are interested in + !! \param jmin [1] minimum col index + !! \param jmax [a\%get_ncols()] maximum col index + !! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:)) + !! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1] + !! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1] + !! ( iren cannot be specified with rscale/cscale) + !! \param append [false] append to ia,ja + !! \param nzin [none] if append, then first new entry should go in entry nzin+1 + !! \param u [none] copy of the complementary triangle + !! + ! + interface + subroutine psb_ld_coo_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) + import + class(psb_ld_coo_sparse_mat), intent(in) :: a + class(psb_ld_coo_sparse_mat), intent(out) :: l + integer(psb_ipk_),intent(out) :: info + integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + class(psb_ld_coo_sparse_mat), optional, intent(out) :: u + end subroutine psb_ld_coo_tril + end interface + + ! + !> Function triu: + !! \memberof psb_d_coo_sparse_mat + !! \brief Copy the upper triangle, i.e. all entries + !! A(I,J) such that DIAG <= J-I + !! default value is DIAG=0, i.e. upper triangle from + !! the main diagonal up. + !! DIAG= 1 means copy the strictly upper triangle + !! DIAG=-1 means copy the upper triangle plus the first diagonal + !! of the lower triangle. + !! Moreover, apply a clipping by copying entries A(I,J) only if + !! IMIN<=I<=IMAX + !! JMIN<=J<=JMAX + !! Optionally copies the lower triangle at the same time + !! + !! \param u the output (sub)matrix + !! \param info return code + !! \param diag [0] the last diagonal (J-I) to be considered. + !! \param imin [1] the minimum row index we are interested in + !! \param imax [a\%get_nrows()] the minimum row index we are interested in + !! \param jmin [1] minimum col index + !! \param jmax [a\%get_ncols()] maximum col index + !! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:)) + !! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1] + !! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1] + !! ( iren cannot be specified with rscale/cscale) + !! \param append [false] append to ia,ja + !! \param nzin [none] if append, then first new entry should go in entry nzin+1 + !! \param l [none] copy of the complementary triangle + !! + ! + interface + subroutine psb_ld_coo_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) + import + class(psb_ld_coo_sparse_mat), intent(in) :: a + class(psb_ld_coo_sparse_mat), intent(out) :: u + integer(psb_ipk_),intent(out) :: info + integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + class(psb_ld_coo_sparse_mat), optional, intent(out) :: l + end subroutine psb_ld_coo_triu + end interface ! diff --git a/base/serial/impl/psb_d_coo_impl.F90 b/base/serial/impl/psb_d_coo_impl.F90 index 99d4e2f63..1f53606aa 100644 --- a/base/serial/impl/psb_d_coo_impl.F90 +++ b/base/serial/impl/psb_d_coo_impl.F90 @@ -3908,7 +3908,7 @@ subroutine psb_d_coo_triu(a,u,info,& i = a%ia(k) j = a%ja(k) if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then - if ((j-i)=diag_) then nzin = nzin + 1 u%ia(nzin) = i @@ -5904,6 +5904,601 @@ subroutine psb_ld_coo_reinit(a,clear) end subroutine psb_ld_coo_reinit +! +! COO implementation of tril/triu +! +subroutine psb_ld_coo_tril(a,l,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale,u) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_tril + implicit none + + class(psb_ld_coo_sparse_mat), intent(in) :: a + class(psb_ld_coo_sparse_mat), intent(out) :: l + integer(psb_ipk_),intent(out) :: info + integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + class(psb_ld_coo_sparse_mat), optional, intent(out) :: u + + integer(psb_ipk_) :: err_act + integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, & + & nzlin, nzuin, nz, nzin, nzout, i, j, k + character(len=20) :: name='tril' + logical :: rscale_, cscale_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + + if (present(diag)) then + diag_ = diag + else + diag_ = 0 + end if + if (present(imin)) then + imin_ = imin + else + imin_ = 1 + end if + if (present(imax)) then + imax_ = imax + else + imax_ = a%get_nrows() + end if + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + end if + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + end if + if (present(rscale)) then + rscale_ = rscale + else + rscale_ = .true. + end if + if (present(cscale)) then + cscale_ = cscale + else + cscale_ = .true. + end if + + if (rscale_) then + mb = imax_ - imin_ +1 + else + mb = imax_ + endif + if (cscale_) then + nb = jmax_ - jmin_ +1 + else + nb = jmax_ + endif +#if defined(PSB_OPENMP) + block + integer(psb_lpk_), allocatable :: lrws(:),urws(:) + integer(psb_lpk_) :: lpnt, upnt, lnz, unz + call psb_realloc(mb,lrws,info) + !$omp workshare + lrws(:) = 0 + !$omp end workshare + nz = a%get_nzeros() + call l%allocate(mb,nb,nz) + !write(0,*) 'Invocation of COO%TRIL', present(u),nz + if (present(u)) then + nzlin = l%get_nzeros() ! At this point it should be 0 + call u%allocate(mb,nb,nz) + nzuin = u%get_nzeros() ! At this point it should be 0 + if (info == 0) call psb_realloc(mb,urws,info) + !$omp workshare + urws(:) = 0 + !$omp end workshare + !write(0,*) 'omp version of COO%TRIL/TRIU' + lnz = 0 + unz = 0 + !$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz) + loop1: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)<=diag_) then + !$omp atomic update + lrws(i-imin_+1) = lrws(i-imin_+1) +1 + !$omp end atomic + lnz = lnz + 1 + else + !$omp atomic update + urws(i-imin_+1) = urws(i-imin_+1) +1 + !$omp end atomic + unz = unz + 1 + end if + end if + end do loop1 + !$omp end parallel do + + call psi_exscan(mb,lrws,info) + call psi_exscan(mb,urws,info) + !write(0,*) lrws(:), urws(:) + !$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a) + loop2: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)<=diag_) then + !$omp atomic capture + lrws(i-imin_+1) = lrws(i-imin_+1) +1 + lpnt = lrws(i-imin_+1) + !$omp end atomic + l%ia(lpnt) = a%ia(k) + l%ja(lpnt) = a%ja(k) + l%val(lpnt) = a%val(k) + else + !$omp atomic capture + urws(i-imin_+1) = urws(i-imin_+1) +1 + upnt = urws(i-imin_+1) + !$omp end atomic + u%ia(upnt) = a%ia(k) + u%ja(upnt) = a%ja(k) + u%val(upnt) = a%val(k) + end if + end if + end do loop2 + !$omp end parallel do + !write(0,*) 'End of copyout',lnz,unz + call l%set_nzeros(lnz) + call l%fix(info) + call u%set_nzeros(unz) + call u%fix(info) + nzout = u%get_nzeros() + if (rscale_) then + !$omp workshare + u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1 + !$omp end workshare + end if + if (cscale_) then + !$omp workshare + u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1 + !$omp end workshare + end if + if ((diag_ >=-1).and.(imin_ == jmin_)) then + call u%set_triangle(.true.) + call u%set_lower(.false.) + end if + else + lnz = 0 + !$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws) reduction(+: lnz) + loop3: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)<=diag_) then + !$omp atomic update + lrws(i-imin_+1) = lrws(i-imin_+1) +1 + !$omp end atomic + lnz = lnz + 1 + end if + end if + end do loop3 + !$omp end parallel do + call psi_exscan(mb,lrws,info) + !$omp parallel do private(i,j,k,lpnt) shared(imin_,imax_,a) + loop4: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)<=diag_) then + !$omp atomic capture + lrws(i-imin_+1) = lrws(i-imin_+1) +1 + lpnt = lrws(i-imin_+1) + !$omp end atomic + l%ia(lpnt) = a%ia(k) + l%ja(lpnt) = a%ja(k) + l%val(lpnt) = a%val(k) + end if + end if + end do loop4 + !$omp end parallel do + call l%set_nzeros(lnz) + call l%fix(info) + end if + nzout = l%get_nzeros() + if (rscale_) then + !$omp workshare + l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 + !$omp end workshare + end if + if (cscale_) then + !$omp workshare + l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 + !$omp end workshare + end if + + if ((diag_ <= 0).and.(imin_ == jmin_)) then + call l%set_triangle(.true.) + call l%set_lower(.true.) + end if + end block + +#else + nz = a%get_nzeros() + call l%allocate(mb,nb,nz) + if (present(u)) then + nzlin = l%get_nzeros() ! At this point it should be 0 + call u%allocate(mb,nb,nz) + nzuin = u%get_nzeros() ! At this point it should be 0 + associate(val =>a%val, ja => a%ja, ia=>a%ia) + loop1: do k=1,nz + i = ia(k) + j = ja(k) + if ((jmin_<=j).and.(j<=jmax_)) then + if ((j-i)<=diag_) then + nzlin = nzlin + 1 + l%ia(nzlin) = i + l%ja(nzlin) = ja(k) + l%val(nzlin) = val(k) + else + nzuin = nzuin + 1 + u%ia(nzuin) = i + u%ja(nzuin) = ja(k) + u%val(nzuin) = val(k) + end if + end if + end do loop1 + end associate + + call l%set_nzeros(nzlin) + call u%set_nzeros(nzuin) + call u%fix(info) + nzout = u%get_nzeros() + if (rscale_) & + & u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1 + if ((diag_ >=-1).and.(imin_ == jmin_)) then + call u%set_triangle(.true.) + call u%set_lower(.false.) + end if + else + nzin = l%get_nzeros() ! At this point it should be 0 + + associate(val =>a%val, ja => a%ja, ia=>a%ia) + loop2: do k=1,nz + i = ia(k) + j = ja(k) + if ((jmin_<=j).and.(j<=jmax_)) then + if ((j-i)<=diag_) then + nzin = nzin + 1 + l%ia(nzin) = i + l%ja(nzin) = j + l%val(nzin) = val(k) + end if + end if + end do loop2 + end associate + call l%set_nzeros(nzin) + end if + call l%fix(info) + nzout = l%get_nzeros() + if (rscale_) & + & l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 + if ((diag_ <= 0).and.(imin_ == jmin_)) then + call l%set_triangle(.true.) + call l%set_lower(.true.) + end if +#endif + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_ld_coo_tril + +subroutine psb_ld_coo_triu(a,u,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale,l) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_triu + implicit none + + class(psb_ld_coo_sparse_mat), intent(in) :: a + class(psb_ld_coo_sparse_mat), intent(out) :: u + integer(psb_ipk_),intent(out) :: info + integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + class(psb_ld_coo_sparse_mat), optional, intent(out) :: l + + integer(psb_ipk_) :: err_act + integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, & + &nzlin, nzuin, nz, nzin, nzout, i, j, k + character(len=20) :: name='triu' + logical :: rscale_, cscale_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + + if (present(diag)) then + diag_ = diag + else + diag_ = 0 + end if + if (present(imin)) then + imin_ = imin + else + imin_ = 1 + end if + if (present(imax)) then + imax_ = imax + else + imax_ = a%get_nrows() + end if + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + end if + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + end if + if (present(rscale)) then + rscale_ = rscale + else + rscale_ = .true. + end if + if (present(cscale)) then + cscale_ = cscale + else + cscale_ = .true. + end if + + if (rscale_) then + mb = imax_ - imin_ +1 + else + mb = imax_ + endif + if (cscale_) then + nb = jmax_ - jmin_ +1 + else + nb = jmax_ + endif + +#if defined(PSB_OPENMP) + block + integer(psb_lpk_), allocatable :: lrws(:),urws(:) + integer(psb_lpk_) :: lpnt, upnt, lnz, unz + call psb_realloc(mb,urws,info) + !$omp workshare + urws(:) = 0 + !$omp end workshare + nz = a%get_nzeros() + call u%allocate(mb,nb,nz) + !write(0,*) 'Invocation of COO%TRIL', present(u),nz + if (present(l)) then + nzuin = u%get_nzeros() ! At this point it should be 0 + call l%allocate(mb,nb,nz) + nzlin = l%get_nzeros() ! At this point it should be 0 + if (info == 0) call psb_realloc(mb,urws,info) + !$omp workshare + lrws(:) = 0 + !$omp end workshare + !write(0,*) 'omp version of COO%TRIL/TRIU' + lnz = 0 + unz = 0 + !$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz) + loop1: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)>=diag_) then + !$omp atomic update + lrws(i-imin_+1) = lrws(i-imin_+1) +1 + !$omp end atomic + lnz = lnz + 1 + else + !$omp atomic update + urws(i-imin_+1) = urws(i-imin_+1) +1 + !$omp end atomic + unz = unz + 1 + end if + end if + end do loop1 + !$omp end parallel do + + call psi_exscan(mb,lrws,info) + call psi_exscan(mb,urws,info) + !write(0,*) lrws(:), urws(:) + !$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a) + loop2: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)>=diag_) then + !$omp atomic capture + lrws(i-imin_+1) = lrws(i-imin_+1) +1 + lpnt = lrws(i-imin_+1) + !$omp end atomic + l%ia(lpnt) = a%ia(k) + l%ja(lpnt) = a%ja(k) + l%val(lpnt) = a%val(k) + else + !$omp atomic capture + urws(i-imin_+1) = urws(i-imin_+1) +1 + upnt = urws(i-imin_+1) + !$omp end atomic + u%ia(upnt) = a%ia(k) + u%ja(upnt) = a%ja(k) + u%val(upnt) = a%val(k) + end if + end if + end do loop2 + !$omp end parallel do + !write(0,*) 'End of copyout',lnz,unz + call l%set_nzeros(lnz) + call l%fix(info) + call u%set_nzeros(unz) + call u%fix(info) + nzout = l%get_nzeros() + if (rscale_) then + !$omp workshare + l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 + !$omp end workshare + end if + if (cscale_) then + !$omp workshare + l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 + !$omp end workshare + end if + if ((diag_ <=-1).and.(imin_ == jmin_)) then + call l%set_triangle(.true.) + call l%set_lower(.false.) + end if + else + unz = 0 + !$omp parallel do private(i,j,k) shared(imin_,imax_,a,urws) reduction(+: unz) + loop3: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)>=diag_) then + !$omp atomic update + urws(i-imin_+1) = urws(i-imin_+1) +1 + !$omp end atomic + unz = unz + 1 + end if + end if + end do loop3 + !$omp end parallel do + call psi_exscan(mb,urws,info) + !$omp parallel do private(i,j,k,upnt) shared(imin_,imax_,a) + loop4: do k=1,nz + i = a%ia(k) + j = a%ja(k) + if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then + if ((j-i)>=diag_) then + !$omp atomic capture + urws(i-imin_+1) = urws(i-imin_+1) +1 + upnt = urws(i-imin_+1) + !$omp end atomic + u%ia(upnt) = a%ia(k) + u%ja(upnt) = a%ja(k) + u%val(upnt) = a%val(k) + end if + end if + end do loop4 + !$omp end parallel do + call u%set_nzeros(unz) + call u%fix(info) + end if + nzout = u%get_nzeros() + if (rscale_) then + !$omp workshare + u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1 + !$omp end workshare + end if + if (cscale_) then + !$omp workshare + u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1 + !$omp end workshare + end if + + if ((diag_ >= 0).and.(imin_ == jmin_)) then + call u%set_triangle(.true.) + call u%set_upper(.true.) + end if + end block + + +#else + nz = a%get_nzeros() + call u%allocate(mb,nb,nz) + if (present(l)) then + nzuin = u%get_nzeros() ! At this point it should be 0 + call l%allocate(mb,nb,nz) + nzlin = l%get_nzeros() ! At this point it should be 0 + associate(val =>a%val, ja => a%ja, ia=>a%ia) + loop1: do k=1,nz + i = ia(k) + j = ja(k) + if ((jmin_<=j).and.(j<=jmax_)) then + if ((j-i)>=diag_) then + nzuin = nzuin + 1 + u%ia(nzuin) = i + u%ja(nzuin) = ja(k) + u%val(nzuin) = val(k) + else + nzlin = nzlin + 1 + l%ia(nzlin) = i + l%ja(nzlin) = ja(k) + l%val(nzlin) = val(k) + end if + end if + end do loop1 + end associate + call u%set_nzeros(nzuin) + call l%set_nzeros(nzlin) + call l%fix(info) + nzout = l%get_nzeros() + if (rscale_) & + & l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 + if ((diag_ <=1).and.(imin_ == jmin_)) then + call l%set_triangle(.true.) + call l%set_lower(.true.) + end if + else + nzin = u%get_nzeros() ! At this point it should be 0 + associate(val =>a%val, ja => a%ja, ia=>a%ia) + loop2: do k=1,nz + i = ia(k) + j = ja(k) + if ((jmin_<=j).and.(j<=jmax_)) then + if ((j-i)>=diag_) then + nzin = nzin + 1 + u%ia(nzin) = i + u%ja(nzin) = ja(k) + u%val(nzin) = val(k) + end if + end if + end do loop2 + end associate + call u%set_nzeros(nzin) + end if + call u%fix(info) + nzout = u%get_nzeros() + if (rscale_) & + & u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1 + + if ((diag_ >= 0).and.(imin_ == jmin_)) then + call u%set_triangle(.true.) + call u%set_lower(.false.) + end if +#endif + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_ld_coo_triu subroutine psb_ld_coo_trim(a) diff --git a/base/tools/psb_dglobtril.F90 b/base/tools/psb_dglobtril.F90 index c4dcd819c..af08081d6 100644 --- a/base/tools/psb_dglobtril.F90 +++ b/base/tools/psb_dglobtril.F90 @@ -50,17 +50,19 @@ Subroutine psb_dglobtril(a,desc_a,b,info,diag,imin,imax,jmin,jmax) ! .. Local Scalars .. integer(psb_ipk_) :: i, j, err_act,m,& & nz - integer(psb_lpk_) :: gidx, lnz + integer(psb_lpk_) :: gidx, lnz, gnr, gnc,lnr,lnc type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np integer(psb_mpk_) :: icomm, minfo type(psb_d_coo_sparse_mat) :: dtcoo - type(psb_ld_coo_sparse_mat) :: ldtcoo + type(psb_ld_coo_sparse_mat) :: ldtcoo, ldcootril type(psb_ldspmat_type) :: ldglob,ldtril integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: ierr(5) character(len=20) :: name, ch_err + character(len=50) :: fname + integer :: iout name='psb_dglobtril' info = psb_success_ @@ -82,23 +84,25 @@ Subroutine psb_dglobtril(a,desc_a,b,info,diag,imin,imax,jmin,jmax) & Write(debug_unit,*) me,' ',trim(name),& & ': start',diag + gnr = desc_a%get_global_rows() + gnc = desc_a%get_global_cols() call a%a%cp_to_lcoo(ldtcoo,info) lnz = ldtcoo%get_nzeros() call desc_a%l2gip(ldtcoo%ia(1:lnz),info,owned=.false.) call desc_a%l2gip(ldtcoo%ja(1:lnz),info,owned=.false.) - call ldglob%mv_from(ldtcoo) - call ldglob%tril(ldtril,info,& - & diag=diag,imin=imin,imax=imax,jmin=jmin,jmax=jmax) - call ldglob%free() - call ldtril%mv_to(ldtcoo) - lnz = ldtcoo%get_nzeros() - call desc_a%g2lip(ldtcoo%ia(1:lnz),info,owned=.false.) - call desc_a%g2lip(ldtcoo%ja(1:lnz),info,owned=.false.) - call b%mv_from_lb(ldtcoo) - + call ldtcoo%tril(ldcootril,info,& + & diag=diag,imax=gnr,jmax=gnc) + lnz = ldcootril%get_nzeros() + call desc_a%g2lip(ldcootril%ia(1:lnz),info,owned=.false.) + call desc_a%g2lip(ldcootril%ja(1:lnz),info,owned=.false.) + lnc = desc_a%get_local_cols() + call ldcootril%set_nrows(lnr) + call ldcootril%set_ncols(lnc) + call ldcootril%fix(info) + call b%mv_from_lb(ldcootril) + if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),': end' - call psb_erractionrestore(err_act) return diff --git a/base/tools/psb_dglobtriu.F90 b/base/tools/psb_dglobtriu.F90 index 34a52da72..1eefbcbcb 100644 --- a/base/tools/psb_dglobtriu.F90 +++ b/base/tools/psb_dglobtriu.F90 @@ -50,13 +50,13 @@ Subroutine psb_dglobtriu(a,desc_a,b,info,diag,imin,imax,jmin,jmax) ! .. Local Scalars .. integer(psb_ipk_) :: i, j, err_act,m,& & nz - integer(psb_lpk_) :: gidx, lnz + integer(psb_lpk_) :: gidx, lnz, lnr, lnc, gnr, gnc type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me, np integer(psb_mpk_) :: icomm, minfo type(psb_d_coo_sparse_mat) :: dtcoo - type(psb_ld_coo_sparse_mat) :: ldtcoo + type(psb_ld_coo_sparse_mat) :: ldtcoo, ldcootriu type(psb_ldspmat_type) :: ldglob,ldtriu integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: ierr(5) @@ -82,20 +82,23 @@ Subroutine psb_dglobtriu(a,desc_a,b,info,diag,imin,imax,jmin,jmax) & Write(debug_unit,*) me,' ',trim(name),& & ': start',diag + gnr = desc_a%get_global_rows() + gnc = desc_a%get_global_cols() call a%a%cp_to_lcoo(ldtcoo,info) lnz = ldtcoo%get_nzeros() call desc_a%l2gip(ldtcoo%ia(1:lnz),info,owned=.false.) call desc_a%l2gip(ldtcoo%ja(1:lnz),info,owned=.false.) - call ldglob%mv_from(ldtcoo) - call ldglob%triu(ldtriu,info,& - & diag=diag,imin=imin,imax=imax,jmin=jmin,jmax=jmax) - call ldglob%free() - call ldtriu%mv_to(ldtcoo) - lnz = ldtcoo%get_nzeros() - call desc_a%g2lip(ldtcoo%ia(1:lnz),info,owned=.false.) - call desc_a%g2lip(ldtcoo%ja(1:lnz),info,owned=.false.) - call b%mv_from_lb(ldtcoo) - + call ldtcoo%triu(ldcootriu,info,& + & diag=diag,imax=gnr,jmax=gnc) + lnz = ldcootriu%get_nzeros() + call desc_a%g2lip(ldcootriu%ia(1:lnz),info,owned=.false.) + call desc_a%g2lip(ldcootriu%ja(1:lnz),info,owned=.false.) + lnc = desc_a%get_local_cols() + call ldcootriu%set_nrows(lnr) + call ldcootriu%set_ncols(lnc) + call ldcootriu%fix(info) + call b%mv_from_lb(ldcootriu) + if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),': end'