fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

New GMRES and MMIO_IMPL improved and templated

Browse Source
maint-3.9.0
sfilippone 22 hours ago
parent 9e1c7b775e
commit 2e40d962e2
9 changed files with 872 additions and 556 deletions
Show Stats Download Patch File Download Diff File

217
linsolve/impl/psb_crgmres.f90
Unescape Escape View File

@ -130,7 +130,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_c_vect_type) :: w, w1, xt type(psb_c_vect_type) :: w, w1, xt
real(psb_spk_) :: tmp real(psb_spk_) :: tmp
complex(psb_spk_) :: scal, gm, rti, rti1 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 & n_row, n_col, nl
integer(psb_lpk_) :: mglob integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. 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 if (present(istop)) then
istop_ = istop istop_ = istop
else else
istop_ = 2 istop_ = psb_get_istop_default()
endif 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_ info=psb_err_invalid_istop_
err=info err=info
call psb_errpush(info,name,i_err=(/istop_/)) call psb_errpush(info,name,i_err=(/istop_/))
goto 9999 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 if (present(itmax)) then
litmax = itmax itmax_ = itmax
else else
litmax = 1000 itmax_ = 1000
endif endif
if (present(itrace)) then 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() & w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info) ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info) bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info) bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(cone,b,czero,v(1),desc_a,info) call psb_geaxpby(cone,b,czero,v(1),desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
@ -273,7 +287,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
r0n2 = psb_genrm2(v(1),desc_a,info) r0n2 = psb_genrm2(v(1),desc_a,info)
endif end select
errnum = czero errnum = czero
errden = cone errden = cone
deps = eps deps = eps
@ -324,27 +339,32 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
! !
! check convergence ! check convergence
! !
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info) rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info) xni = psb_geamax(x,desc_a,info)
errnum = rni errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
else if (istop_ == 2) then case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = bn2 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) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) & 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 prec%apply(v(i),w1,desc_a,info)
call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux) 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) select case(istop_)
call psb_geaxpby(-h(k,i),v(k),cone,w,desc_a,info) case(psb_istop_ani_)
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
! !
! build x and then compute the residual and its infinity norm ! build x and then compute the residual and its infinity norm
! !
rst = rs rst = rs
call w1%set(czero) call psb_geaxpby(cone,x,czero,xt,desc_a,info)
call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rst,size(rst,1)) call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
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,b,czero,w1,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) call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info) 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 errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
! !
case(psb_istop_bn2_)
else if (istop_ == 2) then
! !
! compute the residual 2-norm as byproduct of the solution ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = bn2 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 ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (errnum <= eps*errden) then 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) call psb_geaxpby(cone,xt,czero,x,desc_a,info)
! = x = xt ! = x = xt
else if (istop_ == 2) then 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) !
! build x
! end select
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
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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 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 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
if (itx >= litmax) then 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 >= itmax_) then
if (itrace_ > 0) then if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) & if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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) 9999 call psb_error_handler(err_act)
return 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 end subroutine psb_crgmres_vect

217
linsolve/impl/psb_drgmres.f90
Unescape Escape View File

