From 2e40d962e257d0e27b2db454da6663c77bdb1ad1 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Mon, 8 Jun 2026 13:21:26 +0200 Subject: [PATCH] New GMRES and MMIO_IMPL improved and templated --- linsolve/impl/psb_crgmres.f90 | 219 ++++++++++++++---------- linsolve/impl/psb_drgmres.f90 | 219 ++++++++++++++---------- linsolve/impl/psb_srgmres.f90 | 219 ++++++++++++++---------- linsolve/impl/psb_zrgmres.f90 | 219 ++++++++++++++---------- linsolve/psb_base_linsolve_conv_mod.f90 | 33 +++- util/psb_c_mmio_impl.f90 | 171 +++++++++--------- util/psb_d_mmio_impl.f90 | 53 ++++-- util/psb_s_mmio_impl.f90 | 124 ++++++++++---- util/psb_z_mmio_impl.f90 | 171 +++++++++--------- 9 files changed, 872 insertions(+), 556 deletions(-) diff --git a/linsolve/impl/psb_crgmres.f90 b/linsolve/impl/psb_crgmres.f90 index 673d13994..3c20b1e87 100644 --- a/linsolve/impl/psb_crgmres.f90 +++ b/linsolve/impl/psb_crgmres.f90 @@ -130,7 +130,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& type(psb_c_vect_type) :: w, w1, xt real(psb_spk_) :: tmp complex(psb_spk_) :: scal, gm, rti, rti1 - integer(psb_ipk_) ::litmax, naux, it, k, itrace_,& + integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,& & n_row, n_col, nl integer(psb_lpk_) :: mglob Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. @@ -172,24 +172,35 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then - litmax = itmax + itmax_ = itmax else - litmax = 1000 + itmax_ = 1000 endif if (present(itrace)) then @@ -253,12 +264,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(cone,b,czero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +287,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = czero errden = cone deps = eps @@ -324,27 +339,32 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = sone + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart + if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart if (itrace_ > 0) & & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) @@ -360,42 +380,18 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& call prec%apply(v(i),w1,desc_a,info) call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux) ! + call mgs(i,h,v,w,rs,c,s,desc_a,info) - do k = 1, i - h(k,i) = psb_gedot(v(k),w,desc_a,info) - call psb_geaxpby(-h(k,i),v(k),cone,w,desc_a,info) - end do - h(i+1,i) = psb_genrm2(w,desc_a,info) - scal=cone/h(i+1,i) - call psb_geaxpby(scal,w,czero,v(i+1),desc_a,info) - do k=2,i - call crot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1)) - enddo - - - rti = h(i,i) - rti1 = h(i+1,i) - call crotg(rti,rti1,tmp,s(i)) - c(i) = cmplx(tmp,szero) - call crot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i)) - h(i+1,i) = czero - call crot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! rst = rs - call w1%set(czero) - call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rst,size(rst,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rst(1:i) - do k=1, i - call psb_geaxpby(rst(k),v(k),cone,xt,desc_a,info) - end do - call prec%apply(xt,desc_a,info) - call psb_geaxpby(cone,x,cone,xt,desc_a,info) + call psb_geaxpby(cone,x,czero,xt,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call psb_geaxpby(cone,b,czero,w1,desc_a,info) call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux) rni = psb_geamax(w1,desc_a,info) @@ -403,8 +399,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +407,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = sone + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,28 +422,18 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(cone,xt,czero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:i) - call w1%set(czero) - do k=1, i - call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(cone,w,cone,x,desc_a,info) - end if + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,25 +446,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:nl) - call w1%set(czero) - do k=1, nl - call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(cone,w,cone,x,desc_a,info) - end if + + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + end select - if (itx >= litmax) then + if (itx >= itmax_) then if (itrace_ > 0) then if (mod(itx,itrace_)/=0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -502,6 +484,67 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& 9999 call psb_error_handler(err_act) return +contains + ! + ! Advance one step the Gram-Schmidt process, and apply + ! Givens rotations to keep H triangular + ! + subroutine mgs(n,h,v,w,rs,c,s,desc_a,info) + complex(psb_spk_) :: c(:), s(:), h(:,:), rs(:) + type(psb_c_vect_type) :: v(:), w + type(psb_desc_type) :: desc_a + complex(psb_spk_) :: scal, gm, rti, rti1 + real(psb_spk_) :: tmp + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + ! + + do k = 1, n + h(k,n) = psb_gedot(v(k),w,desc_a,info) + call psb_geaxpby(-h(k,n),v(k),cone,w,desc_a,info) + end do + h(n+1,n) = psb_genrm2(w,desc_a,info) + scal=cone/h(n+1,n) + call psb_geaxpby(scal,w,czero,v(n+1),desc_a,info) + do k=2,n + call crot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1)) + enddo + + rti = h(n,n) + rti1 = h(n+1,n) + call crotg(rti,rti1,tmp,s(n)) + c(n) = cmplx(tmp,szero) + call crot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n)) + call crot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n)) + end subroutine mgs + + ! + ! Rebuild solution X from the space V using the factor + ! stored in R + ! + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + complex(psb_spk_) :: c(:), s(:), rs(:), h(:,:) + type(psb_c_vect_type) :: v(:), w, w1, x + type(psb_desc_type) :: desc_a + class(psb_cprec_type) :: prec + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + + ! + ! + ! build x + ! + call ctrsm('l','u','n','n',n,1,cone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:n) + call w1%zero() + do k=1, n + call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(cone,w,cone,x,desc_a,info) + end subroutine rebuildx end subroutine psb_crgmres_vect diff --git a/linsolve/impl/psb_drgmres.f90 b/linsolve/impl/psb_drgmres.f90 index e03e043f7..7144f8955 100644 --- a/linsolve/impl/psb_drgmres.f90 +++ b/linsolve/impl/psb_drgmres.f90 @@ -130,7 +130,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& type(psb_d_vect_type) :: w, w1, xt real(psb_dpk_) :: tmp real(psb_dpk_) :: scal, gm, rti, rti1 - integer(psb_ipk_) ::litmax, naux, it, k, itrace_,& + integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,& & n_row, n_col, nl integer(psb_lpk_) :: mglob Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. @@ -172,24 +172,35 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then - litmax = itmax + itmax_ = itmax else - litmax = 1000 + itmax_ = 1000 endif if (present(itrace)) then @@ -253,12 +264,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(done,b,dzero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +287,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = dzero errden = done deps = eps @@ -324,27 +339,32 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = done + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart + if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart if (itrace_ > 0) & & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) @@ -360,42 +380,18 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& call prec%apply(v(i),w1,desc_a,info) call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux) ! + call mgs(i,h,v,w,rs,c,s,desc_a,info) - do k = 1, i - h(k,i) = psb_gedot(v(k),w,desc_a,info) - call psb_geaxpby(-h(k,i),v(k),done,w,desc_a,info) - end do - h(i+1,i) = psb_genrm2(w,desc_a,info) - scal=done/h(i+1,i) - call psb_geaxpby(scal,w,dzero,v(i+1),desc_a,info) - do k=2,i - call drot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1)) - enddo - - - rti = h(i,i) - rti1 = h(i+1,i) - call drotg(rti,rti1,tmp,s(i)) - c(i) = cmplx(tmp,dzero) - call drot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i)) - h(i+1,i) = dzero - call drot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! rst = rs - call w1%set(dzero) - call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rst,size(rst,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rst(1:i) - do k=1, i - call psb_geaxpby(rst(k),v(k),done,xt,desc_a,info) - end do - call prec%apply(xt,desc_a,info) - call psb_geaxpby(done,x,done,xt,desc_a,info) + call psb_geaxpby(done,x,dzero,xt,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call psb_geaxpby(done,b,dzero,w1,desc_a,info) call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux) rni = psb_geamax(w1,desc_a,info) @@ -403,8 +399,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +407,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = done + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,28 +422,18 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(done,xt,dzero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:i) - call w1%set(dzero) - do k=1, i - call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(done,w,done,x,desc_a,info) - end if + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,25 +446,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:nl) - call w1%set(dzero) - do k=1, nl - call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(done,w,done,x,desc_a,info) - end if + + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + end select - if (itx >= litmax) then + if (itx >= itmax_) then if (itrace_ > 0) then if (mod(itx,itrace_)/=0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -502,6 +484,67 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& 9999 call psb_error_handler(err_act) return +contains + ! + ! Advance one step the Gram-Schmidt process, and apply + ! Givens rotations to keep H triangular + ! + subroutine mgs(n,h,v,w,rs,c,s,desc_a,info) + real(psb_dpk_) :: c(:), s(:), h(:,:), rs(:) + type(psb_d_vect_type) :: v(:), w + type(psb_desc_type) :: desc_a + real(psb_dpk_) :: scal, gm, rti, rti1 + real(psb_dpk_) :: tmp + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + ! + + do k = 1, n + h(k,n) = psb_gedot(v(k),w,desc_a,info) + call psb_geaxpby(-h(k,n),v(k),done,w,desc_a,info) + end do + h(n+1,n) = psb_genrm2(w,desc_a,info) + scal=done/h(n+1,n) + call psb_geaxpby(scal,w,dzero,v(n+1),desc_a,info) + do k=2,n + call drot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1)) + enddo + + rti = h(n,n) + rti1 = h(n+1,n) + call drotg(rti,rti1,tmp,s(n)) + c(n) = cmplx(tmp,dzero) + call drot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n)) + call drot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n)) + end subroutine mgs + + ! + ! Rebuild solution X from the space V using the factor + ! stored in R + ! + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + real(psb_dpk_) :: c(:), s(:), rs(:), h(:,:) + type(psb_d_vect_type) :: v(:), w, w1, x + type(psb_desc_type) :: desc_a + class(psb_dprec_type) :: prec + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + + ! + ! + ! build x + ! + call dtrsm('l','u','n','n',n,1,done,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:n) + call w1%zero() + do k=1, n + call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(done,w,done,x,desc_a,info) + end subroutine rebuildx end subroutine psb_drgmres_vect diff --git a/linsolve/impl/psb_srgmres.f90 b/linsolve/impl/psb_srgmres.f90 index cf455385e..54c62588d 100644 --- a/linsolve/impl/psb_srgmres.f90 +++ b/linsolve/impl/psb_srgmres.f90 @@ -130,7 +130,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& type(psb_s_vect_type) :: w, w1, xt real(psb_spk_) :: tmp real(psb_spk_) :: scal, gm, rti, rti1 - integer(psb_ipk_) ::litmax, naux, it, k, itrace_,& + integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,& & n_row, n_col, nl integer(psb_lpk_) :: mglob Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. @@ -172,24 +172,35 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then - litmax = itmax + itmax_ = itmax else - litmax = 1000 + itmax_ = 1000 endif if (present(itrace)) then @@ -253,12 +264,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(sone,b,szero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +287,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = szero errden = sone deps = eps @@ -324,27 +339,32 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = sone + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart + if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart if (itrace_ > 0) & & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) @@ -360,42 +380,18 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& call prec%apply(v(i),w1,desc_a,info) call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux) ! + call mgs(i,h,v,w,rs,c,s,desc_a,info) - do k = 1, i - h(k,i) = psb_gedot(v(k),w,desc_a,info) - call psb_geaxpby(-h(k,i),v(k),sone,w,desc_a,info) - end do - h(i+1,i) = psb_genrm2(w,desc_a,info) - scal=sone/h(i+1,i) - call psb_geaxpby(scal,w,szero,v(i+1),desc_a,info) - do k=2,i - call srot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1)) - enddo - - - rti = h(i,i) - rti1 = h(i+1,i) - call srotg(rti,rti1,tmp,s(i)) - c(i) = cmplx(tmp,szero) - call srot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i)) - h(i+1,i) = szero - call srot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! rst = rs - call w1%set(szero) - call strsm('l','u','n','n',i,1,sone,h,size(h,1),rst,size(rst,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rst(1:i) - do k=1, i - call psb_geaxpby(rst(k),v(k),sone,xt,desc_a,info) - end do - call prec%apply(xt,desc_a,info) - call psb_geaxpby(sone,x,sone,xt,desc_a,info) + call psb_geaxpby(sone,x,szero,xt,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call psb_geaxpby(sone,b,szero,w1,desc_a,info) call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux) rni = psb_geamax(w1,desc_a,info) @@ -403,8 +399,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +407,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = sone + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,28 +422,18 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(sone,xt,szero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:i) - call w1%set(szero) - do k=1, i - call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(sone,w,sone,x,desc_a,info) - end if + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,25 +446,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:nl) - call w1%set(szero) - do k=1, nl - call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(sone,w,sone,x,desc_a,info) - end if + + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + end select - if (itx >= litmax) then + if (itx >= itmax_) then if (itrace_ > 0) then if (mod(itx,itrace_)/=0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -502,6 +484,67 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& 9999 call psb_error_handler(err_act) return +contains + ! + ! Advance one step the Gram-Schmidt process, and apply + ! Givens rotations to keep H triangular + ! + subroutine mgs(n,h,v,w,rs,c,s,desc_a,info) + real(psb_spk_) :: c(:), s(:), h(:,:), rs(:) + type(psb_s_vect_type) :: v(:), w + type(psb_desc_type) :: desc_a + real(psb_spk_) :: scal, gm, rti, rti1 + real(psb_spk_) :: tmp + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + ! + + do k = 1, n + h(k,n) = psb_gedot(v(k),w,desc_a,info) + call psb_geaxpby(-h(k,n),v(k),sone,w,desc_a,info) + end do + h(n+1,n) = psb_genrm2(w,desc_a,info) + scal=sone/h(n+1,n) + call psb_geaxpby(scal,w,szero,v(n+1),desc_a,info) + do k=2,n + call srot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1)) + enddo + + rti = h(n,n) + rti1 = h(n+1,n) + call srotg(rti,rti1,tmp,s(n)) + c(n) = cmplx(tmp,szero) + call srot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n)) + call srot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n)) + end subroutine mgs + + ! + ! Rebuild solution X from the space V using the factor + ! stored in R + ! + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + real(psb_spk_) :: c(:), s(:), rs(:), h(:,:) + type(psb_s_vect_type) :: v(:), w, w1, x + type(psb_desc_type) :: desc_a + class(psb_sprec_type) :: prec + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + + ! + ! + ! build x + ! + call strsm('l','u','n','n',n,1,sone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:n) + call w1%zero() + do k=1, n + call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(sone,w,sone,x,desc_a,info) + end subroutine rebuildx end subroutine psb_srgmres_vect diff --git a/linsolve/impl/psb_zrgmres.f90 b/linsolve/impl/psb_zrgmres.f90 index 1d2aa3484..fb7ae52e0 100644 --- a/linsolve/impl/psb_zrgmres.f90 +++ b/linsolve/impl/psb_zrgmres.f90 @@ -130,7 +130,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& type(psb_z_vect_type) :: w, w1, xt real(psb_dpk_) :: tmp complex(psb_dpk_) :: scal, gm, rti, rti1 - integer(psb_ipk_) ::litmax, naux, it, k, itrace_,& + integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,& & n_row, n_col, nl integer(psb_lpk_) :: mglob Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. @@ -172,24 +172,35 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& if (present(istop)) then istop_ = istop else - istop_ = 2 + istop_ = psb_get_istop_default() endif -! -! ISTOP_ = 1: Normwise backward error, infinity norm -! ISTOP_ = 2: ||r||/||b||, 2-norm -! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then + + if (.not.psb_is_valid_istop(istop_)) then info=psb_err_invalid_istop_ err=info call psb_errpush(info,name,i_err=(/istop_/)) goto 9999 - endif + end if + ! + ! istop_ = 1: normwise backward error, infinity norm + ! istop_ = 2: ||r||/||b|| norm 2 + ! + select case(istop_) + case(psb_istop_ani_,psb_istop_bn2_,& + & psb_istop_rn2_abs_,psb_istop_rrn2_) + ! nothing needed + case default + ! should never get here + info=psb_err_internal_error_ + err=info + call psb_errpush(info,name,a_err="invalid istop_") + goto 9999 + end select if (present(itmax)) then - litmax = itmax + itmax_ = itmax else - litmax = 1000 + itmax_ = 1000 endif if (present(itrace)) then @@ -253,12 +264,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& & w%get_nrows(),w1%get_nrows() - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ani = psb_spnrmi(a,desc_a,info) bni = psb_geamax(b,desc_a,info) - else if (istop_ == 2) then - bn2 = psb_genrm2(b,desc_a,info) - else if (istop_ == 3) then + case(psb_istop_bn2_) + bn2 = psb_genrm2(b,desc_a,info) + case(psb_istop_rn2_abs_) + ! do nothing + case(psb_istop_rrn2_) call psb_geaxpby(zone,b,zzero,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ @@ -273,7 +287,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& goto 9999 end if r0n2 = psb_genrm2(v(1),desc_a,info) - endif + end select + errnum = zzero errden = zone deps = eps @@ -324,27 +339,32 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& ! ! check convergence ! - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) rni = psb_geamax(v(1),desc_a,info) xni = psb_geamax(x,desc_a,info) errnum = rni errden = (ani*xni+bni) - else if (istop_ == 2) then + case(psb_istop_bn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = bn2 - else if (istop_ == 3) then + case(psb_istop_rn2_abs_) + rni = psb_genrm2(v(1),desc_a,info) + errnum = rni + errden = done + case(psb_istop_rrn2_) rni = psb_genrm2(v(1),desc_a,info) errnum = rni errden = r0n2 - endif + end select if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) goto 9999 end if - if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart + if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart if (itrace_ > 0) & & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) @@ -360,42 +380,18 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& call prec%apply(v(i),w1,desc_a,info) call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux) ! + call mgs(i,h,v,w,rs,c,s,desc_a,info) - do k = 1, i - h(k,i) = psb_gedot(v(k),w,desc_a,info) - call psb_geaxpby(-h(k,i),v(k),zone,w,desc_a,info) - end do - h(i+1,i) = psb_genrm2(w,desc_a,info) - scal=zone/h(i+1,i) - call psb_geaxpby(scal,w,zzero,v(i+1),desc_a,info) - do k=2,i - call zrot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1)) - enddo - - - rti = h(i,i) - rti1 = h(i+1,i) - call zrotg(rti,rti1,tmp,s(i)) - c(i) = cmplx(tmp,dzero) - call zrot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i)) - h(i+1,i) = zzero - call zrot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i)) - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) ! ! build x and then compute the residual and its infinity norm ! rst = rs - call w1%set(zzero) - call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rst,size(rst,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rst(1:i) - do k=1, i - call psb_geaxpby(rst(k),v(k),zone,xt,desc_a,info) - end do - call prec%apply(xt,desc_a,info) - call psb_geaxpby(zone,x,zone,xt,desc_a,info) + call psb_geaxpby(zone,x,zzero,xt,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call psb_geaxpby(zone,b,zzero,w1,desc_a,info) call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux) rni = psb_geamax(w1,desc_a,info) @@ -403,8 +399,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& errnum = rni errden = (ani*xni+bni) ! - - else if (istop_ == 2) then + case(psb_istop_bn2_) ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -412,7 +407,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = bn2 - else if (istop_ == 3) then + + case(psb_istop_rn2_abs_) + rni = abs(rs(i+1)) + errnum = rni + errden = done + + case(psb_istop_rrn2_) + ! ! compute the residual 2-norm as byproduct of the solution ! procedure of the least-squares problem @@ -420,28 +422,18 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& rni = abs(rs(i+1)) errnum = rni errden = r0n2 - endif + end select if (errnum <= eps*errden) then - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(zone,xt,zzero,x,desc_a,info) ! = x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:i) - call w1%set(zzero) - do k=1, i - call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(zone,w,zone,x,desc_a,info) - end if + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + + end select if (itrace_ > 0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -454,25 +446,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& end do inner - if (istop_ == 1) then + select case(istop_) + case(psb_istop_ani_) call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt - else if (istop_ == 2) then - ! - ! build x - ! - call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Rebuild x-> RS:',rs(1:nl) - call w1%set(zzero) - do k=1, nl - call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info) - end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(zone,w,zone,x,desc_a,info) - end if + + case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + end select - if (itx >= litmax) then + if (itx >= itmax_) then if (itrace_ > 0) then if (mod(itx,itrace_)/=0) & & call log_conv(methdname,me,itx,ione,errnum,errden,deps) @@ -502,6 +484,67 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& 9999 call psb_error_handler(err_act) return +contains + ! + ! Advance one step the Gram-Schmidt process, and apply + ! Givens rotations to keep H triangular + ! + subroutine mgs(n,h,v,w,rs,c,s,desc_a,info) + complex(psb_dpk_) :: c(:), s(:), h(:,:), rs(:) + type(psb_z_vect_type) :: v(:), w + type(psb_desc_type) :: desc_a + complex(psb_dpk_) :: scal, gm, rti, rti1 + real(psb_dpk_) :: tmp + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + ! + + do k = 1, n + h(k,n) = psb_gedot(v(k),w,desc_a,info) + call psb_geaxpby(-h(k,n),v(k),zone,w,desc_a,info) + end do + h(n+1,n) = psb_genrm2(w,desc_a,info) + scal=zone/h(n+1,n) + call psb_geaxpby(scal,w,zzero,v(n+1),desc_a,info) + do k=2,n + call zrot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1)) + enddo + + rti = h(n,n) + rti1 = h(n+1,n) + call zrotg(rti,rti1,tmp,s(n)) + c(n) = cmplx(tmp,dzero) + call zrot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n)) + call zrot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n)) + end subroutine mgs + + ! + ! Rebuild solution X from the space V using the factor + ! stored in R + ! + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + complex(psb_dpk_) :: c(:), s(:), rs(:), h(:,:) + type(psb_z_vect_type) :: v(:), w, w1, x + type(psb_desc_type) :: desc_a + class(psb_zprec_type) :: prec + integer(psb_ipk_) :: info + integer(psb_ipk_) :: k,n + + ! + ! + ! build x + ! + call ztrsm('l','u','n','n',n,1,zone,h,size(h,1),rs,size(rs,1)) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' Rebuild x-> RS:',rs(1:n) + call w1%zero() + do k=1, n + call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info) + end do + call prec%apply(w1,w,desc_a,info) + call psb_geaxpby(zone,w,zone,x,desc_a,info) + end subroutine rebuildx end subroutine psb_zrgmres_vect diff --git a/linsolve/psb_base_linsolve_conv_mod.f90 b/linsolve/psb_base_linsolve_conv_mod.f90 index 8e45614f4..0e45d11e9 100644 --- a/linsolve/psb_base_linsolve_conv_mod.f90 +++ b/linsolve/psb_base_linsolve_conv_mod.f90 @@ -41,10 +41,20 @@ Module psb_base_linsolve_conv_mod module procedure psb_d_end_conv end interface + ! + integer(psb_ipk_), parameter :: psb_istop_min_ = 1 + integer(psb_ipk_), parameter :: psb_istop_ani_ = 1 + integer(psb_ipk_), parameter :: psb_istop_bn2_ = 2 + integer(psb_ipk_), parameter :: psb_istop_rn2_abs_ = 3 + integer(psb_ipk_), parameter :: psb_istop_rrn2_ = 4 + integer(psb_ipk_), parameter :: psb_istop_scbn2_ = 5 + integer(psb_ipk_), parameter :: psb_istop_max_ = 3 + + ! Fields in controls and values integer(psb_ipk_), parameter :: psb_ik_bni_=1, psb_ik_rni_=2, psb_ik_ani_=3 integer(psb_ipk_), parameter :: psb_ik_xni_=4, psb_ik_bn2_=5, psb_ik_r0n2_=6 integer(psb_ipk_), parameter :: psb_ik_xn2_=7, psb_ik_errnum_=8, psb_ik_errden_=9 - integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11 + integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11, psb_ik_rn2_abs_=12 integer(psb_ipk_), parameter :: psb_ik_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3 integer(psb_ipk_), parameter :: psb_ik_ivsz_=16 type psb_itconv_type @@ -52,8 +62,29 @@ Module psb_base_linsolve_conv_mod real(psb_dpk_) :: values(psb_ik_ivsz_) end type psb_itconv_type + integer(psb_ipk_), save :: psb_istop_default = psb_istop_bn2_ + contains + function psb_is_valid_istop(istop) result(res) + integer(psb_ipk_) :: istop + logical :: res + + res = ((psb_istop_min_<=istop).and.(istop<=psb_istop_max_)) + end function psb_is_valid_istop + + function psb_get_istop_default() result(res) + integer(psb_ipk_) :: res + + res = psb_istop_default + end function psb_get_istop_default + + subroutine psb_set_istop_default(val) + integer(psb_ipk_) :: val + if ((psb_istop_min_<=val).and.(val<=psb_istop_max_)) & + & psb_istop_default = val + end subroutine psb_set_istop_default + subroutine log_header(methdname) !use psb_base_mod implicit none diff --git a/util/psb_c_mmio_impl.f90 b/util/psb_c_mmio_impl.f90 index 68690d06c..ad653462a 100644 --- a/util/psb_c_mmio_impl.f90 +++ b/util/psb_c_mmio_impl.f90 @@ -36,12 +36,13 @@ subroutine mm_cvet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_spk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -52,6 +53,7 @@ subroutine mm_cvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -76,15 +78,15 @@ subroutine mm_cvet_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow),stat = ircode) if (ircode /= 0) goto 993 - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i) = cmplx(bre,bim,kind=psb_spk_) + do i=1,nrow + read(infile,fmt=*,end=902) b(i) end do + end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,12 +111,13 @@ subroutine mm_cvet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_spk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -125,6 +128,7 @@ subroutine mm_cvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -152,15 +156,10 @@ subroutine mm_cvet2_read(b, info, iunit, filename) if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow,ncol),stat = ircode) if (ircode /= 0) goto 993 - do j=1, ncol - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i,j) = cmplx(bre,bim,kind=psb_spk_) - end do - end do + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -189,8 +188,10 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -201,6 +202,7 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -209,17 +211,19 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -241,8 +245,10 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -253,6 +259,7 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -262,20 +269,20 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) endif endif - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = 1 write(outfile,*) nrow,ncol - write(frmtv,'(a,i0,a)') '(',2*ncol,'(es26.18,1x))' + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' do i=1,size(b,1) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -298,7 +305,7 @@ subroutine mm_cvect_read(b, info, iunit, filename) complex(psb_spk_), allocatable :: bv(:) call mm_array_read(bv, info, iunit, filename) - call b%bld(bv) + if (info == 0) call b%bld(bv) end subroutine mm_cvect_read @@ -331,9 +338,10 @@ subroutine cmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_c_coo_sparse_mat), allocatable :: acoo - real(psb_spk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -344,6 +352,7 @@ subroutine cmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -369,24 +378,27 @@ subroutine cmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = cone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -398,37 +410,34 @@ subroutine cmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = cone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -456,10 +465,11 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -470,6 +480,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -481,7 +492,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -492,6 +503,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) return end subroutine cmm_mat_write + subroutine lcmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none @@ -501,12 +513,13 @@ subroutine lcmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_lc_coo_sparse_mat), allocatable :: acoo - real(psb_spk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -517,6 +530,7 @@ subroutine lcmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -542,24 +556,27 @@ subroutine lcmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = cone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -571,37 +588,34 @@ subroutine lcmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = cone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -629,10 +643,10 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -643,6 +657,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -654,7 +669,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return diff --git a/util/psb_d_mmio_impl.f90 b/util/psb_d_mmio_impl.f90 index 374eeed92..4fc917a61 100644 --- a/util/psb_d_mmio_impl.f90 +++ b/util/psb_d_mmio_impl.f90 @@ -39,8 +39,10 @@ subroutine mm_dvet_read(b, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -51,6 +53,7 @@ subroutine mm_dvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -83,7 +86,7 @@ subroutine mm_dvet_read(b, info, iunit, filename) end do end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -111,8 +114,10 @@ subroutine mm_dvet2_read(b, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -123,6 +128,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -153,7 +159,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -182,8 +188,10 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -194,6 +202,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -202,7 +211,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' @@ -210,9 +219,11 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) ncol = size(b,2) write(outfile,*) nrow, ncol - write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol) + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -234,8 +245,10 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -246,6 +259,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -268,7 +282,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -324,9 +338,10 @@ subroutine dmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_d_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -337,6 +352,7 @@ subroutine dmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -421,7 +437,7 @@ subroutine dmm_mat_read(a, info, iunit, filename) call a%mv_from(acoo) end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -449,10 +465,11 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -463,6 +480,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -474,7 +492,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -498,9 +516,10 @@ subroutine ldmm_mat_read(a, info, iunit, filename) integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr integer(psb_ipk_) :: ircode, infile type(psb_ld_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -511,6 +530,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -595,7 +615,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename) call a%mv_from(acoo) end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -623,10 +643,10 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -637,6 +657,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -648,7 +669,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -658,5 +679,3 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) write(psb_err_unit,*) 'Error while opening ',filename return end subroutine ldmm_mat_write - - diff --git a/util/psb_s_mmio_impl.f90 b/util/psb_s_mmio_impl.f90 index 94e886849..0ffd9a93d 100644 --- a/util/psb_s_mmio_impl.f90 +++ b/util/psb_s_mmio_impl.f90 @@ -36,11 +36,13 @@ subroutine mm_svet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -51,6 +53,7 @@ subroutine mm_svet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -83,8 +86,7 @@ subroutine mm_svet_read(b, info, iunit, filename) end do end if ! read right hand sides - - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,11 +111,13 @@ subroutine mm_svet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -124,6 +128,7 @@ subroutine mm_svet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -154,8 +159,7 @@ subroutine mm_svet2_read(b, info, iunit, filename) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -184,8 +188,10 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -196,6 +202,7 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -204,17 +211,19 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) outfile=6 endif endif - + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -236,8 +245,10 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -248,6 +259,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -270,7 +282,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -326,9 +338,10 @@ subroutine smm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_s_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -339,6 +352,7 @@ subroutine smm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -370,9 +384,14 @@ subroutine smm_mat_read(a, info, iunit, filename) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + call acoo%set_nzeros(nnzero) else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so @@ -391,17 +410,34 @@ subroutine smm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then + call acoo%allocate(nrow,ncol,2*nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + nzr = nnzero + do i=1,nnzero + if (acoo%ia(i) /= acoo%ja(i)) then + nzr = nzr + 1 + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) + end if + end do + call acoo%set_nzeros(nzr) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -429,10 +465,11 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -443,6 +480,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -454,7 +492,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -465,6 +503,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) return end subroutine smm_mat_write + subroutine lsmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none @@ -474,12 +513,13 @@ subroutine lsmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_ls_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -490,6 +530,7 @@ subroutine lsmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -521,9 +562,14 @@ subroutine lsmm_mat_read(a, info, iunit, filename) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + call acoo%set_nzeros(nnzero) else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so @@ -542,17 +588,34 @@ subroutine lsmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then + call acoo%allocate(nrow,ncol,2*nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + nzr = nnzero + do i=1,nnzero + if (acoo%ia(i) /= acoo%ja(i)) then + nzr = nzr + 1 + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) + end if + end do + call acoo%set_nzeros(nzr) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -580,10 +643,10 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -594,6 +657,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -605,7 +669,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return diff --git a/util/psb_z_mmio_impl.f90 b/util/psb_z_mmio_impl.f90 index 948de2833..aa66f0edd 100644 --- a/util/psb_z_mmio_impl.f90 +++ b/util/psb_z_mmio_impl.f90 @@ -36,12 +36,13 @@ subroutine mm_zvet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_dpk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -52,6 +53,7 @@ subroutine mm_zvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -76,15 +78,15 @@ subroutine mm_zvet_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow),stat = ircode) if (ircode /= 0) goto 993 - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i) = cmplx(bre,bim,kind=psb_dpk_) + do i=1,nrow + read(infile,fmt=*,end=902) b(i) end do + end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,12 +111,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_dpk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -125,6 +128,7 @@ subroutine mm_zvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -149,18 +153,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow,ncol),stat = ircode) if (ircode /= 0) goto 993 - do j=1, ncol - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i,j) = cmplx(bre,bim,kind=psb_dpk_) - end do - end do + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -189,8 +188,10 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -201,6 +202,7 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -209,17 +211,19 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -241,8 +245,10 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -253,6 +259,7 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -269,13 +276,13 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) ncol = 1 write(outfile,*) nrow,ncol - write(frmtv,'(a,i0,a)') '(',2*ncol,'(es26.18,1x))' + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' do i=1,size(b,1) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -331,9 +338,10 @@ subroutine zmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_z_coo_sparse_mat), allocatable :: acoo - real(psb_dpk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -344,6 +352,7 @@ subroutine zmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -369,24 +378,27 @@ subroutine zmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = zone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -398,37 +410,34 @@ subroutine zmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = zone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -446,6 +455,7 @@ subroutine zmm_mat_read(a, info, iunit, filename) return end subroutine zmm_mat_read + subroutine zmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none @@ -455,10 +465,11 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -469,6 +480,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -480,7 +492,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -501,12 +513,13 @@ subroutine lzmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_lz_coo_sparse_mat), allocatable :: acoo - real(psb_dpk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -517,6 +530,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -542,24 +556,27 @@ subroutine lzmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = zone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -571,37 +588,34 @@ subroutine lzmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = zone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -619,6 +633,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename) return end subroutine lzmm_mat_read + subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none @@ -628,10 +643,10 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -642,6 +657,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -653,7 +669,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -663,4 +679,3 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) write(psb_err_unit,*) 'Error while opening ',filename return end subroutine lzmm_mat_write -