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

Fixed small bugs in TRANS handling.

Browse Source
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent 55a22b9c73
commit 9bb64be14d
2 changed files with 258 additions and 256 deletions
Show Stats Download Patch File Download Diff File

327
src/psblas/psb_dspmm.f90
Unescape Escape View File

@ -119,14 +119,14 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow == -1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol /= 1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
endif
ia = 1
@ -134,45 +134,46 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
ix = 1
if (present(jx)) then
ijx = jx
ijx = jx
else
ijx = 1
ijx = 1
endif
iy = 1
if (present(jy)) then
ijy = jy
ijy = jy
else
ijy = 1
ijy = 1
endif
if (present(doswap)) then
idoswap = doswap
idoswap = doswap
else
idoswap = 1
idoswap = 1
endif
if (present(k)) then
ik = min(k,size(x,2)-ijx+1)
ik = min(ik,size(y,2)-ijy+1)
ik = min(k,size(x,2)-ijx+1)
ik = min(ik,size(y,2)-ijy+1)
else
ik = min(size(x,2)-ijx+1,size(y,2)-ijy+1)
ik = min(size(x,2)-ijx+1,size(y,2)-ijy+1)
endif
if (present(trans)) then
if((trans.eq.'N').or.(trans.eq.'T')) then
itrans = trans
else if (trans.eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
if ((trans.eq.'N').or.(trans.eq.'T')&
& .or.(trans.eq.'n').or.(trans.eq.'t')) then
itrans = trans
else if ((trans.eq.'C').or.(trans.eq.'c')) then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
else
itrans = 'N'
itrans = 'N'
endif
m = desc_a%matrix_data(psb_m_)
@ -187,168 +188,168 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
if (a%pr(1) /= 0) liwork = liwork + n * ik
if (a%pl(1) /= 0) liwork = liwork + m * ik
if (present(work)) then
if(size(work).lt.liwork) then
call psb_realloc(liwork,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
if(size(work).lt.liwork) then
call psb_realloc(liwork,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
end if
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
iwork(1)=0.d0
! checking for matrix correctness
call psb_chkmat(m,n,ia,ja,desc_a%matrix_data,info,iia,jja)
if(info.ne.0) then
info=4010
ch_err='psb_chkmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info=4010
ch_err='psb_chkmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (itrans.eq.'N') then
! Matrix is not transposed
if((ja.ne.ix).or.(ia.ne.iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(n,ik,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if(idoswap.lt.0) x(nrow:ncol,1:ik)=0.d0
ib1=min(nb,ik)
xp => x(iix:lldx,jjx:jjx+ib1-1)
if(idoswap.gt.0)&
& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& ib1,dzero,xp,desc_a,iwork,info)
! Matrix is not transposed
if((ja.ne.ix).or.(ia.ne.iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(n,ik,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if(idoswap.lt.0) x(nrow:ncol,1:ik)=0.d0
ib1=min(nb,ik)
xp => x(iix:lldx,jjx:jjx+ib1-1)
if(idoswap.gt.0)&
& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& ib1,dzero,xp,desc_a,iwork,info)
!!$ & call PSI_dSwapData(ior(SWAP_SEND,SWAP_RECV),ib1,&
!!$ & dzero,x(iix,jjx),lldx,desc_a%matrix_data,&
!!$ & desc_a%halo_index,iwork,liwork,info)
blk: do i=1, ik, nb
ib=ib1
ib1 = max(0,min(nb,(ik)-(i-1+ib)))
xp => x(iix:lldx,jjx+i+ib-1:jjx+i+ib+ib1-2)
if((ib1.gt.0).and.(idoswap.gt.0))&
& call psi_swapdata(psb_swap_send_,ib1,&
& dzero,xp,desc_a,iwork,info)
blk: do i=1, ik, nb
ib=ib1
ib1 = max(0,min(nb,(ik)-(i-1+ib)))
xp => x(iix:lldx,jjx+i+ib-1:jjx+i+ib+ib1-2)
if((ib1.gt.0).and.(idoswap.gt.0))&
& call psi_swapdata(psb_swap_send_,ib1,&
& dzero,xp,desc_a,iwork,info)
!!$ & call PSI_dSwapData(SWAP_SEND,ib1,&
!!$ & dzero,x(iix,jjx+i+ib-1),lldx,desc_a%matrix_data,&
!!$ & desc_a%halo_index,iwork,liwork,info)
if(info.ne.0) exit blk
! local Matrix-vector product
call dcsmm(itrans,nrow,ib,ncol,alpha,a%pr,a%fida,&
& a%descra,a%aspk,a%ia1,a%ia2,a%infoa,a%pl,&
& x(iix,jjx+i-1),lldx,beta,y(iiy,jjy+i-1),lldy,&
& iwork,liwork,info)
if(info.ne.0) exit blk
if((ib1.gt.0).and.(idoswap.gt.0))&
& call psi_swapdata(psb_swap_send_,ib1,&
& dzero,xp,desc_a,iwork,info)
if(info.ne.0) exit blk
! local Matrix-vector product
call dcsmm(itrans,nrow,ib,ncol,alpha,a%pr,a%fida,&
& a%descra,a%aspk,a%ia1,a%ia2,a%infoa,a%pl,&
& x(iix,jjx+i-1),lldx,beta,y(iiy,jjy+i-1),lldy,&
& iwork,liwork,info)
if(info.ne.0) exit blk
if((ib1.gt.0).and.(idoswap.gt.0))&
& call psi_swapdata(psb_swap_send_,ib1,&
& dzero,xp,desc_a,iwork,info)
!!$ & call PSI_dSwapData(SWAP_RECV,ib1,&
!!$ & dzero,x(iix,jjx+i+ib-1),lldx,desc_a%matrix_data,&
!!$ & desc_a%halo_index,iwork,liwork,info)
if(info.ne.0) exit blk
end do blk
if(info.ne.0) exit blk
end do blk
if(info.ne.0) then
info = 4011
call psb_errpush(info,name)
goto 9999
end if
if(info.ne.0) then
info = 4011
call psb_errpush(info,name)
goto 9999
end if
else
! Matrix is transposed
if((ja.ne.iy).or.(ia.ne.ix)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
if(desc_a%ovrlap_elem(1).ne.-1) then
info = 3070
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(n,ik,size(y,1),iy,ijy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if(idoswap.lt.0) y(nrow:ncol,1:ik)=0.d0
! local Matrix-vector product
call dcsmm(itrans,ncol,ik,nrow,alpha,a%pr,a%fida,&
& a%descra,a%aspk,a%ia1,a%ia2,a%infoa,a%pl,&
& x(iix,jjx),lldx,beta,y(iiy,jjy),lldy,&
& iwork,liwork,info)
if(info.ne.0) then
info = 4010
ch_err='dcsmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
yp => y(iiy:lldy,jjy:jjy+ik-1)
if(idoswap.gt.0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& ik,done,yp,desc_a,iwork,info)
! Matrix is transposed
if((ja.ne.iy).or.(ia.ne.ix)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
if(desc_a%ovrlap_elem(1).ne.-1) then
info = 3070
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(n,ik,size(y,1),iy,ijy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if(idoswap.lt.0) y(nrow:ncol,1:ik)=0.d0
! local Matrix-vector product
call dcsmm(itrans,ncol,ik,nrow,alpha,a%pr,a%fida,&
& a%descra,a%aspk,a%ia1,a%ia2,a%infoa,a%pl,&
& x(iix,jjx),lldx,beta,y(iiy,jjy),lldy,&
& iwork,liwork,info)
if(info.ne.0) then
info = 4010
ch_err='dcsmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
yp => y(iiy:lldy,jjy:jjy+ik-1)
if(idoswap.gt.0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& ik,done,yp,desc_a,iwork,info)
!!$ & call PSI_dSwapTran(ior(SWAP_SEND,SWAP_RECV),&
!!$ & ik,done,y(iiy,jjy),lldy,desc_a%matrix_data,&
!!$ & desc_a%halo_index,iwork,liwork,info)
if(info.ne.0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(info.ne.0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
@ -362,8 +363,8 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return
end subroutine psb_dspmm

187
src/psblas/psb_dspsm.f90
Unescape Escape View File

@ -117,14 +117,14 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow == -1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol /= 1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
endif
! just this case right now
@ -133,23 +133,23 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
ix = 1
if (present(jx)) then
ijx = jx
ijx = jx
else
ijx = 1
ijx = 1
endif
iy = 1
if (present(jy)) then
ijy = jy
ijy = jy
else
ijy = 1
ijy = 1
endif
if (present(k)) then
ik = min(k,size(x,2)-ijx+1)
ik = min(ik,size(y,2)-ijy+1)
ik = min(k,size(x,2)-ijx+1)
ik = min(ik,size(y,2)-ijy+1)
else
ik = min(size(x,2)-ijx+1,size(y,2)-ijy+1)
ik = min(size(x,2)-ijx+1,size(y,2)-ijy+1)
endif
if (present(choice)) then
@ -165,19 +165,20 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
endif
if (present(trans)) then
if((trans.eq.'N').or.(trans.eq.'T')) then
itrans = trans
else if (trans.eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
if ((trans.eq.'N').or.(trans.eq.'T')&
& .or.(trans.eq.'n').or.(trans.eq.'t')) then
itrans = trans
else if ((trans.eq.'C').or.(trans.eq.'c')) then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
else
itrans = 'N'
itrans = 'N'
endif
m = desc_a%matrix_data(psb_m_)
@ -187,9 +188,9 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
lldy = size(y,1)
if((lldx.lt.ncol).or.(lldy.lt.ncol)) then
info=3010
call psb_errpush(info,name)
goto 9999
info=3010
call psb_errpush(info,name)
goto 9999
end if
! check for presence/size of a work area
@ -197,34 +198,34 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
if (a%pr(1) /= 0) llwork = liwork + m * ik
if (a%pl(1) /= 0) llwork = llwork + m * ik
if (present(work)) then
if(size(work).lt.liwork) then
call psb_realloc(liwork,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
if(size(work).lt.liwork) then
call psb_realloc(liwork,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
end if
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
iwork(1)=0.d0
if(present(d)) then
lld = size(d)
id => d
lld = size(d)
id => d
else
lld=1
allocate(id(1))
id=1.d0
lld=1
allocate(id(1))
id=1.d0
end if
! checking for matrix correctness
@ -233,25 +234,25 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect/mat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info=4010
ch_err='psb_chkvect/mat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(ja.ne.ix) then
! this case is not yet implemented
info = 3030
! this case is not yet implemented
info = 3030
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
! this case is not yet implemented
info = 3040
end if
if(info.ne.0) then
call psb_errpush(info,name)
goto 9999
call psb_errpush(info,name)
goto 9999
end if
! Perform local triangular system solve
@ -260,47 +261,47 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
& a%pl,x(iix,jjx),lldx,beta,y(iiy,jjy),lldy,&
& iwork,liwork,info)
if(info.ne.0) then
info = 4010
ch_err='dcssm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info = 4010
ch_err='dcssm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! update overlap elements
if(lchoice.gt.0) then
yp => y(iiy:lldy,jjy:jjy+ik-1)
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),ik,&
& done,yp,desc_a,iwork,info)
yp => y(iiy:lldy,jjy:jjy+ik-1)
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),ik,&
& done,yp,desc_a,iwork,info)
!!$ call PSI_dSwapData(ior(SWAP_SEND,SWAP_RECV),ik,&
!!$ & done,y,lldy,desc_a%matrix_data,desc_a%ovrlap_index,&
!!$ & iwork,liwork,info)
i=0
! switch on update type
select case (lchoice)
case(psb_square_root_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:)/&
& sqrt(real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_)))
i = i+2
end do
case(psb_avg_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:)/&
& real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_))
i = i+2
end do
case(psb_sum_)
! do nothing
case default
! wrong value for choice argument
info = 70
int_err=(/10,lchoice,0,0,0/)
call psb_errpush(info,name,i_err=int_err)
goto 9999
end select
i=0
! switch on update type
select case (lchoice)
case(psb_square_root_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:)/&
& sqrt(real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_)))
i = i+2
end do
case(psb_avg_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_),:)/&
& real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_))
i = i+2
end do
case(psb_sum_)
! do nothing
case default
! wrong value for choice argument
info = 70
int_err=(/10,lchoice,0,0,0/)
call psb_errpush(info,name,i_err=int_err)
goto 9999
end select
end if
if(.not.present(work)) deallocate(iwork)
@ -313,8 +314,8 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return
end subroutine psb_dspsm

Write Preview
Loading…
Cancel
Save