@ -130,7 +130,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_d_vect_type) :: w, w1, xt type(psb_d_vect_type) :: w, w1, xt
real(psb_dpk_) :: tmp real(psb_dpk_) :: tmp
real(psb_dpk_) :: scal, gm, rti, rti1 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 & n_row, n_col, nl
integer(psb_lpk_) :: mglob integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. 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 if (present(istop)) then
istop_ = istop istop_ = istop
else else
istop_ = 2 istop_ = psb_get_istop_default()
endif 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_ info=psb_err_invalid_istop_
err=info err=info
call psb_errpush(info,name,i_err=(/istop_/)) call psb_errpush(info,name,i_err=(/istop_/))
goto 9999 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 if (present(itmax)) then
litmax = itmax itmax_ = itmax
else else
litmax = 1000 itmax_ = 1000
endif endif
if (present(itrace)) then 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() & w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info) ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info) bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info) bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(done,b,dzero,v(1),desc_a,info) call psb_geaxpby(done,b,dzero,v(1),desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
@ -273,7 +287,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
r0n2 = psb_genrm2(v(1),desc_a,info) r0n2 = psb_genrm2(v(1),desc_a,info)
endif end select
errnum = dzero errnum = dzero
errden = done errden = done
deps = eps deps = eps
@ -324,27 +339,32 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
! !
! check convergence ! check convergence
! !
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info) rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info) xni = psb_geamax(x,desc_a,info)
errnum = rni errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
else if (istop_ == 2) then case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = bn2 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) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) & 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 prec%apply(v(i),w1,desc_a,info)
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux) 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) select case(istop_)
call psb_geaxpby(-h(k,i),v(k),done,w,desc_a,info) case(psb_istop_ani_)
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
! !
! build x and then compute the residual and its infinity norm ! build x and then compute the residual and its infinity norm
! !
rst = rs rst = rs
call w1%set(dzero) call psb_geaxpby(done,x,dzero,xt,desc_a,info)
call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rst,size(rst,1)) call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
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,b,dzero,w1,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) call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info) 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 errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
! !
case(psb_istop_bn2_)
else if (istop_ == 2) then
! !
! compute the residual 2-norm as byproduct of the solution ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = bn2 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 ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (errnum <= eps*errden) then 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) call psb_geaxpby(done,xt,dzero,x,desc_a,info)
! = x = xt ! = x = xt
else if (istop_ == 2) then 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) !
! build x
! end select
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
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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 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 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
if (itx >= litmax) then 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 >= itmax_) then
if (itrace_ > 0) then if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) & if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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) 9999 call psb_error_handler(err_act)
return 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 end subroutine psb_drgmres_vect

217
linsolve/impl/psb_srgmres.f90
Unescape Escape View File

@ -130,7 +130,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_s_vect_type) :: w, w1, xt type(psb_s_vect_type) :: w, w1, xt
real(psb_spk_) :: tmp real(psb_spk_) :: tmp
real(psb_spk_) :: scal, gm, rti, rti1 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 & n_row, n_col, nl
integer(psb_lpk_) :: mglob integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. 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 if (present(istop)) then
istop_ = istop istop_ = istop
else else
istop_ = 2 istop_ = psb_get_istop_default()
endif 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_ info=psb_err_invalid_istop_
err=info err=info
call psb_errpush(info,name,i_err=(/istop_/)) call psb_errpush(info,name,i_err=(/istop_/))
goto 9999 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 if (present(itmax)) then
litmax = itmax itmax_ = itmax
else else
litmax = 1000 itmax_ = 1000
endif endif
if (present(itrace)) then 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() & w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info) ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info) bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info) bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(sone,b,szero,v(1),desc_a,info) call psb_geaxpby(sone,b,szero,v(1),desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
@ -273,7 +287,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
r0n2 = psb_genrm2(v(1),desc_a,info) r0n2 = psb_genrm2(v(1),desc_a,info)
endif end select
errnum = szero errnum = szero
errden = sone errden = sone
deps = eps deps = eps
@ -324,27 +339,32 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
! !
! check convergence ! check convergence
! !
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info) rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info) xni = psb_geamax(x,desc_a,info)
errnum = rni errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
else if (istop_ == 2) then case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = bn2 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) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) & 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 prec%apply(v(i),w1,desc_a,info)
call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux) 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) select case(istop_)
call psb_geaxpby(-h(k,i),v(k),sone,w,desc_a,info) case(psb_istop_ani_)
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
! !
! build x and then compute the residual and its infinity norm ! build x and then compute the residual and its infinity norm
! !
rst = rs rst = rs
call w1%set(szero) call psb_geaxpby(sone,x,szero,xt,desc_a,info)
call strsm('l','u','n','n',i,1,sone,h,size(h,1),rst,size(rst,1)) call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
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,b,szero,w1,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) call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info) 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 errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
! !
case(psb_istop_bn2_)
else if (istop_ == 2) then
! !
! compute the residual 2-norm as byproduct of the solution ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = bn2 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 ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (errnum <= eps*errden) then 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) call psb_geaxpby(sone,xt,szero,x,desc_a,info)
! = x = xt ! = x = xt
else if (istop_ == 2) then 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) !
! build x
! end select
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
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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 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 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
if (itx >= litmax) then 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 >= itmax_) then
if (itrace_ > 0) then if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) & if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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) 9999 call psb_error_handler(err_act)
return 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 end subroutine psb_srgmres_vect

217
linsolve/impl/psb_zrgmres.f90
Unescape Escape View File

