fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-06 22:55:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merged debug infrastructure, internal docs and html headers.

Browse Source
This commit is contained in:
Salvatore Filippone
2007-12-21 13:32:16 +00:00
parent 8eada499ff
commit 1de99a499a
182 changed files with 3181 additions and 2837 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Changelog
+2
View File
@@ -1,6 +1,8 @@
Changelog. A lot less detailed than usual, at least for past
history.
2007/12/21: Merged in debug infrastructure, internal and html docs.
2007/11/14: Fix INTENT(IN) on X vector in preconditioner routines.
2007/10/15: Repackaged the sorting routines in a submodule of their
base/comm/psb_dgather.f90
+2 -2
View File
@@ -38,7 +38,7 @@
! the distributed pieces.
! locx - real,dimension(:,:). The local piece of the distributed
! matrix to be gathered.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
@@ -209,7 +209,7 @@ end subroutine psb_dgatherm
! distributed pieces.
! locx - real,dimension(:). The local piece of the ditributed
! vector to be gathered.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! root - integer. The process that has to own the
! global matrix. If -1 all
base/comm/psb_dhalo.f90
+2 -2
View File
@@ -36,7 +36,7 @@
!
! Arguments:
! x - real,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! alpha - real(optional). Scale factor.
! jx - integer(optional). The starting column of the global matrix.
@@ -252,7 +252,7 @@ end subroutine psb_dhalom
!
! Arguments:
! x - real,dimension(:). The local part of the dense vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! alpha - real(optional). Scale factor.
! work - real(optional). Work area.
base/comm/psb_dovrl.f90
+2 -2
View File
@@ -36,7 +36,7 @@
!
! Arguments:
! x(:,:) - real The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! jx - integer(optional). The starting column of the global matrix
! ik - integer(optional). The number of columns to gather.
@@ -254,7 +254,7 @@ end subroutine psb_dovrlm
!
! Arguments:
! x(:) - real The local part of the dense vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! work - real(optional). A work area.
! update - integer(optional). Type of update:
base/comm/psb_dscatter.F90
+3 -4
View File
@@ -37,10 +37,9 @@
! Arguments:
! globx - real,dimension(:,:). The global matrix to scatter.
! locx - real,dimension(:,:). The local piece of the ditributed matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer(optional). The process that owns the global matrix.
! If -1 all
! iroot - integer(optional). The process that owns the global matrix. If -1 all
! the processes have a copy. Default -1.
!
subroutine psb_dscatterm(globx, locx, desc_a, info, iroot)
@@ -267,7 +266,7 @@ end subroutine psb_dscatterm
! Arguments:
! globx - real,dimension(:). The global vector to scatter.
! locx - real,dimension(:). The local piece of the ditributed vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer(optional). The process that owns the global vector. If -1 all
! the processes have a copy.
base/comm/psb_igather.f90
+2 -2
View File
@@ -38,7 +38,7 @@
! the distributed pieces.
! locx - integer,dimension(:,:). The local piece of the distributed
! matrix to be gathered.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
@@ -209,7 +209,7 @@ end subroutine psb_igatherm
! distributed pieces.
! locx - integer,dimension(:). The local piece of the ditributed
! vector to be gathered.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
base/comm/psb_ihalo.f90
+2 -2
View File
@@ -37,7 +37,7 @@
!
! Arguments:
! x - integer,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! alpha - real(optional). Scale factor.
! jx - integer(optional). The starting column of the global matrix.
@@ -254,7 +254,7 @@ end subroutine psb_ihalom
!
! Arguments:
! x - integer,dimension(:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! alpha - real(optional). Scale factor.
! jx - integer(optional). The starting column of the global matrix.
base/comm/psb_iscatter.F90
+2 -2
View File
@@ -37,7 +37,7 @@
! Arguments:
! globx - integer,dimension(:,:). The global matrix to scatter.
! locx - integer,dimension(:,:). The local piece of the ditributed matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer(optional). The process that owns the global matrix. If -1 all
! the processes have a copy.
@@ -264,7 +264,7 @@ end subroutine psb_iscatterm
! Arguments:
! globx - integer,dimension(:). The global vector to scatter.
! locx - integer,dimension(:). The local piece of the ditributed vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer(optional). The process that owns the global vector. If -1 all
! the processes have a copy.
base/comm/psb_zgather.f90
+2 -2
View File
@@ -38,7 +38,7 @@
! the distributed pieces.
! locx - cplx,dimension(:,:). The local piece of the distributed
! matrix to be gathered.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
@@ -211,7 +211,7 @@ end subroutine psb_zgatherm
! the distributed pieces.
! locx - cplx,dimension(:). The local piece of the distributed
! vector to be gathered.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
base/comm/psb_zhalo.f90
+2 -2
View File
@@ -36,7 +36,7 @@
!
! Arguments:
! x - real,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! alpha - complex(optional). Scale factor.
! jx - integer(optional). The starting column of the global matrix.
@@ -251,7 +251,7 @@ end subroutine psb_zhalom
!
! Arguments:
! x - real,dimension(:). The local part of the dense vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! alpha - complex(optional). Scale factor.
! jx - integer(optional). The starting column of the global matrix.
base/comm/psb_zovrl.f90
+2 -2
View File
@@ -36,7 +36,7 @@
!
! Arguments:
! x(:,:) - complex The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! jx - integer(optional). The starting column of the global matrix
! ik - integer(optional). The number of columns to gather.
@@ -253,7 +253,7 @@ end subroutine psb_zovrlm
!
! Arguments:
! x(:) - complex The local part of the dense vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! work - real(optional). A work area.
! update - integer(optional). Type of update:
base/comm/psb_zscatter.F90
+2 -2
View File
@@ -37,7 +37,7 @@
! Arguments:
! globx - complex,dimension(:,:). The global matrix to scatter.
! locx - complex,dimension(:,:). The local piece of the distributed matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer(optional). The process that owns the global matrix.
! If -1 all the processes have a copy.
@@ -269,7 +269,7 @@ end subroutine psb_zscatterm
! Arguments:
! globx - complex,dimension(:). The global vector to scatter.
! locx - complex,dimension(:). The local piece of the ditributed vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! iroot - integer(optional). The process that owns the global vector. If -1 all
! the processes have a copy.
base/internals/psi_compute_size.f90
+5 -3
View File
@@ -52,11 +52,13 @@ subroutine psi_compute_size(desc_data, index_in, dl_lda, info)
integer, allocatable :: counter_recv(:), counter_dl(:)
! ...parameters
logical, parameter :: debug=.false.
integer :: debug_level, debug_unit
character(len=20) :: name
name='psi_compute_size'
call psb_get_erraction(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = 0
ictxt = desc_data(psb_ctxt_)
@@ -113,8 +115,8 @@ subroutine psi_compute_size(desc_data, index_in, dl_lda, info)
! computing max global value of dl_lda
call psb_amx(ictxt, dl_lda)
if (debug) then
write(0,*) 'psi_compute_size: ',dl_lda
if (debug_level>=psb_debug_inner_) then
write(debug_unit,*) me,' ',trim(name),': ',dl_lda
endif
call psb_erractionrestore(err_act)
base/internals/psi_crea_bnd_elem.f90
+1 -1
View File
@@ -39,7 +39,7 @@
!
! Arguments:
! bndel(:) - integer, allocatable Array containing the output list
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
!
subroutine psi_crea_bnd_elem(bndel,desc_a,info)
base/internals/psi_crea_index.f90
+16 -6
View File
@@ -71,12 +71,14 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
! ...parameters...
integer, allocatable :: dep_list(:,:), length_dl(:)
integer,parameter :: root=psb_root_,no_comm=-1
logical,parameter :: debug=.false.
integer :: debug_level, debug_unit
character(len=20) :: name
info = 0
name='psi_crea_index'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
@@ -99,7 +101,8 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
! ...extract dependence list (ordered list of identifer process
! which every process must communcate with...
if (debug) write(*,*) 'crea_halo: calling extract_dep_list'
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),': calling extract_dep_list'
mode = 1
call psi_extract_dep_list(desc_a%matrix_data,index_in,&
@@ -109,10 +112,12 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
goto 9999
end if
if (debug) write(0,*) 'crea_index: from extract_dep_list',&
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),': from extract_dep_list',&
& me,length_dl(0),index_in(1), ':',dep_list(:length_dl(me),me)
! ...now process root contains dependence list of all processes...
if (debug) write(0,*) 'crea_index: root sorting dep list'
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),': root sorting dep list'
call psi_dl_check(dep_list,max(1,dl_lda),np,length_dl)
@@ -123,12 +128,14 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
goto 9999
end if
if(debug) write(0,*)'in psi_crea_index calling psi_desc_index',&
if(debug_level >= psb_debug_inner_)&
& write(debug_unit,*) me,' ',trim(name),': calling psi_desc_index',&
& size(index_out)
! Do the actual format conversion.
call psi_desc_index(desc_a,index_in,dep_list(1:,me),&
& length_dl(me),nsnd,nrcv, index_out,glob_idx,info)
if(debug) write(0,*)'out of psi_desc_index',&
if(debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),': out of psi_desc_index',&
& size(index_out)
nxch = length_dl(me)
if(info /= 0) then
@@ -137,6 +144,9 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
end if
deallocate(dep_list,length_dl)
if(debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),': done'
call psb_erractionrestore(err_act)
return
base/internals/psi_crea_ovr_elem.f90
+4 -2
View File
@@ -38,7 +38,7 @@
!
! Arguments:
! ovr_elem(:) - integer, allocatable Array containing the output list
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
!
subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
@@ -84,6 +84,8 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
if (usetree) then
!
! This is now here just for historical reasons.
!
! While running through the column indices exchanged with other procs
! we have to record them in overlap_elem. We do this by maintaining
@@ -148,7 +150,7 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
else if (.not.usetree) then
! Simple alternative.
insize = size(desc_overlap)
insize = max(1,(insize+1)/2)
allocate(telem(insize,2),stat=info)
base/internals/psi_desc_index.F90
+15 -12
View File
@@ -131,12 +131,15 @@ subroutine psi_desc_index(desc,index_in,dep_list,&
integer :: ihinsz,ntot,k,err_act,nidx,&
& idxr, idxs, iszs, iszr, nesd, nerv, icomm
logical,parameter :: debug=.false., usempi=.true.
character(len=20) :: name
logical,parameter :: usempi=.true.
integer :: debug_level, debug_unit
character(len=20) :: name
info = 0
name='psi_desc_index'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc)
icomm = psb_cd_get_mpic(desc)
@@ -147,8 +150,8 @@ subroutine psi_desc_index(desc,index_in,dep_list,&
goto 9999
endif
if (debug) then
write(0,*) me,'start desc_index'
if (debug_level >= psb_debug_inner_) then
write(debug_unit,*) me,' ',trim(name),': start'
call psb_barrier(ictxt)
endif
@@ -203,8 +206,8 @@ subroutine psi_desc_index(desc,index_in,dep_list,&
if ((iszs /= idxs).or.(iszr /= idxr)) then
write(0,*) 'strange results?', iszs,idxs,iszr,idxr
end if
if (debug) then
write(0,*) me,'computed sizes ',iszr,iszs
if (debug_level >= psb_debug_inner_) then
write(debug_unit,*) me,' ',trim(name),': computed sizes ',iszr,iszs
call psb_barrier(ictxt)
endif
@@ -223,8 +226,8 @@ subroutine psi_desc_index(desc,index_in,dep_list,&
goto 9999
end if
if (debug) then
write(0,*) me,'computed allocated workspace ',iszr,iszs
if (debug_level >= psb_debug_inner_) then
write(debug_unit,*) me,' ',trim(name),': computed allocated workspace ',iszr,iszs
call psb_barrier(ictxt)
endif
allocate(sndbuf(iszs),rcvbuf(iszr),stat=info)
@@ -264,8 +267,8 @@ subroutine psi_desc_index(desc,index_in,dep_list,&
i = i + nerv + 1
end do
if (debug) then
write(0,*) me,' prepared send buffer '
if (debug_level >= psb_debug_inner_) then
write(debug_unit,*) me,' ',trim(name),': prepared send buffer '
call psb_barrier(ictxt)
endif
!
@@ -317,8 +320,8 @@ subroutine psi_desc_index(desc,index_in,dep_list,&
goto 9999
end if
if (debug) then
write(0,*) me,'end desc_index'
if (debug_level >= psb_debug_inner_) then
write(debug_unit,*) me,' ',trim(name),': done'
call psb_barrier(ictxt)
endif
base/internals/psi_dswapdata.F90
+5 -6
View File
@@ -29,9 +29,9 @@
!!$
!!$
!
! File: psi_Xswapdata.F90
! File: psi_dswapdata.F90
!
! Subroutine: psi_Xswapdatam
! Subroutine: psi_dswapdatam
! Does the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
@@ -68,7 +68,7 @@
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:,:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
@@ -504,9 +504,8 @@ end subroutine psi_dswapdatam
!!$
!!$
!
! File: psi_Xswapdata.F90
!
! Subroutine: psi_Xswapdatav
! Subroutine: psi_dswapdatav
! Does the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
@@ -543,7 +542,7 @@ end subroutine psi_dswapdatam
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
base/internals/psi_dswaptran.F90
+5 -6
View File
@@ -29,9 +29,9 @@
!!$
!!$
!
! File: psi_Xswaptran.F90
! File: psi_dswaptran.F90
!
! Subroutine: psi_Xswaptranm
! Subroutine: psi_dswaptranm
! Does the data exchange among processes. This is similar to Xswapdata, but
! the list is read "in reverse", i.e. indices that are normally SENT are used
! for the RECEIVE part and vice-versa. This is the basic data exchange operation
@@ -72,7 +72,7 @@
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:,:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
@@ -498,9 +498,8 @@ end subroutine psi_dswaptranm
!!$
!!$
!
! File: psi_Xswaptran.F90
!
! Subroutine: psi_Xswaptranv
! Subroutine: psi_dswaptranv
! Does the data exchange among processes. This is similar to Xswapdata, but
! the list is read "in reverse", i.e. indices that are normally SENT are used
! for the RECEIVE part and vice-versa. This is the basic data exchange operation
@@ -541,7 +540,7 @@ end subroutine psi_dswaptranm
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
base/internals/psi_extrct_dl.F90
+22 -9
View File
@@ -141,11 +141,13 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
! .....local scalars...
integer i,me,nprow,pointer_dep_list,proc,j,err_act
integer ictxt, err, icomm
logical, parameter :: debug=.false.
integer :: debug_level, debug_unit
character name*20
name='psi_extrct_dl'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = 0
ictxt = desc_data(psb_ctxt_)
@@ -156,11 +158,13 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
length_dl(i) = 0
enddo
i=1
if (debug) write(0,*) 'extract: info ',info,desc_data(psb_dec_type_)
if (debug_level >= psb_debug_inner_)&
& write(debug_unit,*) me,' ',trim(name),': start ',info,desc_data(psb_dec_type_)
pointer_dep_list=1
if (psb_is_bld_dec(desc_data(psb_dec_type_))) then
do while (desc_str(i) /= -1)
if (debug) write(0,*) me,' extract: looping ',i,&
if (debug_level >= psb_debug_inner_)&
& write(debug_unit,*) me,' ',trim(name),' : looping ',i,&
& desc_str(i),desc_str(i+1),desc_str(i+2)
! ...with different decomposition type we have different
@@ -169,7 +173,8 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
! ..if number of element to be exchanged !=0
proc=desc_str(i)
if ((proc < 0).or.(proc.ge.nprow)) then
if (debug) write(0,*) 'extract error ',i,desc_str(i)
if (debug_level >= psb_debug_inner_)&
& write(debug_unit,*) me,' ',trim(name),': error ',i,desc_str(i)
info = 9999
int_err(1) = i
int_err(2) = desc_str(i)
@@ -203,7 +208,8 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
enddo
else if (psb_is_upd_dec(desc_data(psb_dec_type_))) then
do while (desc_str(i) /= -1)
if (debug) write(0,*) 'extract: looping ',i,desc_str(i)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),': looping ',i,desc_str(i)
! ...with different decomposition type we have different
! structure of indices lists............................
@@ -240,7 +246,6 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
i=i+desc_str(i+1)+2
enddo
else
write(0,*) 'invalid dec_type',desc_data(psb_dec_type_)
info = 2020
goto 9999
endif
@@ -249,7 +254,8 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
! ... check for errors...
998 continue
if (debug) write(0,*) 'extract: info ',info
if (debug_level >= psb_debug_inner_)&
& write(debug_unit,*) me,' ',trim(name),': info ',info
err = info
if (err /= 0) goto 9999
@@ -257,11 +263,18 @@ subroutine psi_extract_dep_list(desc_data,desc_str,dep_list,&
call psb_sum(ictxt,length_dl(0:np))
call psb_get_mpicomm(ictxt,icomm )
allocate(itmp(dl_lda),stat=info)
if (info /= 0) goto 9999
if (info /= 0) then
info=4000
goto 9999
endif
itmp(1:dl_lda) = dep_list(1:dl_lda,me)
call mpi_allgather(itmp,dl_lda,mpi_integer,&
& dep_list,dl_lda,mpi_integer,icomm,info)
deallocate(itmp)
deallocate(itmp,stat=info)
if (info /= 0) then
info=4000
goto 9999
endif
call psb_erractionrestore(err_act)
return
base/internals/psi_fnd_owner.f90
+10 -11
View File
@@ -40,7 +40,7 @@
! idx(:) - integer Required indices on the calling process
! iprc(:) - integer, allocatable Output: process identifiers for the corresponding
! indices
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
!
subroutine psi_fnd_owner(nv,idx,iprc,desc,info)
@@ -64,7 +64,6 @@ subroutine psi_fnd_owner(nv,idx,iprc,desc,info)
integer :: i,n_row,n_col, err_act,ih,icomm,hsize
integer :: ictxt,np,me
logical, parameter :: debug=.false., debugwrt=.false.
character(len=20) :: name
info = 0
@@ -92,20 +91,13 @@ subroutine psi_fnd_owner(nv,idx,iprc,desc,info)
endif
if (.not.(psb_is_ok_desc(desc))) then
write(0,*) 'Invalid input descriptor in psi_fnd_owner'
call psb_errpush(4010,name,a_err='invalid desc')
goto 9999
end if
!
! The basic idea is very simple.
! First we figure out the total number of requests.
! Second we build the aggregate list of requests (with psb_amx)
! Third, we figure out locally whether we own the indices (whoever is
! asking for them) and build our part of the reply (we shift process
! indices so that they're 1-based)
! Fourth, we do a psb_amx on the replies so that we have everybody's answers
! Fifth, we extract the answers for our local query, and shift back the
! process indices to 0-based.
Allocate(hidx(np+1),hsz(np),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
@@ -124,12 +116,16 @@ subroutine psi_fnd_owner(nv,idx,iprc,desc,info)
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
! Second we build the aggregate list of requests (with psb_amx)
helem(:) = 0
ih = hidx(me+1)
do i=1, hsz(me+1)
helem(ih+i-1) = idx(i)
end do
call psb_amx(ictxt,helem,info)
! Third, we figure out locally whether we own the indices (whoever is
! asking for them) and build our part of the reply (we shift process
! indices so that they're 1-based)
call psi_idx_cnv(hsize,helem,desc,info,owned=.true.)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psi_idx_cnv')
@@ -144,8 +140,11 @@ subroutine psi_fnd_owner(nv,idx,iprc,desc,info)
end if
end do
! Fourth, we do a psb_amx on the replies so that we have everybody's answers
call psb_amx(ictxt,hproc,info)
! Fifth, we extract the answers for our local query, and shift back the
! process indices to 0-based.
if (nv > 0) then
call psb_realloc(nv,iprc,info)
ih = hidx(me+1)
base/internals/psi_idx_cnv.f90
+3 -7
View File
@@ -38,7 +38,7 @@
! Arguments:
! nv - integer Number of indices required
! idxin(:) - integer Required indices, overwritten on output.
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! mask(:) - logical, optional Only do the conversion for specific indices.
! owned - logical,optional Restrict to local indices, no halo
@@ -61,7 +61,6 @@ subroutine psi_idx_cnv1(nv,idxin,desc,info,mask,owned)
integer :: np, me
integer :: nrow,ncol, err_act
integer, allocatable :: idxout(:)
logical, parameter :: debug=.false.
integer, parameter :: relocsz=200
character(len=20) :: name
logical, pointer :: mask_(:)
@@ -183,7 +182,6 @@ end subroutine psi_idx_cnv1
!!$
!!$
!
! File: psi_idx_cnv.f90
!
! Subroutine: psi_idx_cnv2
! Converts a bunch of indices from global to local numbering.
@@ -193,7 +191,7 @@ end subroutine psi_idx_cnv1
! nv - integer Number of indices required
! idxin(:) - integer Required indices
! idxout(:) - integer Output values, negative for invalid input.
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! mask(:) - logical, optional Only do the conversion for specific indices.
! owned - logical,optional Restrict to local indices, no halo
@@ -215,7 +213,6 @@ subroutine psi_idx_cnv2(nv,idxin,idxout,desc,info,mask,owned)
integer :: i,ictxt,mglob, nglob
integer :: np, me
integer :: nrow,ncol, ip, err_act,lip
logical, parameter :: debug=.false.
integer, parameter :: relocsz=200
character(len=20) :: name
logical, pointer :: mask_(:)
@@ -405,7 +402,6 @@ end subroutine psi_idx_cnv2
!!$
!!$
!
! File: psi_idx_cnv.f90
!
! Subroutine: psi_idx_cnvs
! Converts an index from global to local numbering.
@@ -414,7 +410,7 @@ end subroutine psi_idx_cnv2
! Arguments:
! idxin - integer Required index
! idxout - integer Output value, negative for invalid input.
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! mask - logical, optional Only do the conversion if true.
! owned - logical,optional Restrict to local indices, no halo
base/internals/psi_idx_ins_cnv.f90
+3 -11
View File
@@ -43,7 +43,7 @@
! nv - integer Number of indices required
! idxin(:) - integer Required indices, overwritten on output
! output is negative for masked entries
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! mask(:) - logical, optional Only do the conversion for specific indices.
!
@@ -64,7 +64,6 @@ subroutine psi_idx_ins_cnv1(nv,idxin,desc,info,mask)
integer :: np, me
integer :: nrow,ncol, err_act
integer, allocatable :: idxout(:)
logical, parameter :: debug=.false.
integer, parameter :: relocsz=200
character(len=20) :: name
logical, pointer :: mask_(:)
@@ -179,7 +178,6 @@ end subroutine psi_idx_ins_cnv1
!!$
!!$
!
! File: psi_idx_ins_cnv.f90
!
! Subroutine: psi_idx_ins_cnv2
! Converts a bunch of indices from global to local numbering.
@@ -193,7 +191,7 @@ end subroutine psi_idx_ins_cnv1
! nv - integer Number of indices required
! idxin(:) - integer Required indices
! idxout(:) - integer Output values (negative for masked entries)
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! mask(:) - logical, optional Only do the conversion for specific indices.
!
@@ -214,7 +212,6 @@ subroutine psi_idx_ins_cnv2(nv,idxin,idxout,desc,info,mask)
integer :: i,ictxt,k,mglob, nglob
integer :: np, me, isize
integer :: pnt_halo,nrow,ncol, nh, ip, err_act,lip,nxt
logical, parameter :: debug=.false.
integer, parameter :: relocsz=200
character(len=20) :: name,ch_err
logical, pointer :: mask_(:)
@@ -293,7 +290,6 @@ subroutine psi_idx_ins_cnv2(nv,idxin,idxout,desc,info,mask)
if (ncol > isize) then
nh = ncol + max(nv,relocsz)
call psb_realloc(nh,desc%loc_to_glob,info,pad=-1)
if (debug) write(0,*) 'done realloc ',nh
if (info /= 0) then
info=1
ch_err='psb_realloc'
@@ -346,9 +342,6 @@ subroutine psi_idx_ins_cnv2(nv,idxin,idxout,desc,info,mask)
if (ncol > isize) then
nh = ncol + max(nv,relocsz)
call psb_realloc(nh,desc%loc_to_glob,info,pad=-1)
if (me==0) then
if (debug) write(0,*) 'done realloc ',nh
end if
if (info /= 0) then
info=3
ch_err='psb_realloc'
@@ -437,7 +430,6 @@ end subroutine psi_idx_ins_cnv2
!!$
!!$
!
! File: psi_idx_ins_cnv.f90
!
! Subroutine: psi_idx_ins_cnvs
! Converts an index from global to local numbering.
@@ -450,7 +442,7 @@ end subroutine psi_idx_ins_cnv2
! Arguments:
! idxin - integer Required index s
! idxout - integer Output value (negative for masked entries)
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! info - integer. return code.
! mask - logical, optional Only do the conversion for specific indices.
!
base/internals/psi_iswapdata.F90
+5 -6
View File
@@ -29,9 +29,9 @@
!!$
!!$
!
! File: psi_Xswapdata.F90
! File: psi_iswapdata.F90
!
! Subroutine: psi_Xswapdatam
! Subroutine: psi_iswapdatam
! Does the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
@@ -68,7 +68,7 @@
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:,:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
@@ -503,9 +503,8 @@ end subroutine psi_iswapdatam
!!$
!!$
!
! File: psi_Xswapdata.F90
!
! Subroutine: psi_Xswapdatav
! Subroutine: psi_iswapdatav
! Does the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
@@ -542,7 +541,7 @@ end subroutine psi_iswapdatam
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
base/internals/psi_iswaptran.F90
+5 -6
View File
@@ -29,9 +29,9 @@
!!$
!!$
!
! File: psi_Xswaptran.F90
! File: psi_iswaptran.F90
!
! Subroutine: psi_Xswaptranm
! Subroutine: psi_iswaptranm
! Does the data exchange among processes. This is similar to Xswapdata, but
! the list is read "in reverse", i.e. indices that are normally SENT are used
! for the RECEIVE part and vice-versa. This is the basic data exchange operation
@@ -72,7 +72,7 @@
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:,:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
@@ -497,9 +497,8 @@ end subroutine psi_iswaptranm
!!$
!!$
!
! File: psi_Xswaptran.F90
!
! Subroutine: psi_Xswaptranv
! Subroutine: psi_iswaptranv
! Does the data exchange among processes. This is similar to Xswapdata, but
! the list is read "in reverse", i.e. indices that are normally SENT are used
! for the RECEIVE part and vice-versa. This is the basic data exchange operation
@@ -540,7 +539,7 @@ end subroutine psi_iswaptranm
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
base/internals/psi_ldsc_pre_halo.f90
+1 -2
View File
@@ -43,7 +43,7 @@
!
!
! Arguments:
! desc - type(<psb_desc_type>). The communication descriptor.
! desc - type(psb_desc_type). The communication descriptor.
! ext_hv - logical Should we work on the halo_index.
! info - integer. return code.
!
@@ -66,7 +66,6 @@ subroutine psi_ldsc_pre_halo(desc,ext_hv,info)
integer :: i,j,np,me,lhalo,nhalo,&
& n_col, err_act, key, ih, nh, idx, nk,icomm
integer :: ictxt,n_row
logical, parameter :: debug=.false., debugwrt=.false.
character(len=20) :: name,ch_err
info = 0
base/internals/psi_sort_dl.f90
+9 -5
View File
@@ -41,20 +41,23 @@ subroutine psi_sort_dl(dep_list,l_dep_list,np,info)
integer :: np,dep_list(:,:), l_dep_list(:)
integer :: idg, iupd, idgp, iedges, iidx, iich,ndgmx, isz, err_act
integer :: i, info
integer, allocatable :: work(:)
logical, parameter :: debug=.false.
character(len=20) :: name
integer, allocatable :: work(:)
integer :: debug_level, debug_unit
character(len=20) :: name
name='psi_sort_dl'
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = 0
ndgmx = 0
do i=1,np
ndgmx = ndgmx + l_dep_list(i)
if (debug) write(0,*) i,l_dep_list(i)
if (debug_level >= psb_debug_inner_)&
& write(debug_unit,*) name,': ',i,l_dep_list(i)
enddo
idg = 1
iupd = idg+np
@@ -63,7 +66,8 @@ subroutine psi_sort_dl(dep_list,l_dep_list,np,info)
iidx = iedges + 2*ndgmx
iich = iidx + ndgmx
isz = iich + ndgmx
if (debug)write(0,*) 'psi_sort_dl: ndgmx ',ndgmx,isz
if (debug_level >= psb_debug_inner_)&
& write(debug_unit,*) name,': ndgmx ',ndgmx,isz
allocate(work(isz))
! call srtlist(dep_list, dl_lda, l_dep_list, np, info)
base/internals/psi_zswapdata.F90
+5 -6
View File
@@ -29,9 +29,9 @@
!!$
!!$
!
! File: psi_Xswapdata.F90
! File: psi_zswapdata.F90
!
! Subroutine: psi_Xswapdatam
! Subroutine: psi_zswapdatam
! Does the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
@@ -68,7 +68,7 @@
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:,:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
@@ -503,9 +503,8 @@ end subroutine psi_zswapdatam
!!$
!!$
!
! File: psi_Xswapdata.F90
!
! Subroutine: psi_Xswapdatav
! Subroutine: psi_zswapdatav
! Does the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
@@ -542,7 +541,7 @@ end subroutine psi_zswapdatam
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
base/internals/psi_zswaptran.F90
+5 -6
View File
@@ -29,9 +29,9 @@
!!$
!!$
!
! File: psi_Xswaptran.F90
! File: psi_zswaptran.F90
!
! Subroutine: psi_Xswaptranm
! Subroutine: psi_zswaptranm
! Does the data exchange among processes. This is similar to Xswapdata, but
! the list is read "in reverse", i.e. indices that are normally SENT are used
! for the RECEIVE part and vice-versa. This is the basic data exchange operation
@@ -72,7 +72,7 @@
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:,:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
@@ -498,9 +498,8 @@ end subroutine psi_zswaptranm
!!$
!!$
!
! File: psi_Xswaptran.F90
!
! Subroutine: psi_Xswaptranv
! Subroutine: psi_zswaptranv
! Does the data exchange among processes. This is similar to Xswapdata, but
! the list is read "in reverse", i.e. indices that are normally SENT are used
! for the RECEIVE part and vice-versa. This is the basic data exchange operation
@@ -541,7 +540,7 @@ end subroutine psi_zswaptranm
! n - integer Number of columns in Y
! beta - X Choose overwrite or sum.
! y(:) - X The data area
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - X Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
base/modules/psb_desc_type.f90
+36 -26
View File
@@ -29,12 +29,10 @@
!!$
!!$
!
! Module to define desc_a,
! structure for coomunications.
!
! Typedef: psb_desc_type
! package: psb_descriptor_type
! Defines a communication descriptor
!
module psb_descriptor_type
use psb_const_mod
@@ -98,7 +96,23 @@ module psb_descriptor_type
!
! DESC data structure.
! type: psb_desc_type
!
! Communication Descriptor data structure.
!| type psb_desc_type
!| integer, allocatable :: matrix_data(:)
!| integer, allocatable :: halo_index(:), ext_index(:)
!| integer, allocatable :: bnd_elem(:)
!| integer, allocatable :: ovrlap_index(:)
!| integer, allocatable :: ovrlap_elem(:)
!| integer, allocatable :: loc_to_glob(:)
!| integer, allocatable :: glob_to_loc (:)
!| integer, allocatable :: hashv(:), glb_lc(:,:), ptree(:)
!| integer, allocatable :: lprm(:)
!| integer, allocatable :: idx_space(:)
!| end type psb_desc_type
!
!
! This is the most important data structure: it holds all the data
! necessary to organize data exchange. The pattern of communication
! among processes depends not only on the allocation of portions of
@@ -174,22 +188,22 @@ module psb_descriptor_type
! risk a deadlock. NOTE: This is the format when the state is PSB_ASB_.
! See below for BLD. The end-of-list is marked with a -1.
!
! notation stored in explanation
! --------------- --------------------------- -----------------------------------
! process_id index_v(p+proc_id_) identifier of process with which
! data is exchanged.
! n_elements_recv index_v(p+n_elem_recv_) number of elements to receive.
! elements_recv index_v(p+elem_recv_+i) indexes of local elements to
! receive. these are stored in the
! array from location p+elem_recv_ to
! location p+elem_recv_+
! index_v(p+n_elem_recv_)-1.
! n_elements_send index_v(p+n_elem_send_) number of elements to send.
! elements_send index_v(p+elem_send_+i) indexes of local elements to
! send. these are stored in the
! array from location p+elem_send_ to
! location p+elem_send_+
! index_v(p+n_elem_send_)-1.
!| notation stored in explanation
!| --------------- --------------------------- -----------------------------------
!| process_id index_v(p+proc_id_) identifier of process with which
!| data is exchanged.
!| n_elements_recv index_v(p+n_elem_recv_) number of elements to receive.
!| elements_recv index_v(p+elem_recv_+i) indexes of local elements to
!| receive. these are stored in the
!| array from location p+elem_recv_ to
!| location p+elem_recv_+
!| index_v(p+n_elem_recv_)-1.
!| n_elements_send index_v(p+n_elem_send_) number of elements to send.
!| elements_send index_v(p+elem_send_+i) indexes of local elements to
!| send. these are stored in the
!| array from location p+elem_send_ to
!| location p+elem_send_+
!| index_v(p+n_elem_send_)-1.
!
! This organization is valid for both halo and overlap indices; overlap entries
! need to be updated to ensure that a variable at a given global index
@@ -234,9 +248,7 @@ module psb_descriptor_type
!
! 10. ovrlap_elem contains a list of overlap indices together with their degree
! of overlap, i.e. the number of processes "owning" them.
!
!
! Yes, it is complex, but it does the following:
! It is complex, but it does the following:
! 1. Allows a purely local matrix/stencil buildup phase, requiring only
! one synch point at the end (CDASB)
! 2. Takes shortcuts when the problem size is not too large (the default threshold
@@ -248,7 +260,6 @@ module psb_descriptor_type
!
!
!
type psb_desc_type
integer, allocatable :: matrix_data(:)
integer, allocatable :: halo_index(:), ext_index(:)
@@ -278,7 +289,6 @@ contains
Type(psb_desc_type), intent(in) :: desc
Integer :: psb_cd_sizeof
!locals
logical, parameter :: debug=.false.
integer :: val
integer, external :: SizeofPairSearchTree
base/modules/psb_error_mod.F90
+75 -12
View File
@@ -31,13 +31,19 @@
module psb_error_mod
integer, parameter, public :: psb_act_ret_=0, psb_act_abort_=1, psb_no_err_=0
integer, parameter, public :: psb_debug_ext_=1, psb_debug_outer_=2
integer, parameter, public :: psb_debug_comp_=3, psb_debug_inner_=4
integer, parameter, public :: psb_debug_serial_=8, psb_debug_serial_comp_=9
!
! Error handling
!
public psb_errpush, psb_error, psb_get_errstatus,&
& psb_get_errverbosity, psb_set_errverbosity,psb_errcomm, &
& psb_erractionsave, psb_erractionrestore, &
& psb_get_erraction, psb_set_erraction
& psb_get_erraction, psb_set_erraction, &
& psb_get_debug_level, psb_set_debug_level,&
& psb_get_debug_unit, psb_set_debug_unit,&
& psb_get_serial_debug_level, psb_set_serial_debug_level
interface psb_error
module procedure psb_serror
@@ -49,29 +55,37 @@ module psb_error_mod
type psb_errstack_node
integer :: err_code=0 ! the error code
character(len=20) :: routine='' ! the name of the routine generating the error
integer,dimension(5) :: i_err_data=0 ! array of integer data to complete the error msg
! the error code
integer :: err_code=0
! the name of the routine generating the error
character(len=20) :: routine=''
! array of integer data to complete the error msg
integer,dimension(5) :: i_err_data=0
! real(kind(1.d0))(dim=10) :: r_err_data=0.d0 ! array of real data to complete the error msg
! complex(dim=10) :: c_err_data=0.c0 ! array of complex data to complete the error msg
character(len=40) :: a_err_data='' ! array of character data to complete the error msg
type(psb_errstack_node), pointer :: next ! pointer to the next element in the stack
! array of character data to complete the error msg
character(len=40) :: a_err_data=''
! pointer to the next element in the stack
type(psb_errstack_node), pointer :: next
end type psb_errstack_node
type psb_errstack
type(psb_errstack_node), pointer :: top => null() ! pointer to the top element of the stack
integer :: n_elems=0 ! number of entries in the stack
! pointer to the top element of the stack
type(psb_errstack_node), pointer :: top => null()
! number of entries in the stack
integer :: n_elems=0
end type psb_errstack
type(psb_errstack),save :: error_stack ! the PSBLAS-2.0 error stack
integer,save :: error_status=0 ! the error status (maybe not here)
integer,save :: verbosity_level=1 ! the verbosity level (maybe not here)
integer,save :: err_action=1
type(psb_errstack), save :: error_stack ! the PSBLAS-2.0 error stack
integer, save :: error_status=0 ! the error status (maybe not here)
integer, save :: verbosity_level=1 ! the verbosity level (maybe not here)
integer, save :: err_action=1
integer, save :: debug_level=0, debug_unit=0, serial_debug_level=0
contains
@@ -105,6 +119,49 @@ contains
end subroutine psb_erractionrestore
function psb_get_debug_level()
integer :: psb_get_debug_level
psb_get_debug_level = debug_level
end function psb_get_debug_level
subroutine psb_set_debug_level(level)
integer, intent(in) :: level
if (level >= 0) then
debug_level = level
else
debug_level = 0
end if
end subroutine psb_set_debug_level
function psb_get_debug_unit()
integer :: psb_get_debug_unit
psb_get_debug_unit = debug_unit
end function psb_get_debug_unit
subroutine psb_set_debug_unit(unit)
integer, intent(in) :: unit
if (unit >= 0) then
debug_unit = unit
else
debug_unit = 0
end if
end subroutine psb_set_debug_unit
function psb_get_serial_debug_level()
integer :: psb_get_serial_debug_level
psb_get_serial_debug_level = serial_debug_level
end function psb_get_serial_debug_level
subroutine psb_set_serial_debug_level(level)
integer, intent(in) :: level
if (level >= 0) then
serial_debug_level = level
else
serial_debug_level = 0
end if
end subroutine psb_set_serial_debug_level
! checks wether an error has occurred on one of the porecesses in the execution pool
subroutine psb_errcomm(ictxt, err)
integer, intent(in) :: ictxt
@@ -452,10 +509,16 @@ contains
write (0,'("local index is: ",i0," and global index is:",i0)')i_e_d(1:2)
case(3110)
write (0,'("Before you call this routine, you must assembly sparse matrix")')
case(3111)
write (0,'("Before you call this routine, you must initialize the preconditioner")')
case(3112)
write (0,'("Before you call this routine, you must build the preconditioner")')
case(3111:3999)
write(0,'("miscellaneus error. code: ",i0)')err_c
case(4000)
write(0,'("Allocation/deallocation error")')
case(4001)
write(0,'("Internal error: ",a)')a_e_d
case(4010)
write (0,'("Error from call to subroutine ",a)')a_e_d
case(4011)
base/modules/psb_penv_mod.F90
+1 -34
View File
@@ -72,7 +72,7 @@ module psb_penv_mod
module procedure psb_ibcasts, psb_ibcastv, psb_ibcastm,&
& psb_dbcasts, psb_dbcastv, psb_dbcastm,&
& psb_zbcasts, psb_zbcastv, psb_zbcastm,&
& psb_hbcasts, psb_lbcasts, psb_lbcastv, psb_hbcastv1
& psb_hbcasts, psb_lbcasts, psb_lbcastv
end interface
@@ -546,39 +546,6 @@ contains
end subroutine psb_hbcasts
subroutine psb_hbcastv1(ictxt,dat,root,length)
#ifdef MPI_H
include 'mpif.h'
#endif
#ifdef MPI_MOD
use mpi
#endif
integer, intent(in) :: ictxt
character(len=1), intent(inout) :: dat(:)
integer, intent(in), optional :: root,length
integer :: iam, np, root_,icomm,length_,info
#if !defined(SERIAL_MPI)
if (present(root)) then
root_ = root
else
root_ = psb_root_
endif
if (present(length)) then
length_ = length
else
length_ = size(dat)
endif
call psb_info(ictxt,iam,np)
call psb_get_mpicomm(ictxt,icomm)
call mpi_bcast(dat,length_,MPI_CHARACTER,root_,icomm,info)
#endif
end subroutine psb_hbcastv1
subroutine psb_lbcasts(ictxt,dat,root)
#ifdef MPI_H
include 'mpif.h'
base/modules/psb_spmat_type.f90
+61 -93
View File
@@ -28,18 +28,16 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Module to define D_SPMAT, structure !!
!! for sparse matrix. !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! package: psb_spmat_type
! Data structure(s) for sparse matrices
!
module psb_spmat_type
use psb_error_mod
use psb_realloc_mod
use psb_const_mod
implicit none
! Typedef: psb_dspmat_type
! Contains a sparse matrix
!
!
! Queries into spmat%info
!
@@ -55,22 +53,40 @@ module psb_spmat_type
integer, parameter :: psb_dupl_=11, psb_upd_=12
integer, parameter :: psb_ifasize_=16
!
!
! Possible matrix states.
! Types: psb_dspmat_type, psb_zspmat_type
!
!| type psb_dspmat_type
!| integer :: m, k ! Row & column size
!| character(len=5) :: fida ! Storage format: CSR,COO etc.
!| character(len=11) :: descra ! Matrix type: encodes general, triang.
!| integer :: infoa(psb_ifasize_) ! Additional integer info
!| real(kind(1.d0)), allocatable :: aspk(:) ! Coefficients
!| integer, allocatable :: ia1(:), ia2(:) ! Row/column indices encoded
!| integer, allocatable :: pl(:), pr(:) ! Row/column permutation
!| end type psb_dspmat_type
!| type psb_zspmat_type
!| integer :: m, k
!| character(len=5) :: fida
!| character(len=11) :: descra
!| integer :: infoa(psb_ifasize_)
!| complex(kind(1.d0)), allocatable :: aspk(:)
!| integer, allocatable :: ia1(:), ia2(:)
!| integer, allocatable :: pl(:), pr(:)
!| end type psb_zspmat_type
!
! A sparse matrix can move between states according to the
! following state transition table.
! In Out Routine
! ----------------------------------
! Null Build psb_sp_all
! Build Build psb_coins
! Build Assembled psb_spcnv
! Assembled Assembled psb_spcnv
! Assembled Update psb_sp_reinit
! Update Update psb_coins
! Update Assembled psb_spcnv
! * unchanged psb_sp_reall
! Assembled Null psb_sp_free
!| In Out Routine
!| ----------------------------------
!| Null Build psb_sp_all
!| Build Build psb_coins
!| Build Assembled psb_spcnv
!| Assembled Assembled psb_spcnv
!| Assembled Update psb_sp_reinit
!| Update Update psb_coins
!| Update Assembled psb_spcnv
!| * unchanged psb_sp_reall
!| Assembled Null psb_sp_free
!
! Note that psb_spcnv is overloaded in two flavours,
! psb_spcnv(a,info)
@@ -81,7 +97,7 @@ module psb_spmat_type
!
integer, parameter :: psb_spmat_null_=0, psb_spmat_bld_=1
integer, parameter :: psb_spmat_asb_=2, psb_spmat_upd_=4
integer, parameter :: psb_ireg_flgs_=10, psb_ip2_=0
integer, parameter :: psb_iflag_=2, psb_ichk_=3
integer, parameter :: psb_nnzt_=4, psb_zero_=5,psb_ipc_=6
@@ -105,35 +121,21 @@ module psb_spmat_type
type psb_dspmat_type
! Rows & columns
integer :: m, k
! Identify the representation method. Es: CSR, JAD, ...
character(len=5) :: fida
! describe some chacteristics of sparse matrix
character(len=11) :: descra
! Contains some additional informations on sparse matrix
integer :: infoa(psb_ifasize_)
! Contains sparse matrix coefficients
real(kind(1.d0)), allocatable :: aspk(:)
! Contains indeces that describes sparse matrix structure
integer, allocatable :: ia1(:), ia2(:)
! Permutations matrix
integer, allocatable :: pl(:), pr(:)
end type psb_dspmat_type
type psb_zspmat_type
! Rows & columns
integer :: m, k
! Identify the representation method. Es: CSR, JAD, ...
character(len=5) :: fida
! describe some chacteristics of sparse matrix
character(len=11) :: descra
! Contains some additional informations on sparse matrix
integer :: infoa(psb_ifasize_)
! Contains sparse matrix coefficients
complex(kind(1.d0)), allocatable :: aspk(:)
! Contains indeces that describes sparse matrix structure
integer, allocatable :: ia1(:), ia2(:)
! Permutations matrix
integer, allocatable :: pl(:), pr(:)
end type psb_zspmat_type
@@ -243,12 +245,12 @@ contains
return
end function psb_get_zsp_ncols
integer function psb_get_dsp_nnzeros(a)
type(psb_dspmat_type), intent(in) :: a
integer :: ires,info
call psb_sp_info(psb_nztotreq_,a,ires,info)
if (info == 0) then
psb_get_dsp_nnzeros = ires
@@ -260,7 +262,7 @@ contains
integer function psb_get_zsp_nnzeros(a)
type(psb_zspmat_type), intent(in) :: a
integer :: ires,info
call psb_sp_info(psb_nztotreq_,a,ires,info)
if (info == 0) then
psb_get_zsp_nnzeros = ires
@@ -272,7 +274,7 @@ contains
integer function psb_get_dsp_nzsize(a)
type(psb_dspmat_type), intent(in) :: a
integer :: ires,info
call psb_sp_info(psb_nzsizereq_,a,ires,info)
if (info == 0) then
psb_get_dsp_nzsize = ires
@@ -284,7 +286,7 @@ contains
integer function psb_get_zsp_nzsize(a)
type(psb_zspmat_type), intent(in) :: a
integer :: ires,info
call psb_sp_info(psb_nzsizereq_,a,ires,info)
if (info == 0) then
psb_get_zsp_nzsize = ires
@@ -298,7 +300,7 @@ contains
integer, intent(in) :: ir
type(psb_dspmat_type), intent(in) :: a
integer :: ires,info
call psb_sp_info(psb_nzrowreq_,a,ires,info,iaux=ir)
if (info == 0) then
psb_get_dsp_nnz_row = ires
@@ -310,7 +312,7 @@ contains
integer, intent(in) :: ir
type(psb_zspmat_type), intent(in) :: a
integer :: ires,info
call psb_sp_info(psb_nzrowreq_,a,ires,info,iaux=ir)
if (info == 0) then
psb_get_zsp_nnz_row = ires
@@ -347,7 +349,7 @@ contains
logical, parameter :: debug=.false.
logical :: clear_
character(len=20) :: name
info = 0
name = 'psb_sp_reinit'
@@ -356,7 +358,7 @@ contains
else
clear_ = .true.
end if
select case(psb_sp_getifld(psb_state_,a,info))
case(psb_spmat_asb_)
@@ -628,7 +630,7 @@ contains
End Subroutine psb_dspclone
! Will be changed to use MOVE_ALLOC
subroutine psb_dsp_transfer(a, b,info)
implicit none
@@ -637,9 +639,6 @@ contains
Type(psb_dspmat_type), intent(inout) :: B
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
@@ -668,16 +667,13 @@ contains
Type(psb_dspmat_type), intent(inout) :: A
Integer, intent(in) :: field,val
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
!!$ call psb_realloc(psb_ifasize_,a%infoa,info)
if (info == 0) &
& call psb_setifield(val,field,a%infoa,psb_ifasize_,info)
Return
@@ -699,7 +695,6 @@ contains
!locals
Integer :: i1, i2, ia
logical, parameter :: debug=.false.
info = 0
call psb_sp_trimsize(a,i1,i2,ia,info)
@@ -718,7 +713,6 @@ contains
Integer, intent(out) :: info
!locals
Integer :: i1, i2, ia
logical, parameter :: debug=.false.
info = 0
call psb_sp_trimsize(a,i1,i2,ia,info)
@@ -739,7 +733,6 @@ contains
!locals
Integer :: nza
logical, parameter :: debug=.false.
info = 0
if (psb_sp_getifld(psb_upd_,a,info) == psb_upd_perm_) then
@@ -775,7 +768,7 @@ contains
Return
End Subroutine psb_dsp_trimsize
function psb_dsp_getifld(field,a,info)
implicit none
!....Parameters...
@@ -786,7 +779,6 @@ contains
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
integer :: val
info = 0
@@ -799,7 +791,7 @@ contains
endif
call psb_getifield(val,field,a%infoa,psb_ifasize_,info)
psb_dsp_getifld = val
Return
@@ -813,11 +805,10 @@ contains
Integer :: psb_dspsizeof
!locals
logical, parameter :: debug=.false.
integer :: val
val = 4*size(a%infoa)
if (allocated(a%aspk)) then
val = val + 8 * size(a%aspk)
endif
@@ -835,7 +826,7 @@ contains
val = val + 4 * size(a%pr)
endif
psb_dspsizeof = val
Return
@@ -848,7 +839,6 @@ contains
Type(psb_dspmat_type), intent(inout) :: A
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
integer :: iret
info = 0
@@ -902,10 +892,9 @@ contains
logical, intent(in), optional :: clear
!locals
logical, parameter :: debug=.false.
logical :: clear_
character(len=20) :: name
info = 0
name = 'psb_sp_reinit'
@@ -914,7 +903,7 @@ contains
else
clear_ = .true.
end if
select case(psb_sp_getifld(psb_state_,a,info))
case(psb_spmat_asb_)
@@ -1078,8 +1067,6 @@ contains
Integer, intent(in) :: ni1,ni2,nz
Integer, intent(inout) :: info
!locals
logical, parameter :: debug=.false.
info = 0
call psb_realloc(nz,a%aspk,info)
@@ -1113,9 +1100,6 @@ contains
Integer, intent(in), optional :: ifc
integer :: ifc_
!locals
logical, parameter :: debug=.false.
info = 0
if (nnz.lt.0) then
info=45
@@ -1154,10 +1138,6 @@ contains
Type(psb_zspmat_type), intent(out) :: B
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
INFO = 0
call psb_nullify_sp(b)
@@ -1191,9 +1171,6 @@ contains
Type(psb_zspmat_type), intent(inout) :: B
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
call psb_transfer( a%aspk, b%aspk , info)
@@ -1221,17 +1198,13 @@ contains
Integer, intent(in) :: field,val
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
!!$ call psb_realloc(psb_ifasize_,a%infoa,info)
if (info == 0) &
& call psb_setifield(val,field,a%infoa,psb_ifasize_,info)
Return
@@ -1247,7 +1220,6 @@ contains
!locals
Integer :: nza
logical, parameter :: debug=.false.
info = 0
if (psb_sp_getifld(psb_upd_,a,info) == psb_upd_perm_) then
@@ -1294,7 +1266,6 @@ contains
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
integer :: val
info = 0
@@ -1307,7 +1278,7 @@ contains
endif
call psb_getifield(val,field,a%infoa,psb_ifasize_,info)
psb_zsp_getifld = val
Return
@@ -1321,11 +1292,10 @@ contains
Integer :: psb_zspsizeof
!locals
logical, parameter :: debug=.false.
integer :: val
val = 4*size(a%infoa)
if (allocated(a%aspk)) then
val = val + 16 * size(a%aspk)
endif
@@ -1342,8 +1312,8 @@ contains
if (allocated(a%pr)) then
val = val + 4 * size(a%pr)
endif
psb_zspsizeof = val
Return
@@ -1357,8 +1327,6 @@ contains
!....Parameters...
Type(psb_zspmat_type), intent(inout) :: A
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
base/modules/psi_mod.f90
+12 -9
View File
@@ -316,11 +316,14 @@ contains
integer, allocatable :: idx_out(:)
! ...parameters
integer :: debug_level, debug_unit
logical, parameter :: debug=.false.
character(len=20) :: name
name='psi_bld_cdesc'
call psb_get_erraction(err_act)
debug_level = psb_get_debug_level()
debug_unit = psb_get_debug_unit()
info = 0
ictxt = cdesc%matrix_data(psb_ctxt_)
@@ -334,7 +337,7 @@ contains
! first the halo index
if (debug) write(0,*) me,'Calling crea_index on halo'
if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on halo'
call psi_crea_index(cdesc,halo_in, idx_out,.false.,nxch,nsnd,nrcv,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
@@ -345,12 +348,12 @@ contains
cdesc%matrix_data(psb_thal_snd_) = nsnd
cdesc%matrix_data(psb_thal_rcv_) = nrcv
if (debug) write(0,*) me,'Done crea_index on halo'
if (debug) write(0,*) me,'Calling crea_index on ext'
if (debug_level>0) write(debug_unit,*) me,'Done crea_index on halo'
if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on ext'
! then ext index
if (debug) write(0,*) me,'Calling crea_index on ext'
if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on ext'
call psi_crea_index(cdesc,ext_in, idx_out,.false.,nxch,nsnd,nrcv,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
@@ -361,8 +364,8 @@ contains
cdesc%matrix_data(psb_text_snd_) = nsnd
cdesc%matrix_data(psb_text_rcv_) = nrcv
if (debug) write(0,*) me,'Done crea_index on ext'
if (debug) write(0,*) me,'Calling crea_index on ovrlap'
if (debug_level>0) write(debug_unit,*) me,'Done crea_index on ext'
if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on ovrlap'
! then the overlap index
@@ -381,10 +384,10 @@ contains
cdesc%matrix_data(psb_tovr_snd_) = nsnd
cdesc%matrix_data(psb_tovr_rcv_) = nrcv
if (debug) write(0,*) me,'Calling crea_ovr_elem'
if (debug_level>0) write(debug_unit,*) me,'Calling crea_ovr_elem'
! next ovrlap_elem
call psi_crea_ovr_elem(cdesc%ovrlap_index,cdesc%ovrlap_elem,info)
if (debug) write(0,*) me,'Done crea_ovr_elem'
if (debug_level>0) write(debug_unit,*) me,'Done crea_ovr_elem'
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_ovr_elem')
goto 9999
@@ -398,7 +401,7 @@ contains
call psb_errpush(4010,name,a_err='psi_crea_bnd_elem')
goto 9999
end if
if (debug) write(0,*) me,'Done crea_bnd_elem'
if (debug_level>0) write(debug_unit,*) me,'Done crea_bnd_elem'
call psb_erractionrestore(err_act)
return
base/psblas/psb_damax.f90
+4 -4
View File
@@ -39,7 +39,7 @@
!
! Arguments:
! x - real,dimension(:,:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
@@ -166,7 +166,7 @@ end function psb_damax
!
! Arguments:
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_damaxv (x,desc_a, info)
@@ -288,7 +288,7 @@ end function psb_damaxv
! Arguments:
! res - real. The result.
! x - real,dimension(:,:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
@@ -410,7 +410,7 @@ end subroutine psb_damaxvs
! Arguments:
! res - real. The result.
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_dmamaxs (res,x,desc_a, info,jx)
base/psblas/psb_dasum.f90
+3 -3
View File
@@ -39,7 +39,7 @@
!
! Arguments:
! x - real,dimension(:,:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
@@ -184,7 +184,7 @@ end function psb_dasum
!
! Arguments:
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_dasumv (x,desc_a, info)
@@ -322,7 +322,7 @@ end function psb_dasumv
! Arguments:
! res - real. The result.
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
base/psblas/psb_daxpby.f90
+2 -3
View File
@@ -44,7 +44,7 @@
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! beta - real. The scalar used to multiply each component of sub( Y ).
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ).
! jy - integer(optional). The column offset for sub( Y ).
@@ -199,7 +199,7 @@ end subroutine psb_daxpby
! x - real,dimension(:). The input vector containing the entries of X.
! beta - real. The scalar used to multiply each component of Y.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
@@ -220,7 +220,6 @@ subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
integer :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, m, iiy, jjy
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='psb_dgeaxpby'
if(psb_get_errstatus().ne.0) return
base/psblas/psb_ddot.f90
+4 -4
View File
@@ -42,7 +42,7 @@
! Arguments:
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ).
! jy - integer(optional). The column offset for sub( Y ).
@@ -199,7 +199,7 @@ end function psb_ddot
! Arguments:
! x - real,dimension(:). The input vector containing the entries of X.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_ddotv(x, y,desc_a, info)
@@ -337,7 +337,7 @@ end function psb_ddotv
! res - real. The result.
! x - real,dimension(:). The input vector containing the entries of X.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_ddotvs(res, x, y,desc_a, info)
@@ -478,7 +478,7 @@ end subroutine psb_ddotvs
! res - real. The result.
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_dmdots(res, x, y, desc_a, info)
base/psblas/psb_dnrm2.f90
+3 -3
View File
@@ -39,7 +39,7 @@
!
! Arguments:
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ).
!
@@ -178,7 +178,7 @@ end function psb_dnrm2
!
! Arguments:
! x - real,dimension(:). The input vector containing the entries of X.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_dnrm2v(x, desc_a, info)
@@ -313,7 +313,7 @@ end function psb_dnrm2v
! Arguments:
! res - real. The result.
! x - real,dimension(:). The input vector containing the entries of X.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_dnrm2vs(res, x, desc_a, info)
base/psblas/psb_dnrmi.f90
+2 -2
View File
@@ -36,8 +36,8 @@
! normi := max(abs(sum(A(i,j))))
!
! Arguments:
! a - type(<psb_dspmat_type>). The sparse matrix containing A.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! a - type(psb_dspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_dnrmi(a,desc_a,info)
base/psblas/psb_dspmm.f90
+24 -16
View File
@@ -64,11 +64,11 @@
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_dspmat_type>). The sparse matrix containing A.
! a - type(psb_dspmat_type). The sparse matrix containing A.
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! beta - real. The scalar beta.
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! k - integer(optional). The number of right-hand sides.
@@ -404,11 +404,11 @@ end subroutine psb_dspmm
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_dspmat_type>). The sparse matrix containing A.
! a - type(psb_dspmat_type). The sparse matrix containing A.
! x - real,dimension(:). The input vector containing the entries of X.
! beta - real. The scalar beta.
! y - real,dimension(:. The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! work - real,dimension(:)(optional). Working area.
@@ -449,15 +449,16 @@ 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.
integer :: debug_level, debug_unit
name='psb_dspmv'
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
if (np == -1) then
info = 2010
@@ -509,7 +510,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
liwork= 2*ncol
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 (size(work) >= liwork) then
aliw =.false.
@@ -520,7 +521,6 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
aliw=.true.
end if
aliw=.true.
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= 0) then
@@ -533,7 +533,8 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
iwork => work
endif
if (debug) write(0,*) me,name,' Allocated work ', info
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' Allocated work ', info
! checking for matrix correctness
call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
if(info /= 0) then
@@ -543,7 +544,8 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
goto 9999
end if
if (debug) write(0,*) me,name,' Checkmat ', info
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' Checkmat ', info
if (itrans == 'N') then
! Matrix is not transposed
if((ja /= ix).or.(ia /= iy)) then
@@ -624,10 +626,12 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
yp => y(iiy:lldy)
yp(nrow+1:ncol)=dzero
if (debug) write(0,*) me,name,' checkvect ', info
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' checkvect ', info
! local Matrix-vector product
call psb_csmm(alpha,a,xp,beta,yp,info,trans=itrans)
if (debug) write(0,*) me,name,' csmm ', info
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info
if(info /= 0) then
info = 4010
ch_err='dcsmm'
@@ -638,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 (debug) write(0,*) me,name,' swaptran ', info
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' swaptran ', info
if(info /= 0) then
info = 4010
ch_err='PSI_dSwapTran'
@@ -649,7 +654,8 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
end if
if (aliw) deallocate(iwork,stat=info)
if (debug) write(0,*) me,name,' deallocat ',aliw, info
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' deallocat ',aliw, info
if(info /= 0) then
info = 4010
ch_err='Deallocate iwork'
@@ -660,8 +666,10 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
nullify(iwork)
call psb_erractionrestore(err_act)
if (debug) call psb_barrier(ictxt)
if (debug) write(0,*) me,name,' Returning '
if (debug_level >= psb_debug_comp_) then
call psb_barrier(ictxt)
write(debug_unit,*) me,' ',trim(name),' Returning '
endif
return
9999 continue
base/psblas/psb_dspsm.f90
+4 -4
View File
@@ -56,11 +56,11 @@
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_dspmat_type>). The sparse matrix containing A.
! a - type(psb_dspmat_type). The sparse matrix containing A.
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! beta - real. The scalar beta.
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with the diagonal matrix D.
@@ -366,11 +366,11 @@ end subroutine psb_dspsm
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_dspmat_type>). The sparse matrix containing A.
! a - type(psb_dspmat_type). The sparse matrix containing A.
! x - real,dimension(:). The input vector containing the entries of X.
! beta - real. The scalar beta.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with the diagonal matrix D.
base/psblas/psb_zamax.f90
+7 -7
View File
@@ -28,9 +28,9 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
! File: psb_damax.f90
! File: psb_zamax.f90
!
! Function: psb_damax
! Function: psb_zamax
! Searches the absolute max of X.
!
! normi := max(abs(sub(X)(i))
@@ -38,8 +38,8 @@
! where sub( X ) denotes X(1:N,JX:).
!
! Arguments:
! x - real,dimension(:,:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! x - complex,dimension(:,:). The input vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
@@ -169,7 +169,7 @@ end function psb_zamax
!
! Arguments:
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_zamaxv (x,desc_a, info)
@@ -296,7 +296,7 @@ end function psb_zamaxv
! Arguments:
! res - real. The result.
! x - real,dimension(:,:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
@@ -422,7 +422,7 @@ end subroutine psb_zamaxvs
! Arguments:
! res - real. The result.
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_zmamaxs (res,x,desc_a, info,jx)
base/psblas/psb_zasum.f90
+3 -3
View File
@@ -39,7 +39,7 @@
!
! Arguments:
! x - real,dimension(:,:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
@@ -189,7 +189,7 @@ end function psb_zasum
!
! Arguments:
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_zasumv (x,desc_a, info)
@@ -333,7 +333,7 @@ end function psb_zasumv
! Arguments:
! res - real. The result.
! x - real,dimension(:). The input vector.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset.
!
base/psblas/psb_zaxpby.f90
+18 -17
View File
@@ -28,9 +28,9 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
! File: psb_daxpby.f90
! File: psb_zaxpby.f90
!
! Subroutine: psb_daxpby
! Subroutine: psb_zaxpby
! Adds one distributed matrix to another,
!
! sub( Y ) := beta * sub( Y ) + alpha * sub( X )
@@ -40,14 +40,14 @@
! sub( Y ) denotes Y(:,JY).
!
! Arguments:
! alpha - real. The scalar used to multiply each component of sub( X ).
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! beta - real. The scalar used to multiply each component of sub( Y ).
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ).
! jy - integer(optional). The column offset for sub( Y ).
! alpha - complex,input The scalar used to multiply each component of X
! x(:,:) - complex,input The input vector containing the entries of X
! beta - real,input The scalar used to multiply each component of Y
! y(:,:) - real,inout The input vector Y
! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer Return code
! jx - integer(optional) The column offset for X
! jy - integer(optional) The column offset for Y
!
subroutine psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
use psb_descriptor_type
@@ -188,18 +188,19 @@ end subroutine psb_zaxpby
!!$
!!$
!
! Subroutine: psb_dgeaxpbyv
! Subroutine: psb_zgeaxpbyv
! Adds one distributed matrix to another,
!
! Y := beta * Y + alpha * X
!
! Arguments:
! alpha - real. The scalar used to multiply each component of X.
! x - real,dimension(:). The input vector containing the entries of X.
! beta - real. The scalar used to multiply each component of Y.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Return code
! alpha - complex,input The scalar used to multiply each component of X
! x(:) - complex,input The input vector containing the entries of X
! beta - real,input The scalar used to multiply each component of Y
! y(:) - real,inout The input vector Y
! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer Return code
!
!
subroutine psb_zaxpbyv(alpha, x, beta,y,desc_a,info)
use psb_descriptor_type
base/psblas/psb_zdot.f90
+4 -4
View File
@@ -42,7 +42,7 @@
! Arguments:
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ).
! jy - integer(optional). The column offset for sub( Y ).
@@ -198,7 +198,7 @@ end function psb_zdot
! Arguments:
! x - real,dimension(:). The input vector containing the entries of X.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_zdotv(x, y,desc_a, info)
@@ -336,7 +336,7 @@ end function psb_zdotv
! res - real. The result.
! x - real,dimension(:). The input vector containing the entries of X.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_zdotvs(res, x, y,desc_a, info)
@@ -476,7 +476,7 @@ end subroutine psb_zdotvs
! res - real. The result.
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_zmdots(res, x, y, desc_a, info)
base/psblas/psb_znrm2.f90
+3 -3
View File
@@ -39,7 +39,7 @@
!
! Arguments:
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ).
!
@@ -177,7 +177,7 @@ end function psb_znrm2
!
! Arguments:
! x - real,dimension(:). The input vector containing the entries of X.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_znrm2v(x, desc_a, info)
@@ -312,7 +312,7 @@ end function psb_znrm2v
! Arguments:
! res - real. The result.
! x - real,dimension(:). The input vector containing the entries of X.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
subroutine psb_znrm2vs(res, x, desc_a, info)
base/psblas/psb_znrmi.f90
+2 -2
View File
@@ -36,8 +36,8 @@
! normi := max(abs(sum(A(i,j))))
!
! Arguments:
! a - type(<psb_dspmat_type>). The sparse matrix containing A.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! a - type(psb_dspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_znrmi(a,desc_a,info)
base/psblas/psb_zspmm.f90
+128 -109
View File
@@ -64,11 +64,11 @@
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_zspmat_type>). The sparse matrix containing A.
! a - type(psb_zspmat_type). The sparse matrix containing A.
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! beta - real. The scalar beta.
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! k - integer(optional). The number of right-hand sides.
@@ -398,11 +398,11 @@ end subroutine psb_zspmm
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_zspmat_type>). The sparse matrix containing A.
! a - type(psb_zspmat_type). The sparse matrix containing A.
! x - real,dimension(:). The input vector containing the entries of X.
! beta - real. The scalar beta.
! y - real,dimension(:. The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! work - real,dimension(:)(optional). Working area.
@@ -443,11 +443,14 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
character :: itrans
character(len=20) :: name, ch_err
logical :: aliw
integer :: debug_level, debug_unit
name='psb_zspmv'
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
@@ -467,9 +470,9 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
ib = 1
if (present(doswap)) then
idoswap = doswap
idoswap = doswap
else
idoswap = 1
idoswap = 1
endif
if (present(trans)) then
@@ -481,7 +484,7 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
goto 9999
end if
else
itrans = 'N'
itrans = 'N'
endif
m = psb_cd_get_global_rows(desc_a)
@@ -496,18 +499,17 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
liwork= 2*ncol
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.
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= 0) then
@@ -517,136 +519,153 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
goto 9999
end if
else
iwork => work
iwork => work
endif
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' Allocated work ', info
! checking for matrix correctness
call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
if(info /= 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 (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' Checkmat ', info
if (itrans == 'N') then
! Matrix is not transposed
if((ja /= ix).or.(ia /= iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! Matrix is not transposed
if((ja /= ix).or.(ia /= 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),ix,jx,desc_a,info,iix,jjx)
if (info == 0) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy)
if(info /= 0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a,info,iix,jjx)
if (info == 0) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy)
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 /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if((iix /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if (idoswap == 0) then
x(nrow+1:ncol)=zzero
else
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& zzero,x,desc_a,iwork,info,data=psb_comm_halo_)
end if
if (idoswap == 0) then
x(nrow+1:ncol)=zzero
else
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& zzero,x,desc_a,iwork,info,data=psb_comm_halo_)
end if
! local Matrix-vector product
call psb_csmm(alpha,a,x(iix:lldx),beta,y(iiy:lldy),info)
! local Matrix-vector product
call psb_csmm(alpha,a,x(iix:lldx),beta,y(iiy:lldy),info)
if(info /= 0) then
info = 4011
call psb_errpush(info,name)
goto 9999
end if
if(info /= 0) then
info = 4011
call psb_errpush(info,name)
goto 9999
end if
else
! Matrix is transposed
if((ja /= iy).or.(ia /= ix)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! Matrix is transposed
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) /= -1) then
info = 3070
call psb_errpush(info,name)
goto 9999
end if
if(desc_a%ovrlap_elem(1) /= -1) then
info = 3070
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx)
if (info == 0)&
& call psb_chkvect(n,ik,size(y),iy,jy,desc_a,info,iiy,jjy)
if(info /= 0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx)
if (info == 0)&
& call psb_chkvect(n,ik,size(y),iy,jy,desc_a,info,iiy,jjy)
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 /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if((iix /= 1).or.(iiy /= 1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
xp => x(iix:lldx)
yp => y(iiy:lldy)
xp => x(iix:lldx)
yp => y(iiy:lldy)
yp(nrow+1:ncol)=zzero
yp(nrow+1:ncol)=zzero
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' checkvect ', info
! local Matrix-vector product
call psb_csmm(alpha,a,xp,beta,yp,info,trans=itrans)
! local Matrix-vector product
call psb_csmm(alpha,a,xp,beta,yp,info,trans=itrans)
if(info /= 0) then
info = 4010
ch_err='zcsmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(info /= 0) then
info = 4010
ch_err='dcsmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(idoswap /= 0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& zone,yp,desc_a,iwork,info)
if(info /= 0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(idoswap /= 0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& zone,yp,desc_a,iwork,info)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' swaptran ', info
if(info /= 0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
if(aliw) deallocate(iwork)
if (aliw) deallocate(iwork,stat=info)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(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_level >= psb_debug_comp_) then
call psb_barrier(ictxt)
write(debug_unit,*) me,' ',trim(name),' Returning '
endif
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
call psb_error(ictxt)
return
end if
return
end subroutine psb_zspmv
base/psblas/psb_zspsm.f90
+4 -4
View File
@@ -56,11 +56,11 @@
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_zspmat_type>). The sparse matrix containing A.
! a - type(psb_zspmat_type). The sparse matrix containing A.
! x - real,dimension(:,:). The input vector containing the entries of sub( X ).
! beta - real. The scalar beta.
! y - real,dimension(:,:). The input vector containing the entries of sub( Y ).
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with the diagonal matrix D.
@@ -369,11 +369,11 @@ end subroutine psb_zspsm
!
! Arguments:
! alpha - real. The scalar alpha.
! a - type(<psb_zspmat_type>). The sparse matrix containing A.
! a - type(psb_zspmat_type). The sparse matrix containing A.
! x - real,dimension(:). The input vector containing the entries of X.
! beta - real. The scalar beta.
! y - real,dimension(:). The input vector containing the entries of Y.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with the diagonal matrix D.
base/serial/aux/dasr.f90
-2
View File
@@ -69,7 +69,6 @@ subroutine dasr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -194,7 +193,6 @@ subroutine dasr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/dasrx.f90
-2
View File
@@ -81,7 +81,6 @@ subroutine dasrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -220,7 +219,6 @@ subroutine dasrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/dsr.f90
-2
View File
@@ -69,7 +69,6 @@ subroutine dsr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -194,7 +193,6 @@ subroutine dsr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/dsrx.f90
-2
View File
@@ -81,7 +81,6 @@ subroutine dsrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -221,7 +220,6 @@ subroutine dsrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/iasr.f90
-2
View File
@@ -69,7 +69,6 @@ subroutine iasr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -194,7 +193,6 @@ subroutine iasr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/iasrx.f90
-2
View File
@@ -81,7 +81,6 @@ subroutine iasrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -220,7 +219,6 @@ subroutine iasrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/imsru.f90
-1
View File
@@ -42,7 +42,6 @@ subroutine imsru(n,x,idir,nout)
nout = 0
if (n<0) then
!!$ write(0,*) 'Error: IMSR: N<0'
return
endif
base/serial/aux/isr.f90
-2
View File
@@ -69,7 +69,6 @@ subroutine isr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -194,7 +193,6 @@ subroutine isr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/isrx.f90
-2
View File
@@ -80,7 +80,6 @@ subroutine isrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -221,7 +220,6 @@ subroutine isrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/zalsr.f90
-2
View File
@@ -70,7 +70,6 @@ subroutine zalsr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -195,7 +194,6 @@ subroutine zalsr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/zalsrx.f90
-2
View File
@@ -82,7 +82,6 @@ subroutine zalsrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -223,7 +222,6 @@ subroutine zalsrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/zasr.f90
-2
View File
@@ -70,7 +70,6 @@ subroutine zasr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -195,7 +194,6 @@ subroutine zasr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/zasrx.f90
-2
View File
@@ -82,7 +82,6 @@ subroutine zasrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -223,7 +222,6 @@ subroutine zasrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/zlsr.f90
-2
View File
@@ -70,7 +70,6 @@ subroutine zlsr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -195,7 +194,6 @@ subroutine zlsr(n,x,dir)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/aux/zlsrx.f90
-2
View File
@@ -82,7 +82,6 @@ subroutine zlsrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
@@ -223,7 +222,6 @@ subroutine zlsrx(n,x,indx,dir,flag)
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
base/serial/coo/Makefile
+1 -1
View File
@@ -3,7 +3,7 @@ include ../../../Make.inc
#
# The object files
#
FOBJS = dcooprt.o dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o\
FOBJS = dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o\
zcoomm.o zcoomv.o zcoonrmi.o zcoorws.o zcoosm.o zcoosv.o
base/serial/coo/dcooprt.f
-80
View File
@@ -1,80 +0,0 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
c
c What if a wrong DESCRA is passed?
c
c
SUBROUTINE DCOOPRT(M,N,DESCRA,AR,IA,JA,INFOA,TITLE,IOUT)
C
C
C .. Scalar Arguments ..
INTEGER M, N, IOUT
C .. Array Arguments ..
DOUBLE PRECISION AR(*)
INTEGER IA(*), JA(*),INFOA(*)
CHARACTER DESCRA*11, TITLE*(*)
C .. Local Scalars ..
INTEGER J
C .. External Subroutines ..
C
C
if ((descra(1:1).eq.'g').or.(descra(1:1).eq.'G')) then
write(iout,fmt=998)
else if ((descra(1:1).eq.'s').or.(descra(1:1).eq.'S')) then
write(iout,fmt=997)
else
write(iout,fmt=998)
endif
nnzero = infoa(1)
write(iout,fmt=992)
write(iout,fmt=996)
write(iout,fmt=996) title
write(iout,fmt=995) 'Number of rows: ',m
write(iout,fmt=995) 'Number of columns: ',n
write(iout,fmt=995) 'Nonzero entries: ',nnzero
write(iout,fmt=996)
write(iout,fmt=992)
write(iout,*) m,n,nnzero
998 format('%%MatrixMarket matrix coordinate real general')
997 format('%%MatrixMarket matrix coordinate real symmetric')
992 format('%======================================== ')
996 format('% ',a)
995 format('% ',a,i9,a,i9,a,i9)
do j=1,nnzero
write(iout,fmt=994) ia(j),ja(j),ar(j)
994 format(i6,1x,i6,1x,e16.8)
enddo
RETURN
END
base/serial/coo/dcoosm.f
+7 -5
View File
@@ -30,9 +30,8 @@ C
C
SUBROUTINE DCOOSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,IA,JA,INFOA,
* B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
use psb_error_mod
IMPLICIT NONE
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
DOUBLE PRECISION ALPHA, BETA
INTEGER LDB, LDC, LWORK, M, N, IERROR
CHARACTER UNITD, TRANST
@@ -43,18 +42,21 @@ C
CHARACTER DIAG, UPLO
INTRINSIC DBLE, IDINT
CHARACTER*20 NAME
integer :: debug_level, debug_unit
NAME = 'DCOOSM\0'
NAME = 'ZCOOSM\0'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IF((ALPHA.NE.1.D0) .OR. (BETA.NE.0.D0))then
IERROR=5
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ENDIF
if (debug) write(*,*) 'DCOOSM ',m
if (debug) write(*,*) 'DCOOSM ',m,ierror
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),':' ,m,ierror
UPLO = '?'
IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') UPLO = 'U'
base/serial/coo/dcoosv.f
+146 -138
View File
@@ -34,10 +34,9 @@ C "right" place, i.e. the last in its row for Lower and the first
C for Upper.
C
SUBROUTINE DCOOSV (UPLO,TRANS,DIAG,N,AS,IA,JA,INFOA,B,X,IERROR)
use psb_error_mod
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D0)
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
INTEGER N,IERROR
CHARACTER DIAG, TRANS, UPLO
DOUBLE PRECISION AS(*), B(*), X(*)
@@ -45,120 +44,129 @@ C
DOUBLE PRECISION ACC
INTEGER I, J, K, NNZ, II
LOGICAL LOW, TRA, UNI
if (debug) write(*,*) 'DCOOSV ',n
if (debug) write(*,*) 'DCOOSV ',n,nnz,diag,trans,uplo
integer :: debug_level, debug_unit
character(len=20) :: name='dcoosv'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
UNI = (DIAG.EQ.'U')
TRA = (TRANS.EQ.'T')
LOW = (UPLO.EQ.'L')
NNZ = INFOA(1)
if (debug) write(*,*) 'DCOOSV ',n,nnz,uni,tra,low,ia(1),ja(1)
IERROR = 0
if (debug) write(*,*) 'DCOOSV ierror ',ierror
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),':',n,nnz,diag,trans,uplo
IF ( .NOT. TRA) THEN
if (debug) write(*,*) 'DCOOSV NOT TRA'
IF (LOW) THEN
if (debug) write(*,*) 'DCOOSV LOW'
IF ( .NOT. UNI) THEN
if (debug) write(*,*) 'DCOOSV NOT UNI'
I = 1
J = I
DO WHILE (I.LE.NNZ)
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J-2
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J-1)
I = J
ENDDO
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT TRA'
IF (LOW) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': LOW'
IF ( .NOT. UNI) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT UNI'
I = 1
J = I
DO WHILE (I.LE.NNZ)
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J-2
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J-1)
I = J
ENDDO
ELSE IF (UNI) THEN
ELSE IF (UNI) THEN
C
C Bug warning: if UNI, some rows might be empty
C
I = 1
if (debug) write(*,*) 'DCOOSV UNILOW ',
+ i,n,nnz,uni,tra,low
DO II = 1, N
if (debug) write(*,*) 'Loop1 COOSV',i,j,ii,x(ii),b(ii)
DO WHILE ((I.LE.NNZ).AND.(IA(I).LT.II))
I = I + 1
c$$$ if (debug) write(*,*) 'Loop2 COOSV',i,ia(i),ii
ENDDO
ACC = ZERO
IF ((I.LE.NNZ).AND.(IA(I).EQ.II)) THEN
J = I + 1
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
DO K = I, J-1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
if (debug) write(*,*) 'Loop COOSV',i,j,ii,x(ii),b(ii)
I = J
I = 1
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': UNILOW ',
+ i,n,nnz,uni,tra,low
DO II = 1, N
DO WHILE ((I.LE.NNZ).AND.(IA(I).LT.II))
I = I + 1
ENDDO
END IF
ELSE IF ( .NOT. LOW) THEN
if (debug) write(*,*) 'DCOOSV NOT LOW'
IF ( .NOT. UNI) THEN
if (debug) write(*,*) 'DCOOSV NOT UNI'
I = NNZ
J = NNZ
DO WHILE (I.GT.0)
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ACC = ZERO
IF ((I.LE.NNZ).AND.(IA(I).EQ.II)) THEN
J = I + 1
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J+2,-1
DO K = I, J-1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J+1)
I = J
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
I = J
ENDDO
END IF
ELSE IF ( .NOT. LOW) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT LOW'
IF ( .NOT. UNI) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT UNI'
I = NNZ
J = NNZ
DO WHILE (I.GT.0)
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J+2,-1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J+1)
I = J
ENDDO
ELSE IF (UNI) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': UNI'
I = NNZ
DO II = N,1,-1
DO WHILE ((I.GT.0).AND.(IA(I).GT.II))
I = I -1
ENDDO
ACC = ZERO
IF ((I.GT.0).AND.(IA(I).EQ.II)) THEN
J = I - 1
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ENDDO
DO K = I, J+1, -1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
I = J
ENDDO
ELSE IF (UNI) THEN
if (debug) write(*,*) 'DCOOSV UNI'
I = NNZ
DO II = N,1,-1
DO WHILE ((I.GT.0).AND.(IA(I).GT.II))
I = I -1
ENDDO
ACC = ZERO
IF ((I.GT.0).AND.(IA(I).EQ.II)) THEN
J = I - 1
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ENDDO
DO K = I, J+1, -1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
if (debug) write(*,*) 'Loop COOSV',i,j,ii,x(ii),b(ii)
I = J
ENDDO
END IF
END IF
END IF
END IF
ELSE IF (TRA) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': TRA'
IERROR = 3010
return
CCCCCCCCCCCCCCCC
@@ -166,48 +174,48 @@ C
C TBF
C
CCCCCCCCCCCCCCCC
DO 180 I = 1, N
X(I) = B(I)
180 CONTINUE
IF (LOW) THEN
IF ( .NOT. UNI) THEN
DO 220 I = N, 1, -1
X(I) = X(I)/AS(IA(I+1)-1)
ACC = X(I)
DO 200 J = IA(I), IA(I+1) - 2
K = JA(J)
X(K) = X(K) - AS(J)*ACC
200 CONTINUE
220 CONTINUE
ELSE IF (UNI) THEN
DO 260 I = N, 1, -1
ACC = X(I)
DO 240 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
240 CONTINUE
260 CONTINUE
END IF
ELSE IF ( .NOT. LOW) THEN
IF ( .NOT. UNI) THEN
DO 300 I = 1, N
X(I) = X(I)/AS(IA(I))
ACC = X(I)
DO 280 J = IA(I) + 1, IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
280 CONTINUE
300 CONTINUE
ELSE IF (UNI) THEN
DO 340 I = 1, N
ACC = X(I)
DO 320 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
320 CONTINUE
340 CONTINUE
END IF
END IF
DO 180 I = 1, N
X(I) = B(I)
180 CONTINUE
IF (LOW) THEN
IF ( .NOT. UNI) THEN
DO 220 I = N, 1, -1
X(I) = X(I)/AS(IA(I+1)-1)
ACC = X(I)
DO 200 J = IA(I), IA(I+1) - 2
K = JA(J)
X(K) = X(K) - AS(J)*ACC
200 CONTINUE
220 CONTINUE
ELSE IF (UNI) THEN
DO 260 I = N, 1, -1
ACC = X(I)
DO 240 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
240 CONTINUE
260 CONTINUE
END IF
ELSE IF ( .NOT. LOW) THEN
IF ( .NOT. UNI) THEN
DO 300 I = 1, N
X(I) = X(I)/AS(IA(I))
ACC = X(I)
DO 280 J = IA(I) + 1, IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
280 CONTINUE
300 CONTINUE
ELSE IF (UNI) THEN
DO 340 I = 1, N
ACC = X(I)
DO 320 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
320 CONTINUE
340 CONTINUE
END IF
END IF
END IF
RETURN
END
base/serial/coo/zcoomm.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
c
c What if a wrong DESCRA is passed?
c
base/serial/coo/zcoomv.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
***********************************************************************
* ZCOOMV. Prolog to be updated. *
* *
base/serial/coo/zcoonrmi.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C ... Compute Infinity norm for sparse matrix in CSR Format ...
DOUBLE PRECISION FUNCTION ZCOONRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
+ INFOA,IERROR)
base/serial/coo/zcoorws.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
SUBROUTINE ZCOORWS(TRANS,M,N,DESCRA,A,IA1,IA2,
& INFOA,ROWSUM,IERROR)
IMPLICIT NONE
base/serial/coo/zcoosm.f
+72 -14
View File
@@ -1,36 +1,77 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
SUBROUTINE ZCOOSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,IA,JA,INFOA,
* B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
COMPLEX*16 ALPHA, BETA
use psb_error_mod
IMPLICIT NONE
COMPLEX*16 ALPHA, BETA
INTEGER LDB, LDC, LWORK, M, N, IERROR
CHARACTER UNITD, TRANST
COMPLEX*16 A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
INTEGER IA(*), JA(*), INFOA(*)
INTEGER IA(*), JA(*), INFOA(*), INT_VAL(5)
CHARACTER DESCRA*11
INTEGER I, K
INTEGER I, K, ERR_ACT
CHARACTER DIAG, UPLO
EXTERNAL XERBLA
INTRINSIC DBLE, IDINT
CHARACTER*20 NAME
integer :: debug_level, debug_unit
NAME = 'DCOOSM\0'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IF((ALPHA.NE.1.D0) .OR. (BETA.NE.0.D0))then
call xerbla('DCSSM ',9)
RETURN
IERROR=5
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ENDIF
if (debug) write(*,*) 'ZCOOSM ',m
if (debug) write(*,*) 'ZCOOSM ',m,ierror
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),':' ,m,ierror
UPLO = '?'
IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') UPLO = 'U'
IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') UPLO = 'L'
IF (UPLO.EQ.'?') THEN
CALL XERBLA('ZCSSM ',10)
RETURN
IERROR=5
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
END IF
IF (DESCRA(3:3).EQ.'N') DIAG = 'N'
IF (DESCRA(3:3).EQ.'U') DIAG = 'U'
IF(UNITD.EQ.'B') THEN
CALL XERBLA('ZCSSM ',11)
RETURN
IERROR=5
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ENDIF
IF (UNITD.EQ.'R') THEN
DO 40 I = 1, N
@@ -44,6 +85,12 @@
CALL ZCOOSV(UPLO,TRANST,DIAG,M,A,IA,JA,INFOA,
+ B(1,I),C(1,I),IERROR)
60 CONTINUE
IF(IERROR.NE.0) THEN
INT_VAL(1)=IERROR
CALL FCPSB_ERRPUSH(4012,NAME,INT_VAL)
GOTO 9999
END IF
IF (UNITD.EQ.'L') THEN
DO 45 I = 1, N
DO 25 K = 1, M
@@ -51,5 +98,16 @@
25 CONTINUE
45 CONTINUE
END IF
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
RETURN
9999 CONTINUE
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
IF ( ERR_ACT .NE. 0 ) THEN
CALL FCPSB_SERROR()
RETURN
ENDIF
RETURN
END
base/serial/coo/zcoosv.f
+175 -138
View File
@@ -1,13 +1,41 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
C Assumption: the triangular matrix has the diagonal element in the
C "right" place, i.e. the last in its row for Lower and the first
C for Upper.
C
SUBROUTINE ZCOOSV (UPLO,TRANS,DIAG,N,AS,IA,JA,INFOA,B,X,IERROR)
use psb_error_mod
COMPLEX*16 ZERO
PARAMETER (ZERO=(0.0D0,0.0D0))
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
INTEGER N,IERROR
CHARACTER DIAG, TRANS, UPLO
COMPLEX*16 AS(*), B(*), X(*)
@@ -15,120 +43,129 @@ C
COMPLEX*16 ACC
INTEGER I, J, K, NNZ, II, JJ
LOGICAL LOW, TRA, UNI
if (debug) write(*,*) 'ZCOOSV ',n
if (debug) write(*,*) 'ZCOOSV ',n,nnz,diag,trans,uplo
integer :: debug_level, debug_unit
character(len=20) :: name='zcoosv'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
UNI = (DIAG.EQ.'U')
TRA = (TRANS.EQ.'T')
LOW = (UPLO.EQ.'L')
NNZ = INFOA(1)
if (debug) write(*,*) 'ZCOOSV ',n,nnz,uni,tra,low,ia(1),ja(1)
IERROR = 0
if (debug) write(*,*) 'ZCOOSV ierror ',ierror
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),':',n,nnz,diag,trans,uplo
IF ( .NOT. TRA) THEN
if (debug) write(*,*) 'ZCOOSV NOT TRA'
IF (LOW) THEN
if (debug) write(*,*) 'ZCOOSV LOW'
IF ( .NOT. UNI) THEN
if (debug) write(*,*) 'ZCOOSV NOT UNI'
I = 1
J = I
DO WHILE (I.LE.NNZ)
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J-2
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J-1)
I = J
ENDDO
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT TRA'
IF (LOW) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': LOW'
IF ( .NOT. UNI) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT UNI'
I = 1
J = I
DO WHILE (I.LE.NNZ)
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J-2
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J-1)
I = J
ENDDO
ELSE IF (UNI) THEN
ELSE IF (UNI) THEN
C
C Bug warning: if UNI, some rows might be empty
C
I = 1
if (debug) write(*,*) 'ZCOOSV UNILOW ',
+ i,n,nnz,uni,tra,low
DO II = 1, N
if (debug) write(*,*) 'Loop1 COOSV',i,j,ii,x(ii),b(ii)
DO WHILE ((I.LE.NNZ).AND.(IA(I).LT.II))
I = I + 1
c$$$ if (debug) write(*,*) 'Loop2 COOSV',i,ia(i),ii
ENDDO
ACC = ZERO
IF ((I.LE.NNZ).AND.(IA(I).EQ.II)) THEN
J = I + 1
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
DO K = I, J-1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
if (debug) write(*,*) 'Loop COOSV',i,j,ii,x(ii),b(ii)
I = J
I = 1
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': UNILOW ',
+ i,n,nnz,uni,tra,low
DO II = 1, N
DO WHILE ((I.LE.NNZ).AND.(IA(I).LT.II))
I = I + 1
ENDDO
END IF
ELSE IF ( .NOT. LOW) THEN
if (debug) write(*,*) 'ZCOOSV NOT LOW'
IF ( .NOT. UNI) THEN
if (debug) write(*,*) 'ZCOOSV NOT UNI'
I = NNZ
J = NNZ
DO WHILE (I.GT.0)
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ACC = ZERO
IF ((I.LE.NNZ).AND.(IA(I).EQ.II)) THEN
J = I + 1
DO WHILE ((J.LE.NNZ).AND.(IA(J).EQ.IA(I)))
J = J+1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J+2,-1
DO K = I, J-1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J+1)
I = J
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
I = J
ENDDO
END IF
ELSE IF ( .NOT. LOW) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT LOW'
IF ( .NOT. UNI) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NOT UNI'
I = NNZ
J = NNZ
DO WHILE (I.GT.0)
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ENDDO
ACC = ZERO
IR = IA(I)
DO K = I, J+2,-1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
X(IR) = (B(IR)-ACC)/AS(J+1)
I = J
ENDDO
ELSE IF (UNI) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': UNI'
I = NNZ
DO II = N,1,-1
DO WHILE ((I.GT.0).AND.(IA(I).GT.II))
I = I -1
ENDDO
ACC = ZERO
IF ((I.GT.0).AND.(IA(I).EQ.II)) THEN
J = I - 1
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ENDDO
DO K = I, J+1, -1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
I = J
ENDDO
ELSE IF (UNI) THEN
if (debug) write(*,*) 'ZCOOSV UNI'
I = NNZ
DO II = N,1,-1
DO WHILE ((I.GT.0).AND.(IA(I).GT.II))
I = I -1
ENDDO
ACC = ZERO
IF ((I.GT.0).AND.(IA(I).EQ.II)) THEN
J = I - 1
DO WHILE ((J.GT.0).AND.(IA(J).EQ.IA(I)))
J = J-1
ENDDO
DO K = I, J+1, -1
JC = JA(K)
ACC = ACC + AS(K)*X(JC)
ENDDO
ELSE
J = I
ENDIF
X(II) = (B(II)-ACC)
if (debug) write(*,*) 'Loop COOSV',i,j,ii,x(ii),b(ii)
I = J
ENDDO
END IF
END IF
END IF
END IF
ELSE IF (TRA) THEN
if (debug_level>=psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': TRA'
IERROR = 3010
return
CCCCCCCCCCCCCCCC
@@ -136,48 +173,48 @@ C
C TBF
C
CCCCCCCCCCCCCCCC
DO 180 I = 1, N
X(I) = B(I)
180 CONTINUE
IF (LOW) THEN
IF ( .NOT. UNI) THEN
DO 220 I = N, 1, -1
X(I) = X(I)/AS(IA(I+1)-1)
ACC = X(I)
DO 200 J = IA(I), IA(I+1) - 2
K = JA(J)
X(K) = X(K) - AS(J)*ACC
200 CONTINUE
220 CONTINUE
ELSE IF (UNI) THEN
DO 260 I = N, 1, -1
ACC = X(I)
DO 240 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
240 CONTINUE
260 CONTINUE
END IF
ELSE IF ( .NOT. LOW) THEN
IF ( .NOT. UNI) THEN
DO 300 I = 1, N
X(I) = X(I)/AS(IA(I))
ACC = X(I)
DO 280 J = IA(I) + 1, IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
280 CONTINUE
300 CONTINUE
ELSE IF (UNI) THEN
DO 340 I = 1, N
ACC = X(I)
DO 320 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
320 CONTINUE
340 CONTINUE
END IF
END IF
DO 180 I = 1, N
X(I) = B(I)
180 CONTINUE
IF (LOW) THEN
IF ( .NOT. UNI) THEN
DO 220 I = N, 1, -1
X(I) = X(I)/AS(IA(I+1)-1)
ACC = X(I)
DO 200 J = IA(I), IA(I+1) - 2
K = JA(J)
X(K) = X(K) - AS(J)*ACC
200 CONTINUE
220 CONTINUE
ELSE IF (UNI) THEN
DO 260 I = N, 1, -1
ACC = X(I)
DO 240 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
240 CONTINUE
260 CONTINUE
END IF
ELSE IF ( .NOT. LOW) THEN
IF ( .NOT. UNI) THEN
DO 300 I = 1, N
X(I) = X(I)/AS(IA(I))
ACC = X(I)
DO 280 J = IA(I) + 1, IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
280 CONTINUE
300 CONTINUE
ELSE IF (UNI) THEN
DO 340 I = 1, N
ACC = X(I)
DO 320 J = IA(I), IA(I+1) - 1
K = JA(J)
X(K) = X(K) - AS(J)*ACC
320 CONTINUE
340 CONTINUE
END IF
END IF
END IF
RETURN
END
base/serial/csr/Makefile
+1 -1
View File
@@ -5,7 +5,7 @@ include ../../../Make.inc
#
FOBJS = dcsrck.o dcsrmm.o dcsrsm.o dcsrmv.o dcsrsv.o dcrnrmi.o \
dcsrprt.o dcsrmv4.o dcsrmv2.o dcsrmv3.o dcsrrws.o\
dcsrmv4.o dcsrmv2.o dcsrmv3.o dcsrrws.o\
zcrnrmi.o zcsrmm.o zcsrrws.o zcsrsm.o zsrmv.o zsrsv.o zcsrck.o
OBJS=$(FOBJS)
base/serial/csr/dcsrprt.f
-84
View File
@@ -1,84 +0,0 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
c
c What if a wrong DESCRA is passed?
c
c
*
*
SUBROUTINE DCSRPRT(M,N,DESCRA,AR,JA,IA,TITLE,IOUT)
C
C
C .. Scalar Arguments ..
INTEGER M, N, IOUT
C .. Array Arguments ..
DOUBLE PRECISION AR(*)
INTEGER IA(*), JA(*)
CHARACTER DESCRA*11, TITLE*(*)
C .. Local Scalars ..
INTEGER I, J, nnzero
C .. External Subroutines ..
C
C
if ((descra(1:1).eq.'g').or.(descra(1:1).eq.'G')) then
write(iout,fmt=998)
else if ((descra(1:1).eq.'s').or.(descra(1:1).eq.'S')) then
write(iout,fmt=997)
else
write(iout,fmt=998)
endif
nnzero = ia(m+1) -1
write(iout,fmt=992)
write(iout,fmt=996)
write(iout,fmt=996) title
write(iout,fmt=995) 'Number of rows: ',m
write(iout,fmt=995) 'Number of columns: ',n
write(iout,fmt=995) 'Nonzero entries: ',nnzero
write(iout,fmt=996)
write(iout,fmt=992)
write(iout,*) m,n,nnzero
998 format('%%MatrixMarket matrix coordinate real general')
997 format('%%MatrixMarket matrix coordinate real symmetric')
992 format('%======================================== ')
996 format('% ',a)
995 format('% ',a,i9,a,i9,a,i9)
do i=1, m
do j=ia(i),ia(i+1)-1
write(iout,fmt=994) i,ja(j),ar(j)
994 format(i6,1x,i6,1x,e16.8)
enddo
enddo
RETURN
END
base/serial/csr/dcsrsm.f
-9
View File
@@ -38,8 +38,6 @@ C
CHARACTER DESCRA*11
INTEGER I, K
CHARACTER DIAG, UPLO
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
C .. Local Arrays ..
CHARACTER*20 NAME
INTEGER INT_VAL(5), err_Act
@@ -97,13 +95,6 @@ C .. Local Arrays ..
GOTO 9999
END IF
if (debug) then
write(0,*) 'Check from DCSRSM'
do k=1,m
write(0,*) k, b(k,1),c(k,1)
enddo
endif
IF (UNITD.EQ.'L') THEN
DO 45 I = 1, N
DO 25 K = 1, M
base/serial/csr/zcrnrmi.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
DOUBLE PRECISION FUNCTION ZCRNRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
& INFOA,IERROR)
IMPLICIT NONE
base/serial/csr/zcsrmm.f
+29
View File
@@ -1,4 +1,33 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
C SUBROUTINE ZCSRMM(TRANSA,M,K,N,ALPHA,DESCRA,AR,
C * JA,IA,B,LDB,BETA,C,LDC,WORK,LWORK)
C
base/serial/csr/zcsrrws.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
SUBROUTINE ZCSRRWS(TRANS,M,N,DESCRA,A,IA1,IA2,
& INFOA,ROWSUM,IERROR)
IMPLICIT NONE
base/serial/csr/zsrmv.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
***********************************************************************
* ZSRMV modified for SPARKER *
* *
base/serial/csr/zsrsv.f
+29
View File
@@ -1,3 +1,32 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
SUBROUTINE ZSRSV (UPLO,TRANS,DIAG,N,AS,JA,IA,B,X)
COMPLEX*16 ZERO
PARAMETER (ZERO = (0.0D0, 0.0D0))
base/serial/dp/dcoco.f
+19 -12
View File
@@ -33,7 +33,7 @@ c
subroutine dcoco(trans,m,n,unitd,d,descra,ar,ia1,ia2,info,
* p1,descrn,arn,ia1n,ia2n,infon,p2,larn,lia1n,
* lia2n,aux,laux,ierror)
use psb_error_mod
use psb_const_mod
use psb_spmat_type
use psb_string_mod
@@ -55,8 +55,7 @@ c .. local scalars ..
integer elem_in, elem_out
logical scale
integer max_nnzero
logical debug
parameter (debug=.false.)
integer :: debug_level, debug_unit
c .. local arrays ..
character*20 name
integer int_val(5)
@@ -69,9 +68,11 @@ c .. external subroutines ..
c .. executable statements ..
c
name = 'dcoco\0'
name = 'dcoco'
ierror = 0
call fcpsb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror)
@@ -81,8 +82,8 @@ c
p2(1) = 0
call psb_getifield(nnz,psb_nnz_,info,psb_ifasize_,ierror)
if (debug) then
write(*,*) 'on entry to dcoco: nnz laux ',
if (debug_level >= psb_debug_serial_) then
write(debug_unit,*) trim(name),': on entry nnz laux ',
+ nnz,laux,larn,lia1n,lia2n
endif
if (laux.lt.nnz+2) then
@@ -119,14 +120,16 @@ c
c
c sort COO data structure
c
if (debug) write(*,*)'first sort',nnz
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': first sort',nnz
do k=1, nnz
arn(k) = ar(k)
ia1n(k) = ia1(k)
ia2n(k) = ia2(k)
enddo
if (debug) write(*,*)'second sort'
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': second sort'
if ((lia2n.ge.(2*nnz+psb_ireg_flgs_+1))
+ .and.(laux.ge.2*(2+nnz))) then
@@ -145,7 +148,9 @@ c
ia2n(ip1+psb_nnzt_) = nnz
ia2n(ip1+psb_nnz_) = 0
ia2n(ip1+psb_ichk_) = nnz+check_flag
if (debug) write(0,*) 'build check :',ia2n(ip1+psb_nnzt_)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': build check :',ia2n(ip1+psb_nnzt_)
c .... order with key ia1n ...
call msort_up(nnz,ia1n,aux,iret)
@@ -183,7 +188,7 @@ c ... error, there are duplicated elements ...
c ... insert only the first duplicated element ...
ia2n(ip2+aux(ipx+elem_in-1)-1) = elem_out
else if (check_flag.eq.psb_dupl_add_) then
c ... add the duplicated element ...
c ... sum the duplicated element ...
arn(elem_out) = arn(elem_out) + arn(elem_in)
ia2n(ip2+aux(ipx+elem_in-1)-1) = elem_out
end if
@@ -230,7 +235,7 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the first duplicated element ...
else if (check_flag.eq.psb_dupl_add_) then
c ... add the duplicated element ...
c ... sum the duplicated element ...
arn(elem_out) = arn(elem_out) + arn(elem_in)
end if
else
@@ -244,7 +249,9 @@ c ... add the duplicated element ...
infon(psb_nnz_) = elem_out
infon(psb_srtd_) = psb_isrtdcoo_
if (debug) write(*,*)'done rebuild COO',infon(1)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': done rebuild COO',infon(1)
else if (toupper(descra(1:1)).eq.'S' .and.
+ toupper(descra(2:2)).eq.'U') then
base/serial/dp/dcocr.f
+41 -51
View File
@@ -34,8 +34,9 @@ C
SUBROUTINE DCOCR(TRANS,M,N,UNITD,D,DESCRA,AR,JA,IA,INFO,
* P1,DESCRN,ARN,IAN1,IAN2,INFON,P2,LARN,LIAN1,
* LIAN2,AUX,LAUX,IERROR)
use psb_const_mod
use psb_error_mod
use psb_spmat_type
use psb_string_mod
IMPLICIT NONE
@@ -57,12 +58,11 @@ C .. Local Scalars ..
integer elem, elem_csr,regen_flag
logical scale
integer max_nnzero
logical debug
parameter (debug=.false.)
integer, allocatable :: itmp(:)
c .. local arrays ..
character*20 name
integer int_val(5)
integer :: debug_level, debug_unit
c
C ...Common variables...
@@ -72,9 +72,11 @@ C .. External Subroutines ..
C .. Executable Statements ..
C
NAME = 'DCOCR\0'
NAME = 'DCOCR'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror)
call psb_getifield(regen_flag,psb_upd_,infon,psb_ifasize_,ierror)
@@ -84,10 +86,9 @@ C
SCALE = (toupper(UNITD).EQ.'L') ! meaningless
P1(1) = 0
P2(1) = 0
nnz = info(1)
if (debug) then
write(0,*) 'On entry to DCOCR: NNZ LAUX ',
if (debug_level >= psb_debug_serial_) then
write(debug_unit,*) trim(name),': On entry NNZ LAUX ',
+ nnz,laux,larn,lian1,lian2
endif
IF (LAUX.LT.NNZ+2) THEN
@@ -139,10 +140,9 @@ C
C
C Sort COO data structure
C
if (debug) write(0,*)'First sort',nnz
c$$$ write(0,*) 'DCOCR Sizes ',lian2,((m+1)+nnz+psb_ireg_flgs_+1),
c$$$ + m+1,nnz,psb_ireg_flgs_,
c$$$ + laux,2*(2+nnz)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': First sort',nnz
if ((regen_flag==psb_upd_perm_).and.
+ (lian2.ge.((m+1)+nnz+psb_ireg_flgs_+1))
+ .and.(laux.ge.2*(2+nnz))) then
@@ -162,17 +162,22 @@ c
ian2(ip1+psb_nnz_) = 0
ian2(ip1+psb_ichk_) = nnz+check_flag
if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Build check :',ian2(ip1+psb_nnzt_)
C .... Order with key IA ...
call msort_up(nnz,itmp,aux,iret)
if (iret.eq.0) call reordvn3(nnz,arn,itmp,ian1,aux(ipx),aux)
if (debug) then
if (debug_level >= psb_debug_serial_) then
do i=1, nnz-1
if (itmp(i).gt.itmp(i+1)) then
write(0,*) 'Sorting error:',i,itmp(i),itmp(i+1)
write(debug_unit,*) trim(name),
+ 'Sorting error:',i,itmp(i),itmp(i+1)
endif
enddo
write(0,*) 'nnz :',m,nnz,itmp(nnz),ian1(nnz)
write(debug_unit,*) trim(name),
+ 'nnz :',m,nnz,itmp(nnz),ian1(nnz)
endif
C .... Order with key JA ...
@@ -200,7 +205,10 @@ c ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr
@@ -233,14 +241,10 @@ c ... error, there are duplicated elements ...
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
if (debug) write(0,*) 'duplicated overwrite perm ',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
if (debug) write(0,*) 'duplicated add perm ',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -280,7 +284,10 @@ C ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
elem = elem+1
@@ -309,13 +316,9 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
if (debug) write(0,*) 'Duplicated overwrite srch',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
if (debug) write(0,*) 'Duplicated add srch',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -324,13 +327,9 @@ c ... sum the duplicated element ...
enddo
endif
if (debug) write(0,*)'Done Rebuild CSR',
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': Done Rebuild CSR',
+ ian2(m+1),ia(elem)
c$$$ if (debug) then
c$$$ do i=ian2(m+1), nnz
c$$$ write(0,*) 'Overflow check :',ia(i),ja(i),ar(i)
c$$$ enddo
c$$$ endif
ELSE IF (toupper(DESCRA(1:1)).EQ.'S' .AND.
+ toupper(DESCRA(2:2)).EQ.'U') THEN
@@ -342,7 +341,6 @@ c$$$ endif
else if (toupper(DESCRA(1:1)).EQ.'T' .AND.
+ toupper(DESCRA(2:2)).EQ.'U') THEN
call msort_up(nnz,itmp,aux,iret)
if (iret.eq.0) call reordvn(nnz,arn,itmp,ian1,aux)
C .... Order with key JA ...
@@ -371,7 +369,10 @@ C ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
elem = elem+1
@@ -400,13 +401,9 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
if (debug) write(0,*) 'Duplicated overwrite srch',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
if (debug) write(0,*) 'Duplicated add srch',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -436,8 +433,6 @@ C .... Order with key JA ...
i = j
enddo
C ... Construct CSR Representation...
elem = 1
elem_csr = 1
@@ -445,7 +440,10 @@ C ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
elem = elem+1
@@ -474,13 +472,9 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
if (debug) write(0,*) 'Duplicated overwrite srch',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
if (debug) write(0,*) 'Duplicated add srch',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -488,13 +482,9 @@ c ... sum the duplicated element ...
ian2(row+1) = elem_csr
enddo
if (debug) write(0,*)'Done Rebuild CSR',
+ ian2(m+1),ia(elem)
if (debug) then
do i=ian2(m+1), nnz
write(0,*) 'Overflow check :',ia(i),ja(i),ar(i)
enddo
endif
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': Done Rebuild CSR',
+ ian2(m+1),ia(elem)
end if
base/serial/dp/dcrco.f
+13 -5
View File
@@ -35,6 +35,7 @@ C
use psb_const_mod
use psb_spmat_type
use psb_string_mod
use psb_error_mod
IMPLICIT NONE
C
@@ -54,23 +55,27 @@ C .. Local Scalars ..
c .. Local Arrays ..
CHARACTER*20 NAME
INTEGER INT_VAL(5)
logical, parameter :: debug=.false.
integer :: debug_level, debug_unit
C .. External Subroutines ..
EXTERNAL MAX_NNZERO
C .. Executable Statements ..
C
NAME = 'DCRCO\0'
NAME = 'DCRCO'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IF (toupper(TRANS).EQ.'N') THEN
SCALE = (toupper(UNITD).EQ.'L') ! meaningless
IP1(1) = 0
IP2(1) = 0
NNZ = IA2(M+1)-1
if (debug) write(0,*) 'CRCO: ',m,n,nnz,' : ',
+ descra,' : ',descrn
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': entry',m,n,nnz,
+ ' : ',descra,' : ',descrn
IF (LARN.LT.NNZ) THEN
IERROR = 60
INT_VAL(1) = 18
@@ -106,7 +111,10 @@ C ... Construct COO Representation...
ENDDO
ENDDO
INFON(psb_nnz_) = elem
if (debug) write(0,*) 'CRCO endloop',m,elem
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': endloop',m,elem
ELSE IF (toupper(DESCRA(1:1)).EQ.'S' .AND.
+ toupper(DESCRA(2:2)).EQ.'U') THEN
base/serial/dp/dcrcr.f
+7 -5
View File
@@ -188,18 +188,20 @@ C .. Intrinsic Functions ..
C .. Executable Statements ..
C
EXIT=.FALSE.
NAME = 'DCOCO\0'
NAME = 'DCRCR'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
C
C Check for argument errors
C
idescra=toupper(descra)
IF(((idescra(1:1) .EQ. 'S' .OR. idescra(1:1) .EQ. 'H' .OR.
& idescra(1:1) .EQ. 'A') .AND. (toupper(unitd) .NE. 'B')) .OR.
IF (((idescra(1:1) .EQ. 'S' .OR. idescra(1:1) .EQ. 'H' .OR.
& idescra(1:1) .EQ. 'A') .AND. (toupper(unitd) .NE. 'B'))
+ .OR.
& (.NOT.((idescra(3:3).EQ.'N').OR.(idescra(3:3).EQ.'L').OR.
+ (idescra(3:3).EQ.'U'))) .OR.
+ toupper(TRANS).NE.'N') THEN
+ (idescra(3:3).EQ.'U')))
+ .OR.
+ toupper(TRANS).NE.'N') THEN
IERROR = 20
ENDIF
IF(LAN.LT.(IA2(M+1)-1)) THEN
base/serial/dp/dcrjd.f
-2
View File
@@ -78,8 +78,6 @@ C .. Local Scalars ..
INTEGER IOFF, ISTROW, NJA, NZ, PIA,
+ PJA, PNG, K, MAX_NG, NG, LJA, ERR_ACT
LOGICAL SCALE
logical debug
parameter (debug=.false.)
CHARACTER UPLO
INTEGER MAX_NNZERO
c .. Local Arrays ..
base/serial/dp/djadrp1.f
-3
View File
@@ -114,8 +114,6 @@ C .. Local Scalars ..
INTEGER I, K, IPG, ERR_ACT
C .. Intrinsic Functions ..
INTRINSIC DBLE
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
C .. Local Arrays ..
CHARACTER*20 NAME
INTEGER INT_VAL(5)
@@ -127,7 +125,6 @@ C
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
IF(toupper(TRANS).EQ.'N') THEN
IF (DEBUG) WRITE(0,*)'DJADRP1:',NG
DO IPG = 1, NG
DO K = IA(2,IPG), IA(3,IPG)-1
DO I = JA(K), JA(K+1) - 1
base/serial/dp/djdco.f
+7 -4
View File
@@ -32,6 +32,7 @@ C
* IP1,DESCRN,ARN,IA1N,IA2N,INFON,IP2,LARN,LIA1N,
* LIA2N,AUX,LAUX,IERROR)
use psb_const_mod
use psb_error_mod
IMPLICIT NONE
C
C .. Scalar Arguments ..
@@ -45,21 +46,23 @@ C .. Array Arguments ..
CHARACTER DESCRA*11, DESCRN*11
C .. Local Scalars ..
INTEGER PIA, PJA, PNG, ERR_ACT
logical debug
parameter (debug=.false.)
integer :: debug_level, debug_unit
c .. Local Arrays ..
CHARACTER*20 NAME
INTEGER INT_VAL(5)
NAME = 'DJDCO\0'
NAME = 'DJDCO'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
PNG = IA2(1)
PIA = IA2(2)
PJA = IA2(3)
if(debug) write(*,*) 'On entry to DJDCO: NNZ LAUX ',
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': On entry NNZ LAUX ',
+ info(1),laux,larn,lia1n,lia2n
CALL DJDCOX(TRANS,M,N,DESCRA,AR,IA2(PIA),IA2(PJA),
base/serial/dp/djdcox.f
+7 -5
View File
@@ -38,6 +38,7 @@ C
use psb_const_mod
use psb_string_mod
use psb_error_mod
IMPLICIT NONE
C
@@ -54,8 +55,7 @@ C .. Local Scalars ..
INTEGER IPX, IPG, NNZ, K, ROW,
* I, J, NZL, IRET, ERR_ACT
LOGICAL SCALE
logical debug
parameter (debug=.false.)
integer :: debug_level, debug_unit
c .. Local Arrays ..
CHARACTER*20 NAME
INTEGER INT_VAL(5)
@@ -63,9 +63,11 @@ c .. Local Arrays ..
C
C .. Executable Statements ..
C
NAME = 'DJDCOX\0'
NAME = 'DJDCOX'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IF (toupper(TRANS).EQ.'N') THEN
C SCALE = (UNITD.EQ.'L') ! meaningless
@@ -80,8 +82,8 @@ C SCALE = (UNITD.EQ.'L') ! meaningless
NNZ = JA(IA(2,NG+1)-1 +1)-1
if (debug) then
write(0,*) 'On entry to DJDCOX: NNZ LAUX ',
if (debug_level >= psb_debug_serial_) then
write(debug_unit,*) trim(name),': On entry NNZ LAUX ',
+ nnz,laux,larn,lia1n,lia2n
endif
IF (LAUX.LT.NNZ+2) THEN
base/serial/dp/dvtfg.f
+11 -5
View File
@@ -34,6 +34,7 @@ C
SUBROUTINE DVTFG (UPLO,M,JA,IA,NG,IPA,IPAT,KLEN,IWORK1,IWORK2,
* IWORK3)
use psb_string_mod
use psb_error_mod
implicit none
C .. Scalar Arguments ..
INTEGER M, NG
@@ -45,10 +46,12 @@ C .. Local Scalars ..
INTEGER I, J, L, L0, L1, LEV, NP, iret
C .. Intrinsic Functions ..
INTRINSIC MAX
logical debug
parameter (debug=.false.)
integer :: debug_level, debug_unit
character(len=20) :: name='DVTFG'
C .. Executable Statements ..
C
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
NG = 0
C
C CHECK ON THE NUMBERS OF THE ELEMENTS OF THE MATRIX
@@ -152,7 +155,9 @@ C
DO 260 L = 1, L1
IPA(L) = IWORK3(L)
260 CONTINUE
if (debug) write(0,*) 'DVTFG: Group ',1,':',(ipa(l),l=1,l1)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Group ',1,':',(ipa(l),l=1,l1)
DO 360 LEV = 2, NG
C
C LOOP ON GROUPS
@@ -177,8 +182,9 @@ C
IPA(L0+L) = IWORK3(L0+L)
320 CONTINUE
ENDIF
if (debug) write(0,*) 'DVTFG: Group ',lev,
+ ':',(ipa(l0+l),l=1,l1)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Group ',lev,':',(ipa(l0+l),l=1,l1)
360 CONTINUE
C
C IPAT = IPA-1
base/serial/dp/zcoco.f
+17 -9
View File
@@ -33,7 +33,7 @@ c
subroutine zcoco(trans,m,n,unitd,d,descra,ar,ia1,ia2,info,
* p1,descrn,arn,ia1n,ia2n,infon,p2,larn,lia1n,
* lia2n,aux,laux,ierror)
use psb_error_mod
use psb_const_mod
use psb_spmat_type
use psb_string_mod
@@ -55,8 +55,7 @@ c .. local scalars ..
integer elem_in, elem_out
logical scale
integer max_nnzero
logical debug
parameter (debug=.false.)
integer :: debug_level, debug_unit
c .. local arrays ..
character*20 name
integer int_val(5)
@@ -72,16 +71,19 @@ c
name = 'zcoco'
ierror = 0
call fcpsb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror)
if (toupper(trans).eq.'N') then
scale = (toupper(unitd).eq.'L') ! meaningless
p1(1) = 0
p2(1) = 0
call psb_getifield(nnz,psb_nnz_,info,psb_ifasize_,ierror)
if (debug) then
write(*,*) 'on entry to dcoco: nnz laux ',
if (debug_level >= psb_debug_serial_) then
write(debug_unit,*) trim(name),': on entry nnz laux ',
+ nnz,laux,larn,lia1n,lia2n
endif
if (laux.lt.nnz+2) then
@@ -118,14 +120,16 @@ c
c
c sort COO data structure
c
if (debug) write(*,*)'first sort',nnz
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': first sort',nnz
do k=1, nnz
arn(k) = ar(k)
ia1n(k) = ia1(k)
ia2n(k) = ia2(k)
enddo
if (debug) write(*,*)'second sort'
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': second sort'
if ((lia2n.ge.(2*nnz+psb_ireg_flgs_+1))
+ .and.(laux.ge.2*(2+nnz))) then
@@ -144,7 +148,9 @@ c
ia2n(ip1+psb_nnzt_) = nnz
ia2n(ip1+psb_nnz_) = 0
ia2n(ip1+psb_ichk_) = nnz+check_flag
if (debug) write(0,*) 'build check :',ia2n(ip1+psb_nnzt_)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': build check :',ia2n(ip1+psb_nnzt_)
c .... order with key ia1n ...
call msort_up(nnz,ia1n,aux,iret)
@@ -244,7 +250,9 @@ c ... sum the duplicated element ...
infon(psb_nnz_) = elem_out
infon(psb_srtd_) = psb_isrtdcoo_
if (debug) write(*,*)'done rebuild COO',infon(1)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': done rebuild COO',infon(1)
else if (toupper(descra(1:1)).eq.'S' .and.
+ toupper(descra(2:2)).eq.'U') then
base/serial/dp/zcocr.f
+43 -57
View File
@@ -36,6 +36,7 @@ C
* LIAN2,AUX,LAUX,IERROR)
use psb_const_mod
use psb_error_mod
use psb_spmat_type
use psb_string_mod
IMPLICIT NONE
@@ -56,13 +57,12 @@ C .. Local Scalars ..
integer ipx, ip1, ip2, check_flag, err_act
integer elem, elem_csr,regen_flag
logical scale
integer max_nnzero
logical debug
parameter (debug=.false.)
integer max_nnzero
integer, allocatable :: itmp(:)
c .. local arrays ..
character*20 name
integer int_val(5)
integer :: debug_level, debug_unit
C
C ...Common variables...
@@ -71,10 +71,11 @@ C .. External Subroutines ..
EXTERNAL MAX_NNZERO
C .. Executable Statements ..
C
NAME = 'ZCOCR\0'
NAME = 'ZCOCR'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror)
call psb_getifield(regen_flag,psb_upd_,infon,psb_ifasize_,ierror)
@@ -84,10 +85,9 @@ C
SCALE = (toupper(UNITD).EQ.'L') ! meaningless
P1(1) = 0
P2(1) = 0
nnz = info(1)
if (debug) then
write(0,*) 'On entry to ZCOCR: NNZ LAUX ',
if (debug_level >= psb_debug_serial_) then
write(debug_unit,*) trim(name),': On entry NNZ LAUX ',
+ nnz,laux,larn,lian1,lian2
endif
IF (LAUX.LT.NNZ+2) THEN
@@ -131,15 +131,17 @@ C
ian1(k) = ja(k)
itmp(k) = ia(k)
enddo
! Mark as unavailable by default.
infon(psb_upd_pnt_) = 0
IF (toupper(descra(1:1)).EQ.'G') THEN
C
C Sort COO data structure
C
if (debug) write(0,*)'First sort',nnz
c$$$ write(0,*) 'ZCOCR Sizes ',lian2,((m+1)+nnz+psb_ireg_flgs_+1),
c$$$ + m+1,nnz,psb_ireg_flgs_,
c$$$ + laux,2*(2+nnz)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': First sort',nnz
if ((regen_flag==psb_upd_perm_).and.
+ (lian2.ge.((m+1)+nnz+psb_ireg_flgs_+1))
+ .and.(laux.ge.2*(2+nnz))) then
@@ -159,24 +161,23 @@ c
ian2(ip1+psb_nnz_) = 0
ian2(ip1+psb_ichk_) = nnz+check_flag
c$$$ write(0,*)'ZCOCR Check: ',ip2,ian2(ip1+psb_iflag_),
c$$$ + ian2(ip1+psb_nnzt_), ian2(ip1+psb_nnz_),
c$$$ + ian2(ip1+psb_ichk_)
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Build check :',ian2(ip1+psb_nnzt_)
c$$$ + ip1,ip2,nnz,ian2(ip1+nnzt_)
if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_)
C .... Order with key IA ...
call msort_up(nnz,itmp,aux,iret)
if (iret.eq.0)
+ call zreordvn3(nnz,arn,itmp,ian1,aux(ipx),aux)
if (debug) then
if (debug_level >= psb_debug_serial_) then
do i=1, nnz-1
if (itmp(i).gt.itmp(i+1)) then
write(0,*) 'Sorting error:',i,itmp(i),itmp(i+1)
write(debug_unit,*) trim(name),
+ 'Sorting error:',i,itmp(i),itmp(i+1)
endif
enddo
write(0,*) 'nnz :',m,nnz,itmp(nnz),ian1(nnz)
write(debug_unit,*) trim(name),
+ 'nnz :',m,nnz,itmp(nnz),ian1(nnz)
endif
C .... Order with key JA ...
@@ -204,7 +205,10 @@ c ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr
@@ -237,14 +241,10 @@ c ... error, there are duplicated elements ...
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
if (debug) write(0,*) 'duplicated overwrite perm ',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
if (debug) write(0,*) 'duplicated add perm ',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -284,7 +284,10 @@ C ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
elem = elem+1
@@ -313,13 +316,9 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
if (debug) write(0,*) 'Duplicated overwrite srch',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
if (debug) write(0,*) 'Duplicated add srch',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -328,13 +327,9 @@ c ... sum the duplicated element ...
enddo
endif
if (debug) write(0,*)'Done Rebuild CSR',
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': Done Rebuild CSR',
+ ian2(m+1),ia(elem)
c$$$ if (debug) then
c$$$ do i=ian2(m+1), nnz
c$$$ write(0,*) 'Overflow check :',ia(i),ja(i),ar(i)
c$$$ enddo
c$$$ endif
ELSE IF (toupper(DESCRA(1:1)).EQ.'S' .AND.
+ toupper(DESCRA(2:2)).EQ.'U') THEN
@@ -346,7 +341,6 @@ c$$$ endif
else if (toupper(DESCRA(1:1)).EQ.'T' .AND.
+ toupper(DESCRA(2:2)).EQ.'U') THEN
call msort_up(nnz,itmp,aux,iret)
if (iret.eq.0) call zreordvn(nnz,arn,itmp,ian1,aux)
C .... Order with key JA ...
@@ -375,7 +369,10 @@ C ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
elem = elem+1
@@ -404,13 +401,9 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
if (debug) write(0,*) 'Duplicated overwrite srch',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
if (debug) write(0,*) 'Duplicated add srch',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -440,8 +433,6 @@ C .... Order with key JA ...
i = j
enddo
C ... Construct CSR Representation...
elem = 1
elem_csr = 1
@@ -449,7 +440,10 @@ C ... Insert first element ...
do row = 1, itmp(1)
ian2(row) = 1
enddo
if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Rebuild CSR',ia(1),elem_csr
ian1(elem_csr) = ian1(elem)
arn(elem_csr) = arn(elem)
elem = elem+1
@@ -478,13 +472,9 @@ c ... error, there are duplicated elements ...
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the last duplicated element ...
arn(elem_csr-1) = arn(elem)
if (debug) write(0,*) 'Duplicated overwrite srch',
+ elem_csr-1,elem
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_csr-1) = arn(elem_csr-1) + arn(elem)
if (debug) write(0,*) 'Duplicated add srch',
+ elem_csr-1,elem
end if
endif
elem = elem + 1
@@ -492,13 +482,9 @@ c ... sum the duplicated element ...
ian2(row+1) = elem_csr
enddo
if (debug) write(0,*)'Done Rebuild CSR',
+ ian2(m+1),ia(elem)
if (debug) then
do i=ian2(m+1), nnz
write(0,*) 'Overflow check :',ia(i),ja(i),ar(i)
enddo
endif
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': Done Rebuild CSR',
+ ian2(m+1),ia(elem)
end if
base/serial/dp/zcrco.f
+15 -5
View File
@@ -35,6 +35,7 @@ C
use psb_const_mod
use psb_spmat_type
use psb_string_mod
use psb_error_mod
IMPLICIT NONE
C
@@ -54,21 +55,27 @@ C .. Local Scalars ..
c .. Local Arrays ..
CHARACTER*20 NAME
INTEGER INT_VAL(5)
integer :: debug_level, debug_unit
C .. External Subroutines ..
EXTERNAL MAX_NNZERO
C .. Executable Statements ..
C
NAME = 'ZCRCO\0'
NAME = 'ZCRCO'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IF (toupper(TRANS).EQ.'N') THEN
SCALE = (toupper(UNITD).EQ.'L') ! meaningless
IP1(1) = 0
IP2(1) = 0
NNZ = IA2(M+1)-1
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': entry',m,n,nnz,
+ ' : ',descra,' : ',descrn
IF (LARN.LT.NNZ) THEN
IERROR = 60
INT_VAL(1) = 18
@@ -93,17 +100,20 @@ C
IF (toupper(DESCRA(1:1)).EQ.'G') THEN
C ... Construct COO Representation...
ELEM = 1
ELEM = 0
DO ROW = 1, M
DO J = IA2(ROW), IA2(ROW+1)-1
ELEM = ELEM + 1
IAN1(ELEM) = ROW
IAN2(ELEM) = IA1(J)
ARN(ELEM) = AR(J)
ELEM = ELEM + 1
ENDDO
ENDDO
INFON(psb_nnz_) = IA2(M+1)-1
INFON(psb_nnz_) = elem
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),': endloop',m,elem
ELSE IF (toupper(DESCRA(1:1)).EQ.'S' .AND.
+ toupper(DESCRA(2:2)).EQ.'U') THEN
base/serial/dp/zcrcr.f
+1 -1
View File
@@ -188,7 +188,7 @@ C .. Intrinsic Functions ..
C .. Executable Statements ..
C
EXIT=.FALSE.
NAME = 'DCOCO\0'
NAME = 'ZCRCR'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
C
base/serial/dp/zcrjd.f
-2
View File
@@ -78,8 +78,6 @@ C .. Local Scalars ..
INTEGER IOFF, ISTROW, NJA, NZ, PIA,
+ PJA, PNG, K, MAX_NG, NG, LJA, ERR_ACT
LOGICAL SCALE
logical debug
parameter (debug=.false.)
CHARACTER UPLO
INTEGER MAX_NNZERO
c .. Local Arrays ..
base/serial/f77/Makefile
+1 -1
View File
@@ -4,7 +4,7 @@ include ../../../Make.inc
# The object files
#
FOBJS = daxpby.o dcsmm.o dcsnmi.o dcsrp.o dcssm.o \
dgelp.o dlpupd.o dswmm.o dswprt.o \
dgelp.o dlpupd.o dswmm.o \
dswsm.o smmp.o dcsrws.o \
zcsnmi.o zaxpby.o zcsmm.o zcssm.o zswmm.o zswsm.o\
zcsrws.o zgelp.o zlpupd.o
base/serial/f77/dcsrp.f
+8 -6
View File
@@ -118,6 +118,7 @@ C
C
SUBROUTINE DCSRP(TRANS,M,N,FIDA,DESCRA,IA1,IA2,INFOA,
+ P,WORK,LWORK,IERROR)
use psb_error_mod
IMPLICIT NONE
C .. Scalar Arguments ..
INTEGER LWORK, M, N, IERROR
@@ -130,9 +131,7 @@ C .. Array Arguments ..
CHARACTER DESCRA*11, FIDA*5
C .. External Subroutines ..
EXTERNAL DCSRRP
logical debug
parameter (debug=.false.)
integer :: debug_level, debug_unit
CHARACTER*20 NAME
C
C .. Executable Statements ..
@@ -140,7 +139,9 @@ C
C
C Check on M, N, TRANS
C
NAME = 'DCSRP\0'
NAME = 'DCSRP'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
@@ -167,14 +168,15 @@ C
C
C Switching on FIDA
C
c$$$ write(0,*) 'DCSRP FORMAT: ',fida
IF (FIDA(1:3).EQ.'CSR') THEN
C
C Permuting CSR structure
C
CALL DCSRRP(TRANS,M,N,DESCRA,IA1,IA2,P,WORK,LWORK)
ELSE IF (FIDA(1:3).EQ.'JAD') THEN
if (debug) write(0,*) 'Calling djadrp',m,p(1),lwork
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),
+ ': Calling djadrp',m,p(1),lwork
CALL DJADRP(TRANS,M,N,DESCRA,IA1,IA2,P,WORK,LWORK)
ELSE
C
base/serial/f77/dcssm.f
+18 -9
View File
@@ -196,6 +196,7 @@ C
+ PL,FIDT,DESCRT,T,IT1,IT2,INFOT,PR,
+ B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
C .. Scalar Arguments ..
use psb_error_mod
IMPLICIT NONE
DOUBLE PRECISION ALPHA, BETA
INTEGER N, LDB, LDC, M, LWORK, IERROR
@@ -211,11 +212,11 @@ C .. Local Scalars ..
C .. Local Array..
INTEGER INT_VAL(5)
CHARACTER*30 STRINGS(2)
CHARACTER NAME*30
CHARACTER NAME*20
integer :: debug_level, debug_unit
C .. Parameters ..
PARAMETER (ZERO=0.D0)
LOGICAL DEBUG
PARAMETER (DEBUG=.FALSE.)
C .. External Subroutines ..
EXTERNAL DSWSM, DLPUPD
C .. Intrinsic Functions ..
@@ -225,10 +226,12 @@ C
C Check for argument errors
C
IERROR = 0
NAME = 'DCSSM\0'
NAME = 'DCSSM'
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
IF (M.LT.0) THEN
IF (M.LT.0) THEN
IERROR = 10
INT_VAL(1) = 2
INT_VAL(2) = M
@@ -299,6 +302,9 @@ C
C
C Both right and left permutations required
C
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': RP LP ',m,n,ierror
CALL DLPUPD(M,N,PR,B,LDB,BETA,WORK,M)
CALL DSWSM(TRANS,M,N,ALPHA,UNITD,D,FIDT,DESCRT,T,IT1,IT2,
& INFOT,WORK,M,ZERO,WORK(P),M,WORK(P+LWORKB),LWORK,IERROR)
@@ -314,7 +320,9 @@ C
C
C Only right permutation required
C
c$$$ write(0,*) 'DCSSM: RP'
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': RP NLP ',m,n,ierror
CALL DLPUPD(M,N,PR,B,LDB,BETA,WORK,M)
CALL DSWSM(TRANS,M,N,ALPHA,UNITD,D,FIDT,DESCRT,T,IT1,IT2,
& INFOT,WORK,M,ZERO,C,LDC,WORK(P),LWORK,IERROR)
@@ -330,7 +338,8 @@ c$$$ write(0,*) 'DCSSM: RP'
C
C Only left permutation required
C
c$$$ write(0,*) 'DCSSM: LP'
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NRP LP ',m,n,ierror
CALL DSWSM(TRANS,M,N,ALPHA,UNITD,D,FIDT,DESCRT,T,IT1,IT2,
& INFOT,B,LDB,BETA,WORK,M,WORK(P),LWORK,IERROR)
LWORKS = IDINT(WORK(P))
@@ -345,8 +354,8 @@ c$$$ write(0,*) 'DCSSM: LP'
C
C Only triangular systems solver required
C
if (debug) write(*,*) 'DCSSM ',m,n
if (debug) write(*,*) 'DCSSM ',m,n,ierror
if (debug_level >= psb_debug_serial_comp_)
+ write(debug_unit,*) trim(name),': NRP NLP ',m,n,ierror
CALL DSWSM(TRANS,M,N,ALPHA,UNITD,D,FIDT,DESCRT,T,IT1,IT2,
& INFOT,B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
LWORKS = IDINT(WORK(1))

Some files were not shown because too many files have changed in this diff Show More