Fixed bug in swaptran. Added psb_krylov interface.

psblas3-type-indexed
Salvatore Filippone 19 years ago
parent f861c89418
commit 182f53486b

@ -118,20 +118,20 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
allocate(sdsz(0:nprow-1), rvsz(0:nprow-1), bsdidx(0:nprow-1),&
& brvidx(0:nprow-1), rvhd(0:nprow-1), prcid(0:nprow-1),&
& ptp(0:nprow-1), stat=info)
if(info.ne.0) then
if(info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
swap_mpi = iand(flag,psb_swap_mpi_).ne.0
swap_sync = iand(flag,psb_swap_sync_).ne.0
swap_send = iand(flag,psb_swap_send_).ne.0
swap_recv = iand(flag,psb_swap_recv_).ne.0
swap_mpi = iand(flag,psb_swap_mpi_) /= 0
swap_sync = iand(flag,psb_swap_sync_) /= 0
swap_send = iand(flag,psb_swap_send_) /= 0
swap_recv = iand(flag,psb_swap_recv_) /= 0
if(present(data)) then
if(data.eq.psb_comm_halo_) then
if(data == psb_comm_halo_) then
d_idx => desc_a%halo_index
else if(data.eq.psb_comm_ovr_) then
else if(data == psb_comm_ovr_) then
d_idx => desc_a%ovrlap_index
else
d_idx => desc_a%halo_index
@ -140,16 +140,16 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
d_idx => desc_a%halo_index
end if
idxs = 0
idxr = 0
idxs = 1
idxr = 1
totxch = 0
point_to_proc = 1
rvhd(:) = mpi_request_null
! prepare info for communications
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm.ne.-1)
if(proc_to_comm .ne. myrow) totxch = totxch+1
do while (proc_to_comm /= -1)
if(proc_to_comm /= myrow) totxch = totxch+1
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -174,7 +174,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
all=.false.
else
allocate(sndbuf(idxs),rcvbuf(idxr), stat=info)
if(info.ne.0) then
if(info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
@ -187,7 +187,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
! gather elements into sendbuffer for swapping
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -205,7 +205,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
call mpi_alltoallv(rcvbuf,rvsz,brvidx,&
& mpi_double_precision,sndbuf,sdsz,&
& bsdidx,mpi_double_precision,icomm,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -215,7 +215,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
! scatter elements from receivebuffer after swapping
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -232,7 +232,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if (proc_to_comm .lt. myrow) then
@ -255,7 +255,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
call psi_gth(nerv,n,d_idx(idx_pt:idx_pt+nerv-1),&
& y,rcvbuf(rcv_pt:rcv_pt+nerv*n-1))
call dgesd2d(ictxt,nerv,n,rcvbuf(rcv_pt),nerv,proc_to_comm,0)
else if (proc_to_comm .eq. myrow) then
else if (proc_to_comm == myrow) then
! I send to myself
idx_pt = point_to_proc+psb_elem_recv_
rcv_pt = bsdidx(proc_to_comm)
@ -269,11 +269,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
idx_pt = point_to_proc+nerv+psb_elem_send_
snd_pt = bsdidx(proc_to_comm)
call psi_sct(nesd,n,d_idx(idx_pt:idx_pt+nesd-1),&
@ -294,17 +294,17 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
! First I post all the non blocking receives
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
p2ptag = krecvid(ictxt,proc_to_comm,myrow)
snd_pt = brvidx(proc_to_comm)
call mpi_irecv(sndbuf(rcv_pt),sdsz(proc_to_comm),&
& mpi_double_precision,prcid(proc_to_comm),&
& p2ptag, icomm,rvhd(proc_to_comm),iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -319,7 +319,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
! Then I post all the blocking sends
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -328,12 +328,12 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
call psi_gth(nerv,n,d_idx(idx_pt:idx_pt+nerv-1),&
& y,rcvbuf(rcv_pt:rcv_pt+nerv*n-1))
if(proc_to_comm .ne. myrow) then
if(proc_to_comm /= myrow) then
p2ptag=ksendid(ictxt,proc_to_comm,myrow)
call mpi_send(rcvbuf(rcv_pt),rvsz(proc_to_comm),&
& mpi_double_precision,prcid(proc_to_comm),&
& p2ptag,icomm,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -347,14 +347,14 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
if(.false.) then
do i=1, totxch
call mpi_waitany(nprow,rvhd,ixrec,p2pstat,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (ixrec .ne. mpi_undefined) then
if (ixrec /= mpi_undefined) then
ixrec=ixrec-1 ! mpi_waitany returns an 1 to nprow index
point_to_proc = ptp(ixrec)
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
@ -375,11 +375,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm .eq. myrow) then
if(proc_to_comm == myrow) then
idx_pt = point_to_proc+nerv+psb_elem_send_
rcv_pt = brvidx(proc_to_comm)
call psi_sct(nesd,n,d_idx(idx_pt:idx_pt+nesd-1),&
@ -393,13 +393,13 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
call mpi_wait(rvhd(proc_to_comm),p2pstat,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -428,7 +428,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -446,11 +446,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
snd_pt = bsdidx(proc_to_comm)
call dgerv2d(ictxt,nesd,n,sndbuf(snd_pt),nesd,proc_to_comm,0)
idx_pt = point_to_proc+nerv+psb_elem_send_
@ -473,7 +473,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
& brvidx,rvhd,prcid,&
& ptp,stat=info)
if(all) deallocate(sndbuf,rcvbuf,stat=info)
if(info.ne.0) then
if(info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
@ -483,7 +483,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
if (err_act == act_abort) then
call psb_error(ictxt)
return
end if
@ -616,20 +616,20 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
allocate(sdsz(0:nprow-1), rvsz(0:nprow-1), bsdidx(0:nprow-1),&
& brvidx(0:nprow-1), rvhd(0:nprow-1), prcid(0:nprow-1),&
& ptp(0:nprow-1), stat=info)
if(info.ne.0) then
if(info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
swap_mpi = iand(flag,psb_swap_mpi_).ne.0
swap_sync = iand(flag,psb_swap_sync_).ne.0
swap_send = iand(flag,psb_swap_send_).ne.0
swap_recv = iand(flag,psb_swap_recv_).ne.0
swap_mpi = iand(flag,psb_swap_mpi_) /= 0
swap_sync = iand(flag,psb_swap_sync_) /= 0
swap_send = iand(flag,psb_swap_send_) /= 0
swap_recv = iand(flag,psb_swap_recv_) /= 0
if(present(data)) then
if(data.eq.psb_comm_halo_) then
if(data == psb_comm_halo_) then
d_idx => desc_a%halo_index
else if(data.eq.psb_comm_ovr_) then
else if(data == psb_comm_ovr_) then
d_idx => desc_a%ovrlap_index
else
d_idx => desc_a%halo_index
@ -638,8 +638,8 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
d_idx => desc_a%halo_index
end if
idxs = 0
idxr = 0
idxs = 1
idxr = 1
totxch = 0
point_to_proc = 1
rvhd(:) = mpi_request_null
@ -647,8 +647,8 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
! prepare info for communications
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm.ne.-1)
if(proc_to_comm .ne. myrow) totxch = totxch+1
do while (proc_to_comm /= -1)
if(proc_to_comm /= myrow) totxch = totxch+1
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -673,7 +673,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
all=.false.
else
allocate(sndbuf(idxs),rcvbuf(idxr), stat=info)
if(info.ne.0) then
if(info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
@ -686,7 +686,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
! gather elements into sendbuffer for swapping
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -703,7 +703,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
call mpi_alltoallv(rcvbuf,rvsz,brvidx,&
& mpi_double_precision,sndbuf,sdsz,&
& bsdidx,mpi_double_precision,icomm,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -713,7 +713,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
! scatter elements from receivebuffer after swapping
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -730,7 +730,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if (proc_to_comm .lt. myrow) then
@ -753,7 +753,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
call psi_gth(nerv,d_idx(idx_pt:idx_pt+nerv-1),&
& y,rcvbuf(rcv_pt:rcv_pt+nerv-1))
call dgesd2d(ictxt,nerv,1,rcvbuf(rcv_pt),nerv,proc_to_comm,0)
else if (proc_to_comm .eq. myrow) then
else if (proc_to_comm == myrow) then
! I send to myself
idx_pt = point_to_proc+psb_elem_recv_
rcv_pt = bsdidx(proc_to_comm)
@ -767,11 +767,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
idx_pt = point_to_proc+nerv+psb_elem_send_
snd_pt = bsdidx(proc_to_comm)
call psi_sct(nesd,d_idx(idx_pt:idx_pt+nesd-1),&
@ -792,17 +792,18 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
! First I post all the non blocking receives
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if (proc_to_comm /= myrow) then
p2ptag = krecvid(ictxt,proc_to_comm,myrow)
snd_pt = brvidx(proc_to_comm)
!!$ write(0,*) myrow,'Swaptran ',snd_pt,' to',proc_to_comm
call mpi_irecv(sndbuf(snd_pt),sdsz(proc_to_comm),&
& mpi_double_precision,prcid(proc_to_comm),&
& p2ptag, icomm,rvhd(proc_to_comm),iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -817,7 +818,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
! Then I post all the blocking sends
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -826,12 +827,12 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
call psi_gth(nerv,d_idx(idx_pt:idx_pt+nerv-1),&
& y,rcvbuf(rcv_pt:rcv_pt+nerv-1))
if(proc_to_comm .ne. myrow) then
if(proc_to_comm /= myrow) then
p2ptag=ksendid(ictxt,proc_to_comm,myrow)
call mpi_send(rcvbuf(rcv_pt),rvsz(proc_to_comm),&
& mpi_double_precision,prcid(proc_to_comm),&
& p2ptag,icomm,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -845,14 +846,14 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
if(.false.) then
do i=1, totxch
call mpi_waitany(nprow,rvhd,ixrec,p2pstat,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (ixrec .ne. mpi_undefined) then
if (ixrec /= mpi_undefined) then
ixrec=ixrec-1 ! mpi_waitany returns an 1 to nprow index
point_to_proc = ptp(ixrec)
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
@ -873,11 +874,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm .eq. myrow) then
if(proc_to_comm == myrow) then
idx_pt = point_to_proc+nerv+psb_elem_send_
rcv_pt = brvidx(proc_to_comm)
call psi_sct(nesd,d_idx(idx_pt:idx_pt+nesd-1),&
@ -892,13 +893,13 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
call mpi_wait(rvhd(proc_to_comm),p2pstat,iret)
if(iret.ne.mpi_success) then
if(iret /= mpi_success) then
int_err(1) = iret
info=400
call psb_errpush(info,name,i_err=int_err)
@ -926,7 +927,7 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
@ -944,11 +945,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
point_to_proc = 1
proc_to_comm = d_idx(point_to_proc+psb_proc_id_)
do while (proc_to_comm .ne. -1)
do while (proc_to_comm /= -1)
nerv = d_idx(point_to_proc+psb_n_elem_recv_)
nesd = d_idx(point_to_proc+nerv+psb_n_elem_send_)
if(proc_to_comm.ne.myrow) then
if(proc_to_comm /= myrow) then
snd_pt = bsdidx(proc_to_comm)
call dgerv2d(ictxt,nesd,1,sndbuf(snd_pt),nesd,proc_to_comm,0)
idx_pt = point_to_proc+psb_elem_recv_
@ -971,8 +972,14 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
deallocate(sdsz,rvsz,bsdidx,&
& brvidx,rvhd,prcid,&
& ptp,stat=info)
if(all) deallocate(sndbuf,rcvbuf,stat=info)
if(info.ne.0) then
if ((info==0).and. (all)) then
write(0,*) myrow, 'Deallocating sndbuf? '
deallocate(sndbuf,rcvbuf,stat=info)
end if
if(info /= 0) then
write(0,*) 'Error on deallocate ',info
call psb_errpush(4000,name)
goto 9999
end if
@ -982,9 +989,10 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
if (err_act == act_abort) then
call psb_error(ictxt)
return
end if
return
end subroutine psi_dswaptranv

@ -140,8 +140,8 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data)
d_idx => desc_a%halo_index
end if
idxs = 0
idxr = 0
idxs = 1
idxr = 1
totxch = 0
point_to_proc = 1
rvhd(:) = mpi_request_null
@ -635,8 +635,8 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data)
d_idx => desc_a%halo_index
end if
idxs = 0
idxr = 0
idxs = 1
idxr = 1
totxch = 0
point_to_proc = 1
rvhd(:) = mpi_request_null

@ -156,8 +156,8 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data)
d_idx => desc_a%halo_index
end if
idxs = 0
idxr = 0
idxs = 1
idxr = 1
totxch = 0
point_to_proc = 1
rvhd(:) = mpi_request_null
@ -670,8 +670,8 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data)
d_idx => desc_a%halo_index
end if
idxs = 0
idxr = 0
idxs = 1
idxr = 1
totxch = 0
point_to_proc = 1
rvhd(:) = mpi_request_null

@ -216,8 +216,9 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,&
if (itx.ge.litmax) exit restart
it = 0
call psb_geaxpby(done,b,dzero,r,desc_a,info)
call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
call psb_geaxpby(done,r,dzero,rt,desc_a,info)
if (info == 0) call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
if (debug) write(0,*) me,' Done spmm',info
if (info == 0) call psb_geaxpby(done,r,dzero,rt,desc_a,info)
if(info.ne.0) then
info=4011
call psb_errpush(info,name)

@ -30,6 +30,12 @@
!!$
Module psb_methd_mod
interface psb_krylov
module procedure psb_dkrylov, psb_zkrylov
end interface
interface psb_cg
subroutine psb_dcg(a,prec,b,x,eps,&
& desc_a,info,itmax,iter,err,itrace,istop)
@ -177,6 +183,174 @@ Module psb_methd_mod
real(kind(1.d0)), optional, intent(out) :: err
end subroutine psb_zcgs
end interface
contains
Subroutine psb_dkrylov(method,a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace,irst,istop)
use psb_serial_mod
use psb_descriptor_type
Use psb_prec_type
use psb_string_mod
!!$ parameters
character(len=*) :: method
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
type(psb_dprec_type), intent(in) :: prec
Real(Kind(1.d0)), Intent(in) :: b(:)
Real(Kind(1.d0)), Intent(inout) :: x(:)
Real(Kind(1.d0)), Intent(in) :: eps
integer, intent(out) :: info
Integer, Optional, Intent(in) :: itmax, itrace, irst,istop
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
integer :: itmax_, itrace_, irst_, istop_, iter_
real(kind(1.d0)) :: err_
if (present(itmax)) then
itmax_ = itmax
else
itmax_ = 1000
end if
if (present(itrace)) then
itrace_ = itrace
else
itrace_ = -1
end if
if (present(irst)) then
irst_ = irst
else
irst_ = 1
end if
if (present(istop)) then
istop_ = istop
else
istop_ = 1
end if
select case(toupper(method))
case('CG')
call psb_cg(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
case('CGS')
call psb_cgs(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
case('BICG')
call psb_bicg(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
case('BICGSTAB')
call psb_bicgstab(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
case('RGMRES')
call psb_rgmres(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,irst_,istop_)
case('BICGSTABL')
call psb_bicgstabl(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,irst_,istop_)
case default
call psb_bicgstab(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
end select
if (present(err)) then
err = err_
endif
if (present(iter)) then
iter = iter_
endif
end subroutine psb_dkrylov
Subroutine psb_zkrylov(method,a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace,irst,istop)
use psb_serial_mod
use psb_descriptor_type
Use psb_prec_type
use psb_string_mod
!!$ parameters
character(len=*) :: method
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
type(psb_zprec_type), intent(in) :: prec
complex(Kind(1.d0)), Intent(in) :: b(:)
complex(Kind(1.d0)), Intent(inout) :: x(:)
Real(Kind(1.d0)), Intent(in) :: eps
integer, intent(out) :: info
Integer, Optional, Intent(in) :: itmax, itrace, irst,istop
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
integer :: itmax_, itrace_, irst_, istop_, iter_
real(kind(1.d0)) :: err_
if (present(itmax)) then
itmax_ = itmax
else
itmax_ = 1000
end if
if (present(itrace)) then
itrace_ = itrace
else
itrace_ = -1
end if
if (present(irst)) then
irst_ = irst
else
irst_ = 1
end if
if (present(istop)) then
istop_ = istop
else
istop_ = 1
end if
select case(toupper(method))
!!$ case('CG')
!!$ call psb_cg(a,prec,b,x,eps,desc_a,info,&
!!$ &itmax_,iter_,err_,itrace_,istop_)
case('CGS')
call psb_cgs(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
!!$ case('BICG')
!!$ call psb_bicg(a,prec,b,x,eps,desc_a,info,&
!!$ &itmax_,iter_,err_,itrace_,istop_)
case('BICGSTAB')
call psb_bicgstab(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
!!$ case('RGMRES')
!!$ call psb_rgmres(a,prec,b,x,eps,desc_a,info,&
!!$ &itmax_,iter_,err_,itrace_,irst_,istop_)
!!$ case('BICGSTABL')
!!$ call psb_bicgstabl(a,prec,b,x,eps,desc_a,info,&
!!$ &itmax_,iter_,err_,itrace_,irst_,istop_)
case default
call psb_bicgstab(a,prec,b,x,eps,desc_a,info,&
&itmax_,iter_,err_,itrace_,istop_)
end select
if (present(err)) then
err = err_
endif
if (present(iter)) then
iter = iter_
endif
end subroutine psb_zkrylov
end module psb_methd_mod

@ -112,7 +112,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
logical :: aliw
name='psb_dspmm'
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
@ -162,9 +162,9 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
endif
if (present(trans)) then
if ( (toupper(trans).eq.'N').or.(toupper(trans).eq.'T')) then
if ( (toupper(trans) == 'N').or.(toupper(trans) == 'T')) then
itrans = toupper(trans)
else if (toupper(trans).eq.'C') then
else if (toupper(trans) == 'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
@ -200,7 +200,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
if (aliw) then
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
@ -214,7 +214,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
! checking for matrix correctness
call psb_chkmat(m,n,ia,ja,desc_a%matrix_data,info,iia,jja)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_chkmat'
call psb_errpush(info,name,a_err=ch_err)
@ -222,9 +222,9 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
end if
if (itrans.eq.'N') then
if (itrans == 'N') then
! Matrix is not transposed
if((ja.ne.ix).or.(ia.ne.iy)) then
if((ja /= ix).or.(ia /= iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
@ -234,14 +234,14 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
! 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
if(info /= 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
if((iix /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
@ -251,7 +251,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
ib1=min(nb,ik)
xp => x(iix:lldx,jjx:jjx+ib1-1)
if(idoswap.gt.0)&
if(idoswap > 0)&
& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& ib1,dzero,xp,desc_a,iwork,info)
@ -260,26 +260,26 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
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))&
if((ib1 > 0).and.(idoswap > 0))&
& call psi_swapdata(psb_swap_send_,ib1,&
& dzero,xp,desc_a,iwork,info)
if(info.ne.0) exit blk
if(info /= 0) exit blk
! local Matrix-vector product
call psb_csmm(alpha,a,x(iix:lldx,jjx:jjx+ib-1),&
& beta,y(iiy:lldy,jjy:jjy+ib-1),info,trans=itrans)
if(info.ne.0) exit blk
if(info /= 0) exit blk
if((ib1.gt.0).and.(idoswap.gt.0))&
if((ib1 > 0).and.(idoswap > 0))&
& call psi_swapdata(psb_swap_send_,ib1,&
& dzero,xp,desc_a,iwork,info)
if(info.ne.0) exit blk
if(info /= 0) exit blk
end do blk
if(info.ne.0) then
if(info /= 0) then
info = 4011
call psb_errpush(info,name)
goto 9999
@ -287,14 +287,14 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
else
! Matrix is transposed
if((ja.ne.iy).or.(ia.ne.ix)) then
if((ja /= iy).or.(ia /= 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
if(desc_a%ovrlap_elem(1) /= -1) then
info = 3070
call psb_errpush(info,name)
goto 9999
@ -303,14 +303,14 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
! 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
if(info /= 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
if((iix /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
@ -324,7 +324,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
call psb_csmm(alpha,a,x(iix:lldx,jjx:jjx+ik-1),&
& beta,y(iiy:lldy,jjy:jjy+ik-1),info,trans=itrans)
if(info.ne.0) then
if(info /= 0) then
info = 4010
ch_err='csmm'
call psb_errpush(info,name,a_err=ch_err)
@ -336,7 +336,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& ik,done,yp,desc_a,iwork,info)
if(info.ne.0) then
if(info /= 0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
@ -354,7 +354,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
if (err_act == act_abort) then
call psb_error(ictxt)
return
end if
@ -450,9 +450,10 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
character :: itrans
character(len=20) :: name, ch_err
logical :: aliw
logical, parameter :: debug=.false.
name='psb_dspmv'
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
@ -487,9 +488,9 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
endif
if (present(trans)) then
if ( (toupper(trans).eq.'N').or.(toupper(trans).eq.'T')) then
if ( (toupper(trans) == 'N').or.(toupper(trans) == 'T')) then
itrans = toupper(trans)
else if (toupper(trans).eq.'C') then
else if (toupper(trans) == 'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
@ -515,43 +516,43 @@ subroutine psb_dspmv(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
! write(0,*)'---->>>',work(1)
if (present(work)) then
if (present(work)) then
if (size(work) >= liwork) then
aliw =.false.
aliw =.false.
else
aliw=.true.
aliw=.true.
endif
else
aliw=.true.
aliw=.true.
end if
aliw=.true.
aliw=.true.
if (aliw) then
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
allocate(iwork(liwork),stat=info)
if(info /= 0) then
info=4010
ch_err='Allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
else
iwork => work
iwork => work
endif
if (debug) write(0,*) myrow,name,' Allocated work ', info
! checking for matrix correctness
call psb_chkmat(m,n,ia,ja,desc_a%matrix_data,info,iia,jja)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_chkmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (itrans.eq.'N') then
if (debug) write(0,*) myrow,name,' Checkmat ', info
if (itrans == 'N') then
! Matrix is not transposed
if((ja.ne.ix).or.(ia.ne.iy)) then
if((ja /= ix).or.(ia /= iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
@ -561,14 +562,14 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
if(info /= 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
if((iix /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
@ -585,7 +586,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
! local Matrix-vector product
call psb_csmm(alpha,a,x(iix:lldx),beta,y(iiy:lldy),info)
if(info.ne.0) then
if(info /= 0) then
info = 4011
call psb_errpush(info,name)
goto 9999
@ -593,14 +594,14 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
else
! Matrix is transposed
if((ja.ne.iy).or.(ia.ne.ix)) then
if((ja /= iy).or.(ia /= 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
if(desc_a%ovrlap_elem(1) /= -1) then
info = 3070
call psb_errpush(info,name)
goto 9999
@ -609,14 +610,14 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(n,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
if(info /= 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
if((iix /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
@ -627,11 +628,11 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
yp => y(iiy:lldy)
yp(nrow+1:ncol)=dzero
if (debug) write(0,*) myrow,name,' checkvect ', info
! local Matrix-vector product
call psb_csmm(alpha,a,xp,beta,yp,info,trans=itrans)
if(info.ne.0) then
if (debug) write(0,*) myrow,name,' csmm ', info
if(info /= 0) then
info = 4010
ch_err='dcsmm'
call psb_errpush(info,name,a_err=ch_err)
@ -641,7 +642,8 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
if(idoswap /= 0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& done,yp,desc_a,iwork,info)
if(info.ne.0) then
if (debug) write(0,*) myrow,name,' swaptran ', info
if(info /= 0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
@ -650,16 +652,26 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
end if
if(aliw) deallocate(iwork)
if (aliw) deallocate(iwork,stat=info)
if (debug) write(0,*) myrow,name,' deallocat ',aliw, info
if(info /= 0) then
info = 4010
ch_err='Deallocate iwork'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
nullify(iwork)
call psb_erractionrestore(err_act)
if (debug) call blacs_barrier(ictxt,'A')
if (debug) write(0,*) myrow,name,' Returning '
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
if (err_act == act_abort) then
call psb_error(ictxt)
return
end if

@ -67,22 +67,46 @@ subroutine psb_dcsmm(alpha,a,b,beta,c,info,trans)
lb = size(b,1)
lc = size(c,1)
iwsz = 2*m*n
allocate(work(iwsz))
allocate(work(iwsz),stat=info)
if (info /= 0) then
info = 4010
ch_err='allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call dcsmm(trans_,m,n,k,alpha,&
& a%pl,a%fida,a%descra,a%aspk,a%ia1,a%ia2,a%infoa,a%pr,&
& b,lb,beta,c,lc,work,iwsz,info)
deallocate(work)
call psb_erractionrestore(err_act)
if(info.ne.0) then
if (err_act.eq.act_abort) then
call psb_error()
return
end if
if (info.ne.0) then
info = 4010
ch_err='Serial csmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
goto 9999
end if
deallocate(work,stat=info)
if (info.ne.0) then
info = 4010
ch_err='Deallocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == act_abort) then
call psb_error()
return
end if
return
end subroutine psb_dcsmm

@ -155,6 +155,9 @@ program df_sample
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
call random_seed()
call random_number(b_col_glob(1:m_problem))
b_col_glob(1:m_problem) = 2.0d0 * b_col_glob(1:m_problem) - 1.0d0
endif
call gebs2d(ictxt,'a',b_col_glob(1:m_problem))
else
@ -272,22 +275,8 @@ program df_sample
iparm = 0
call blacs_barrier(ictxt,'all')
t1 = mpi_wtime()
if (cmethd.eq.'BICGSTAB') then
call psb_bicgstab(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace,istop=istopc)
else if (cmethd.eq.'BICG') then
call psb_bicg(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd.eq.'CGS') then
call psb_cgs(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd.eq.'CG') then
call psb_cg(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd.eq.'BICGSTABL') then
call psb_bicgstabl(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,ierr,itrace,ml)
endif
call psb_krylov(cmethd,a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace,istop=istopc,irst=ml)
call blacs_barrier(ictxt,'all')
t2 = mpi_wtime() - t1
call psb_amx(ictxt,t2)
@ -299,7 +288,7 @@ program df_sample
if (amroot) then
call psb_prec_descr(6,pre)
write(*,'("Matrix: ",a)')mtrx_file
write(*,'("Computed solution on ",i4," processors")')np
write(*,'("Computed solution on ",i8," processors")')np
write(*,'("Iterations to convergence: ",i6)')iter
write(*,'("Error indicator on exit: ",f7.2)')err
write(*,'("Time to buil prec. : ",es10.4)')tprec

@ -272,13 +272,8 @@ program zf_sample
iparm = 0
call blacs_barrier(ictxt,'all')
t1 = mpi_wtime()
if (cmethd.eq.'BICGSTAB') then
call psb_bicgstab(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace,istop=istopc)
else if (cmethd.eq.'CGS') then
call psb_cgs(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace)
endif
call psb_krylov(cmethd,a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace,istop=istopc,irst=ml)
call blacs_barrier(ictxt,'all')
t2 = mpi_wtime() - t1
call psb_amx(ictxt,t2)

@ -201,21 +201,8 @@ program pde90
call psb_barrier(ictxt)
t1 = mpi_wtime()
eps = 1.d-9
if (cmethd == 'BICGSTAB') then
call psb_bicgstab(a,pre,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd == 'CGS') then
call psb_cgs(a,pre,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd == 'CG') then
call psb_cg(a,pre,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd == 'BICGSTABL') then
call psb_bicgstabl(a,pre,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,ml)
else
write(0,*) 'unknown method ',cmethd
end if
call psb_krylov(cmethd,a,pre,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop=istopc,irst=ml)
if(info.ne.0) then
info=4010

Loading…
Cancel
Save