@ -130,7 +130,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_z_vect_type) :: w, w1, xt type(psb_z_vect_type) :: w, w1, xt
real(psb_dpk_) :: tmp real(psb_dpk_) :: tmp
complex(psb_dpk_) :: scal, gm, rti, rti1 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 & n_row, n_col, nl
integer(psb_lpk_) :: mglob integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. 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 if (present(istop)) then
istop_ = istop istop_ = istop
else else
istop_ = 2 istop_ = psb_get_istop_default()
endif 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_ info=psb_err_invalid_istop_
err=info err=info
call psb_errpush(info,name,i_err=(/istop_/)) call psb_errpush(info,name,i_err=(/istop_/))
goto 9999 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 if (present(itmax)) then
litmax = itmax itmax_ = itmax
else else
litmax = 1000 itmax_ = 1000
endif endif
if (present(itrace)) then 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() & w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info) ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info) bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info) bn2 = psb_genrm2(b,desc_a,info)
else if (istop_ == 3) then case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(zone,b,zzero,v(1),desc_a,info) call psb_geaxpby(zone,b,zzero,v(1),desc_a,info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
@ -273,7 +287,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999 goto 9999
end if end if
r0n2 = psb_genrm2(v(1),desc_a,info) r0n2 = psb_genrm2(v(1),desc_a,info)
endif end select
errnum = zzero errnum = zzero
errden = zone errden = zone
deps = eps deps = eps
@ -324,27 +339,32 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
! !
! check convergence ! check convergence
! !
if (istop_ == 1) then select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info) rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info) xni = psb_geamax(x,desc_a,info)
errnum = rni errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
else if (istop_ == 2) then case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = bn2 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) rni = psb_genrm2(v(1),desc_a,info)
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_non_ info=psb_err_from_subroutine_non_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps) & 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 prec%apply(v(i),w1,desc_a,info)
call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux) 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) select case(istop_)
call psb_geaxpby(-h(k,i),v(k),zone,w,desc_a,info) case(psb_istop_ani_)
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
! !
! build x and then compute the residual and its infinity norm ! build x and then compute the residual and its infinity norm
! !
rst = rs rst = rs
call w1%set(zzero) call psb_geaxpby(zone,x,zzero,xt,desc_a,info)
call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rst,size(rst,1)) call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
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,b,zzero,w1,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) call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info) 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 errnum = rni
errden = (ani*xni+bni) errden = (ani*xni+bni)
! !
case(psb_istop_bn2_)
else if (istop_ == 2) then
! !
! compute the residual 2-norm as byproduct of the solution ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = bn2 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 ! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem ! 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)) rni = abs(rs(i+1))
errnum = rni errnum = rni
errden = r0n2 errden = r0n2
endif end select
if (errnum <= eps*errden) then 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) call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
! = x = xt ! = x = xt
else if (istop_ == 2) then 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) !
! build x
! end select
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
if (itrace_ > 0) & if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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 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 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
if (itx >= litmax) then 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 >= itmax_) then
if (itrace_ > 0) then if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) & if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps) & 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) 9999 call psb_error_handler(err_act)
return 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 end subroutine psb_zrgmres_vect

33
linsolve/psb_base_linsolve_conv_mod.f90
Unescape Escape View File

@ -41,10 +41,20 @@ Module psb_base_linsolve_conv_mod
module procedure psb_d_end_conv module procedure psb_d_end_conv
end interface 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_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_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_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_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3
integer(psb_ipk_), parameter :: psb_ik_ivsz_=16 integer(psb_ipk_), parameter :: psb_ik_ivsz_=16
type psb_itconv_type type psb_itconv_type
@ -52,8 +62,29 @@ Module psb_base_linsolve_conv_mod
real(psb_dpk_) :: values(psb_ik_ivsz_) real(psb_dpk_) :: values(psb_ik_ivsz_)
end type psb_itconv_type end type psb_itconv_type
integer(psb_ipk_), save :: psb_istop_default = psb_istop_bn2_
contains 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) subroutine log_header(methdname)
!use psb_base_mod !use psb_base_mod
implicit none implicit none

169
util/psb_c_mmio_impl.f90
Unescape Escape View File

