diff --git a/src/psblas/psb_dspmm.f90 b/src/psblas/psb_dspmm.f90 index 7e975095..d61ff0ef 100644 --- a/src/psblas/psb_dspmm.f90 +++ b/src/psblas/psb_dspmm.f90 @@ -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 diff --git a/src/psblas/psb_dspsm.f90 b/src/psblas/psb_dspsm.f90 index f3a34e5c..a00e1fbc 100644 --- a/src/psblas/psb_dspsm.f90 +++ b/src/psblas/psb_dspsm.f90 @@ -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