From f24d39ec18f630a318528425f0f5cde7e02e8da5 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 29 Jan 2008 09:06:17 +0000 Subject: [PATCH] psblas: base/modules/psb_tools_mod.f90 base/tools/psb_dallc.f90 base/tools/psb_ialloc.f90 base/tools/psb_zallc.f90 krylov/psb_dcgstabl.f90 Fixed stupid bug in BiCGSTAB(L): the column index of a distributed matrix might start from something else than 0. Changed allocation routines. --- base/modules/psb_tools_mod.f90 | 16 ++--- base/tools/psb_dallc.f90 | 97 +++++++++++++------------------ base/tools/psb_ialloc.f90 | 103 ++++++++++++++------------------- base/tools/psb_zallc.f90 | 97 +++++++++++++------------------ krylov/psb_dcgstabl.f90 | 6 +- 5 files changed, 132 insertions(+), 187 deletions(-) diff --git a/base/modules/psb_tools_mod.f90 b/base/modules/psb_tools_mod.f90 index b52cb7df..31357d3d 100644 --- a/base/modules/psb_tools_mod.f90 +++ b/base/modules/psb_tools_mod.f90 @@ -35,13 +35,13 @@ Module psb_tools_mod interface psb_geall ! 2-D double precision version - subroutine psb_dalloc(x, desc_a, info, n) + subroutine psb_dalloc(x, desc_a, info, n, lb) use psb_descriptor_type implicit none real(kind(1.d0)), allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer,intent(out) :: info - integer, optional, intent(in) :: n + integer, optional, intent(in) :: n, lb end subroutine psb_dalloc ! 1-D double precision version subroutine psb_dallocv(x, desc_a,info,n) @@ -52,12 +52,12 @@ Module psb_tools_mod integer, optional, intent(in) :: n end subroutine psb_dallocv ! 2-D integer version - subroutine psb_ialloc(x, desc_a, info,n) + subroutine psb_ialloc(x, desc_a, info,n, lb) use psb_descriptor_type integer, allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, optional, intent(in) :: n + integer, optional, intent(in) :: n, lb end subroutine psb_ialloc subroutine psb_iallocv(x, desc_a,info,n) use psb_descriptor_type @@ -66,16 +66,16 @@ Module psb_tools_mod integer :: info integer, optional, intent(in) :: n end subroutine psb_iallocv - ! 2-D double precision version - subroutine psb_zalloc(x, desc_a, info, n) + ! 2-D double complex version + subroutine psb_zalloc(x, desc_a, info, n, lb) use psb_descriptor_type implicit none complex(kind(1.d0)), allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer :: info - integer, optional, intent(in) :: n + integer, optional, intent(in) :: n, lb end subroutine psb_zalloc - ! 1-D double precision version + ! 1-D double complex version subroutine psb_zallocv(x, desc_a,info,n) use psb_descriptor_type complex(kind(1.d0)), allocatable, intent(out) :: x(:) diff --git a/base/tools/psb_dallc.f90 b/base/tools/psb_dallc.f90 index 12da806f..f1f3bcdd 100644 --- a/base/tools/psb_dallc.f90 +++ b/base/tools/psb_dallc.f90 @@ -40,11 +40,13 @@ ! desc_a - the communication descriptor. ! info - Return code ! n - optional number of columns. -subroutine psb_dalloc(x, desc_a, info, n) +! lb - optional lower bound on column indices +subroutine psb_dalloc(x, desc_a, info, n, lb) !....allocate dense matrix for psblas routines..... use psb_descriptor_type use psb_const_mod use psb_error_mod + use psb_realloc_mod use psb_penv_mod implicit none @@ -53,15 +55,15 @@ subroutine psb_dalloc(x, desc_a, info, n) real(kind(1.d0)), allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer,intent(out) :: info - integer, optional, intent(in) :: n + integer, optional, intent(in) :: n, lb !locals - integer :: np,me,err,n_col,n_row,i,j,err_act + integer :: np,me,err,nr,i,j,err_act integer :: ictxt,n_ integer :: int_err(5), exch(3) character(len=20) :: name - name='psb_dgeall_m' + name='psb_geall' if(psb_get_errstatus() /= 0) return info=0 err=0 @@ -105,34 +107,24 @@ subroutine psb_dalloc(x, desc_a, info, n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col,n_),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_col*n_ - call psb_errpush(info,name,int_err,a_err='real(kind(1.d0))') - goto 9999 - endif - do j=1,n_ - do i=1,n_col - x(i,j) = 0.0d0 - end do - end do + nr = max(1,psb_cd_get_local_cols(desc_a)) else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row,n_),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_row*n_ - call psb_errpush(info,name,int_err,a_err='real(kind(1.d0))') - goto 9999 - endif - do j = 1, n_ - do i=1,n_row - x(i,j) = 0.0d0 - end do - end do + nr = max(1,psb_cd_get_local_rows(desc_a)) + else + info = 4001 + call psb_errpush(info,name,int_err,a_err='Invalid desc_a') + goto 9999 + endif + + call psb_realloc(nr,n_,x,info,lb2=lb) + if (info /= 0) then + info=4025 + int_err(1)=nr*n_ + call psb_errpush(info,name,int_err,a_err='real(kind(1.d0))') + goto 9999 endif + + x(:,:) = dzero call psb_erractionrestore(err_act) return @@ -184,7 +176,7 @@ end subroutine psb_dalloc ! The descriptor may be in either the build or assembled state. ! ! Arguments: -! x - the matrix to be allocated. +! x(:) - the matrix to be allocated. ! desc_a - the communication descriptor. ! info - return code subroutine psb_dallocv(x, desc_a,info,n) @@ -211,7 +203,7 @@ subroutine psb_dallocv(x, desc_a,info,n) if(psb_get_errstatus() /= 0) return info=0 - name='psb_dgeall_v' + name='psb_geall' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() @@ -237,35 +229,24 @@ subroutine psb_dallocv(x, desc_a,info,n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_col - call psb_errpush(info,name,int_err,a_err='real(kind(1.d0))') - goto 9999 - endif - do i=1,n_col - x(i) = 0.0d0 - end do - + nr = max(1,psb_cd_get_local_cols(desc_a)) else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_row - call psb_errpush(info,name,int_err,a_err='real(kind(1.d0))') - goto 9999 - endif - do i=1,n_row - x(i) = 0.0d0 - end do + nr = max(1,psb_cd_get_local_rows(desc_a)) else - if (debug_level > psb_debug_ext_) & - & write(debug_unit,*) me,name,& - & ': Did not allocate anything because of dectype',psb_cd_get_dectype(desc_a) + info = 4001 + call psb_errpush(info,name,int_err,a_err='Invalid desc_a') + goto 9999 + endif + + call psb_realloc(nr,x,info) + if (info /= 0) then + info=4025 + int_err(1)=nr + call psb_errpush(info,name,int_err,a_err='real(kind(1.d0))') + goto 9999 endif + + x(:) = dzero call psb_erractionrestore(err_act) return diff --git a/base/tools/psb_ialloc.f90 b/base/tools/psb_ialloc.f90 index 9ffc6e21..08df4f91 100644 --- a/base/tools/psb_ialloc.f90 +++ b/base/tools/psb_ialloc.f90 @@ -39,11 +39,13 @@ ! desc_a - the communication descriptor. ! info - possibly returns an error code ! n - optional number of columns. -subroutine psb_ialloc(x, desc_a, info, n) +! lb - optional lower bound on column indices +subroutine psb_ialloc(x, desc_a, info, n, lb) !....allocate dense matrix for psblas routines..... use psb_descriptor_type use psb_const_mod use psb_error_mod + use psb_realloc_mod use psb_penv_mod implicit none @@ -52,18 +54,19 @@ subroutine psb_ialloc(x, desc_a, info, n) integer, allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, optional, intent(in) :: n - + integer, optional, intent(in) :: n, lb !locals - integer :: np,me,n_col,n_row,i,j,err_act + integer :: np,me,err,nr,i,j,err_act integer :: ictxt,n_ integer :: int_err(5), exch(3) character(len=20) :: name + name='psb_geall' if(psb_get_errstatus() /= 0) return info=0 - name='psb_ialloc' + err=0 + int_err(1)=0 call psb_erractionsave(err_act) ictxt=psb_cd_get_context(desc_a) @@ -103,34 +106,24 @@ subroutine psb_ialloc(x, desc_a, info, n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col,n_),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_col*n_ - call psb_errpush(info,name,int_err,a_err='integer') - goto 9999 - endif - do j=1,n_ - do i=1,n_col - x(i,j) = 0 - end do - end do + nr = max(1,psb_cd_get_local_cols(desc_a)) else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row,n_),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_row*n_ - call psb_errpush(info,name,int_err,a_err='integer') - goto 9999 - endif - do j = 1, n_ - do i=1,n_row - x(i,j) = 0 - end do - end do + nr = max(1,psb_cd_get_local_rows(desc_a)) + else + info = 4001 + call psb_errpush(info,name,int_err,a_err='Invalid desc_a') + goto 9999 + endif + + call psb_realloc(nr,n_,x,info,lb2=lb) + if (info /= 0) then + info=4025 + int_err(1)=nr*n_ + call psb_errpush(info,name,int_err,a_err='integer') + goto 9999 endif + + x(:,:) = izero call psb_erractionrestore(err_act) return @@ -182,10 +175,9 @@ end subroutine psb_ialloc ! The descriptor may be in either the build or assembled state. ! ! Arguments: -! m - integer. The number of rows. -! x - integer,dimension(:). The matrix to be allocated. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Eventually returns an error code +! x(:) - the matrix to be allocated. +! desc_a - the communication descriptor. +! info - return code subroutine psb_iallocv(x, desc_a, info,n) !....allocate sparse matrix structure for psblas routines..... use psb_descriptor_type @@ -210,7 +202,7 @@ subroutine psb_iallocv(x, desc_a, info,n) if(psb_get_errstatus() /= 0) return info=0 - name='psb_igeall_v' + name='psb_geall' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() @@ -236,33 +228,24 @@ subroutine psb_iallocv(x, desc_a, info,n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_col - call psb_errpush(info,name,int_err,a_err='integer') - goto 9999 - endif + nr = max(1,psb_cd_get_local_cols(desc_a)) else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_row - call psb_errpush(info,name,int_err,a_err='integer') - goto 9999 - endif - do i=1,n_row - x(i) = 0.0d0 - end do + nr = max(1,psb_cd_get_local_rows(desc_a)) else - if (debug_level > psb_debug_ext_) & - & write(debug_unit,*) me,name,& - & ': Did not allocate anything because of dectype',psb_cd_get_dectype(desc_a) + info = 4001 + call psb_errpush(info,name,int_err,a_err='Invalid desc_a') + goto 9999 endif - - x = 0 + + call psb_realloc(nr,x,info) + if (info /= 0) then + info=4025 + int_err(1)=nr + call psb_errpush(info,name,int_err,a_err='integer') + goto 9999 + endif + + x(:) = izero call psb_erractionrestore(err_act) return diff --git a/base/tools/psb_zallc.f90 b/base/tools/psb_zallc.f90 index 2ca1edf1..1fa763f5 100644 --- a/base/tools/psb_zallc.f90 +++ b/base/tools/psb_zallc.f90 @@ -40,11 +40,13 @@ ! desc_a - the communication descriptor. ! info - Return code ! n - optional number of columns. -subroutine psb_zalloc(x, desc_a, info, n) +! lb - optional lower bound on column indices +subroutine psb_zalloc(x, desc_a, info, n, lb) !....allocate dense matrix for psblas routines..... use psb_descriptor_type use psb_const_mod use psb_error_mod + use psb_realloc_mod use psb_penv_mod implicit none @@ -53,15 +55,15 @@ subroutine psb_zalloc(x, desc_a, info, n) complex(kind(1.d0)), allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer,intent(out) :: info - integer, optional, intent(in) :: n + integer, optional, intent(in) :: n, lb !locals - integer :: np,me,err,n_col,n_row,i,j,err_act + integer :: np,me,err,nr,i,j,err_act integer :: ictxt,n_ integer :: int_err(5),exch(3) character(len=20) :: name - name='psb_zgeall_m' + name='psb_geall' if(psb_get_errstatus() /= 0) return info=0 err=0 @@ -105,34 +107,24 @@ subroutine psb_zalloc(x, desc_a, info, n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col,n_),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_col*n_ - call psb_errpush(info,name,int_err,a_err='complex(kind(1.d0))') - goto 9999 - endif - do j=1,n_ - do i=1,n_col - x(i,j) = 0.0d0 - end do - end do + nr = max(1,psb_cd_get_local_cols(desc_a)) else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row,n_),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_row*n_ - call psb_errpush(info,name,int_err,a_err='complex(kind(1.d0))') - goto 9999 - endif - do j = 1, n_ - do i=1,n_row - x(i,j) = 0.0d0 - end do - end do + nr = max(1,psb_cd_get_local_rows(desc_a)) + else + info = 4001 + call psb_errpush(info,name,int_err,a_err='Invalid desc_a') + goto 9999 + endif + + call psb_realloc(nr,n_,x,info,lb2=lb) + if (info /= 0) then + info=4025 + int_err(1)=nr*n_ + call psb_errpush(info,name,int_err,a_err='complex(kind(1.d0))') + goto 9999 endif + + x(:,:) = zzero call psb_erractionrestore(err_act) return @@ -183,7 +175,7 @@ end subroutine psb_zalloc ! The descriptor may be in either the build or assembled state. ! ! Arguments: -! x - the matrix to be allocated. +! x(:) - the matrix to be allocated. ! desc_a - the communication descriptor. ! info - return code subroutine psb_zallocv(x, desc_a,info,n) @@ -210,7 +202,7 @@ subroutine psb_zallocv(x, desc_a,info,n) if(psb_get_errstatus() /= 0) return info=0 - name='psb_zgeall_v' + name='psb_geall' call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() @@ -236,35 +228,24 @@ subroutine psb_zallocv(x, desc_a,info,n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_col - call psb_errpush(info,name,int_err,a_err='complex(kind(1.d0))') - goto 9999 - endif - do i=1,n_col - x(i) = 0.0d0 - end do - + nr = max(1,psb_cd_get_local_cols(desc_a)) else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row),stat=info) - if (info /= 0) then - info=4025 - int_err(1)=n_row - call psb_errpush(info,name,int_err,a_err='complex(kind(1.d0))') - goto 9999 - endif - do i=1,n_row - x(i) = 0.0d0 - end do + nr = max(1,psb_cd_get_local_rows(desc_a)) else - if (debug_level > psb_debug_ext_) & - & write(debug_unit,*) me,name,& - & ': Did not allocate anything because of dectype',psb_cd_get_dectype(desc_a) + info = 4001 + call psb_errpush(info,name,int_err,a_err='Invalid desc_a') + goto 9999 + endif + + call psb_realloc(nr,x,info) + if (info /= 0) then + info=4025 + int_err(1)=nr + call psb_errpush(info,name,int_err,a_err='complex(kind(1.d0))') + goto 9999 endif + + x(:) = zzero call psb_erractionrestore(err_act) return diff --git a/krylov/psb_dcgstabl.f90 b/krylov/psb_dcgstabl.f90 index f912abde..94cdc4af 100644 --- a/krylov/psb_dcgstabl.f90 +++ b/krylov/psb_dcgstabl.f90 @@ -208,8 +208,8 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is goto 9999 end if if (info == 0) Call psb_geall(wwrk,desc_a,info,n=10) - if (info == 0) Call psb_geall(uh,desc_a,info,n=nl+1) - if (info == 0) Call psb_geall(rh,desc_a,info,n=nl+1) + if (info == 0) Call psb_geall(uh,desc_a,info,n=nl+1,lb=0) + if (info == 0) Call psb_geall(rh,desc_a,info,n=nl+1,lb=0) if (info == 0) Call psb_geasb(wwrk,desc_a,info) if (info == 0) Call psb_geasb(uh,desc_a,info) if (info == 0) Call psb_geasb(rh,desc_a,info) @@ -385,7 +385,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is if (info == 0) call psb_gefree(wwrk,desc_a,info) if (info == 0) call psb_gefree(uh,desc_a,info) if (info == 0) call psb_gefree(rh,desc_a,info) - if(info/=0) then + if (info/=0) then call psb_errpush(info,name) goto 9999 end if