@ -36,12 +36,13 @@ subroutine mm_cvet_read(b, info, iunit, filename)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename 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
real(psb_spk_) :: bre, bim
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -52,6 +53,7 @@ subroutine mm_cvet_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -76,15 +78,15 @@ subroutine mm_cvet_read(b, info, iunit, filename)
read(line,fmt=*)nrow,ncol 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) allocate(b(nrow),stat = ircode)
if (ircode /= 0) goto 993 if (ircode /= 0) goto 993
do i=1, nrow do i=1,nrow
read(infile,fmt=*,end=902) bre,bim read(infile,fmt=*,end=902) b(i)
b(i) = cmplx(bre,bim,kind=psb_spk_)
end do end do
end if ! read right hand sides end if ! read right hand sides
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! open failed
@ -109,12 +111,13 @@ subroutine mm_cvet2_read(b, info, iunit, filename)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename 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
real(psb_spk_) :: bre, bim
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -125,6 +128,7 @@ subroutine mm_cvet2_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then 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 if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then
allocate(b(nrow,ncol),stat = ircode) allocate(b(nrow,ncol),stat = ircode)
if (ircode /= 0) goto 993 if (ircode /= 0) goto 993
do j=1, ncol read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),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
end if ! read right hand sides end if ! read right hand sides
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -201,6 +202,7 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -210,16 +212,18 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename)
endif endif
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)') '% '//trim(header)
write(outfile,'(a)') '% ' write(outfile,'(a)') '% '
nrow = size(b,1) nrow = size(b,1)
ncol = size(b,2) ncol = size(b,2)
write(outfile,*) nrow,ncol write(outfile,*) nrow, ncol
write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))'
if (outfile /= 6) close(outfile) write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol)
if (opened) close(outfile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -253,6 +259,7 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -262,20 +269,20 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename)
endif endif
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)') '% '//trim(header)
write(outfile,'(a)') '% ' write(outfile,'(a)') '% '
nrow = size(b,1) nrow = size(b,1)
ncol = 1 ncol = 1
write(outfile,*) nrow,ncol 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) do i=1,size(b,1)
write(outfile,frmtv) b(i) write(outfile,frmtv) b(i)
end do end do
if (outfile /= 6) close(outfile) if (opened) close(outfile)
return return
! open failed ! open failed
@ -298,7 +305,7 @@ subroutine mm_cvect_read(b, info, iunit, filename)
complex(psb_spk_), allocatable :: bv(:) complex(psb_spk_), allocatable :: bv(:)
call mm_array_read(bv, info, iunit, filename) call mm_array_read(bv, info, iunit, filename)
call b%bld(bv) if (info == 0) call b%bld(bv)
end subroutine mm_cvect_read 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_) :: nrow, ncol, nnzero
integer(psb_ipk_) :: ircode, i,nzr,infile integer(psb_ipk_) :: ircode, i,nzr,infile
type(psb_c_coo_sparse_mat), allocatable :: acoo type(psb_c_coo_sparse_mat), allocatable :: acoo
real(psb_spk_) :: are, aim logical :: opened
info = psb_success_
info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -344,6 +352,7 @@ subroutine cmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -369,24 +378,27 @@ subroutine cmm_mat_read(a, info, iunit, filename)
allocate(acoo, stat=ircode) allocate(acoo, stat=ircode)
if (ircode /= 0) goto 993 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) call acoo%allocate(nrow,ncol,nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
end do end do
call acoo%set_nzeros(nnzero) 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 are generally working with non-symmetric matrices, so
! we de-symmetrize what we are about to read ! we de-symmetrize what we are about to read
call acoo%allocate(nrow,ncol,2*nnzero) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
end do end do
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
@ -398,37 +410,34 @@ subroutine cmm_mat_read(a, info, iunit, filename)
end if end if
end do end do
call acoo%set_nzeros(nzr) 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 else if ((psb_tolower(type) == 'pattern').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) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i)
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
end do end do
acoo%val(:) = cone
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
if (acoo%ia(i) /= acoo%ja(i)) then if (acoo%ia(i) /= acoo%ja(i)) then
nzr = nzr + 1 nzr = nzr + 1
acoo%val(nzr) = conjg(acoo%val(i))
acoo%ia(nzr) = acoo%ja(i) acoo%ia(nzr) = acoo%ja(i)
acoo%ja(nzr) = acoo%ia(i) acoo%ja(nzr) = acoo%ia(i)
end if end if
end do end do
call acoo%set_nzeros(nzr) call acoo%set_nzeros(nzr)
call acoo%fix(info)
call a%mv_from(acoo)
else else
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
info=904 info=904
end if 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 return
! open failed ! open failed
@ -456,10 +465,11 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -470,6 +480,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -481,7 +492,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return
@ -492,6 +503,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename)
return return
end subroutine cmm_mat_write end subroutine cmm_mat_write
subroutine lcmm_mat_read(a, info, iunit, filename) subroutine lcmm_mat_read(a, info, iunit, filename)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -501,12 +513,13 @@ subroutine lcmm_mat_read(a, info, iunit, filename)
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character :: mmheader*15, fmt*15, object*10, type*10, sym*15
character(1024) :: line character(1024) :: line
integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
integer(psb_ipk_) :: ircode,infile integer(psb_ipk_) :: ircode, infile
type(psb_lc_coo_sparse_mat), allocatable :: acoo type(psb_lc_coo_sparse_mat), allocatable :: acoo
real(psb_spk_) :: are, aim logical :: opened
info = psb_success_
info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -517,6 +530,7 @@ subroutine lcmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -542,24 +556,27 @@ subroutine lcmm_mat_read(a, info, iunit, filename)
allocate(acoo, stat=ircode) allocate(acoo, stat=ircode)
if (ircode /= 0) goto 993 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) call acoo%allocate(nrow,ncol,nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
end do end do
call acoo%set_nzeros(nnzero) 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 are generally working with non-symmetric matrices, so
! we de-symmetrize what we are about to read ! we de-symmetrize what we are about to read
call acoo%allocate(nrow,ncol,2*nnzero) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
end do end do
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
@ -571,37 +588,34 @@ subroutine lcmm_mat_read(a, info, iunit, filename)
end if end if
end do end do
call acoo%set_nzeros(nzr) 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 else if ((psb_tolower(type) == 'pattern').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) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i)
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
end do end do
acoo%val(:) = cone
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
if (acoo%ia(i) /= acoo%ja(i)) then if (acoo%ia(i) /= acoo%ja(i)) then
nzr = nzr + 1 nzr = nzr + 1
acoo%val(nzr) = conjg(acoo%val(i))
acoo%ia(nzr) = acoo%ja(i) acoo%ia(nzr) = acoo%ja(i)
acoo%ja(nzr) = acoo%ia(i) acoo%ja(nzr) = acoo%ia(i)
end if end if
end do end do
call acoo%set_nzeros(nzr) call acoo%set_nzeros(nzr)
call acoo%fix(info)
call a%mv_from(acoo)
else else
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
info=904 info=904
end if 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 return
! open failed ! open failed
@ -629,10 +643,10 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -643,6 +657,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -654,7 +669,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return

