From 182f53486b432454eaf7affb74885a4941bf4f3d Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 6 Jul 2006 12:14:56 +0000 Subject: [PATCH] Fixed bug in swaptran. Added psb_krylov interface. --- src/internals/psi_dswaptran.f90 | 158 +++++++++++++++-------------- src/internals/psi_iswaptran.f90 | 8 +- src/internals/psi_zswaptran.f90 | 8 +- src/methd/psb_dbicg.f90 | 5 +- src/modules/psb_methd_mod.f90 | 174 ++++++++++++++++++++++++++++++++ src/psblas/psb_dspmm.f90 | 126 ++++++++++++----------- src/serial/psb_dcsmm.f90 | 44 ++++++-- test/Fileread/df_sample.f90 | 23 ++--- test/Fileread/zf_sample.f90 | 9 +- test/pargen/ppde90.f90 | 17 +--- 10 files changed, 381 insertions(+), 191 deletions(-) diff --git a/src/internals/psi_dswaptran.f90 b/src/internals/psi_dswaptran.f90 index 60ed5e9a..9729f667 100644 --- a/src/internals/psi_dswaptran.f90 +++ b/src/internals/psi_dswaptran.f90 @@ -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 diff --git a/src/internals/psi_iswaptran.f90 b/src/internals/psi_iswaptran.f90 index 599149b4..0ffd2ce7 100644 --- a/src/internals/psi_iswaptran.f90 +++ b/src/internals/psi_iswaptran.f90 @@ -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 diff --git a/src/internals/psi_zswaptran.f90 b/src/internals/psi_zswaptran.f90 index a70e478a..0b7b6567 100644 --- a/src/internals/psi_zswaptran.f90 +++ b/src/internals/psi_zswaptran.f90 @@ -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 diff --git a/src/methd/psb_dbicg.f90 b/src/methd/psb_dbicg.f90 index 0d9983e6..b374c815 100644 --- a/src/methd/psb_dbicg.f90 +++ b/src/methd/psb_dbicg.f90 @@ -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) diff --git a/src/modules/psb_methd_mod.f90 b/src/modules/psb_methd_mod.f90 index 121a1170..4cdb62a0 100644 --- a/src/modules/psb_methd_mod.f90 +++ b/src/modules/psb_methd_mod.f90 @@ -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 diff --git a/src/psblas/psb_dspmm.f90 b/src/psblas/psb_dspmm.f90 index 341577d3..f7cc3f80 100644 --- a/src/psblas/psb_dspmm.f90 +++ b/src/psblas/psb_dspmm.f90 @@ -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 diff --git a/src/serial/psb_dcsmm.f90 b/src/serial/psb_dcsmm.f90 index f31f4394..5e255f69 100644 --- a/src/serial/psb_dcsmm.f90 +++ b/src/serial/psb_dcsmm.f90 @@ -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 diff --git a/test/Fileread/df_sample.f90 b/test/Fileread/df_sample.f90 index 211c6505..5ab45b06 100644 --- a/test/Fileread/df_sample.f90 +++ b/test/Fileread/df_sample.f90 @@ -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 diff --git a/test/Fileread/zf_sample.f90 b/test/Fileread/zf_sample.f90 index 05e6b9f6..6f9e7747 100644 --- a/test/Fileread/zf_sample.f90 +++ b/test/Fileread/zf_sample.f90 @@ -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) diff --git a/test/pargen/ppde90.f90 b/test/pargen/ppde90.f90 index 236faaf6..6093230a 100644 --- a/test/pargen/ppde90.f90 +++ b/test/pargen/ppde90.f90 @@ -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