51
util/psb_d_mmio_impl.f90
Unescape Escape View File

@ -39,8 +39,10 @@ subroutine mm_dvet_read(b, info, iunit, 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,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -51,6 +53,7 @@ subroutine mm_dvet_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -83,7 +86,7 @@ subroutine mm_dvet_read(b, info, iunit, filename)
end do end do
end if ! read right hand sides end if ! read right hand sides
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -123,6 +128,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then 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) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol)
end if ! read right hand sides end if ! read right hand sides
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -194,6 +202,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -210,9 +219,11 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
ncol = size(b,2) ncol = size(b,2)
write(outfile,*) nrow, ncol 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 return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -246,6 +259,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -268,7 +282,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
write(outfile,frmtv) b(i) write(outfile,frmtv) b(i)
end do end do
if (outfile /= 6) close(outfile) if (opened) close(outfile)
return return
! open failed ! open failed
@ -324,9 +338,10 @@ subroutine dmm_mat_read(a, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: nrow, ncol, nnzero
integer(psb_ipk_) :: ircode, i,nzr,infile integer(psb_ipk_) :: ircode, i,nzr,infile
type(psb_d_coo_sparse_mat), allocatable :: acoo type(psb_d_coo_sparse_mat), allocatable :: acoo
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -337,6 +352,7 @@ subroutine dmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -421,7 +437,7 @@ subroutine dmm_mat_read(a, info, iunit, filename)
call a%mv_from(acoo) call a%mv_from(acoo)
end if end if
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! open failed
@ -449,10 +465,11 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -463,6 +480,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -474,7 +492,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return
@ -498,9 +516,10 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
integer(psb_ipk_) :: ircode, infile integer(psb_ipk_) :: ircode, infile
type(psb_ld_coo_sparse_mat), allocatable :: acoo type(psb_ld_coo_sparse_mat), allocatable :: acoo
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -511,6 +530,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -595,7 +615,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
call a%mv_from(acoo) call a%mv_from(acoo)
end if end if
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! open failed
@ -623,10 +643,10 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -637,6 +657,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -648,7 +669,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return
@ -658,5 +679,3 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
write(psb_err_unit,*) 'Error while opening ',filename write(psb_err_unit,*) 'Error while opening ',filename
return return
end subroutine ldmm_mat_write end subroutine ldmm_mat_write

122
util/psb_s_mmio_impl.f90
Unescape Escape View File

@ -36,11 +36,13 @@ subroutine mm_svet_read(b, info, iunit, filename)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename 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,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -51,6 +53,7 @@ subroutine mm_svet_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -83,8 +86,7 @@ subroutine mm_svet_read(b, info, iunit, filename)
end do end do
end if ! read right hand sides end if ! read right hand sides
if (opened) close(infile)
if (infile /= 5) close(infile)
return return
! open failed ! open failed
@ -109,11 +111,13 @@ subroutine mm_svet2_read(b, info, iunit, filename)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename 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,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -124,6 +128,7 @@ subroutine mm_svet2_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then 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) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol)
end if ! read right hand sides end if ! read right hand sides
if (opened) close(infile)
if (infile /= 5) close(infile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -196,6 +202,7 @@ subroutine mm_svet2_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -210,11 +217,13 @@ subroutine mm_svet2_write(b, header, info, iunit, filename)
write(outfile,'(a)') '% ' write(outfile,'(a)') '% '
nrow = size(b,1) nrow = size(b,1)
ncol = size(b,2) 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 return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -248,6 +259,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -270,7 +282,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename)
write(outfile,frmtv) b(i) write(outfile,frmtv) b(i)
end do end do
if (outfile /= 6) close(outfile) if (opened) close(outfile)
return return
! open failed ! open failed
@ -326,9 +338,10 @@ subroutine smm_mat_read(a, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: nrow, ncol, nnzero
integer(psb_ipk_) :: ircode, i,nzr,infile integer(psb_ipk_) :: ircode, i,nzr,infile
type(psb_s_coo_sparse_mat), allocatable :: acoo type(psb_s_coo_sparse_mat), allocatable :: acoo
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -339,6 +352,7 @@ subroutine smm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then 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) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
end do end do
call acoo%set_nzeros(nnzero) 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 else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then
! we are generally working with non-symmetric matrices, so ! we are generally working with non-symmetric matrices, so
@ -391,17 +410,34 @@ subroutine smm_mat_read(a, info, iunit, filename)
end if end if
end do end do
call acoo%set_nzeros(nzr) 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 else
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
info=904 info=904
end if 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 return
! open failed ! open failed
@ -429,10 +465,11 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -443,6 +480,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -454,7 +492,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return
@ -465,6 +503,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename)
return return
end subroutine smm_mat_write end subroutine smm_mat_write
subroutine lsmm_mat_read(a, info, iunit, filename) subroutine lsmm_mat_read(a, info, iunit, filename)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -474,12 +513,13 @@ subroutine lsmm_mat_read(a, info, iunit, filename)
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character :: mmheader*15, fmt*15, object*10, type*10, sym*15
character(1024) :: line character(1024) :: line
integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
integer(psb_ipk_) :: ircode,infile integer(psb_ipk_) :: ircode, infile
type(psb_ls_coo_sparse_mat), allocatable :: acoo type(psb_ls_coo_sparse_mat), allocatable :: acoo
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -490,6 +530,7 @@ subroutine lsmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then 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) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
end do end do
call acoo%set_nzeros(nnzero) 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 else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then
! we are generally working with non-symmetric matrices, so ! we are generally working with non-symmetric matrices, so
@ -542,17 +588,34 @@ subroutine lsmm_mat_read(a, info, iunit, filename)
end if end if
end do end do
call acoo%set_nzeros(nzr) 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 else
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
info=904 info=904
end if 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 return
! open failed ! open failed
@ -580,10 +643,10 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -594,6 +657,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -605,7 +669,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return

169
util/psb_z_mmio_impl.f90
Unescape Escape View File

@ -36,12 +36,13 @@ subroutine mm_zvet_read(b, info, iunit, filename)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename 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
real(psb_dpk_) :: bre, bim
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -52,6 +53,7 @@ subroutine mm_zvet_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -76,15 +78,15 @@ subroutine mm_zvet_read(b, info, iunit, filename)
read(line,fmt=*)nrow,ncol 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) allocate(b(nrow),stat = ircode)
if (ircode /= 0) goto 993 if (ircode /= 0) goto 993
do i=1, nrow do i=1,nrow
read(infile,fmt=*,end=902) bre,bim read(infile,fmt=*,end=902) b(i)
b(i) = cmplx(bre,bim,kind=psb_dpk_)
end do end do
end if ! read right hand sides end if ! read right hand sides
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! open failed
@ -109,12 +111,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename 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
real(psb_dpk_) :: bre, bim
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024 & line*1024
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -125,6 +128,7 @@ subroutine mm_zvet2_read(b, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -149,18 +153,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename)
read(line,fmt=*)nrow,ncol 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) allocate(b(nrow,ncol),stat = ircode)
if (ircode /= 0) goto 993 if (ircode /= 0) goto 993
do j=1, ncol read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),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
end if ! read right hand sides end if ! read right hand sides
if (infile /= 5) close(infile) if (opened) close(infile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -201,6 +202,7 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -210,16 +212,18 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename)
endif endif
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)') '% '//trim(header)
write(outfile,'(a)') '% ' write(outfile,'(a)') '% '
nrow = size(b,1) nrow = size(b,1)
ncol = size(b,2) ncol = size(b,2)
write(outfile,*) nrow,ncol write(outfile,*) nrow, ncol
write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))'
if (outfile /= 6) close(outfile) write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol)
if (opened) close(outfile)
return return
! open failed ! 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 integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile
character(len=80) :: frmtv character(len=80) :: frmtv
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
outfile=6 outfile=6
@ -253,6 +259,7 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename)
outfile=99 outfile=99
endif endif
open(outfile,file=filename, err=901, action='WRITE') open(outfile,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -269,13 +276,13 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename)
ncol = 1 ncol = 1
write(outfile,*) nrow,ncol 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) do i=1,size(b,1)
write(outfile,frmtv) b(i) write(outfile,frmtv) b(i)
end do end do
if (outfile /= 6) close(outfile) if (opened) close(outfile)
return return
! open failed ! open failed
@ -331,9 +338,10 @@ subroutine zmm_mat_read(a, info, iunit, filename)
integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: nrow, ncol, nnzero
integer(psb_ipk_) :: ircode, i,nzr,infile integer(psb_ipk_) :: ircode, i,nzr,infile
type(psb_z_coo_sparse_mat), allocatable :: acoo type(psb_z_coo_sparse_mat), allocatable :: acoo
real(psb_dpk_) :: are, aim logical :: opened
info = psb_success_
info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -344,6 +352,7 @@ subroutine zmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -369,24 +378,27 @@ subroutine zmm_mat_read(a, info, iunit, filename)
allocate(acoo, stat=ircode) allocate(acoo, stat=ircode)
if (ircode /= 0) goto 993 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) call acoo%allocate(nrow,ncol,nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
end do end do
call acoo%set_nzeros(nnzero) 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 are generally working with non-symmetric matrices, so
! we de-symmetrize what we are about to read ! we de-symmetrize what we are about to read
call acoo%allocate(nrow,ncol,2*nnzero) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
end do end do
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
@ -398,37 +410,34 @@ subroutine zmm_mat_read(a, info, iunit, filename)
end if end if
end do end do
call acoo%set_nzeros(nzr) 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
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
call acoo%allocate(nrow,ncol,2*nnzero) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i)
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
end do end do
acoo%val(:) = zone
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
if (acoo%ia(i) /= acoo%ja(i)) then if (acoo%ia(i) /= acoo%ja(i)) then
nzr = nzr + 1 nzr = nzr + 1
acoo%val(nzr) = conjg(acoo%val(i))
acoo%ia(nzr) = acoo%ja(i) acoo%ia(nzr) = acoo%ja(i)
acoo%ja(nzr) = acoo%ia(i) acoo%ja(nzr) = acoo%ia(i)
end if end if
end do end do
call acoo%set_nzeros(nzr) call acoo%set_nzeros(nzr)
call acoo%fix(info)
call a%mv_from(acoo)
else else
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
info=904 info=904
end if 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 return
! open failed ! open failed
@ -446,6 +455,7 @@ subroutine zmm_mat_read(a, info, iunit, filename)
return return
end subroutine zmm_mat_read end subroutine zmm_mat_read
subroutine zmm_mat_write(a,mtitle,info,iunit,filename) subroutine zmm_mat_write(a,mtitle,info,iunit,filename)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -455,10 +465,11 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -469,6 +480,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -480,7 +492,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return
@ -501,12 +513,13 @@ subroutine lzmm_mat_read(a, info, iunit, filename)
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character :: mmheader*15, fmt*15, object*10, type*10, sym*15
character(1024) :: line character(1024) :: line
integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
integer(psb_ipk_) :: ircode,infile integer(psb_ipk_) :: ircode, infile
type(psb_lz_coo_sparse_mat), allocatable :: acoo type(psb_lz_coo_sparse_mat), allocatable :: acoo
real(psb_dpk_) :: are, aim logical :: opened
info = psb_success_
info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
infile=5 infile=5
@ -517,6 +530,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename)
infile=99 infile=99
endif endif
open(infile,file=filename, status='OLD', err=901, action='READ') open(infile,file=filename, status='OLD', err=901, action='READ')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -542,24 +556,27 @@ subroutine lzmm_mat_read(a, info, iunit, filename)
allocate(acoo, stat=ircode) allocate(acoo, stat=ircode)
if (ircode /= 0) goto 993 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) call acoo%allocate(nrow,ncol,nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
end do end do
call acoo%set_nzeros(nnzero) 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 are generally working with non-symmetric matrices, so
! we de-symmetrize what we are about to read ! we de-symmetrize what we are about to read
call acoo%allocate(nrow,ncol,2*nnzero) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
end do end do
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
@ -571,37 +588,34 @@ subroutine lzmm_mat_read(a, info, iunit, filename)
end if end if
end do end do
call acoo%set_nzeros(nzr) 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
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
call acoo%allocate(nrow,ncol,2*nnzero) call acoo%allocate(nrow,ncol,2*nnzero)
do i=1,nnzero do i=1,nnzero
read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i)
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
end do end do
acoo%val(:) = zone
nzr = nnzero nzr = nnzero
do i=1,nnzero do i=1,nnzero
if (acoo%ia(i) /= acoo%ja(i)) then if (acoo%ia(i) /= acoo%ja(i)) then
nzr = nzr + 1 nzr = nzr + 1
acoo%val(nzr) = conjg(acoo%val(i))
acoo%ia(nzr) = acoo%ja(i) acoo%ia(nzr) = acoo%ja(i)
acoo%ja(nzr) = acoo%ia(i) acoo%ja(nzr) = acoo%ia(i)
end if end if
end do end do
call acoo%set_nzeros(nzr) call acoo%set_nzeros(nzr)
call acoo%fix(info)
call a%mv_from(acoo)
else else
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
info=904 info=904
end if 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 return
! open failed ! open failed
@ -619,6 +633,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename)
return return
end subroutine lzmm_mat_read end subroutine lzmm_mat_read
subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -628,10 +643,10 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
integer(psb_ipk_), optional, intent(in) :: iunit integer(psb_ipk_), optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename character(len=*), optional, intent(in) :: filename
integer(psb_ipk_) :: iout integer(psb_ipk_) :: iout
logical :: opened
info = psb_success_ info = psb_success_
opened = .false.
if (present(filename)) then if (present(filename)) then
if (filename == '-') then if (filename == '-') then
iout=6 iout=6
@ -642,6 +657,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
iout=99 iout=99
endif endif
open(iout,file=filename, err=901, action='WRITE') open(iout,file=filename, err=901, action='WRITE')
opened = .true.
endif endif
else else
if (present(iunit)) then if (present(iunit)) then
@ -653,7 +669,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
call a%print(iout,head=mtitle) call a%print(iout,head=mtitle)
if (iout /= 6) close(iout) if (opened) close(iout)
return return
@ -663,4 +679,3 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
write(psb_err_unit,*) 'Error while opening ',filename write(psb_err_unit,*) 'Error while opening ',filename
return return
end subroutine lzmm_mat_write end subroutine lzmm_mat_write

Write Preview
Loading…
Cancel
Save