fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merged in the version with ALLOCATABLE components. We are now moving

Browse Source 
towards version 2.1.
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent a675cdd316
commit 9620471951
121 changed files with 1728 additions and 1815 deletions
Show Stats Download Patch File Download Diff File

8
Changelog
Unescape Escape View File

@ -1,6 +1,14 @@
Changelog. A lot less detailed than usual, at least for past
history.
2006/11/09: The allocatable version works, but under gcc42 there is a
strange problem when using -fbounds-check. Be careful!
2006/11/08: Merged the allocatable version; hope everything works!
2006/11/08: Branched version psblas2-2-0-maint, and defined tag
2.0.2.6
2006/11/01: Merged changes in the handling of data exchange.
2006/10/03: Merged in the multilevel preconditioner stuff. This is

4
Make.inc
Unescape Escape View File

@ -10,8 +10,8 @@ F90=/usr/local/gcc42/bin/gfortran
FC=/usr/local/gcc42/bin/gfortran
F77=$(FC)
CC=/usr/local/gcc42/bin/gcc
F90COPT= -O3 -march=pentium4 -msse2 -mfpmath=sse -ggdb -fbounds-check
FCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse -ggdb -fbounds-check
F90COPT= -O3 -march=pentium4 -msse2 -mfpmath=sse -ggdb
FCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse -ggdb
CCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse -ggdb
####################### Section 2 #######################

2
Make.inc.g95
Unescape Escape View File

@ -1,6 +1,6 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
.SUFFIXES: .f90 $(.mod) .F90
####################### Section 1 #######################

6
Make.inc.gcc42
Unescape Escape View File

@ -1,6 +1,6 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
.SUFFIXES: .f90 $(.mod) .F90
####################### Section 1 #######################
@ -11,8 +11,8 @@ FC=/usr/local/gcc42/bin/gfortran
F77=$(FC)
CC=/usr/local/gcc42/bin/gcc
F90COPT= -O3 -march=pentium4 -msse2 -mfpmath=sse
FCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse
CCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse
FCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse
CCOPT=-O3 -march=pentium4 -msse2 -mfpmath=sse
####################### Section 2 #######################
# Define your linker and linker flags here #

86
Make.inc.ifc71
Unescape Escape View File

@ -1,86 +0,0 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
####################### Section 1 #######################
# Define your compilers and compiler flags here #
##########################################################
F90=/opt/intel/compiler70/ia32/bin/ifc
FC=/opt/intel/compiler70/ia32/bin/ifc
F77=$(FC)
F90COPT=-O3
FCOPT=-O3
CC=gcc
CCOPT=-O3
####################### Section 2 #######################
# Define your linker and linker flags here #
##########################################################
F90LINK=/usr/local/mpich-ifc71/bin/mpif90
FLINK=/usr/local/mpich-ifc71/bin/mpif77
MPF90=/usr/local/mpich-ifc71/bin/mpif90
MPF77=/usr/local/mpich-ifc71/bin/mpif77
MPCC=/usr/local/mpich-ifc71/bin/mpicc
####################### Section 3 #######################
# Specify paths to libraries #
##########################################################
BLAS=-lblasintel -L$(HOME)/LIB
BLACS=-lmpiblacsifc71 -L$(HOME)/LIB
EXTRA_BLACS_ENV_OBJS=extra_env.o
####################### Section 4 #######################
# Other useful tools&defines #
##########################################################
SLUDIR=/usr/local/SuperLU_3.0
SLU=-lslu_lx_ifc8 -L$(SLUDIR)
SLUDEF=-DHave_SLU_ -I$(SLUDIR)
UMFDIR=$(HOME)/LIB/Umfpack_gcc41
UMF=-lumfpack -lamd -L$(UMFDIR)
UMFDEF=-DHave_UMF_ -I$(UMFDIR)
CDEFINES=-DAdd_ $(SLUDEF) $(UMFDEF)
AR=ar -cur
RANLIB=ranlib
####################### Section 5 #######################
# Do not edit this #
##########################################################
LIBDIR = lib
LIBNAME = libpsblas.a
TYPEMODS = PSB_SPMAT_TYPE$(.mod) PSB_DESCRIPTOR_TYPE$(.mod) PSB_PREC_TYPE$(.mod) PSB_REALLOC_MOD$(.mod)
CONSTMODS = PSB_CONST_MOD$(.mod)
BLASMODS = $(TYPEMODS) PSB_PSBLAS_MOD$(.mod) PSB_COMM_MOD$(.mod)
METHDMODS = PSB_METHD_MOD$(.mod)
TOOLSMODS = $(CONSTMODS) PSI_MOD$(.mod) PSB_TOOLS_MOD$(.mod) PSB_SERIAL_MOD$(.mod)
PRECMODS = PSB_PREC_MOD$(.mod)
ERRORMODS = PSB_ERROR_MOD$(.mod)
F90MODS= $(BLASMODS) $(PRECMODS) $(METHDMODS) $(TOOLSMODS) $(ERRORMODS) STRING$(.mod)
MODS=$(LIBDIR)/PSB_CONST_MOD$(.mod) $(LIBDIR)/PSB_SPMAT_TYPE$(.mod) $(LIBDIR)/PSB_REALLOC_MOD$(.mod) \
$(LIBDIR)/PSB_DESCRIPTOR_TYPE$(.mod) $(LIBDIR)/PSB_PREC_TYPE$(.mod) $(LIBDIR)/parts.fh \
$(LIBDIR)/PSB_SERIAL_MOD$(.mod) $(LIBDIR)/PSB_COMM_MOD$(.mod) $(LIBDIR)/PSB_ERROR_MOD$(.mod)
# Under Linux/gmake there is a rule interpreting .mod as Modula source!
$(.mod).o:
.f.o:
$(FC) $(FCOPT) $(INCDIRS) -c $<
.c.o:
$(CC) $(CCOPT) $(INCDIRS) $(CDEFINES) -c $<
.f$(.mod):
$(F90) $(FCOPT) $(INCDIRS) -c $<
.f90$(.mod):
$(F90) $(F90COPT) $(INCDIRS) -c $<
.f90.o:
$(F90) $(F90COPT) $(INCDIRS) -c $<

89
Make.inc.ifc8
Unescape Escape View File

@ -1,89 +0,0 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
####################### Section 1 #######################
# Define your compilers and compiler flags here #
##########################################################
IFC8=/opt/intel_fc_80
F90=$(IFC8)/bin/ifort
FC=$(IFC8)/bin/ifort
F77=$(FC)
CC=gcc
F90COPT= -O3
FCOPT=-O3
CCOPT=-O3
####################### Section 2 #######################
# Define your linker and linker flags here #
##########################################################
F90LINK=/usr/local/mpich-ifc80/bin/mpif90 -g -CB -no_cpprt
FLINK=/usr/local/mpich-ifc80/bin/mpif77 -g -CB -no_cpprt
MPF90=/usr/local/mpich-ifc80/bin/mpif90 -g -CB -no_cpprt
MPF77=/usr/local/mpich-ifc80/bin/mpif77 -g -CB -no_cpprt
MPCC=/usr/local/mpich-ifc80/bin/mpicc -g -CB -no_cpprt
####################### Section 3 #######################
# Specify paths to libraries #
##########################################################
BLAS=-lblas-intel -L$(HOME)/NUMERICAL/LIB
BLACS=-lmpiblacs-intel -L$(HOME)/NUMERICAL/LIB
EXTRA_BLACS_ENV_OBJS=extra_env.o
####################### Section 4 #######################
# Other useful tools&defines #
##########################################################
SLUDIR=/usr/local/SuperLU_3.0
SLU=-lslu_lx_ifc8 -L$(SLUDIR)
SLUDEF=-DHave_SLU_ -I$(SLUDIR)
UMFDIR=$(HOME)/LIB/Umfpack_gcc41
UMF=-lumfpack -lamd -L$(UMFDIR)
UMFDEF=-DHave_UMF_ -I$(UMFDIR)
# Add -DLargeFptr for 64-bit addresses
CDEFINES=-DAdd_ $(SLUDEF) $(UMFDEF)
AR=ar -cur
RANLIB=ranlib
####################### Section 5 #######################
# Do not edit this #
##########################################################
LIBDIR = lib
LIBNAME = libpsblas.a
TYPEMODS = psb_spmat_type$(.mod) psb_descriptor_type$(.mod) psb_prec_type$(.mod) psb_realloc_mod$(.mod)
CONSTMODS = psb_const_mod$(.mod)
BLASMODS = $(TYPEMODS) psb_psblas_mod$(.mod) psb_comm_mod$(.mod)
METHDMODS = psb_methd_mod$(.mod)
TOOLSMODS = $(CONSTMODS) psi_mod$(.mod) psb_tools_mod$(.mod) psb_serial_mod$(.mod)
PRECMODS = psb_prec_mod$(.mod)
ERRORMODS = psb_error_mod$(.mod)
F90MODS= $(BLASMODS) $(PRECMODS) $(METHDMODS) $(TOOLSMODS) $(ERRORMODS) string$(.mod)
MODS=$(LIBDIR)/psb_const_mod$(.mod) $(LIBDIR)/psb_spmat_type$(.mod) $(LIBDIR)/psb_realloc_mod$(.mod) \
$(LIBDIR)/psb_descriptor_type$(.mod) $(LIBDIR)/psb_prec_type$(.mod) $(LIBDIR)/parts.fh \
$(LIBDIR)/psb_serial_mod$(.mod) $(LIBDIR)/psb_comm_mod$(.mod) $(LIBDIR)/psb_error_mod$(.mod)
# Under Linux/gmake there is a rule interpreting .mod as Modula source!
$(.mod).o:
.f.o:
$(FC) $(FCOPT) $(INCDIRS) -c $<
.c.o:
$(CC) $(CCOPT) $(INCDIRS) $(CDEFINES) -c $<
.f$(.mod):
$(F90) $(FCOPT) $(INCDIRS) -c $<
.f90$(.mod):
$(F90) $(F90COPT) $(INCDIRS) -c $<
.f90.o:
$(F90) $(F90COPT) $(INCDIRS) -c $<

2
Make.inc.ifc9
Unescape Escape View File

@ -1,6 +1,6 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
.SUFFIXES: .f90 $(.mod) .F90
####################### Section 1 #######################

86
Make.inc.pgf90
Unescape Escape View File

@ -1,86 +0,0 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
####################### Section 1 #######################
# Define your compilers and compiler flags here #
##########################################################
F90=pgf90
FC=pgf90
F77=$(FC)
F90COPT=-fast -g77libs
FCOPT=-fast -g77libs
CC=gcc
CCOPT=-O2 -g -ggdb -pg
####################### Section 2 #######################
# Define your linker and linker flags here #
##########################################################
F90LINK=/usr/local/mpich-pgi/bin/mpif90
FLINK=/usr/local/mpich-pgi/bin/mpif77
MPF90=/usr/local/mpich-pgi/bin/mpif90
MPF77=/usr/local/mpich-pgi/bin/mpif77
MPCC=/usr/local/mpich-pgi/bin/mpicc
####################### Section 3 #######################
# Specify paths to libraries #
##########################################################
BLAS=-lblas-pgi -L$(HOME)/LIB
BLACS=-lmpiblacs-pgi -L$(HOME)/LIB
EXTRA_BLACS_ENV_OBJS=extra_env.o
####################### Section 4 #######################
# Other useful tools&defines #
##########################################################
#SLUDIR=/usr/local/SuperLU_3.0
#SLU=-lslu_lx_gfort -L$(SLUDIR)
#SLUDEF=-DHave_SLU_ -I$(SLUDIR)
#UMFDIR=$(HOME)/LIB/Umfpack_gcc41
#UMF=-lumfpack -lamd -L$(UMFDIR)
#UMFDEF=-DHave_UMF_ -I$(UMFDIR)
# Add -DLargeFptr for 64-bit addresses
CDEFINES=-DAdd_ $(SLUDEF) $(UMFDEF)
AR=ar -cur
RANLIB=ranlib
####################### Section 5 #######################
# Do not edit this #
##########################################################
LIBDIR = lib
LIBNAME = libpsblas.a
TYPEMODS = psb_spmat_type$(.mod) psb_descriptor_type$(.mod) psb_prec_type$(.mod) psb_realloc_mod$(.mod)
CONSTMODS = psb_const_mod$(.mod)
BLASMODS = $(TYPEMODS) psb_psblas_mod$(.mod) psb_comm_mod$(.mod)
METHDMODS = psb_methd_mod$(.mod)
TOOLSMODS = $(CONSTMODS) psi_mod$(.mod) psb_tools_mod$(.mod) psb_serial_mod$(.mod)
PRECMODS = psb_prec_mod$(.mod)
ERRORMODS = psb_error_mod$(.mod)
F90MODS= $(BLASMODS) $(PRECMODS) $(METHDMODS) $(TOOLSMODS) $(ERRORMODS) string$(.mod)
MODS=$(LIBDIR)/psb_const_mod$(.mod) $(LIBDIR)/psb_spmat_type$(.mod) $(LIBDIR)/psb_realloc_mod$(.mod) \
$(LIBDIR)/psb_descriptor_type$(.mod) $(LIBDIR)/psb_prec_type$(.mod) $(LIBDIR)/parts.fh \
$(LIBDIR)/psb_serial_mod$(.mod) $(LIBDIR)/psb_comm_mod$(.mod) $(LIBDIR)/psb_error_mod$(.mod)
# Under Linux/gmake there is a rule interpreting .mod as Modula source!
$(.mod).o:
.f.o:
$(FC) $(FCOPT) $(INCDIRS) -c $<
.c.o:
$(CC) $(CCOPT) $(INCDIRS) $(CDEFINES) -c $<
.f$(.mod):
$(F90) $(FCOPT) $(INCDIRS) -c $<
.f90$(.mod):
$(F90) $(F90COPT) $(INCDIRS) -c $<
.f90.o:
$(F90) $(F90COPT) $(INCDIRS) -c $<

2
Make.inc.rs6k
Unescape Escape View File

@ -1,6 +1,6 @@
.mod=.mod
.fh=.fh
.SUFFIXES: .f90 $(.mod)
.SUFFIXES: .f90 $(.mod) .F90
####################### Section 1 #######################

21
README
Unescape Escape View File

@ -1,4 +1,4 @@
This directory contains the PSBLAS library, version 2.0.2.
This directory contains the PSBLAS library, version 2.1.0
Version 1.0 of the library was described in:
@ -19,24 +19,17 @@ LINUX:
There finally exist a GNU Fortran 95 implementation: we are using the
development snapshots from GCC 3.5.0, later 4.1 and 4.2 since July
2004, and it appears to work. The 4.2 version of GNU compilers is now
our reference platform.
our reference platform; the official 4.2.0 release is due pretty
soon. It now includes support for ALLOCATABLES.
For the PGI compilers, we used them in conjunction with gcc, NOT
pgcc. Note that with pgi 3.6 we have horrible performance, due to
spurious array copies when calling Fortran 77 codes from Fortran 90;
this is fixed in version 4 and later.
The Lahey version we got access to (6.0 and 6.1) seems to suffer from
the same extra copies problem; this is most apparent in the matrix
spurious extra copies problem; this is most apparent in the matrix
build process.
For the Intel compilers, we used ifc versions 7, 8 and 9; with version 6.0
you need to change the way modules are handled, but we recommend to migrate
to the new version anyway. Moreover, with versions prior to 7.1, there
is a strange error in pargen/ppde90: the compiler did not like the
INTERFACE for the dummy argument subroutine PARTS, it wanted an
EXTERNAL specification. Again, please move to 7.1.
For the Intel compilers, we recommend moving to version 9; previous
versions of the library have been compiled with version 7 and 8 of
ifort.
IBM SP.
The library has been tested on an IBM SP2, SP4 and SP5, with XLC and XLF

7
src/comm/psb_dgather.f90
Unescape Escape View File

@ -95,6 +95,8 @@ subroutine psb_dgatherm(globx, locx, desc_a, info, iroot,&
end if
if (root==-1) then
iiroot=0
else
iiroot = root
endif
if (present(iiglobx)) then
@ -258,7 +260,7 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,&
! locals
integer :: int_err(5), ictxt, np, me, &
& err_act, n, root, iiroot, ilocx, iglobx, jlocx,&
& err_act, n, root, ilocx, iglobx, jlocx,&
& jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
character(len=20) :: name, ch_err
@ -289,9 +291,6 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,&
else
root = -1
end if
if (root==-1) then
iiroot=0
endif
jglobx=1
if (present(iiglobx)) then

30
src/comm/psb_dhalo.f90
Unescape Escape View File

@ -31,7 +31,8 @@
! File: psb_dhalo.f90
!
! Subroutine: psb_dhalom
! This subroutine performs the exchange of the halo elements in a distributed dense matrix between all the processes.
! This subroutine performs the exchange of the halo elements in a
! distributed dense matrix between all the processes.
!
! Parameters:
! x - real,dimension(:,:). The local part of the dense matrix.
@ -69,6 +70,7 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
real(kind(1.d0)),pointer :: iwork(:), xp(:,:)
character :: ltran
character(len=20) :: name, ch_err
logical :: aliw
name='psb_dhalom'
if(psb_get_errstatus().ne.0) return
@ -148,8 +150,10 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
aliw=.false.
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -158,7 +162,9 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
!!$ write(0,*) 'halom ',liwork
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -189,8 +195,7 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
if (aliw) deallocate(iwork)
call psb_erractionrestore(err_act)
return
@ -240,7 +245,8 @@ end subroutine psb_dhalom
!!$
!
! Subroutine: psb_dhalov
! This subroutine performs the exchange of the halo elements in a distributed dense vector between all the processes.
! This subroutine performs the exchange of the halo elements in a
! distributed dense vector between all the processes.
!
! Parameters:
! x - real,dimension(:). The local part of the dense vector.
@ -276,6 +282,7 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode)
real(kind(1.d0)),pointer :: iwork(:)
character :: ltran
character(len=20) :: name, ch_err
logical :: aliw
name='psb_dhalov'
if(psb_get_errstatus().ne.0) return
@ -337,8 +344,10 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode)
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
aliw=.false.
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -347,7 +356,7 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode)
end if
end if
else
call psb_realloc(liwork,iwork,info)
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -366,13 +375,12 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode)
end if
if(info.ne.0) then
ch_err='PSI_dSwap...'
ch_err='PSI_swapdata'
call psb_errpush(4010,name,a_err=ch_err)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
if (aliw) deallocate(iwork)
call psb_erractionrestore(err_act)
return

40
src/comm/psb_dovrl.f90
Unescape Escape View File

@ -31,15 +31,16 @@
! File: psb_dovrl.f90
!
! Subroutine: psb_dovrlm
! This subroutine performs the exchange of the overlap elements in a distributed dense matrix between all the processes.
! This subroutine performs the exchange of the overlap elements in a
! distributed dense matrix between all the processes.
!
! Parameters:
! x - real,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code.
! jx - integer(optional). The starting column of the global matrix.
! ik - integer(optional). The number of columns to gather.
! work - real(optional). A working area.
! x - real,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. A return code.
! jx - integer(optional). The starting column of the global matrix
! ik - integer(optional). The number of columns to gather.
! work - real(optional). A work area.
! update - integer(optional). ???.
!
subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,update)
@ -65,6 +66,7 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,update)
real(kind(1.d0)),pointer :: iwork(:), xp(:,:)
logical :: do_update
character(len=20) :: name, ch_err
logical :: aliw
name='psb_dovrlm'
if(psb_get_errstatus().ne.0) return
@ -136,8 +138,10 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,update)
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
aliw=.false.
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -146,7 +150,8 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,update)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -154,7 +159,7 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,update)
goto 9999
end if
end if
! exchange overlap elements
if(do_update) then
xp => x(iix:size(x,1),jjx:jjx+k-1)
@ -194,7 +199,7 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,update)
goto 9999
end select
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
@ -246,7 +251,8 @@ end subroutine psb_dovrlm
!!$
!!$
! Subroutine: psb_dovrlv
! This subroutine performs the exchange of the overlap elements in a distributed dense vector between all the processes.
! This subroutine performs the exchange of the overlap elements in a
! distributed dense vector between all the processes.
!
! Parameters:
! x - real,dimension(:). The local part of the dense vector.
@ -278,6 +284,7 @@ subroutine psb_dovrlv(x,desc_a,info,work,update)
real(kind(1.d0)),pointer :: iwork(:)
logical :: do_update
character(len=20) :: name, ch_err
logical :: aliw
name='psb_dovrlv'
if(psb_get_errstatus().ne.0) return
@ -335,8 +342,10 @@ subroutine psb_dovrlv(x,desc_a,info,work,update)
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
aliw=.false.
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -345,7 +354,8 @@ subroutine psb_dovrlv(x,desc_a,info,work,update)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -392,7 +402,7 @@ subroutine psb_dovrlv(x,desc_a,info,work,update)
goto 9999
end select
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)

108
src/comm/psb_ihalo.f90
Unescape Escape View File

@ -32,7 +32,8 @@
! File: psb_ihalo.f90
!
! Subroutine: psb_ihalom
! This subroutine performs the exchange of the halo elements in a distributed dense matrix between all the processes.
! This subroutine performs the exchange of the halo elements in a
! distributed dense matrix between all the processes.
!
! Parameters:
! x - integer,dimension(:,:). The local part of the dense matrix.
@ -70,6 +71,7 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
integer, pointer :: xp(:,:), iwork(:)
character :: ltran
character(len=20) :: name, ch_err
logical :: aliw
name='psb_ihalom'
if(psb_get_errstatus().ne.0) return
@ -88,49 +90,49 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
ix = 1
if (present(jx)) then
ijx = jx
ijx = jx
else
ijx = 1
ijx = 1
endif
m = desc_a%matrix_data(psb_m_)
n = desc_a%matrix_data(psb_n_)
nrow = desc_a%matrix_data(psb_n_row_)
maxk=size(x,2)-ijx+1
if(present(ik)) then
if(ik.gt.maxk) then
k=maxk
else
k=ik
end if
if(ik.gt.maxk) then
k=maxk
else
k=ik
end if
else
k = maxk
k = maxk
end if
if (present(tran)) then
ltran = tran
ltran = tran
else
ltran = 'N'
ltran = 'N'
endif
if (present(mode)) then
imode = mode
imode = mode
else
imode = IOR(psb_swap_send_,psb_swap_recv_)
imode = IOR(psb_swap_send_,psb_swap_recv_)
endif
! check vector correctness
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
end if
if (iix.ne.1) then
info=3040
call psb_errpush(info,name)
info=3040
call psb_errpush(info,name)
end if
err=info
@ -149,43 +151,46 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
liwork=nrow
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
if(size(work).ge.liwork) then
aliw=.false.
iwork => work
else
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
end if
else
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
xp => x(iix:size(x,1),jjx:jjx+k-1)
! exchange halo elements
if(ltran.eq.'N') then
call psi_swapdata(imode,k,0,xp,&
& desc_a,iwork,info)
call psi_swapdata(imode,k,0,xp,&
& desc_a,iwork,info)
else if((ltran.eq.'T').or.(ltran.eq.'H')) then
call psi_swaptran(imode,k,1,xp,&
& desc_a,iwork,info)
call psi_swaptran(imode,k,1,xp,&
& desc_a,iwork,info)
end if
if(info.ne.0) then
call psb_errpush(4010,name,a_err='PSI_iSwap...')
goto 9999
call psb_errpush(4010,name,a_err='PSI_iSwap...')
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
@ -195,8 +200,8 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(ictxt)
return
call psb_error(ictxt)
return
end if
return
end subroutine psb_ihalom
@ -236,7 +241,8 @@ end subroutine psb_ihalom
!!$
! Subroutine: psb_ihalov
! This subroutine performs the exchange of the halo elements in a distributed dense matrix between all the processes.
! This subroutine performs the exchange of the halo elements in a
! distributed dense matrix between all the processes.
!
! Parameters:
! x - integer,dimension(:). The local part of the dense matrix.
@ -272,6 +278,7 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode)
integer,pointer :: iwork(:)
character :: ltran
character(len=20) :: name, ch_err
logical :: aliw
name='psb_ihalov'
if(psb_get_errstatus().ne.0) return
@ -334,9 +341,11 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode)
liwork=nrow
if (present(work)) then
if(size(work).ge.liwork) then
aliw=.false.
iwork => work
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -345,7 +354,8 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -364,11 +374,11 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode)
end if
if(info.ne.0) then
call psb_errpush(4010,name,a_err='PSI_iSwap...')
call psb_errpush(4010,name,a_err='PSI_iswapdata')
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)

2
src/comm/psb_zgather.f90
Unescape Escape View File

@ -96,6 +96,8 @@ subroutine psb_zgatherm(globx, locx, desc_a, info, iroot,&
end if
if (root==-1) then
iiroot=0
else
iiroot = root
endif
if (present(iiglobx)) then

29
src/comm/psb_zhalo.f90
Unescape Escape View File

@ -31,7 +31,8 @@
! File: psb_zhalo.f90
!
! Subroutine: psb_zhalom
! This subroutine performs the exchange of the halo elements in a distributed dense matrix between all the processes.
! This subroutine performs the exchange of the halo elements in a
! distributed dense matrix between all the processes.
!
! Parameters:
! x - real,dimension(:,:). The local part of the dense matrix.
@ -69,6 +70,7 @@ subroutine psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
complex(kind(1.d0)),pointer :: iwork(:), xp(:,:)
character :: ltran
character(len=20) :: name, ch_err
logical :: aliw
name='psb_zhalom'
if(psb_get_errstatus().ne.0) return
@ -147,9 +149,11 @@ subroutine psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
liwork=nrow
if (present(work)) then
if(size(work).ge.liwork) then
aliw=.false.
iwork => work
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -158,7 +162,9 @@ subroutine psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -178,12 +184,12 @@ subroutine psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
end if
if(info.ne.0) then
ch_err='PSI_dSwap...'
ch_err='PSI_zswapdata'
call psb_errpush(4010,name,a_err=ch_err)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
@ -234,7 +240,8 @@ end subroutine psb_zhalom
!!$
!
! Subroutine: psb_zhalov
! This subroutine performs the exchange of the halo elements in a distributed dense vector between all the processes.
! This subroutine performs the exchange of the halo elements in a
! distributed dense vector between all the processes.
!
! Parameters:
! x - real,dimension(:). The local part of the dense vector.
@ -269,6 +276,7 @@ subroutine psb_zhalov(x,desc_a,info,alpha,work,tran,mode)
complex(kind(1.d0)),pointer :: iwork(:)
character :: ltran
character(len=20) :: name, ch_err
logical :: aliw
name='psb_zhalov'
if(psb_get_errstatus().ne.0) return
@ -329,9 +337,11 @@ subroutine psb_zhalov(x,desc_a,info,alpha,work,tran,mode)
liwork=nrow
if (present(work)) then
if(size(work).ge.liwork) then
aliw=.false.
iwork => work
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -340,7 +350,8 @@ subroutine psb_zhalov(x,desc_a,info,alpha,work,tran,mode)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -364,7 +375,7 @@ subroutine psb_zhalov(x,desc_a,info,alpha,work,tran,mode)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)

38
src/comm/psb_zovrl.f90
Unescape Escape View File

@ -31,15 +31,16 @@
! File: psb_zovrl.f90
!
! Subroutine: psb_zovrlm
! This subroutine performs the exchange of the overlap elements in a distributed dense matrix between all the processes.
! This subroutine performs the exchange of the overlap elements in a
! distributed dense matrix between all the processes.
!
! Parameters:
! x - real,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code.
! jx - integer(optional). The starting column of the global matrix.
! ik - integer(optional). The number of columns to gather.
! work - real(optional). A working area.
! x - real,dimension(:,:). The local part of the dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Returns an output code.
! jx - integer(optional). The starting column of the global matrix.
! ik - integer(optional). The number of columns to gather.
! work - real(optional). A working area.
! update - integer(optional). ???.
!
subroutine psb_zovrlm(x,desc_a,info,jx,ik,work,update)
@ -65,6 +66,7 @@ subroutine psb_zovrlm(x,desc_a,info,jx,ik,work,update)
complex(kind(1.d0)),pointer :: iwork(:), xp(:,:)
logical :: do_update
character(len=20) :: name, ch_err
logical :: aliw
name='psb_zovrlm'
if(psb_get_errstatus().ne.0) return
@ -135,9 +137,11 @@ subroutine psb_zovrlm(x,desc_a,info,jx,ik,work,update)
liwork=ncol
if (present(work)) then
if(size(work).ge.liwork) then
aliw=.false.
iwork => work
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -146,7 +150,8 @@ subroutine psb_zovrlm(x,desc_a,info,jx,ik,work,update)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -194,7 +199,7 @@ subroutine psb_zovrlm(x,desc_a,info,jx,ik,work,update)
goto 9999
end select
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
@ -246,7 +251,8 @@ end subroutine psb_zovrlm
!!$
!!$
! Subroutine: psb_zovrlv
! This subroutine performs the exchange of the overlap elements in a distributed dense vector between all the processes.
! This subroutine performs the exchange of the overlap elements in a
! distributed dense vector between all the processes.
!
! Parameters:
! x - real,dimension(:). The local part of the dense vector.
@ -278,6 +284,7 @@ subroutine psb_zovrlv(x,desc_a,info,work,update)
complex(kind(1.d0)),pointer :: iwork(:)
logical :: do_update
character(len=20) :: name, ch_err
logical :: aliw
name='psb_zovrlv'
if(psb_get_errstatus().ne.0) return
@ -334,9 +341,11 @@ subroutine psb_zovrlv(x,desc_a,info,work,update)
liwork=ncol
if (present(work)) then
if(size(work).ge.liwork) then
aliw=.false.
iwork => work
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -345,7 +354,8 @@ subroutine psb_zovrlv(x,desc_a,info,work,update)
end if
end if
else
call psb_realloc(liwork,iwork,info)
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -393,7 +403,7 @@ subroutine psb_zovrlv(x,desc_a,info,work,update)
goto 9999
end select
if(.not.present(work)) deallocate(iwork)
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)

10
src/internals/psi_crea_bnd_elem.f90
Unescape Escape View File

@ -34,11 +34,11 @@ subroutine psi_crea_bnd_elem(bndel,desc_a,info)
use psb_error_mod
implicit none
integer, pointer :: bndel(:)
integer, allocatable :: bndel(:)
type(psb_desc_type) :: desc_a
integer, intent(out) :: info
integer, pointer :: work(:)
integer, allocatable :: work(:)
integer :: i, j, nr, ns, k, irv, err_act
character(len=20) :: name
@ -84,19 +84,19 @@ subroutine psi_crea_bnd_elem(bndel,desc_a,info)
if (.true.) then
if (j>0) then
allocate(bndel(j),stat=info)
call psb_realloc(j,bndel,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
bndel(1:j) = work(1:j)
else
if (associated(bndel)) then
if (allocated(bndel)) then
deallocate(bndel)
end if
end if
else
allocate(bndel(j+1),stat=info)
call psb_realloc(j+1,bndel,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999

6
src/internals/psi_crea_index.f90
Unescape Escape View File

@ -39,7 +39,7 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info,nxch,nsnd,nrcv
integer, intent(in) :: index_in(:)
integer, pointer :: index_out(:)
integer, allocatable :: index_out(:)
logical :: glob_idx
! ....local scalars...
@ -70,8 +70,8 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info
& isglob_in,info)
integer :: desc_data(:),index_in(:),dep_list(:)
integer :: loc_to_glob(:),glob_to_loc(:)
integer,pointer :: desc_index(:)
integer :: length_dl,nsnd,nrcv, info
integer, allocatable :: desc_index(:)
integer :: length_dl,nsnd,nrcv,info
logical :: isglob_in
end subroutine psi_desc_index
end interface

53
src/internals/psi_crea_ovr_elem.f90
Unescape Escape View File

@ -28,16 +28,20 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
use psb_realloc_mod
use psb_error_mod
use psb_penv_mod
implicit none
! ...parameter arrays....
integer :: desc_overlap(:)
integer, pointer :: ovr_elem(:)
integer, allocatable, intent(inout) :: ovr_elem(:)
integer, intent(out) :: info
! ...local scalars...
integer :: i,pnt_new_elem,ret,j, info
integer :: i,pnt_new_elem,ret,j,iret
integer :: dim_ovr_elem
integer :: pairtree(2)
@ -45,15 +49,21 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
integer :: psi_exist_ovr_elem,dim
external :: psi_exist_ovr_elem
integer :: nel, ip, ix, iel, insize
integer :: nel, ip, ix, iel, insize, err_act
integer, allocatable :: telem(:,:)
logical, parameter :: usetree=.false.
character(len=20) :: name
info = 0
name='psi_crea_ovr_elem'
if (associated(ovr_elem)) then
dim_ovr_elem=size(ovr_elem)
if (allocated(ovr_elem)) then
dim_ovr_elem = size(ovr_elem)
else
dim_ovr_elem = 0
dim_ovr_elem = 0
endif
@ -76,7 +86,6 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
call initpairsearchtree(pairtree,info)
do while (desc_overlap(i).ne.-1)
! ...loop over all procs of desc_overlap list....
i=i+1
do j=1,desc_overlap(i)
! ....loop over all overlap indices referred to act proc.....
@ -92,8 +101,10 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
dim_ovr_elem=max(((3*dim_ovr_elem)/2+2),pnt_new_elem+100)
call psb_realloc(dim_ovr_elem,ovr_elem,info)
if (info /= 0) then
write(0,*) 'Error in CREA_OVR_ELEM'
endif
info = 4000
call psb_errpush(info,name)
goto 9999
end if
endif
ovr_elem(pnt_new_elem)=desc_overlap(i+j)
ovr_elem(pnt_new_elem+1)=2
@ -113,17 +124,22 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
dim_ovr_elem=pnt_new_elem
call psb_realloc(dim_ovr_elem,ovr_elem,info)
if (info /= 0) then
write(0,*) 'Error in CREA_OVR_ELEM'
endif
info = 4000
call psb_errpush(info,name)
goto 9999
end if
ovr_elem(pnt_new_elem)=-1
call freepairsearchtree(pairtree)
else
insize = size(desc_overlap)
insize = max(1,(insize+1)/2)
allocate(telem(insize,2),stat=info)
if (info /= 0) then
write(0,*) 'Error in CREA_OVR_ELEM'
info = 4000
call psb_errpush(info,name)
goto 9999
endif
i = 1
nel = 0
@ -169,4 +185,15 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
ovr_elem(iel) = -1
deallocate(telem)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == act_abort) then
call psb_error()
return
end if
return
end subroutine psi_crea_ovr_elem

16
src/internals/psi_desc_index.f90
Unescape Escape View File

@ -42,8 +42,8 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
! ...array parameters.....
integer :: desc_data(:),index_in(:),dep_list(:)
integer :: loc_to_glob(:),glob_to_loc(:)
integer,pointer :: desc_index(:)
integer :: length_dl, nsnd,nrcv,info
integer,allocatable :: desc_index(:)
integer :: length_dl,nsnd,nrcv,info
logical :: isglob_in
! ....local scalars...
integer :: j,me,np,i,proc
@ -55,7 +55,7 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
integer,allocatable :: brvindx(:),rvsz(:),&
& bsdindx(:),sdsz(:), sndbuf(:), rcvbuf(:)
integer :: ihinsz,ntot,k,err_act, l_di, &
integer :: ihinsz,ntot,k,err_act,nidx,&
& idxr, idxs, iszs, iszr, nesd, nerv, icomm
logical,parameter :: debug=.false., usempi=.true.
@ -136,13 +136,13 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
endif
ntot = (3*(count((sdsz>0).or.(rvsz>0)))+ iszs + iszr) + 1
if (associated(desc_index)) then
l_di = size(desc_index)
if (allocated(desc_index)) then
nidx = size(desc_index)
else
l_di = 0
nidx = 0
endif
if (l_di < ntot) then
if (nidx < ntot) then
call psb_realloc(ntot,desc_index,info)
endif
if (info /= 0) then

2
src/internals/psi_gthsct.f90
Unescape Escape View File

@ -43,6 +43,7 @@ module psi_gthsct_mod
end interface
contains
subroutine psi_dgthm(n,k,idx,x,y)
implicit none
@ -117,7 +118,6 @@ contains
end if
end subroutine psi_dsctm
subroutine psi_dsctv(n,idx,x,beta,y)
implicit none

2
src/internals/psi_sort_dl.f90
Unescape Escape View File

@ -40,7 +40,7 @@ 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, pointer :: work(:)
integer, allocatable :: work(:)
logical, parameter :: debug=.false.
character(len=20) :: name

2
src/methd/psb_dbicg.f90
Unescape Escape View File

@ -98,7 +98,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,&
integer, optional, intent(out) :: iter
real(kind(1.d0)), optional, intent(out) :: err
!!$ local data
real(kind(1.d0)), pointer :: aux(:),wwrk(:,:)
real(kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
real(kind(1.d0)), pointer :: ww(:), q(:),&
& r(:), p(:), zt(:), pt(:), z(:), rt(:),qt(:)
integer :: int_err(5)

4
src/methd/psb_dcg.f90
Unescape Escape View File

@ -98,8 +98,8 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ Local data
real(kind(1.d0)), pointer :: aux(:), q(:), p(:),&
& r(:), z(:), w(:), wwrk(:,:)
real(kind(1.d0)), allocatable, target :: aux(:), wwrk(:,:)
real(kind(1.d0)), pointer :: q(:), p(:), r(:), z(:), w(:)
real(kind(1.d0)) ::rerr
real(kind(1.d0)) ::alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,&
& sigma

2
src/methd/psb_dcgs.f90
Unescape Escape View File

@ -96,7 +96,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ local data
Real(Kind(1.d0)), Pointer :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), Pointer :: ww(:), q(:),&
& r(:), p(:), v(:), s(:), z(:), f(:), rt(:),qt(:),uv(:)
Real(Kind(1.d0)) :: rerr

9
src/methd/psb_dcgstab.f90
Unescape Escape View File

@ -96,7 +96,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ Local data
Real(Kind(1.d0)), Pointer :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), Pointer :: q(:),&
& r(:), p(:), v(:), s(:), t(:), z(:), f(:)
Real(Kind(1.d0)) :: rerr
@ -278,6 +278,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,&
rho_old = rho
rho = psb_gedot(q,r,desc_a,info)
!!$ write(0,'(i2," rho old ",2(f,2x))')me,rho,rho_old
If (debug) Write(0,*) 'Bi-CGSTAB RHO:',rho
If (rho==dzero) Then
If (debug) Write(0,*) 'Bi-CGSTAB Itxation breakdown R',rho
Exit iteration
@ -301,7 +302,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,&
If (debug) Write(0,*) 'Bi-CGSTAB Iteration breakdown S1', sigma
Exit iteration
Endif
If (debug) Write(0,*) 'Bi-CGSTAB SIGMA:',sigma
alpha = rho/sigma
Call psb_geaxpby(done,r,dzero,s,desc_a,info)
if(info.ne.0) then
@ -336,7 +337,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,&
tau = psb_gedot(t,s,desc_a,info)
omega = tau/sigma
If (debug) Write(0,*) 'Bi-CGSTAB OMEGA:',omega
If (omega==dzero) Then
If (debug) Write(0,*) 'BI-CGSTAB ITERATION BREAKDOWN O',omega
Exit iteration
@ -364,7 +365,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,&
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0)) &
& write(*,'(a,i4,3(2x,es10.4)))') &
& write(*,'(a,i4,3(2x,es10.4))') &
& 'bicgstab: ',itx,rerr
Endif

5
src/methd/psb_dcgstabl.f90
Unescape Escape View File

@ -103,10 +103,9 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ local data
Real(Kind(1.d0)), Pointer :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:),uh(:,:), rh(:,:)
Real(Kind(1.d0)), Pointer :: ww(:), q(:), r(:), rt0(:), p(:), v(:), &
& s(:), t(:), z(:), f(:), uh(:,:), rh(:,:), &
& gamma(:), gamma1(:), gamma2(:), taum(:,:), sigma(:)
& s(:), t(:), z(:), f(:), gamma(:), gamma1(:), gamma2(:), taum(:,:), sigma(:)
Real(Kind(1.d0)) :: rerr
Integer :: litmax, naux, m, mglob, it, itrace_,&

5
src/methd/psb_dgmresr.f90
Unescape Escape View File

@ -105,9 +105,8 @@ Subroutine psb_dgmresr(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ local data
Real(Kind(1.d0)), Pointer :: aux(:)
Real(Kind(1.d0)), Pointer :: w(:), v(:,:), &
& c(:),s(:), h(:,:), rs(:), rr(:,:)
Real(Kind(1.d0)), allocatable, target :: aux(:),w(:), v(:,:)
Real(Kind(1.d0)), allocatable :: c(:),s(:), h(:,:), rs(:), rr(:,:)
Real(Kind(1.d0)) :: rerr, scal, gm
Integer ::litmax, liter, naux, m, mglob, it,k, itrace_,&
& np,me, n_row, n_col, nl, int_err(5)

2
src/methd/psb_zcgs.f90
Unescape Escape View File

@ -96,7 +96,7 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ local data
Complex(Kind(1.d0)), Pointer :: aux(:),wwrk(:,:)
Complex(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Complex(Kind(1.d0)), Pointer :: ww(:), q(:),&
& r(:), p(:), v(:), s(:), z(:), f(:), rt(:),qt(:),uv(:)
Real(Kind(1.d0)) :: rerr

2
src/methd/psb_zcgstab.f90
Unescape Escape View File

@ -96,7 +96,7 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,&
Integer, Optional, Intent(out) :: iter
Real(Kind(1.d0)), Optional, Intent(out) :: err
!!$ Local data
Complex(Kind(1.d0)), Pointer :: aux(:),wwrk(:,:)
Complex(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Complex(Kind(1.d0)), Pointer :: q(:),&
& r(:), p(:), v(:), s(:), t(:), z(:), f(:)
Real(Kind(1.d0)) :: rerr

49
src/modules/psb_desc_type.f90
Unescape Escape View File

@ -42,37 +42,36 @@ module psb_descriptor_type
! desc_type contains data for communications.
type psb_desc_type
! contain decomposition informations
integer, pointer :: matrix_data(:)=>null()
! contain index of halo elements to send/receive
integer, pointer :: halo_index(:)=>null()
! contain indices of boundary elements
integer, pointer :: bnd_elem(:)=>null()
! contain index of overlap elements to send/receive
integer, pointer :: ovrlap_index(:)=>null()
! contain for each local overlap element, the number of times
! that is duplicated
integer, pointer :: ovrlap_elem(:)=>null()
! contain for each local element the corresponding global index
integer, pointer :: loc_to_glob(:)=>null()
! contain for each global element the corresponding local index,
! if exist.
integer, pointer :: glob_to_loc (:)=>null()
! local renumbering induced by sparse matrix storage.
integer, pointer :: lprm(:)=>null()
! index space in case it is not just the contiguous range 1:n
integer, pointer :: idx_space(:)=>null()
! contain decomposition informations
integer, allocatable :: matrix_data(:)
! contain index of halo elements to send/receive
integer, allocatable :: halo_index(:)
! contain indices of boundary elements
integer, allocatable :: bnd_elem(:)
! contain index of overlap elements to send/receive
integer, allocatable :: ovrlap_index(:)
! contain for each local overlap element, the number of times
! that is duplicated
integer, allocatable :: ovrlap_elem(:)
! contain for each local element the corresponding global index
integer, allocatable :: loc_to_glob(:)
! contain for each global element the corresponding local index,
! if exist.
integer, allocatable :: glob_to_loc (:)
! local renumbering induced by sparse matrix storage.
integer, allocatable :: lprm(:)
! index space in case it is not just the contiguous range 1:n
integer, allocatable :: idx_space(:)
end type psb_desc_type
contains
subroutine psb_nullify_desc(desc)
type(psb_desc_type), intent(inout) :: desc
nullify(desc%matrix_data,desc%loc_to_glob,desc%glob_to_loc,&
& desc%halo_index,desc%bnd_elem,desc%ovrlap_elem,&
& desc%ovrlap_index, desc%lprm, desc%idx_space)!,&
! & desc%halo_pt,desc%ovrlap_pt)
!!$ nullify(desc%matrix_data,desc%loc_to_glob,desc%glob_to_loc,&
!!$ &desc%halo_index,desc%bnd_elem,desc%ovrlap_elem,&
!!$ &desc%ovrlap_index, desc%lprm, desc%idx_space)
end subroutine psb_nullify_desc

3
src/modules/psb_prec_mod.f90
Unescape Escape View File

@ -70,7 +70,7 @@ module psb_prec_mod
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
end subroutine psb_dprecset
subroutine psb_zprecset(prec,ptype,info,iv,rs,rv)
subroutine psb_zprecset(prec,ptype,info,iv,rs,rv,ilev,nlev)
use psb_serial_mod
use psb_descriptor_type
use psb_prec_type
@ -81,6 +81,7 @@ module psb_prec_mod
integer, optional, intent(in) :: iv(:)
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
integer, optional, intent(in) :: nlev,ilev
end subroutine psb_zprecset
end interface

166
src/modules/psb_prec_type.f90
Unescape Escape View File

@ -84,16 +84,16 @@ module psb_prec_type
type psb_dbaseprc_type
type(psb_dspmat_type), pointer :: av(:) => null() !
real(kind(1.d0)), pointer :: d(:) => null()
type(psb_desc_type), pointer :: desc_data => null(), desc_ac=>null()! !
integer, pointer :: iprcparm(:) => null() !
real(kind(1.d0)), pointer :: dprcparm(:) => null() !
integer, pointer :: perm(:) => null(), invperm(:) => null()
integer, pointer :: mlia(:) => null(), nlaggr(:) => null() !
type(psb_dspmat_type), pointer :: base_a => null() !
type(psb_desc_type), pointer :: base_desc => null() !
real(kind(1.d0)), pointer :: dorig(:) => null() !
type(psb_dspmat_type), allocatable :: av(:)
real(kind(1.d0)), allocatable :: d(:)
type(psb_desc_type) :: desc_data , desc_ac
integer, allocatable :: iprcparm(:)
real(kind(1.d0)), allocatable :: dprcparm(:)
integer, allocatable :: perm(:), invperm(:)
integer, allocatable :: mlia(:), nlaggr(:)
type(psb_dspmat_type), pointer :: base_a => null() !
type(psb_desc_type), pointer :: base_desc=> null() !
real(kind(1.d0)), allocatable :: dorig(:)
end type psb_dbaseprc_type
@ -141,28 +141,28 @@ module psb_prec_type
! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs.
!
type psb_dprec_type
type(psb_dbaseprc_type), pointer :: baseprecv(:) => null()
type(psb_dbaseprc_type), allocatable :: baseprecv(:)
! contain type of preconditioning to be performed
integer :: prec, base_prec
end type psb_dprec_type
type psb_zbaseprc_type
type(psb_zspmat_type), pointer :: av(:) => null() !
complex(kind(1.d0)), pointer :: d(:) => null()
type(psb_desc_type), pointer :: desc_data => null() , desc_ac=>null()!
integer, pointer :: iprcparm(:) => null() !
real(kind(1.d0)), pointer :: dprcparm(:) => null() !
integer, pointer :: perm(:) => null(), invperm(:) => null()
integer, pointer :: mlia(:) => null(), nlaggr(:) => null() !
type(psb_zspmat_type), pointer :: base_a => null() !
type(psb_desc_type), pointer :: base_desc => null() !
complex(kind(1.d0)), pointer :: dorig(:) => null() !
type(psb_zspmat_type), allocatable :: av(:)
complex(kind(1.d0)), allocatable :: d(:)
type(psb_desc_type) :: desc_data , desc_ac
integer, allocatable :: iprcparm(:)
real(kind(1.d0)), allocatable :: dprcparm(:)
integer, allocatable :: perm(:), invperm(:)
integer, allocatable :: mlia(:), nlaggr(:)
type(psb_zspmat_type), pointer :: base_a => null() !
type(psb_desc_type), pointer :: base_desc => null() !
complex(kind(1.d0)), allocatable :: dorig(:)
end type psb_zbaseprc_type
type psb_zprec_type
type(psb_zbaseprc_type), pointer :: baseprecv(:) => null()
type(psb_zbaseprc_type), allocatable :: baseprecv(:)
! contain type of preconditioning to be performed
integer :: prec, base_prec
end type psb_zprec_type
@ -229,7 +229,7 @@ contains
integer :: ilev
write(iout,*) 'Preconditioner description'
if (associated(p%baseprecv)) then
if (allocated(p%baseprecv)) then
if (size(p%baseprecv)>=1) then
write(iout,*) 'Base preconditioner'
select case(p%baseprecv(1)%iprcparm(p_type_))
@ -253,7 +253,7 @@ contains
end if
if (size(p%baseprecv)>=2) then
do ilev = 2, size(p%baseprecv)
if (.not.associated(p%baseprecv(ilev)%iprcparm)) then
if (.not.allocated(p%baseprecv(ilev)%iprcparm)) then
write(iout,*) 'Inconsistent MLPREC part!'
return
endif
@ -267,7 +267,8 @@ contains
write(iout,*) 'Smoother: ', &
& smooth_kinds(p%baseprecv(ilev)%iprcparm(smth_kind_))
if (p%baseprecv(ilev)%iprcparm(smth_kind_) /= no_smth_) then
write(iout,*) 'Smoothing omega: ', p%baseprecv(ilev)%dprcparm(smooth_omega_)
write(iout,*) 'Smoothing omega: ', &
& p%baseprecv(ilev)%dprcparm(smooth_omega_)
write(iout,*) 'Smoothing position: ',&
& smooth_names(p%baseprecv(ilev)%iprcparm(smth_pos_))
end if
@ -372,7 +373,7 @@ contains
type(psb_zprec_type), intent(in) :: p
write(iout,*) 'Preconditioner description'
if (associated(p%baseprecv)) then
if (allocated(p%baseprecv)) then
if (size(p%baseprecv)>=1) then
write(iout,*) 'Base preconditioner'
select case(p%baseprecv(1)%iprcparm(p_type_))
@ -395,7 +396,7 @@ contains
end select
end if
if (size(p%baseprecv)>=2) then
if (.not.associated(p%baseprecv(2)%iprcparm)) then
if (.not.allocated(p%baseprecv(2)%iprcparm)) then
write(iout,*) 'Inconsistent MLPREC part!'
return
endif
@ -660,11 +661,14 @@ contains
info = 0
if (associated(p%d)) then
! Actually we migh just deallocate the top level array, except
! for the inner UMFPACK or SLU stuff
if (allocated(p%d)) then
deallocate(p%d,stat=info)
end if
if (associated(p%av)) then
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
@ -674,53 +678,40 @@ contains
end if
enddo
deallocate(p%av,stat=info)
p%av => null()
end if
if (associated(p%desc_data)) then
if (associated(p%desc_data%matrix_data)) then
call psb_cdfree(p%desc_data,info)
end if
deallocate(p%desc_data)
endif
if (associated(p%desc_ac)) then
if (associated(p%desc_ac%matrix_data)) then
call psb_cdfree(p%desc_ac,info)
end if
deallocate(p%desc_ac)
endif
if (associated(p%dprcparm)) then
! Do we really need the two below? Probably not.
! call psb_cdfree(p%desc_data,info)
! call psb_cdfree(p%desc_ac,info)
if (allocated(p%dprcparm)) then
deallocate(p%dprcparm,stat=info)
end if
if (associated(p%base_a)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_a)
endif
if (associated(p%base_desc)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_desc)
endif
if (associated(p%dorig)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_a)
! This is a pointer to something else, must not free it here.
nullify(p%base_desc)
if (allocated(p%dorig)) then
deallocate(p%dorig,stat=info)
endif
if (associated(p%mlia)) then
if (allocated(p%mlia)) then
deallocate(p%mlia,stat=info)
endif
if (associated(p%nlaggr)) then
if (allocated(p%nlaggr)) then
deallocate(p%nlaggr,stat=info)
endif
if (associated(p%perm)) then
if (allocated(p%perm)) then
deallocate(p%perm,stat=info)
endif
if (associated(p%invperm)) then
if (allocated(p%invperm)) then
deallocate(p%invperm,stat=info)
endif
if (associated(p%iprcparm)) then
if (allocated(p%iprcparm)) then
if (p%iprcparm(f_type_)==f_slu_) then
call psb_dslu_free(p%iprcparm(slu_ptr_),info)
end if
@ -737,8 +728,10 @@ contains
use psb_descriptor_type
type(psb_dbaseprc_type), intent(inout) :: p
nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,&
& p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac)
nullify(p%base_a)
nullify(p%base_desc)
!!$ nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,&
!!$ & p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac)
end subroutine psb_nullify_dbaseprec
@ -752,11 +745,11 @@ contains
info = 0
if (associated(p%d)) then
if (allocated(p%d)) then
deallocate(p%d,stat=info)
end if
if (associated(p%av)) then
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
@ -766,52 +759,40 @@ contains
end if
enddo
deallocate(p%av,stat=info)
p%av => null()
end if
if (associated(p%desc_data)) then
if (associated(p%desc_data%matrix_data)) then
call psb_cdfree(p%desc_data,info)
end if
deallocate(p%desc_data)
endif
if (associated(p%desc_ac)) then
if (associated(p%desc_ac%matrix_data)) then
call psb_cdfree(p%desc_ac,info)
end if
deallocate(p%desc_ac)
endif
if (associated(p%dprcparm)) then
! call psb_cdfree(p%desc_data,info)
! call psb_cdfree(p%desc_ac,info)
if (allocated(p%dprcparm)) then
deallocate(p%dprcparm,stat=info)
end if
if (associated(p%base_a)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_a)
endif
if (associated(p%base_desc)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_desc)
endif
if (associated(p%dorig)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_a)
! This is a pointer to something else, must not free it here.
nullify(p%base_desc)
if (allocated(p%dorig)) then
deallocate(p%dorig,stat=info)
endif
if (associated(p%mlia)) then
if (allocated(p%mlia)) then
deallocate(p%mlia,stat=info)
endif
if (associated(p%nlaggr)) then
if (allocated(p%nlaggr)) then
deallocate(p%nlaggr,stat=info)
endif
if (associated(p%perm)) then
if (allocated(p%perm)) then
deallocate(p%perm,stat=info)
endif
if (associated(p%invperm)) then
if (allocated(p%invperm)) then
deallocate(p%invperm,stat=info)
endif
if (associated(p%iprcparm)) then
if (allocated(p%iprcparm)) then
if (p%iprcparm(f_type_)==f_slu_) then
call psb_zslu_free(p%iprcparm(slu_ptr_),info)
end if
@ -828,8 +809,9 @@ contains
use psb_descriptor_type
type(psb_zbaseprc_type), intent(inout) :: p
nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,&
& p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data,p%desc_ac)
nullify(p%base_a)
nullify(p%base_desc)
end subroutine psb_nullify_zbaseprec

8
src/modules/psb_psblas_mod.f90
Unescape Escape View File

@ -79,6 +79,14 @@ module psb_psblas_mod
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
end subroutine psb_dmdots
subroutine psb_ddot2v(res, x, y,w,z,desc_a, info)
use psb_descriptor_type
real(kind(1.d0)), intent(in) :: x(:), y(:),w(:), z(:)
real(kind(1.d0)), intent(out) :: res(:)
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
end subroutine psb_ddot2v
subroutine psb_zdotvs(res,x, y, desc_a, info)
use psb_descriptor_type
complex(kind(1.d0)), intent(out) :: res

456
src/modules/psb_realloc_mod.f90
Unescape Escape View File

@ -44,6 +44,12 @@ module psb_realloc_mod
module procedure psb_dreallocatez2
end Interface
interface psb_transfer
module procedure psb_dtransfer1d
module procedure psb_itransfer1d
module procedure psb_ztransfer1d
end interface
Interface psb_safe_cpy
module procedure psb_icpy1d,psb_icpy2d, &
& psb_dcpy1d, psb_dcpy2d, psb_zcpy1d, psb_zcpy2d
@ -54,94 +60,15 @@ module psb_realloc_mod
& psb_dsize1d, psb_dsize2d,&
& psb_zsize1d, psb_zsize2d
end interface
contains
function psb_isize1d(vin)
integer :: psb_isize1d
integer, pointer :: vin(:)
if (.not.associated(vin)) then
psb_isize1d = 0
else
psb_isize1d = size(vin)
end if
end function psb_isize1d
function psb_isize2d(vin,dim)
integer :: psb_isize2d
integer, pointer :: vin(:,:)
integer, optional :: dim
if (.not.associated(vin)) then
psb_isize2d = 0
else
if (present(dim)) then
psb_isize2d = size(vin,dim=dim)
else
psb_isize2d = size(vin)
end if
end if
end function psb_isize2d
function psb_dsize1d(vin)
integer :: psb_dsize1d
real(kind(1.d0)), pointer :: vin(:)
if (.not.associated(vin)) then
psb_dsize1d = 0
else
psb_dsize1d = size(vin)
end if
end function psb_dsize1d
function psb_dsize2d(vin,dim)
integer :: psb_dsize2d
real(kind(1.d0)), pointer :: vin(:,:)
integer, optional :: dim
if (.not.associated(vin)) then
psb_dsize2d = 0
else
if (present(dim)) then
psb_dsize2d = size(vin,dim=dim)
else
psb_dsize2d = size(vin)
end if
end if
end function psb_dsize2d
function psb_zsize1d(vin)
integer :: psb_zsize1d
complex(kind(1.d0)), pointer :: vin(:)
if (.not.associated(vin)) then
psb_zsize1d = 0
else
psb_zsize1d = size(vin)
end if
end function psb_zsize1d
function psb_zsize2d(vin,dim)
integer :: psb_zsize2d
complex(kind(1.d0)), pointer :: vin(:,:)
integer, optional :: dim
if (.not.associated(vin)) then
psb_zsize2d = 0
else
if (present(dim)) then
psb_zsize2d = size(vin,dim=dim)
else
psb_zsize2d = size(vin)
end if
end if
end function psb_zsize2d
Contains
subroutine psb_icpy1d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
Integer,pointer :: vin(:),vout(:)
Integer,allocatable :: vin(:),vout(:)
integer :: info
! ...Local Variables
@ -154,7 +81,7 @@ contains
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
if (allocated(vin)) then
isz = size(vin)
call psb_realloc(isz,vout,info)
if (info /= 0) then
@ -186,7 +113,7 @@ contains
use psb_error_mod
! ...Subroutine Arguments
Integer,pointer :: vin(:,:),vout(:,:)
Integer,allocatable :: vin(:,:),vout(:,:)
integer :: info
! ...Local Variables
@ -199,7 +126,7 @@ contains
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
if (allocated(vin)) then
isz1 = size(vin,1)
isz2 = size(vin,2)
call psb_realloc(isz1,isz2,vout,info)
@ -232,7 +159,7 @@ contains
use psb_error_mod
! ...Subroutine Arguments
real(kind(1.d0)),pointer :: vin(:),vout(:)
real(kind(1.d0)),allocatable :: vin(:),vout(:)
integer :: info
! ...Local Variables
@ -245,7 +172,7 @@ contains
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
if (allocated(vin)) then
isz = size(vin)
call psb_realloc(isz,vout,info)
if (info /= 0) then
@ -277,7 +204,7 @@ contains
use psb_error_mod
! ...Subroutine Arguments
real(kind(1.d0)),pointer :: vin(:,:),vout(:,:)
real(kind(1.d0)),allocatable :: vin(:,:),vout(:,:)
integer :: info
! ...Local Variables
@ -290,7 +217,7 @@ contains
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
if (allocated(vin)) then
isz1 = size(vin,1)
isz2 = size(vin,2)
call psb_realloc(isz1,isz2,vout,info)
@ -323,7 +250,7 @@ contains
use psb_error_mod
! ...Subroutine Arguments
complex(kind(1.d0)),pointer :: vin(:),vout(:)
complex(kind(1.d0)),allocatable :: vin(:),vout(:)
integer :: info
! ...Local Variables
@ -336,7 +263,7 @@ contains
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
if (allocated(vin)) then
isz = size(vin)
call psb_realloc(isz,vout,info)
if (info /= 0) then
@ -368,7 +295,7 @@ contains
use psb_error_mod
! ...Subroutine Arguments
complex(kind(1.d0)),pointer :: vin(:,:),vout(:,:)
complex(kind(1.d0)),allocatable :: vin(:,:),vout(:,:)
integer :: info
! ...Local Variables
@ -381,7 +308,7 @@ contains
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
if (allocated(vin)) then
isz1 = size(vin,1)
isz2 = size(vin,2)
call psb_realloc(isz1,isz2,vout,info)
@ -410,17 +337,96 @@ contains
end subroutine psb_zcpy2d
function psb_isize1d(vin)
integer :: psb_isize1d
integer, pointer :: vin(:)
if (.not.associated(vin)) then
psb_isize1d = 0
else
psb_isize1d = size(vin)
end if
end function psb_isize1d
function psb_isize2d(vin,dim)
integer :: psb_isize2d
integer, pointer :: vin(:,:)
integer, optional :: dim
if (.not.associated(vin)) then
psb_isize2d = 0
else
if (present(dim)) then
psb_isize2d = size(vin,dim=dim)
else
psb_isize2d = size(vin)
end if
end if
end function psb_isize2d
function psb_dsize1d(vin)
integer :: psb_dsize1d
real(kind(1.d0)), pointer :: vin(:)
if (.not.associated(vin)) then
psb_dsize1d = 0
else
psb_dsize1d = size(vin)
end if
end function psb_dsize1d
function psb_dsize2d(vin,dim)
integer :: psb_dsize2d
real(kind(1.d0)), pointer :: vin(:,:)
integer, optional :: dim
if (.not.associated(vin)) then
psb_dsize2d = 0
else
if (present(dim)) then
psb_dsize2d = size(vin,dim=dim)
else
psb_dsize2d = size(vin)
end if
end if
end function psb_dsize2d
function psb_zsize1d(vin)
integer :: psb_zsize1d
complex(kind(1.d0)), pointer :: vin(:)
if (.not.associated(vin)) then
psb_zsize1d = 0
else
psb_zsize1d = size(vin)
end if
end function psb_zsize1d
function psb_zsize2d(vin,dim)
integer :: psb_zsize2d
complex(kind(1.d0)), pointer :: vin(:,:)
integer, optional :: dim
if (.not.associated(vin)) then
psb_zsize2d = 0
else
if (present(dim)) then
psb_zsize2d = size(vin,dim=dim)
else
psb_zsize2d = size(vin)
end if
end if
end function psb_zsize2d
Subroutine psb_dreallocate1i(len,rrax,info,pad)
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len
Integer,pointer :: rrax(:)
Integer,Intent(in) :: len
Integer,allocatable, intent(inout) :: rrax(:)
integer :: info
integer, optional, intent(in) :: pad
! ...Local Variables
Integer,Pointer :: tmp(:)
Integer,allocatable :: tmp(:)
Integer :: dim, err_act, err,i
character(len=20) :: name
logical, parameter :: debug=.false.
@ -431,7 +437,7 @@ contains
if(psb_get_errstatus().ne.0) return
info=0
if (debug) write(0,*) 'reallocate I',len
if (associated(rrax)) then
if (allocated(rrax)) then
dim=size(rrax)
If (dim /= len) Then
Allocate(tmp(len),stat=info)
@ -440,20 +446,10 @@ contains
call psb_errpush(err,name)
goto 9999
end if
if (.true.) then
do i=1, min(len,dim)
tmp(i)=rrax(i)
end do
else
tmp(1:min(len,dim))=rrax(1:min(len,dim))
end if
deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
tmp(1:min(len,dim))=rrax(1:min(len,dim))
call move_alloc(tmp,rrax)
end if
else
dim = 0
@ -490,12 +486,12 @@ contains
! ...Subroutine Arguments
Integer,Intent(in) :: len
Real(kind(1.d0)),pointer :: rrax(:)
Real(kind(1.d0)),allocatable, intent(inout) :: rrax(:)
integer :: info
real(kind(1.d0)), optional, intent(in) :: pad
! ...Local Variables
Real(kind(1.d0)),Pointer :: tmp(:)
Real(kind(1.d0)),allocatable :: tmp(:)
Integer :: dim,err_act,err,i, m
character(len=20) :: name
logical, parameter :: debug=.false.
@ -505,7 +501,7 @@ contains
info = 0
if (debug) write(0,*) 'reallocate D',len
if (associated(rrax)) then
if (allocated(rrax)) then
dim=size(rrax)
If (dim /= len) Then
@ -516,20 +512,10 @@ contains
goto 9999
end if
m = min(dim,len)
if (.true.) then
do i=1,m
tmp(i) = rrax(i)
end do
else
tmp(1:m) = rrax(1:m)
end if
Deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
tmp(1:m) = rrax(1:m)
call move_alloc(tmp,rrax)
End If
else
dim = 0
@ -564,12 +550,12 @@ contains
! ...Subroutine Arguments
Integer,Intent(in) :: len
complex(kind(1.d0)),pointer :: rrax(:)
complex(kind(1.d0)),allocatable, intent(inout):: rrax(:)
integer :: info
complex(kind(1.d0)), optional, intent(in) :: pad
! ...Local Variables
complex(kind(1.d0)),Pointer :: tmp(:)
complex(kind(1.d0)),allocatable :: tmp(:)
Integer :: dim,err_act,err,i, m
character(len=20) :: name
logical, parameter :: debug=.false.
@ -579,7 +565,7 @@ contains
info = 0
if (debug) write(0,*) 'reallocate Z',len
if (associated(rrax)) then
if (allocated(rrax)) then
dim=size(rrax)
If (dim /= len) Then
@ -590,22 +576,10 @@ contains
goto 9999
end if
m = min(dim,len)
!!$ write(0,*) 'DA: copying ',min(len,dim)
if (.true.) then
do i=1,m
tmp(i) = rrax(i)
end do
else
tmp(1:m) = rrax(1:m)
end if
!!$ write(0,*) 'DA: copying done ',m
Deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
tmp(1:m) = rrax(1:m)
call move_alloc(tmp,rrax)
End If
else
dim = 0
@ -640,12 +614,13 @@ contains
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len1,len2
Real(kind(1.d0)),pointer :: rrax(:,:)
Real(kind(1.d0)),allocatable :: rrax(:,:)
integer :: info
real(kind(1.d0)), optional, intent(in) :: pad
! ...Local Variables
Real(kind(1.d0)),Pointer :: tmp(:,:)
Real(kind(1.d0)),allocatable :: tmp(:,:)
Integer :: dim,err_act,err,i, m, dim2
character(len=20) :: name
@ -653,7 +628,7 @@ contains
call psb_erractionsave(err_act)
info = 0
if (associated(rrax)) then
if (allocated(rrax)) then
dim=size(rrax,1)
dim2=size(rrax,2)
@ -665,22 +640,10 @@ contains
goto 9999
end if
m = min(dim,len1)
!!$ write(0,*) 'DA: copying ',min(len,dim)
if (.true.) then
do i=1,m
tmp(i,1:min(len2,dim2)) = rrax(i,1:min(len2,dim2))
end do
else
tmp(1:m,1:min(len2,dim2)) = rrax(1:m,1:min(len2,dim2))
end if
!!$ write(0,*) 'DA: copying done ',m
Deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
tmp(1:m,1:min(len2,dim2)) = rrax(1:m,1:min(len2,dim2))
call move_alloc(tmp,rrax)
End If
else
dim = 0
@ -711,16 +674,19 @@ contains
End Subroutine psb_dreallocated2
Subroutine psb_dreallocatez2(len1,len2,rrax,info,pad)
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len1,len2
complex(kind(1.d0)),pointer :: rrax(:,:)
complex(kind(1.d0)),allocatable :: rrax(:,:)
integer :: info
complex(kind(1.d0)), optional, intent(in) :: pad
! ...Local Variables
complex(kind(1.d0)),Pointer :: tmp(:,:)
complex(kind(1.d0)),allocatable :: tmp(:,:)
Integer :: dim,err_act,err,i, m, dim2
character(len=20) :: name
@ -728,7 +694,7 @@ contains
call psb_erractionsave(err_act)
info = 0
if (associated(rrax)) then
if (allocated(rrax)) then
dim=size(rrax,1)
dim2=size(rrax,2)
@ -740,26 +706,13 @@ contains
goto 9999
end if
m = min(dim,len1)
!!$ write(0,*) 'DA: copying ',min(len,dim)
if (.true.) then
do i=1,m
tmp(i,1:min(len2,dim2)) = rrax(i,1:min(len2,dim2))
end do
else
tmp(1:m,1:min(len2,dim2)) = rrax(1:m,1:min(len2,dim2))
end if
!!$ write(0,*) 'DA: copying done ',m
Deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
tmp(1:m,1:min(len2,dim2)) = rrax(1:m,1:min(len2,dim2))
call move_alloc(tmp,rrax)
End If
else
dim = 0
dim2 = 0
dim = 0
Allocate(rrax(len1,len2),stat=info)
if (info /= 0) then
err=4000
@ -786,16 +739,17 @@ contains
End Subroutine psb_dreallocatez2
Subroutine psb_dreallocatei2(len1,len2,rrax,info,pad)
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len1,len2
integer,pointer :: rrax(:,:)
integer,allocatable :: rrax(:,:)
integer :: info
integer, optional, intent(in) :: pad
! ...Local Variables
integer,Pointer :: tmp(:,:)
integer,allocatable :: tmp(:,:)
Integer :: dim,err_act,err,i, m, dim2
character(len=20) :: name
@ -803,7 +757,7 @@ contains
call psb_erractionsave(err_act)
info = 0
if (associated(rrax)) then
if (allocated(rrax)) then
dim=size(rrax,1)
dim2=size(rrax,2)
@ -815,22 +769,10 @@ contains
goto 9999
end if
m = min(dim,len1)
!!$ write(0,*) 'DA: copying ',min(len,dim)
if (.true.) then
do i=1,m
tmp(i,1:min(len2,dim2)) = rrax(i,1:min(len2,dim2))
end do
else
tmp(1:m,1:min(len2,dim2)) = rrax(1:m,1:min(len2,dim2))
end if
!!$ write(0,*) 'DA: copying done ',m
Deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
tmp(1:m,1:min(len2,dim2)) = rrax(1:m,1:min(len2,dim2))
call move_alloc(tmp,rrax)
End If
else
dim = 0
@ -861,13 +803,12 @@ contains
End Subroutine psb_dreallocatei2
Subroutine psb_dreallocate2i(len,rrax,y,info,pad)
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len
Integer,pointer :: rrax(:),y(:)
Integer,allocatable, intent(inout) :: rrax(:),y(:)
integer :: info
integer, optional, intent(in) :: pad
character(len=20) :: name
@ -912,8 +853,8 @@ contains
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len
Integer,pointer :: rrax(:),y(:)
Real(Kind(1.d0)),pointer :: z(:)
Integer,allocatable, intent(inout) :: rrax(:),y(:)
Real(Kind(1.d0)),allocatable, intent(inout) :: z(:)
integer :: info
character(len=20) :: name
integer :: err_act, err
@ -961,8 +902,8 @@ contains
use psb_error_mod
! ...Subroutine Arguments
Integer,Intent(in) :: len
Integer,pointer :: rrax(:),y(:)
complex(Kind(1.d0)),pointer :: z(:)
Integer,allocatable, intent(inout) :: rrax(:),y(:)
complex(Kind(1.d0)),allocatable, intent(inout) :: z(:)
integer :: info
character(len=20) :: name
integer :: err_act, err
@ -1005,4 +946,97 @@ contains
End Subroutine psb_dreallocate2i1z
Subroutine psb_dtransfer1d(vin,vout,info)
use psb_error_mod
real(kind(1.d0)), allocatable, intent(inout) :: vin(:),vout(:)
integer, intent(out) :: info
!
! To be reimplemented with MOVE_ALLOC
!
info = 0
call move_alloc(vin,vout)
!!$
!!$ if (.not.allocated(vin) ) then
!!$ if (allocated(vout)) then
!!$ deallocate(vout,stat=info)
!!$ end if
!!$ else if (allocated(vin)) then
!!$ if (.not.allocated(vout)) then
!!$ allocate(vout(size(vin)),stat=info)
!!$ if (info /= 0) return
!!$ else
!!$ if (size(vout) /= size(vin)) then
!!$ deallocate(vout,stat=info)
!!$ if (info /= 0) return
!!$ allocate(vout(size(vin)),stat=info)
!!$ if (info /= 0) return
!!$ end if
!!$ end if
!!$ vout = vin
!!$ deallocate(vin,stat=info)
!!$ end if
end Subroutine psb_dtransfer1d
Subroutine psb_ztransfer1d(vin,vout,info)
use psb_error_mod
complex(kind(1.d0)), allocatable, intent(inout) :: vin(:),vout(:)
integer, intent(out) :: info
!
! To be reimplemented with MOVE_ALLOC
!
info = 0
call move_alloc(vin,vout)
!!$ if (.not.allocated(vin) ) then
!!$ if (allocated(vout)) then
!!$ deallocate(vout,stat=info)
!!$ end if
!!$ else if (allocated(vin)) then
!!$ if (.not.allocated(vout)) then
!!$ allocate(vout(size(vin)),stat=info)
!!$ if (info /= 0) return
!!$ else
!!$ if (size(vout) /= size(vin)) then
!!$ deallocate(vout,stat=info)
!!$ if (info /= 0) return
!!$ allocate(vout(size(vin)),stat=info)
!!$ if (info /= 0) return
!!$ end if
!!$ end if
!!$ vout = vin
!!$ deallocate(vin,stat=info)
!!$ end if
end Subroutine psb_ztransfer1d
Subroutine psb_itransfer1d(vin,vout,info)
use psb_error_mod
integer, allocatable, intent(inout) :: vin(:),vout(:)
integer, intent(out) :: info
!
! To be reimplemented with MOVE_ALLOC
!
info = 0
call move_alloc(vin,vout)
!!$ if (.not.allocated(vin) ) then
!!$ if (allocated(vout)) then
!!$ deallocate(vout,stat=info)
!!$ end if
!!$ else if (allocated(vin)) then
!!$ if (.not.allocated(vout)) then
!!$ allocate(vout(size(vin)),stat=info)
!!$ if (info /= 0) return
!!$ else
!!$ if (size(vout) /= size(vin)) then
!!$ deallocate(vout,stat=info)
!!$ if (info /= 0) return
!!$ allocate(vout(size(vin)),stat=info)
!!$ if (info /= 0) return
!!$ end if
!!$ end if
!!$ vout = vin
!!$ deallocate(vin,stat=info)
!!$ end if
end Subroutine psb_itransfer1d
end module psb_realloc_mod

8
src/modules/psb_serial_mod.f90
Unescape Escape View File

@ -55,14 +55,14 @@ module psb_serial_mod
subroutine psb_dcsrws(rw,a,info,trans)
use psb_spmat_type
type(psb_dspmat_type) :: a
real(kind(1.d0)), pointer :: rw(:)
real(kind(1.d0)), allocatable :: rw(:)
integer :: info
character, optional :: trans
end subroutine psb_dcsrws
subroutine psb_zcsrws(rw,a,info,trans)
use psb_spmat_type
type(psb_zspmat_type) :: a
complex(kind(1.d0)), pointer :: rw(:)
complex(kind(1.d0)), allocatable :: rw(:)
integer :: info
character, optional :: trans
end subroutine psb_zcsrws
@ -319,7 +319,7 @@ module psb_serial_mod
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: idx
integer, intent(out) :: n
integer, pointer :: neigh(:)
integer, allocatable :: neigh(:)
integer, intent(out) :: info
integer, optional, intent(in) :: lev
end subroutine psb_dneigh
@ -328,7 +328,7 @@ module psb_serial_mod
type(psb_zspmat_type), intent(in) :: a
integer, intent(in) :: idx
integer, intent(out) :: n
integer, pointer :: neigh(:)
integer, allocatable :: neigh(:)
integer, intent(out) :: info
integer, optional, intent(in) :: lev
end subroutine psb_zneigh

154
src/modules/psb_spmat_type.f90
Unescape Escape View File

@ -49,11 +49,11 @@ module psb_spmat_type
! Contains some additional informations on sparse matrix
integer :: infoa(psb_ifasize_)
! Contains sparse matrix coefficients
real(kind(1.d0)), pointer :: aspk(:)=>null()
real(kind(1.d0)), allocatable :: aspk(:)
! Contains indeces that describes sparse matrix structure
integer, pointer :: ia1(:)=>null(), ia2(:)=>null()
integer, allocatable :: ia1(:), ia2(:)
! Permutations matrix
integer, pointer :: pl(:)=>null(), pr(:)=>null()
integer, allocatable :: pl(:), pr(:)
end type psb_dspmat_type
type psb_zspmat_type
! Rows & columns
@ -65,11 +65,11 @@ module psb_spmat_type
! Contains some additional informations on sparse matrix
integer :: infoa(psb_ifasize_)
! Contains sparse matrix coefficients
complex(kind(1.d0)), pointer :: aspk(:)=>null()
complex(kind(1.d0)), allocatable :: aspk(:)
! Contains indeces that describes sparse matrix structure
integer, pointer :: ia1(:)=>null(), ia2(:)=>null()
integer, allocatable :: ia1(:), ia2(:)
! Permutations matrix
integer, pointer :: pl(:)=>null(), pr(:)=>null()
integer, allocatable :: pl(:), pr(:)
end type psb_zspmat_type
interface psb_nullify_sp
@ -124,8 +124,9 @@ contains
implicit none
type(psb_dspmat_type), intent(inout) :: mat
nullify(mat%aspk,mat%ia1,mat%ia2,mat%pl,mat%pr)
mat%infoa(:) = 0
!!$ nullify(mat%aspk,mat%ia1,mat%ia2,mat%pl,mat%pr)
mat%infoa(:)=0
mat%m=0
mat%k=0
mat%fida=''
@ -264,7 +265,7 @@ contains
a%m=max(0,m)
a%k=max(0,k)
call psb_sp_reall(a,nnz,info)
if (debug) write(0,*) 'Check in ALLOCATE ',info,associated(a%pl),associated(a%pr)
if (debug) write(0,*) 'Check in ALLOCATE ',info,allocated(a%pl),allocated(a%pr)
a%pl(1)=0
a%pr(1)=0
! set infoa fields
@ -348,8 +349,8 @@ contains
call psb_realloc(max(1,a%m),a%pl,info)
if (info /= 0) return
call psb_realloc(max(1,a%k),a%pr,info)
if (debug) write(0,*) associated(a%ia1),associated(a%ia2),&
& associated(a%aspk),associated(a%pl),associated(a%pr),info
if (debug) write(0,*) allocated(a%ia1),allocated(a%ia2),&
& allocated(a%aspk),allocated(a%pl),allocated(a%pr),info
if (info /= 0) return
Return
@ -416,9 +417,8 @@ contains
End Subroutine psb_dspclone
! This is done with pointer assignments, but it
! will be feasible with MOVE_ALLOC when we move
! to ALLOCATABLE components.
! Will be changed to use MOVE_ALLOC
subroutine psb_dsp_transfer(a, b,info)
implicit none
!....Parameters...
@ -431,28 +431,12 @@ contains
info = 0
if (associated(b%pr)) then
deallocate(b%pr,stat=info)
end if
if (associated(b%pl)) then
deallocate(b%pl,stat=info)
end if
if (associated(b%ia2)) then
deallocate(b%ia2,stat=info)
end if
if (associated(b%ia1)) then
deallocate(b%ia1,stat=info)
endif
if (associated(b%aspk)) then
deallocate(b%aspk,stat=info)
endif
b%aspk => a%aspk
b%ia1 => a%ia1
b%ia2 => a%ia2
b%pl => a%pl
b%pr => a%pr
call psb_transfer( a%aspk, b%aspk , info)
call psb_transfer( a%ia1 , b%ia1 , info)
call psb_transfer( a%ia2 , b%ia2 , info)
call psb_transfer( a%pl , b%pl , info)
call psb_transfer( a%pr , b%pr , info)
b%infoa(:) = a%infoa(:)
b%fida = a%fida
b%descra = a%descra
@ -557,11 +541,6 @@ contains
return
endif
!!$ if (.not.associated(a%infoa)) then
!!$ info = -2
!!$ return
!!$ endif
call psb_getifield(val,field,a%infoa,psb_ifasize_,info)
psb_dsp_getifld = val
@ -582,20 +561,20 @@ contains
val = 4*size(a%infoa)
if (associated(a%aspk)) then
if (allocated(a%aspk)) then
val = val + 8 * size(a%aspk)
endif
if (associated(a%ia1)) then
if (allocated(a%ia1)) then
val = val + 4 * size(a%ia1)
endif
if (associated(a%ia2)) then
if (allocated(a%ia2)) then
val = val + 4 * size(a%ia2)
endif
if (associated(a%pl)) then
if (allocated(a%pl)) then
val = val + 4 * size(a%pl)
endif
if (associated(a%pr)) then
if (allocated(a%pr)) then
val = val + 4 * size(a%pr)
endif
@ -613,25 +592,33 @@ contains
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
integer :: iret
info = 0
if (associated(a%aspk)) then
deallocate(a%aspk,STAT=INFO)
if (allocated(a%aspk)) then
!!$ write(0,*) 'Deallocating aspk'
deallocate(a%aspk,STAT=IRET)
!!$ write(0,*) 'Deallocated aspk',iret
if (iret /= 0) info = max(info,1)
endif
if ((info == 0) .and. associated(a%ia1)) then
deallocate(a%ia1,STAT=INFO)
if (allocated(a%ia1)) then
deallocate(a%ia1,STAT=IRET)
if (iret /= 0) info = max(info,2)
endif
if ((info == 0) .and. associated(a%ia2)) then
deallocate(a%ia2,STAT=INFO)
if (allocated(a%ia2)) then
deallocate(a%ia2,STAT=IRET)
if (iret /= 0) info = max(info,3)
endif
if ((info == 0) .and. associated(a%pr)) then
deallocate(a%pr,STAT=INFO)
if (allocated(a%pr)) then
deallocate(a%pr,STAT=IRET)
if (iret /= 0) info = max(info,4)
endif
if ((info == 0) .and. associated(a%pl)) then
deallocate(a%pl,STAT=INFO)
if (allocated(a%pl)) then
deallocate(a%pl,STAT=IRET)
if (iret /= 0) info = max(info,5)
endif
call psb_nullify_sp(a)
!!$ write(0,*) 'End of sp_free ',info
Return
End Subroutine psb_dsp_free
@ -640,8 +627,7 @@ contains
implicit none
type(psb_zspmat_type), intent(inout) :: mat
nullify(mat%aspk,mat%ia1,mat%ia2,mat%pl,mat%pr)
mat%infoa(:) = 0
mat%infoa(:)=0
mat%m=0
mat%k=0
mat%fida=''
@ -945,28 +931,11 @@ contains
info = 0
if (associated(b%pr)) then
deallocate(b%pr,stat=info)
end if
if (associated(b%pl)) then
deallocate(b%pl,stat=info)
end if
if (associated(b%ia2)) then
deallocate(b%ia2,stat=info)
end if
if (associated(b%ia1)) then
deallocate(b%ia1,stat=info)
endif
if (associated(b%aspk)) then
deallocate(b%aspk,stat=info)
endif
b%aspk => a%aspk
b%ia1 => a%ia1
b%ia2 => a%ia2
b%pl => a%pl
b%pr => a%pr
call psb_transfer( a%aspk, b%aspk , info)
call psb_transfer( a%ia1 , b%ia1 , info)
call psb_transfer( a%ia2 , b%ia2 , info)
call psb_transfer( a%pl , b%pl , info)
call psb_transfer( a%pr , b%pr , info)
b%infoa(:) = a%infoa(:)
b%fida = a%fida
b%descra = a%descra
@ -1072,11 +1041,6 @@ contains
return
endif
!!$ if (.not.associated(a%infoa)) then
!!$ info = -2
!!$ return
!!$ endif
call psb_getifield(val,field,a%infoa,psb_ifasize_,info)
psb_zsp_getifld = val
@ -1097,20 +1061,20 @@ contains
val = 4*size(a%infoa)
if (associated(a%aspk)) then
if (allocated(a%aspk)) then
val = val + 16 * size(a%aspk)
endif
if (associated(a%ia1)) then
if (allocated(a%ia1)) then
val = val + 4 * size(a%ia1)
endif
if (associated(a%ia2)) then
if (allocated(a%ia2)) then
val = val + 4 * size(a%ia2)
endif
if (associated(a%pl)) then
if (allocated(a%pl)) then
val = val + 4 * size(a%pl)
endif
if (associated(a%pr)) then
if (allocated(a%pr)) then
val = val + 4 * size(a%pr)
endif
@ -1133,19 +1097,19 @@ contains
info = 0
if (associated(a%aspk)) then
if (allocated(a%aspk)) then
deallocate(a%aspk,STAT=INFO)
endif
if ((info == 0) .and. associated(a%ia1)) then
if (allocated(a%ia1)) then
deallocate(a%ia1,STAT=INFO)
endif
if ((info == 0) .and. associated(a%ia2)) then
if ( allocated(a%ia2)) then
deallocate(a%ia2,STAT=INFO)
endif
if ((info == 0) .and. associated(a%pr)) then
if ( allocated(a%pr)) then
deallocate(a%pr,STAT=INFO)
endif
if ((info == 0) .and. associated(a%pl)) then
if ( allocated(a%pl)) then
deallocate(a%pl,STAT=INFO)
endif
call psb_nullify_sp(a)

52
src/modules/psb_tools_mod.f90
Unescape Escape View File

@ -36,7 +36,7 @@ Module psb_tools_mod
subroutine psb_dalloc(x, desc_a, info, n)
use psb_descriptor_type
implicit none
real(kind(1.d0)), pointer :: x(:,:)
real(kind(1.d0)), allocatable, intent(out) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n
@ -44,7 +44,7 @@ Module psb_tools_mod
! 1-D double precision version
subroutine psb_dallocv(x, desc_a,info,n)
use psb_descriptor_type
real(kind(1.d0)), pointer :: x(:)
real(kind(1.d0)), allocatable, intent(out) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n
@ -52,14 +52,14 @@ Module psb_tools_mod
! 2-D integer version
subroutine psb_ialloc(x, desc_a, info,n)
use psb_descriptor_type
integer, pointer :: x(:,:)
integer, allocatable, intent(out) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer, optional, intent(in) :: n
end subroutine psb_ialloc
subroutine psb_iallocv(x, desc_a,info,n)
use psb_descriptor_type
integer, pointer :: x(:)
integer, allocatable, intent(out) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n
@ -68,7 +68,7 @@ Module psb_tools_mod
subroutine psb_zalloc(x, desc_a, info, n)
use psb_descriptor_type
implicit none
complex(kind(1.d0)), pointer :: x(:,:)
complex(kind(1.d0)), allocatable, intent(out) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n
@ -76,7 +76,7 @@ Module psb_tools_mod
! 1-D double precision version
subroutine psb_zallocv(x, desc_a,info,n)
use psb_descriptor_type
complex(kind(1.d0)), pointer :: x(:)
complex(kind(1.d0)), allocatable, intent(out) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n
@ -89,42 +89,42 @@ Module psb_tools_mod
subroutine psb_dasb(x, desc_a, info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
real(kind(1.d0)), pointer :: x(:,:)
real(kind(1.d0)), allocatable, intent(inout) :: x(:,:)
integer, intent(out) :: info
end subroutine psb_dasb
! 1-D double precision version
subroutine psb_dasbv(x, desc_a, info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
real(kind(1.d0)), pointer :: x(:)
real(kind(1.d0)), allocatable, intent(inout) :: x(:)
integer, intent(out) :: info
end subroutine psb_dasbv
! 2-D integer version
subroutine psb_iasb(x, desc_a, info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
integer, pointer :: x(:,:)
integer, allocatable, intent(inout) :: x(:,:)
integer, intent(out) :: info
end subroutine psb_iasb
! 1-D integer version
subroutine psb_iasbv(x, desc_a, info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
integer, pointer :: x(:)
integer, allocatable, intent(inout) :: x(:)
integer, intent(out) :: info
end subroutine psb_iasbv
! 2-D double precision version
subroutine psb_zasb(x, desc_a, info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
complex(kind(1.d0)), pointer :: x(:,:)
complex(kind(1.d0)), allocatable, intent(inout) :: x(:,:)
integer, intent(out) :: info
end subroutine psb_zasb
! 1-D double precision version
subroutine psb_zasbv(x, desc_a, info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
complex(kind(1.d0)), pointer :: x(:)
complex(kind(1.d0)), allocatable, intent(inout) :: x(:)
integer, intent(out) :: info
end subroutine psb_zasbv
end interface
@ -209,42 +209,42 @@ Module psb_tools_mod
! 2-D double precision version
subroutine psb_dfree(x, desc_a, info)
use psb_descriptor_type
real(kind(1.d0)),pointer :: x(:,:)
real(kind(1.d0)),allocatable, intent(inout) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
end subroutine psb_dfree
! 1-D double precision version
subroutine psb_dfreev(x, desc_a, info)
use psb_descriptor_type
real(kind(1.d0)),pointer :: x(:)
real(kind(1.d0)),allocatable, intent(inout) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
end subroutine psb_dfreev
! 2-D integer version
subroutine psb_ifree(x, desc_a, info)
use psb_descriptor_type
integer,pointer :: x(:,:)
integer,allocatable, intent(inout) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
end subroutine psb_ifree
! 1-D integer version
subroutine psb_ifreev(x, desc_a, info)
use psb_descriptor_type
integer, pointer :: x(:)
integer, allocatable, intent(inout) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
end subroutine psb_ifreev
! 2-D double precision version
subroutine psb_zfree(x, desc_a, info)
use psb_descriptor_type
complex(kind(1.d0)),pointer :: x(:,:)
complex(kind(1.d0)),allocatable, intent(inout) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
end subroutine psb_zfree
! 1-D double precision version
subroutine psb_zfreev(x, desc_a, info)
use psb_descriptor_type
complex(kind(1.d0)),pointer :: x(:)
complex(kind(1.d0)),allocatable, intent(inout) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
end subroutine psb_zfreev
@ -293,7 +293,7 @@ Module psb_tools_mod
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
real(kind(1.d0)),pointer :: x(:,:)
real(kind(1.d0)),intent(inout) :: x(:,:)
integer, intent(in) :: irw(:)
real(kind(1.d0)), intent(in) :: val(:,:)
integer, intent(out) :: info
@ -304,7 +304,7 @@ Module psb_tools_mod
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
real(kind(1.d0)),pointer :: x(:)
real(kind(1.d0)),intent(inout) :: x(:)
integer, intent(in) :: irw(:)
real(kind(1.d0)), intent(in) :: val(:)
integer, intent(out) :: info
@ -315,7 +315,7 @@ Module psb_tools_mod
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
integer,pointer :: x(:,:)
integer,intent(inout) :: x(:,:)
integer, intent(in) :: irw(:)
integer, intent(in) :: val(:,:)
integer, intent(out) :: info
@ -326,7 +326,7 @@ Module psb_tools_mod
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
integer,pointer :: x(:)
integer,intent(inout) :: x(:)
integer, intent(in) :: irw(:)
integer, intent(in) :: val(:)
integer, intent(out) :: info
@ -337,7 +337,7 @@ Module psb_tools_mod
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
complex(kind(1.d0)),pointer :: x(:,:)
complex(kind(1.d0)),intent(inout) :: x(:,:)
integer, intent(in) :: irw(:)
complex(kind(1.d0)), intent(in) :: val(:,:)
integer, intent(out) :: info
@ -348,7 +348,7 @@ Module psb_tools_mod
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
complex(kind(1.d0)),pointer :: x(:)
complex(kind(1.d0)),intent(inout) :: x(:)
integer, intent(in) :: irw(:)
complex(kind(1.d0)), intent(in) :: val(:)
integer, intent(out) :: info
@ -640,7 +640,7 @@ Module psb_tools_mod
subroutine psb_get_ovrlap(ovrel,desc,info)
use psb_descriptor_type
implicit none
integer, pointer :: ovrel(:)
integer, allocatable :: ovrel(:)
type(psb_desc_type), intent(in) :: desc
integer, intent(out) :: info
end subroutine psb_get_ovrlap
@ -654,7 +654,7 @@ contains
use psb_descriptor_type
use psi_mod
implicit none
integer, pointer :: bndel(:)
integer, allocatable :: bndel(:)
type(psb_desc_type), intent(in) :: desc
integer, intent(out) :: info

111
src/modules/psi_mod.f90
Unescape Escape View File

@ -45,7 +45,7 @@ module psi_mod
interface
subroutine psi_crea_bnd_elem(bndel,desc_a,info)
use psb_descriptor_type
integer, pointer :: bndel(:)
integer, allocatable :: bndel(:)
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
end subroutine psi_crea_bnd_elem
@ -54,18 +54,19 @@ module psi_mod
interface
subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info)
use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info,nxch,nsnd,nrcv
integer, intent(in) :: index_in(:)
integer, pointer :: index_out(:)
logical :: glob_idx
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info,nxch,nsnd,nrcv
integer, intent(in) :: index_in(:)
integer, allocatable, intent(inout) :: index_out(:)
logical :: glob_idx
end subroutine psi_crea_index
end interface
interface
subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
integer :: desc_overlap(:)
integer, pointer :: ovr_elem(:)
integer, allocatable, intent(inout) :: ovr_elem(:)
integer, intent(out) :: info
end subroutine psi_crea_ovr_elem
end interface
@ -75,8 +76,8 @@ module psi_mod
& isglob_in,info)
integer :: desc_data(:),index_in(:),dep_list(:)
integer :: loc_to_glob(:),glob_to_loc(:)
integer,pointer :: desc_index(:)
integer :: length_dl,nsnd,nrcv, info
integer,allocatable, intent(inout) :: desc_index(:)
integer :: length_dl,nsnd,nrcv,info
logical :: isglob_in
end subroutine psi_desc_index
end interface
@ -92,48 +93,54 @@ module psi_mod
use psb_descriptor_type
integer, intent(in) :: flag, n
integer, intent(out) :: info
real(kind(1.d0)) :: y(:,:), beta, work(:)
type(psb_desc_type) :: desc_a
real(kind(1.d0)) :: y(:,:), beta
real(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
integer, optional :: data
end subroutine psi_dswapdatam
subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data)
use psb_descriptor_type
integer, intent(in) :: flag
integer, intent(out) :: info
real(kind(1.d0)) :: y(:), beta, work(:)
type(psb_desc_type) :: desc_a
real(kind(1.d0)) :: y(:), beta
real(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
integer, optional :: data
end subroutine psi_dswapdatav
subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data)
use psb_descriptor_type
integer, intent(in) :: flag, n
integer, intent(out) :: info
integer :: y(:,:), beta, work(:)
type(psb_desc_type) :: desc_a
integer :: y(:,:), beta
integer, target :: work(:)
type(psb_desc_type), target :: desc_a
integer, optional :: data
end subroutine psi_iswapdatam
subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data)
use psb_descriptor_type
integer, intent(in) :: flag
integer, intent(out) :: info
integer :: y(:), beta, work(:)
type(psb_desc_type) :: desc_a
integer :: y(:), beta
integer, target :: work(:)
type(psb_desc_type), target :: desc_a
integer, optional :: data
end subroutine psi_iswapdatav
subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data)
use psb_descriptor_type
integer, intent(in) :: flag, n
integer, intent(out) :: info
complex(kind(1.d0)) :: y(:,:), beta, work(:)
type(psb_desc_type) :: desc_a
complex(kind(1.d0)) :: y(:,:), beta
complex(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
integer, optional :: data
end subroutine psi_zswapdatam
subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data)
use psb_descriptor_type
integer, intent(in) :: flag
integer, intent(out) :: info
complex(kind(1.d0)) :: y(:), beta, work(:)
type(psb_desc_type) :: desc_a
complex(kind(1.d0)) :: y(:), beta
complex(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
integer, optional :: data
end subroutine psi_zswapdatav
end interface
@ -144,43 +151,49 @@ module psi_mod
use psb_descriptor_type
integer, intent(in) :: flag, n
integer, intent(out) :: info
real(kind(1.d0)) :: y(:,:), beta, work(:)
type(psb_desc_type) :: desc_a
real(kind(1.d0)) :: y(:,:), beta
real(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
end subroutine psi_dswaptranm
subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info)
use psb_descriptor_type
integer, intent(in) :: flag
integer, intent(out) :: info
real(kind(1.d0)) :: y(:), beta, work(:)
type(psb_desc_type) :: desc_a
real(kind(1.d0)) :: y(:), beta
real(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
end subroutine psi_dswaptranv
subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info)
use psb_descriptor_type
integer, intent(in) :: flag, n
integer, intent(out) :: info
integer :: y(:,:), beta, work(:)
type(psb_desc_type) :: desc_a
integer :: y(:,:), beta
integer,target :: work(:)
type(psb_desc_type), target :: desc_a
end subroutine psi_iswaptranm
subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info)
use psb_descriptor_type
integer, intent(in) :: flag
integer, intent(out) :: info
integer :: y(:), beta, work(:)
type(psb_desc_type) :: desc_a
integer :: y(:), beta
integer,target :: work(:)
type(psb_desc_type), target :: desc_a
end subroutine psi_iswaptranv
subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info)
use psb_descriptor_type
integer, intent(in) :: flag, n
integer, intent(out) :: info
complex(kind(1.d0)) :: y(:,:), beta, work(:)
type(psb_desc_type) :: desc_a
complex(kind(1.d0)) :: y(:,:), beta
complex(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
end subroutine psi_zswaptranm
subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info)
use psb_descriptor_type
integer, intent(in) :: flag
integer, intent(out) :: info
complex(kind(1.d0)) :: y(:), beta, work(:)
type(psb_desc_type) :: desc_a
complex(kind(1.d0)) :: y(:), beta
complex(kind(1.d0)),target :: work(:)
type(psb_desc_type), target :: desc_a
end subroutine psi_zswaptranv
end interface
@ -252,6 +265,7 @@ contains
use psb_error_mod
use psb_penv_mod
use psb_descriptor_type
use psb_realloc_mod
implicit none
! ....scalars parameters....
@ -264,7 +278,7 @@ contains
integer :: ictxt, err_act,nxch,nsnd,nrcv
! ...local array...
integer :: int_err(5)
integer, pointer :: idx_out(:)
integer, allocatable :: idx_out(:)
! ...parameters
logical, parameter :: debug=.false.
@ -286,13 +300,12 @@ contains
! first the halo index
if (debug) write(0,*) me,'Calling crea_index on halo'
idx_out => null()
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')
goto 9999
end if
cdesc%halo_index => idx_out
call psb_transfer(idx_out,cdesc%halo_index,info)
cdesc%matrix_data(psb_thal_xch_) = nxch
cdesc%matrix_data(psb_thal_snd_) = nsnd
cdesc%matrix_data(psb_thal_rcv_) = nrcv
@ -301,30 +314,35 @@ contains
if (debug) write(0,*) me,'Calling crea_index on ovrlap'
! then the overlap index
idx_out => null()
call psi_crea_index(cdesc,ovrlap_in, idx_out,.true.,nxch,nsnd,nrcv,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
cdesc%ovrlap_index => idx_out
call psb_transfer(idx_out,cdesc%ovrlap_index,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_transfer')
goto 9999
end if
cdesc%matrix_data(psb_tovr_xch_) = nxch
cdesc%matrix_data(psb_tovr_snd_) = nsnd
cdesc%matrix_data(psb_tovr_rcv_) = nrcv
if (debug) write(0,*) me,'Calling crea_ovr_elem'
! next ovrlap_elem
call psi_crea_ovr_elem(cdesc%ovrlap_index,cdesc%ovrlap_elem)
call psi_crea_ovr_elem(cdesc%ovrlap_index,cdesc%ovrlap_elem,info)
if (debug) write(0,*) me,'Done crea_ovr_elem'
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_ovr_elem')
goto 9999
end if
! finally bnd_elem
idx_out => null()
call psi_crea_bnd_elem(idx_out,cdesc,info)
if (associated(idx_out)) then
cdesc%bnd_elem => idx_out
else
cdesc%bnd_elem => null()
endif
if (info == 0) call psb_transfer(idx_out,cdesc%bnd_elem,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_bnd_elem')
goto 9999
@ -344,5 +362,4 @@ contains
end subroutine psi_cnv_dsc
end module psi_mod

10
src/prec/psb_dbaseprc_bld.f90
Unescape Escape View File

@ -158,11 +158,11 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd)
call psb_check_def(p%iprcparm(p_type_),'base_prec',&
& diagsc_,is_legal_base_prec)
allocate(p%desc_data,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
!!$ allocate(p%desc_data,stat=info)
!!$ if (info /= 0) then
!!$ call psb_errpush(4010,name,a_err='Allocate')
!!$ goto 9999
!!$ end if
call psb_nullify_desc(p%desc_data)

265
src/prec/psb_dbldaggrmat.f90
Unescape Escape View File

@ -34,7 +34,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dbldaggrmat(a,desc_a,ac,p,desc_p,info)
subroutine psb_dbldaggrmat(a,desc_a,ac,desc_ac,p,info)
use psb_serial_mod
use psb_penv_mod
use psb_prec_type
@ -46,10 +46,10 @@ subroutine psb_dbldaggrmat(a,desc_a,ac,p,desc_p,info)
implicit none
type(psb_dspmat_type), intent(in), target :: a
type(psb_dbaseprc_type), intent(inout) :: p
type(psb_dspmat_type), intent(out), target :: ac
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_p
type(psb_desc_type), intent(inout) :: desc_ac
type(psb_dbaseprc_type), intent(inout), target :: p
integer, intent(out) :: info
logical, parameter :: aggr_dump=.false.
@ -112,10 +112,9 @@ contains
include 'mpif.h'
integer, intent(out) :: info
type(psb_dspmat_type), pointer :: bg
type(psb_dspmat_type) :: b, tmp
integer, pointer :: nzbr(:), idisp(:)
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzbg, ip, ndx,&
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, np, me, nzt,irs,jl,nzl,nlr,&
& icomm,naggrm1, i, j, k, err_act
@ -124,7 +123,6 @@ contains
info=0
call psb_erractionsave(err_act)
bg => ac
call psb_nullify_sp(b)
ictxt = desc_a%matrix_data(psb_ctxt_)
@ -175,7 +173,7 @@ contains
b%fida = 'COO'
b%m=a%m
b%k=a%k
if (.false.) then
if (.true.) then
call psb_csdp(a,b,info)
if(info /= 0) then
info=4010
@ -196,38 +194,38 @@ contains
enddo
else
! Ok, this is extremely dirty because we use pointers from
! one sparse matrix into another. But it gives us something
! in term of performance
jl = 0
do i=1,a%m,50
nlr = min(a%m-i+1,50)
call psb_spgtblk(i,a,b,info,append=.true.,iren=p%mlia,lrw=i+nlr-1)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spgtblk')
goto 9999
end if
call psb_spinfo(psb_nztotreq_,b,nzl,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spinfo')
goto 9999
end if
nzl = nzl - jl
tmp%fida = 'COO'
tmp%infoa(psb_nnz_) = nzl
tmp%aspk => b%aspk(jl+1:jl+nzl)
tmp%ia1 => b%ia1(jl+1:jl+nzl)
tmp%ia2 => b%ia2(jl+1:jl+nzl)
call psb_fixcoo(tmp,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_fixcoo')
goto 9999
end if
nzl = tmp%infoa(psb_nnz_)
b%infoa(psb_nnz_) = jl+nzl
jl = jl + nzl
enddo
!!$ ! Ok, this is extremely dirty because we use pointers from
!!$ ! one sparse matrix into another. But it gives us something
!!$ ! in term of performance
!!$ jl = 0
!!$ do i=1,a%m,50
!!$ nlr = min(a%m-i+1,50)
!!$ call psb_spgtblk(i,a,b,info,append=.true.,iren=p%mlia,lrw=i+nlr-1)
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='spgtblk')
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_spinfo(psb_nztotreq_,b,nzl,info)
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='spinfo')
!!$ goto 9999
!!$ end if
!!$ nzl = nzl - jl
!!$ tmp%fida = 'COO'
!!$ tmp%infoa(psb_nnz_) = nzl
!!$ tmp%aspk => b%aspk(jl+1:jl+nzl)
!!$ tmp%ia1 => b%ia1(jl+1:jl+nzl)
!!$ tmp%ia2 => b%ia2(jl+1:jl+nzl)
!!$ call psb_fixcoo(tmp,info)
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='psb_fixcoo')
!!$ goto 9999
!!$ end if
!!$ nzl = tmp%infoa(psb_nnz_)
!!$ b%infoa(psb_nnz_) = jl+nzl
!!$ jl = jl + nzl
!!$ enddo
end if
@ -242,7 +240,7 @@ contains
if (p%iprcparm(coarse_mat_) == mat_repl_) then
call psb_cdrep(ntaggr,ictxt,desc_p,info)
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
@ -251,8 +249,8 @@ contains
nzbr(:) = 0
nzbr(me+1) = irs
call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,bg,nzbg,info)
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
@ -264,11 +262,11 @@ contains
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,bg%aspk,nzbr,idisp,&
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,bg%ia1,nzbr,idisp,&
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,bg%ia2,nzbr,idisp,&
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
@ -276,13 +274,12 @@ contains
goto 9999
end if
bg%m = ntaggr
bg%k = ntaggr
bg%infoa(psb_nnz_) = nzbg
bg%fida='COO'
bg%descra='G'
call psb_sp_free(b,info)
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
call psb_fixcoo(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
goto 9999
@ -290,14 +287,13 @@ contains
else if (p%iprcparm(coarse_mat_) == mat_distr_) then
call psb_cddec(naggr,ictxt,desc_p,info)
call psb_cddec(naggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cddec')
goto 9999
end if
call psb_sp_clone(b,bg,info)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
@ -310,23 +306,23 @@ contains
!if(.not.associated(p%av(ap_nd_)%aspk)) p%iprcparm(jac_sweeps_) = 1
!------------------------------------------------------------------
! Split BG=M+N N off-diagonal part
call psb_sp_all(bg%m,bg%k,p%av(ap_nd_),nzl,info)
! Split AC=M+N N off-diagonal part
call psb_sp_all(ac%m,ac%k,p%av(ap_nd_),nzl,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
end if
if(.not.associated(p%av(ap_nd_)%aspk)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
if(.not.associated(p%av(ap_nd_)%ia1)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
if(.not.allocated(p%av(ap_nd_)%aspk)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
if(.not.allocated(p%av(ap_nd_)%ia1)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
!write(0,*) 'ok line 238'
k=0
do i=1,nzl
if (bg%ia2(i)>bg%m) then
if (ac%ia2(i)>ac%m) then
k = k + 1
p%av(ap_nd_)%aspk(k) = bg%aspk(i)
p%av(ap_nd_)%ia1(k) = bg%ia1(i)
p%av(ap_nd_)%ia2(k) = bg%ia2(i)
p%av(ap_nd_)%aspk(k) = ac%aspk(i)
p%av(ap_nd_)%ia1(k) = ac%ia1(i)
p%av(ap_nd_)%ia2(k) = ac%ia2(i)
endif
enddo
p%av(ap_nd_)%infoa(psb_nnz_) = k
@ -352,19 +348,17 @@ contains
goto 9999
end if
else
else
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_)
end if
call psb_ipcoo2csr(bg,info)
call psb_ipcoo2csr(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='ipcoo2csr')
goto 9999
end if
deallocate(nzbr,idisp)
call psb_erractionrestore(err_act)
@ -394,10 +388,9 @@ contains
integer, intent(out) :: info
type(psb_dspmat_type), pointer :: bg
type(psb_dspmat_type) :: b
integer, pointer :: nzbr(:), idisp(:), ivall(:)
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzbg, ip, ndx,&
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, np, me, &
& icomm, naggrm1,naggrp1,i,j,err_act,k,nzl
type(psb_dspmat_type), pointer :: am1,am2
@ -418,7 +411,6 @@ contains
ictxt = desc_a%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
bg => ac
call psb_nullify_sp(b)
call psb_nullify_sp(am3)
call psb_nullify_sp(am4)
@ -616,6 +608,15 @@ contains
goto 9999
end if
if (test_dump) then
open(30+me)
write(30+me,*) 'OMEGA: ',omega
do i=1,size(p%dorig)
write(30+me,*) p%dorig(i)
end do
close(30+me)
end if
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,&
@ -633,9 +634,6 @@ contains
call psb_errpush(4010,name,a_err='symbmm 1')
goto 9999
end if
!!$ am1%aspk(:) = 0.d0
if (test_dump) &
& call psb_csprt(50+me,am1,head='% (I-wDA)Pt symbmm ')
call psb_numbmm(am3,am4,am1)
@ -785,10 +783,10 @@ contains
case(mat_distr_)
call psb_sp_clone(b,bg,info)
call psb_sp_clone(b,ac,info)
if(info /= 0) goto 9999
nzbg = bg%infoa(psb_nnz_)
nzl = bg%infoa(psb_nnz_)
nzac = ac%infoa(psb_nnz_)
nzl = ac%infoa(psb_nnz_)
allocate(ivall(ntaggr),stat=info)
if (info /= 0) then
@ -803,21 +801,22 @@ contains
i = i + 1
end do
end do
call psb_cdall(ntaggr,ivall,ictxt,desc_p,info,flag=1)
call psb_cdall(ntaggr,ivall,ictxt,desc_ac,info,flag=1)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdins(nzl,bg%ia1,bg%ia2,desc_p,info)
call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdins')
goto 9999
end if
if (debug) write(0,*) me,'Created aux descr. distr.'
call psb_cdasb(desc_p,info)
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
@ -826,24 +825,24 @@ contains
if (debug) write(0,*) me,'Asmbld aux descr. distr.'
call psb_glob_to_loc(bg%ia1(1:nzl),desc_p,info,iact='I')
call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
goto 9999
end if
call psb_glob_to_loc(bg%ia2(1:nzl),desc_p,info,iact='I')
call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
goto 9999
end if
bg%m=desc_p%matrix_data(psb_n_row_)
bg%k=desc_p%matrix_data(psb_n_col_)
bg%fida='COO'
bg%descra='G'
ac%m=desc_ac%matrix_data(psb_n_row_)
ac%k=desc_ac%matrix_data(psb_n_col_)
ac%fida='COO'
ac%descra='G'
call psb_sp_free(b,info)
if(info /= 0) then
@ -854,8 +853,8 @@ contains
deallocate(ivall,nzbr,idisp)
! Split BG=M+N N off-diagonal part
call psb_sp_all(bg%m,bg%k,p%av(ap_nd_),nzl,info)
! Split AC=M+N N off-diagonal part
call psb_sp_all(ac%m,ac%k,p%av(ap_nd_),nzl,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
@ -863,11 +862,11 @@ contains
k=0
do i=1,nzl
if (bg%ia2(i)>bg%m) then
if (ac%ia2(i)>ac%m) then
k = k + 1
p%av(ap_nd_)%aspk(k) = bg%aspk(i)
p%av(ap_nd_)%ia1(k) = bg%ia1(i)
p%av(ap_nd_)%ia2(k) = bg%ia2(i)
p%av(ap_nd_)%aspk(k) = ac%aspk(i)
p%av(ap_nd_)%ia1(k) = ac%ia1(i)
p%av(ap_nd_)%ia2(k) = ac%ia2(i)
endif
enddo
p%av(ap_nd_)%infoa(psb_nnz_) = k
@ -889,13 +888,13 @@ contains
if (np>1) then
call psb_spinfo(psb_nztotreq_,am1,nzl,info)
call psb_glob_to_loc(am1%ia1(1:nzl),desc_p,info,'I')
call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
end if
endif
am1%k=desc_p%matrix_data(psb_n_col_)
am1%k=desc_ac%matrix_data(psb_n_col_)
if (np>1) then
call psb_ipcsr2coo(am2,info)
@ -905,7 +904,7 @@ contains
end if
nzl = am2%infoa(psb_nnz_)
call psb_glob_to_loc(am2%ia1(1:nzl),desc_p,info,'I')
call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
@ -917,23 +916,24 @@ contains
goto 9999
end if
end if
am2%m=desc_p%matrix_data(psb_n_col_)
am2%m=desc_ac%matrix_data(psb_n_col_)
if (debug) write(0,*) me,'Done ac '
case(mat_repl_)
!
!
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_cdrep(ntaggr,ictxt,desc_p,info)
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,bg,nzbg,info)
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) goto 9999
@ -943,26 +943,26 @@ contains
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,bg%aspk,nzbr,idisp,&
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,bg%ia1,nzbr,idisp,&
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,bg%ia2,nzbr,idisp,&
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) goto 9999
bg%m = ntaggr
bg%k = ntaggr
bg%infoa(psb_nnz_) = nzbg
bg%fida='COO'
bg%descra='G'
call psb_fixcoo(bg,info)
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
call psb_fixcoo(ac,info)
if(info /= 0) goto 9999
call psb_sp_free(b,info)
if(info /= 0) goto 9999
if (me==0) then
if (test_dump) call psb_csprt(80+me,bg,head='% Smoothed aggregate AC.')
if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.')
endif
deallocate(nzbr,idisp)
@ -978,12 +978,12 @@ contains
case(mat_distr_)
call psb_sp_clone(b,bg,info)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
end if
call psb_cddec(naggr,ictxt,desc_p,info)
call psb_cddec(naggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cddec')
goto 9999
@ -999,19 +999,18 @@ contains
case(mat_repl_)
!
!
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_cdrep(ntaggr,ictxt,desc_p,info)
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,bg,nzbg,info)
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
@ -1023,11 +1022,11 @@ contains
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,bg%aspk,nzbr,idisp,&
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,bg%ia1,nzbr,idisp,&
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,bg%ia2,nzbr,idisp,&
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
@ -1036,12 +1035,12 @@ contains
end if
bg%m = ntaggr
bg%k = ntaggr
bg%infoa(psb_nnz_) = nzbg
bg%fida='COO'
bg%descra='G'
call psb_fixcoo(bg,info)
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
call psb_fixcoo(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_fixcoo')
goto 9999
@ -1056,8 +1055,14 @@ contains
deallocate(nzbr,idisp)
end select
call psb_ipcoo2csr(bg,info)
call psb_ipcoo2csr(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_ipcoo2csr')
goto 9999
end if
if (debug) write(0,*) me,'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return
@ -1073,6 +1078,4 @@ contains
end subroutine smooth_aggregate
end subroutine psb_dbldaggrmat

2
src/prec/psb_ddiagsc_bld.f90
Unescape Escape View File

@ -58,7 +58,7 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
! Local scalars
Integer :: err, n_row, n_col,I,j,k,ictxt,&
& me,np,mglob,lw, err_act
real(kind(1.d0)),pointer :: gd(:), work(:)
real(kind(1.d0)),allocatable :: gd(:), work(:)
integer :: int_err(5)
character :: iupd

4
src/prec/psb_dgenaggrmap.f90
Unescape Escape View File

@ -44,10 +44,10 @@ subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info)
integer, intent(in) :: aggr_type
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, pointer :: ilaggr(:),nlaggr(:)
integer, allocatable :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Locals
integer, pointer :: ils(:), neigh(:)
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
logical :: recovery

10
src/prec/psb_dilu_bld.f90
Unescape Escape View File

@ -58,6 +58,7 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info)
use psb_tools_mod
use psb_psblas_mod
use psb_error_mod
use psb_realloc_mod
use psb_penv_mod
implicit none
!
@ -170,7 +171,7 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info)
if (debug) write(0,*)me,': out of psb_asmatbld'
if (debug) call psb_barrier(ictxt)
if (associated(p%av)) then
if (allocated(p%av)) then
if (size(p%av) < bp_ilu_avsz) then
call psb_errpush(4010,name,a_err='Insufficient av size')
goto 9999
@ -179,6 +180,7 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info)
call psb_errpush(4010,name,a_err='AV not associated')
goto 9999
endif
!!$ call psb_csprt(50+me,a,head='% (A)')
nrow_a = desc_a%matrix_data(psb_n_row_)
call psb_spinfo(psb_nztotreq_,a,nztota,info)
@ -208,12 +210,12 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (associated(p%d)) then
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
if (.not.associated(p%d)) then
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
@ -335,6 +337,8 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info)
close(80+me)
endif
!!$ call psb_csprt(60+me,a,head='% (A)')
! ierr = MPE_Log_event( ifcte, 0, "st SIMPLE" )
t6 = mpi_wtime()

1
src/prec/psb_dilu_fct.f90
Unescape Escape View File

@ -87,6 +87,7 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
blck_%m=0
endif
!!$ write(0,*) 'ilu_fct: ',size(l%ia2),size(u%ia2),a%m,blck_%m
call psb_dilu_fctint(m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if(info.ne.0) then

49
src/prec/psb_dmlprc_aply.f90
Unescape Escape View File

@ -105,8 +105,6 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,n_col
real(kind(1.d0)), allocatable :: tx(:),ty(:),x2l(:),y2l(:),&
& b2l(:),tty(:)
character ::diagl, diagu
integer :: ictxt,np,me,i, isz, nrg,nr2l,err_act, iptype, int_err(5)
real(kind(1.d0)) :: omega
@ -117,11 +115,9 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name, ch_err
type psb_mlprec_wrk_type
real(kind(1.d0)), pointer :: tx(:)=>null(),ty(:)=>null(),&
& x2l(:)=>null(),y2l(:)=>null(),&
& b2l(:)=>null(),tty(:)=>null()
real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type psb_mlprec_wrk_type
type(psb_mlprec_wrk_type), pointer :: mlprec_wrk(:)
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
interface psb_baseprc_aply
subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
@ -455,8 +451,6 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
mlprec_wrk(1)%y2l(:) = dzero
mlprec_wrk(1)%x2l(:) = x
call psb_baseprc_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,&
@ -744,7 +738,6 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
case default
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
@ -760,8 +753,6 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end select
call mlprec_wrk_free(mlprec_wrk)
deallocate(mlprec_wrk)
call psb_erractionrestore(err_act)
@ -776,24 +767,24 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
return
contains
subroutine mlprec_wrk_free(wrk)
type(psb_mlprec_wrk_type) :: wrk(:)
! This will not be needed when we have allocatables, as
! it is sufficient to deallocate the container, and
! the compiler is supposed to recursively deallocate the
! various components.
integer i
do i=1, size(wrk)
if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx)
if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty)
if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l)
if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l)
if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l)
if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty)
end do
end subroutine mlprec_wrk_free
!!$contains
!!$ subroutine mlprec_wrk_free(wrk)
!!$ type(psb_mlprec_wrk_type) :: wrk(:)
!!$ ! This will not be needed when we have allocatables, as
!!$ ! it is sufficient to deallocate the container, and
!!$ ! the compiler is supposed to recursively deallocate the
!!$ ! various components.
!!$ integer i
!!$
!!$ do i=1, size(wrk)
!!$ if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx)
!!$ if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty)
!!$ if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l)
!!$ if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l)
!!$ if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l)
!!$ if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty)
!!$ end do
!!$ end subroutine mlprec_wrk_free
end subroutine psb_dmlprc_aply

38
src/prec/psb_dmlprc_bld.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
type(psb_dbaseprc_type), intent(inout),target :: p
integer, intent(out) :: info
type(psb_desc_type), pointer :: desc_p
type(psb_desc_type), pointer :: desc_ac
integer :: i, nrg, nzg, err_act,k
character(len=20) :: name, ch_err
@ -77,21 +77,21 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
integer, intent(in) :: aggr_type
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, pointer :: ilaggr(:),nlaggr(:)
integer, allocatable :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine psb_dgenaggrmap
end interface
interface psb_bldaggrmat
subroutine psb_dbldaggrmat(a,desc_a,ac,p,desc_p,info)
subroutine psb_dbldaggrmat(a,desc_a,ac,desc_ac,p,info)
use psb_prec_type
use psb_descriptor_type
use psb_spmat_type
type(psb_dspmat_type), intent(in), target :: a
type(psb_dbaseprc_type), intent(inout),target :: p
type(psb_dspmat_type), intent(out),target :: ac
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_p
type(psb_dspmat_type), intent(out),target :: ac
type(psb_desc_type), intent(inout) :: desc_ac
type(psb_dbaseprc_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine psb_dbldaggrmat
end interface
@ -105,7 +105,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
call psb_nullify_sp(ac)
if (.not.associated(p%iprcparm)) then
if (.not.allocated(p%iprcparm)) then
info = 2222
call psb_errpush(info,name)
goto 9999
@ -122,7 +122,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
& pre_smooth_,is_legal_ml_smooth_pos)
nullify(p%desc_data)
!!$ nullify(p%desc_data)
select case(p%iprcparm(f_type_))
case(f_ilu_n_)
call psb_check_def(p%iprcparm(ilu_fill_in_),'Level',0,is_legal_ml_lev)
@ -134,10 +134,6 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
& 1,is_legal_jac_sweeps)
nullify(p%d)
! Currently this is ignored by gen_aggrmap, but it could be
! changed in the future. Need to package nlaggr & mlia in a
! private data structure?
@ -150,21 +146,21 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
end if
if (debug) write(0,*) 'Out from genaggrmap',p%nlaggr
nullify(desc_p)
allocate(desc_p)
call psb_nullify_desc(desc_p)
call psb_bldaggrmat(a,desc_a,ac,p,desc_p,info)
nullify(desc_ac)
allocate(desc_ac)
call psb_nullify_desc(desc_ac)
call psb_bldaggrmat(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
info=4010
ch_err='psb_bld_aggrmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) 'Out from bldaggrmat',desc_p%matrix_data(:)
if (debug) write(0,*) 'Out from bldaggrmat',desc_ac%matrix_data(:)
call psb_baseprc_bld(ac,desc_p,p,info)
call psb_baseprc_bld(ac,desc_ac,p,info)
if (debug) write(0,*) 'Out from baseprcbld',info
if(info /= 0) then
info=4010
@ -182,9 +178,9 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
! Hence a separate AC and a TRANSFER function at the end.
!
call psb_sp_transfer(ac,p%av(ac_),info)
p%base_a => p%av(ac_)
p%desc_ac => desc_p
nullify(desc_p)
p%base_a => p%av(ac_)
call psb_cdtransfer(desc_ac,p%desc_ac,info)
p%base_desc => p%desc_ac
call psb_erractionrestore(err_act)
return

2
src/prec/psb_dprc_aply.f90
Unescape Escape View File

@ -112,7 +112,7 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
end if
if (.not.(associated(prec%baseprecv))) then
if (.not.(allocated(prec%baseprecv))) then
write(0,*) 'Inconsistent preconditioner: neither SMTH nor BASE?'
end if
if (size(prec%baseprecv) >1) then

10
src/prec/psb_dprecbld.f90
Unescape Escape View File

@ -83,7 +83,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
end interface
! Local scalars
Integer :: err, I,j,k,ictxt, me,np,lw, err_act
Integer :: err,i,j,k,ictxt, me,np,lw, err_act
integer :: int_err(5)
character :: iupd
@ -116,10 +116,10 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
iupd='F'
endif
if (.not.associated(p%baseprecv)) then
if (.not.allocated(p%baseprecv)) then
!! Error 1: should call precset
info=4010
ch_err='unassociated bpv'
ch_err='unallocated bpv'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
@ -157,7 +157,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
call psb_mlprc_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,&
& p%baseprecv(i),info)
if (info /= 0) then
@ -189,7 +189,7 @@ contains
subroutine init_baseprc_av(p,info)
type(psb_dbaseprc_type), intent(inout) :: p
integer :: info
if (associated(p%av)) then
if (allocated(p%av)) then
! Have not decided what to do yet
end if
allocate(p%av(max_avsz),stat=info)

3
src/prec/psb_dprecfree.f90
Unescape Escape View File

@ -59,12 +59,11 @@ subroutine psb_dprecfree(p,info)
me=-1
if (associated(p%baseprecv)) then
if (allocated(p%baseprecv)) then
do i=1,size(p%baseprecv)
call psb_base_precfree(p%baseprecv(i),info)
end do
deallocate(p%baseprecv)
p%baseprecv => null()
end if
call psb_erractionrestore(err_act)
return

35
src/prec/psb_dprecset.f90
Unescape Escape View File

@ -49,7 +49,6 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
type(psb_dbaseprc_type), pointer :: bpv(:)=>null()
character(len=len(ptype)) :: typeup
integer :: isz, err, nlev_, ilev_, i
@ -61,18 +60,16 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
ilev_ = 1
end if
if (present(nlev)) then
if (associated(p%baseprecv)) then
write(0,*) 'Warning: NLEV is ignored when P is already associated'
if (allocated(p%baseprecv)) then
write(0,*) 'Warning: NLEV is ignored when P is already allocated'
end if
nlev_ = max(1, nlev)
else
nlev_ = 1
end if
if (.not.associated(p%baseprecv)) then
if (.not.allocated(p%baseprecv)) then
allocate(p%baseprecv(nlev_),stat=err)
do i=1, nlev_
call psb_nullify_baseprec(p%baseprecv(i))
end do
else
nlev_ = size(p%baseprecv)
endif
@ -82,14 +79,12 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
info = -1
return
endif
if (.not.associated(p%baseprecv(ilev_)%iprcparm)) then
allocate(p%baseprecv(ilev_)%iprcparm(ifpsz),&
& p%baseprecv(ilev_)%dprcparm(dfpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure',err
endif
end if
call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
select case(toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
@ -147,17 +142,8 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
case ('ML', '2L', '2LEV')
!!$ allocate(p%baseprecv(ilev_)%iprcparm(ifpsz),stat=err)
!!$ if (err/=0) then
!!$ write(0,*)'Precset Memory Failure 2l:2',err
!!$ endif
!!$ allocate(p%baseprecv(ilev_)%dprcparm(dfpsz),stat=err)
!!$ if (err/=0) then
!!$ write(0,*)'Precset Memory Failure 2l:3',err
!!$ endif
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bja_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
@ -175,7 +161,6 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(ml_type_) = iv(1)

4
src/prec/psb_dsp_renum.f90
Unescape Escape View File

@ -58,8 +58,8 @@ subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info)
integer nztota, nztotb, nztmp, nzl, nnr, ir, mglob, mtype, n_row, &
& nrow_a,n_col, nhalo,lovr, ind, iind, pi,nr,ns,i,j,jj,k,kk
integer ::ictxt,np,me, err_act
integer, pointer :: itmp(:), itmp2(:)
real(kind(1.d0)), pointer :: rtmp(:)
integer, allocatable :: itmp(:), itmp2(:)
real(kind(1.d0)), allocatable :: rtmp(:)
real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,mpi_wtime, t7, t8
external mpi_wtime

14
src/prec/psb_zbaseprc_bld.f90
Unescape Escape View File

@ -158,13 +158,13 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd)
call psb_check_def(p%iprcparm(p_type_),'base_prec',&
& diagsc_,is_legal_base_prec)
allocate(p%desc_data,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
call psb_nullify_desc(p%desc_data)
!!$ allocate(p%desc_data,stat=info)
!!$ if (info /= 0) then
!!$ call psb_errpush(4010,name,a_err='Allocate')
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(p_type_))
case (noprec_)

289
src/prec/psb_zbldaggrmat.f90
Unescape Escape View File

@ -34,7 +34,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zbldaggrmat(a,desc_a,ac,p,desc_p,info)
subroutine psb_zbldaggrmat(a,desc_a,ac,desc_ac,p,info)
use psb_serial_mod
use psb_penv_mod
use psb_prec_type
@ -46,10 +46,10 @@ subroutine psb_zbldaggrmat(a,desc_a,ac,p,desc_p,info)
implicit none
type(psb_zspmat_type), intent(in), target :: a
type(psb_zbaseprc_type), intent(inout) :: p
type(psb_zbaseprc_type), intent(inout),target :: p
type(psb_zspmat_type), intent(out), target :: ac
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout),target :: desc_p
type(psb_desc_type), intent(inout) :: desc_ac
integer, intent(out) :: info
logical, parameter :: aggr_dump=.false.
@ -111,18 +111,17 @@ contains
include 'mpif.h'
integer, intent(out) :: info
type(psb_zspmat_type), pointer :: bg
type(psb_zspmat_type) :: b, tmp
integer, pointer :: nzbr(:), idisp(:)
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzbg, ip, ndx,&
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, np, me, nzt,irs,jl,nzl,nlr,&
& icomm,naggrm1, i, j, k, err_act
name='raw_aggregate'
if(psb_get_errstatus().ne.0) return
info=0
call psb_erractionsave(err_act)
bg => ac
call psb_nullify_sp(b)
ictxt = desc_a%matrix_data(psb_ctxt_)
@ -173,7 +172,7 @@ contains
b%fida = 'COO'
b%m=a%m
b%k=a%k
if (.false.) then
if (.true.) then
call psb_csdp(a,b,info)
if(info /= 0) then
info=4010
@ -194,38 +193,38 @@ contains
enddo
else
! Ok, this is extremely dirty because we use pointers from
! one sparse matrix into another. But it gives us something
! in term of performance
jl = 0
do i=1,a%m,50
nlr = min(a%m-i+1,50)
call psb_spgtblk(i,a,b,info,append=.true.,iren=p%mlia,lrw=i+nlr-1)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spgtblk')
goto 9999
end if
call psb_spinfo(psb_nztotreq_,b,nzl,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spinfo')
goto 9999
end if
nzl = nzl - jl
tmp%fida = 'COO'
tmp%infoa(psb_nnz_) = nzl
tmp%aspk => b%aspk(jl+1:jl+nzl)
tmp%ia1 => b%ia1(jl+1:jl+nzl)
tmp%ia2 => b%ia2(jl+1:jl+nzl)
call psb_fixcoo(tmp,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_fixcoo')
goto 9999
end if
nzl = tmp%infoa(psb_nnz_)
b%infoa(psb_nnz_) = jl+nzl
jl = jl + nzl
enddo
!!$ ! Ok, this is extremely dirty because we use pointers from
!!$ ! one sparse matrix into another. But it gives us something
!!$ ! in term of performance
!!$ jl = 0
!!$ do i=1,a%m,50
!!$ nlr = min(a%m-i+1,50)
!!$ call psb_spgtblk(i,a,b,info,append=.true.,iren=p%mlia,lrw=i+nlr-1)
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='spgtblk')
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_spinfo(psb_nztotreq_,b,nzl,info)
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='spinfo')
!!$ goto 9999
!!$ end if
!!$ nzl = nzl - jl
!!$ tmp%fida = 'COO'
!!$ tmp%infoa(psb_nnz_) = nzl
!!$ tmp%aspk => b%aspk(jl+1:jl+nzl)
!!$ tmp%ia1 => b%ia1(jl+1:jl+nzl)
!!$ tmp%ia2 => b%ia2(jl+1:jl+nzl)
!!$ call psb_fixcoo(tmp,info)
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='psb_fixcoo')
!!$ goto 9999
!!$ end if
!!$ nzl = tmp%infoa(psb_nnz_)
!!$ b%infoa(psb_nnz_) = jl+nzl
!!$ jl = jl + nzl
!!$ enddo
end if
call psb_fixcoo(b,info)
@ -245,13 +244,17 @@ contains
if (p%iprcparm(coarse_mat_) == mat_repl_) then
call psb_cdrep(ntaggr,ictxt,desc_p,info)
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = irs
call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,bg,nzbg,info)
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
@ -263,11 +266,11 @@ contains
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,bg%aspk,nzbr,idisp,&
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,&
& mpi_double_complex,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,bg%ia1,nzbr,idisp,&
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,bg%ia2,nzbr,idisp,&
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
@ -275,13 +278,12 @@ contains
goto 9999
end if
bg%m = ntaggr
bg%k = ntaggr
bg%infoa(psb_nnz_) = nzbg
bg%fida='COO'
bg%descra='G'
call psb_sp_free(b,info)
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
call psb_fixcoo(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
goto 9999
@ -289,8 +291,12 @@ contains
else if (p%iprcparm(coarse_mat_) == mat_distr_) then
call psb_cddec(naggr,ictxt,desc_p,info)
call psb_sp_clone(b,bg,info)
call psb_cddec(naggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cddec')
goto 9999
end if
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
@ -303,23 +309,23 @@ contains
!if(.not.associated(p%av(ap_nd_)%aspk)) p%iprcparm(jac_sweeps_) = 1
!------------------------------------------------------------------
! Split BG=M+N N off-diagonal part
call psb_sp_all(bg%m,bg%k,p%av(ap_nd_),nzl,info)
! Split AC=M+N N off-diagonal part
call psb_sp_all(ac%m,ac%k,p%av(ap_nd_),nzl,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
end if
if(.not.associated(p%av(ap_nd_)%aspk)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
if(.not.associated(p%av(ap_nd_)%ia1)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
if(.not.allocated(p%av(ap_nd_)%aspk)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
if(.not.allocated(p%av(ap_nd_)%ia1)) write(0,*) '.not.associated(p%av(ap_nd_)%ia1)'
!write(0,*) 'ok line 238'
k=0
do i=1,nzl
if (bg%ia2(i)>bg%m) then
if (ac%ia2(i)>ac%m) then
k = k + 1
p%av(ap_nd_)%aspk(k) = bg%aspk(i)
p%av(ap_nd_)%ia1(k) = bg%ia1(i)
p%av(ap_nd_)%ia2(k) = bg%ia2(i)
p%av(ap_nd_)%aspk(k) = ac%aspk(i)
p%av(ap_nd_)%ia1(k) = ac%ia1(i)
p%av(ap_nd_)%ia2(k) = ac%ia2(i)
endif
enddo
p%av(ap_nd_)%infoa(psb_nnz_) = k
@ -345,19 +351,17 @@ contains
goto 9999
end if
else
else
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_)
end if
call psb_ipcoo2csr(bg,info)
call psb_ipcoo2csr(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='ipcoo2csr')
goto 9999
end if
deallocate(nzbr,idisp)
call psb_erractionrestore(err_act)
@ -387,11 +391,10 @@ contains
integer, intent(out) :: info
type(psb_zspmat_type), pointer :: bg
type(psb_zspmat_type) :: b
integer, pointer :: nzbr(:), idisp(:), ivall(:)
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzbg, ip, ndx,&
& naggr, np, me,&
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, np, me, &
& icomm, naggrm1,naggrp1,i,j,err_act,k,nzl
type(psb_zspmat_type), pointer :: am1,am2
type(psb_zspmat_type) :: am3,am4
@ -411,7 +414,6 @@ contains
ictxt = desc_a%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
bg => ac
call psb_nullify_sp(b)
call psb_nullify_sp(am3)
call psb_nullify_sp(am4)
@ -523,8 +525,6 @@ contains
endif
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
call psb_ipcoo2csr(am4,info)
@ -542,7 +542,7 @@ contains
!
! WARNING: the cycles below assume that AM3 does have
! its diagonal elements stored explicitly!!!
! Should we swicth to something safer?
! Should we switch to something safer?
!
call psb_spscal(am3,p%dorig,info)
if(info /= 0) goto 9999
@ -604,12 +604,24 @@ contains
am3%aspk(j) = done - omega*am3%aspk(j)
endif
end do
call psb_ipcoo2csr(am3,info)
else
write(0,*) 'Missing implementation of I sum'
call psb_errpush(4010,name)
goto 9999
end if
if (test_dump) then
open(30+me)
write(30+me,*) 'OMEGA: ',omega
do i=1,size(p%dorig)
write(30+me,*) p%dorig(i)
end do
close(30+me)
end if
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,&
& ivc=desc_a%loc_to_glob)
if (debug) write(0,*) me,'Done gather, going for SYMBMM 1'
@ -622,7 +634,7 @@ contains
!
call psb_symbmm(am3,am4,am1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm')
call psb_errpush(4010,name,a_err='symbmm 1')
goto 9999
end if
@ -674,7 +686,7 @@ contains
call psb_symbmm(a,am1,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm')
call psb_errpush(4010,name,a_err='symbmm 2')
goto 9999
end if
@ -736,7 +748,7 @@ contains
if (debug) write(0,*) me,'starting symbmm 3'
call psb_symbmm(am2,am3,b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm')
call psb_errpush(4010,name,a_err='symbmm 3')
goto 9999
end if
@ -774,10 +786,10 @@ contains
case(mat_distr_)
call psb_sp_clone(b,bg,info)
call psb_sp_clone(b,ac,info)
if(info /= 0) goto 9999
nzbg = bg%infoa(psb_nnz_)
nzl = bg%infoa(psb_nnz_)
nzac = ac%infoa(psb_nnz_)
nzl = ac%infoa(psb_nnz_)
allocate(ivall(ntaggr),stat=info)
if (info /= 0) then
@ -792,21 +804,22 @@ contains
i = i + 1
end do
end do
call psb_cdall(ntaggr,ivall,ictxt,desc_p,info,flag=1)
call psb_cdall(ntaggr,ivall,ictxt,desc_ac,info,flag=1)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdins(nzl,bg%ia1,bg%ia2,desc_p,info)
call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdins')
goto 9999
end if
if (debug) write(0,*) me,'Created aux descr. distr.'
call psb_cdasb(desc_p,info)
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
@ -815,24 +828,24 @@ contains
if (debug) write(0,*) me,'Asmbld aux descr. distr.'
call psb_glob_to_loc(bg%ia1(1:nzl),desc_p,info,iact='I')
call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
goto 9999
end if
call psb_glob_to_loc(bg%ia2(1:nzl),desc_p,info,iact='I')
call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
goto 9999
end if
bg%m=desc_p%matrix_data(psb_n_row_)
bg%k=desc_p%matrix_data(psb_n_col_)
bg%fida='COO'
bg%descra='G'
ac%m=desc_ac%matrix_data(psb_n_row_)
ac%k=desc_ac%matrix_data(psb_n_col_)
ac%fida='COO'
ac%descra='G'
call psb_sp_free(b,info)
if(info /= 0) then
@ -843,8 +856,8 @@ contains
deallocate(ivall,nzbr,idisp)
! Split BG=M+N N off-diagonal part
call psb_sp_all(bg%m,bg%k,p%av(ap_nd_),nzl,info)
! Split AC=M+N N off-diagonal part
call psb_sp_all(ac%m,ac%k,p%av(ap_nd_),nzl,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
@ -852,11 +865,11 @@ contains
k=0
do i=1,nzl
if (bg%ia2(i)>bg%m) then
if (ac%ia2(i)>ac%m) then
k = k + 1
p%av(ap_nd_)%aspk(k) = bg%aspk(i)
p%av(ap_nd_)%ia1(k) = bg%ia1(i)
p%av(ap_nd_)%ia2(k) = bg%ia2(i)
p%av(ap_nd_)%aspk(k) = ac%aspk(i)
p%av(ap_nd_)%ia1(k) = ac%ia1(i)
p%av(ap_nd_)%ia2(k) = ac%ia2(i)
endif
enddo
p%av(ap_nd_)%infoa(psb_nnz_) = k
@ -878,13 +891,13 @@ contains
if (np>1) then
call psb_spinfo(psb_nztotreq_,am1,nzl,info)
call psb_glob_to_loc(am1%ia1(1:nzl),desc_p,info,'I')
call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
end if
endif
am1%k=desc_p%matrix_data(psb_n_col_)
am1%k=desc_ac%matrix_data(psb_n_col_)
if (np>1) then
call psb_ipcsr2coo(am2,info)
@ -894,7 +907,7 @@ contains
end if
nzl = am2%infoa(psb_nnz_)
call psb_glob_to_loc(am2%ia1(1:nzl),desc_p,info,'I')
call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
@ -906,19 +919,23 @@ contains
goto 9999
end if
end if
am2%m=desc_p%matrix_data(psb_n_col_)
am2%m=desc_ac%matrix_data(psb_n_col_)
case(mat_repl_)
!
!
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_cdrep(ntaggr,ictxt,desc_p,info)
call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,bg,nzbg,info)
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) goto 9999
@ -928,26 +945,26 @@ contains
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,bg%aspk,nzbr,idisp,&
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,&
& mpi_double_complex,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,bg%ia1,nzbr,idisp,&
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,bg%ia2,nzbr,idisp,&
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) goto 9999
bg%m = ntaggr
bg%k = ntaggr
bg%infoa(psb_nnz_) = nzbg
bg%fida='COO'
bg%descra='G'
call psb_fixcoo(bg,info)
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
call psb_fixcoo(ac,info)
if(info /= 0) goto 9999
call psb_sp_free(b,info)
if(info /= 0) goto 9999
if (me==0) then
if (test_dump) call psb_csprt(80+me,bg,head='% Smoothed aggregate AC.')
if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.')
endif
deallocate(nzbr,idisp)
@ -963,12 +980,16 @@ contains
case(mat_distr_)
call psb_sp_clone(b,bg,info)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
end if
call psb_cddec(naggr,ictxt,desc_p,info)
call psb_cddec(naggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cddec')
goto 9999
end if
call psb_sp_free(b,info)
if(info /= 0) then
@ -980,15 +1001,18 @@ contains
case(mat_repl_)
!
!
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_cdrep(ntaggr,ictxt,desc_p,info)
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,bg,nzbg,info)
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
@ -1000,11 +1024,11 @@ contains
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,bg%aspk,nzbr,idisp,&
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,&
& mpi_double_complex,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,bg%ia1,nzbr,idisp,&
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,bg%ia2,nzbr,idisp,&
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
@ -1013,12 +1037,12 @@ contains
end if
bg%m = ntaggr
bg%k = ntaggr
bg%infoa(psb_nnz_) = nzbg
bg%fida='COO'
bg%descra='G'
call psb_fixcoo(bg,info)
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
call psb_fixcoo(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_fixcoo')
goto 9999
@ -1034,6 +1058,13 @@ contains
end select
call psb_ipcoo2csr(ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_ipcoo2csr')
goto 9999
end if
if (debug) write(0,*) me,'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return

4
src/prec/psb_zgenaggrmap.f90
Unescape Escape View File

@ -44,10 +44,10 @@ subroutine psb_zgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info)
integer, intent(in) :: aggr_type
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, pointer :: ilaggr(:),nlaggr(:)
integer, allocatable :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Locals
integer, pointer :: ils(:), neigh(:)
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
logical :: recovery

6
src/prec/psb_zilu_bld.f90
Unescape Escape View File

@ -175,7 +175,7 @@ subroutine psb_zilu_bld(a,desc_a,p,upd,info)
if (debug) write(0,*)me,': out of psb_asmatbld'
if (debug) call psb_barrier(ictxt)
if (associated(p%av)) then
if (allocated(p%av)) then
if (size(p%av) < bp_ilu_avsz) then
call psb_errpush(4010,name,a_err='Insufficient av size')
goto 9999
@ -213,12 +213,12 @@ subroutine psb_zilu_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (associated(p%d)) then
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
if (.not.associated(p%d)) then
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')

49
src/prec/psb_zmlprc_aply.f90
Unescape Escape View File

@ -104,8 +104,6 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,n_col
complex(kind(1.d0)), allocatable :: tx(:),ty(:),t2l(:),w2l(:),&
& x2l(:),b2l(:),tz(:),tty(:)
character ::diagl, diagu
integer :: ictxt,np,me,i, isz, nrg,nr2l,err_act, iptype, int_err(5)
real(kind(1.d0)) :: omega
@ -116,11 +114,9 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name, ch_err
type psb_mlprec_wrk_type
complex(kind(1.d0)), pointer :: tx(:)=>null(),ty(:)=>null(),&
& x2l(:)=>null(),y2l(:)=>null(),&
& b2l(:)=>null(),tty(:)=>null()
complex(kind(1.d0)), allocatable :: tx(:),ty(:),x2l(:),y2l(:)
end type psb_mlprec_wrk_type
type(psb_mlprec_wrk_type), pointer :: mlprec_wrk(:)
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
interface psb_baseprc_aply
subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
@ -454,8 +450,6 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
mlprec_wrk(1)%y2l(:) = zzero
mlprec_wrk(1)%x2l(:) = x
call psb_baseprc_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,&
@ -743,7 +737,6 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
case default
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
@ -759,8 +752,6 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end select
call mlprec_wrk_free(mlprec_wrk)
deallocate(mlprec_wrk)
call psb_erractionrestore(err_act)
@ -775,24 +766,24 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
return
contains
subroutine mlprec_wrk_free(wrk)
type(psb_mlprec_wrk_type) :: wrk(:)
! This will not be needed when we have allocatables, as
! it is sufficient to deallocate the container, and
! the compiler is supposed to recursively deallocate the
! various components.
integer i
do i=1, size(wrk)
if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx)
if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty)
if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l)
if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l)
if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l)
if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty)
end do
end subroutine mlprec_wrk_free
!!$contains
!!$ subroutine mlprec_wrk_free(wrk)
!!$ type(psb_mlprec_wrk_type) :: wrk(:)
!!$ ! This will not be needed when we have allocatables, as
!!$ ! it is sufficient to deallocate the container, and
!!$ ! the compiler is supposed to recursively deallocate the
!!$ ! various components.
!!$ integer i
!!$
!!$ do i=1, size(wrk)
!!$ if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx)
!!$ if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty)
!!$ if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l)
!!$ if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l)
!!$ if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l)
!!$ if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty)
!!$ end do
!!$ end subroutine mlprec_wrk_free
end subroutine psb_zmlprc_aply

36
src/prec/psb_zmlprc_bld.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
type(psb_zbaseprc_type), intent(inout),target :: p
integer, intent(out) :: info
type(psb_desc_type), pointer :: desc_p
type(psb_desc_type), pointer :: desc_ac
integer :: i, nrg, nzg, err_act,k
character(len=20) :: name, ch_err
@ -77,13 +77,13 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
integer, intent(in) :: aggr_type
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, pointer :: ilaggr(:),nlaggr(:)
integer, allocatable :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine psb_zgenaggrmap
end interface
interface psb_bldaggrmat
subroutine psb_zbldaggrmat(a,desc_a,ac,p,desc_p,info)
subroutine psb_zbldaggrmat(a,desc_a,ac,desc_ac,p,info)
use psb_prec_type
use psb_descriptor_type
use psb_spmat_type
@ -91,7 +91,7 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
type(psb_zbaseprc_type), intent(inout),target :: p
type(psb_zspmat_type), intent(out),target :: ac
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_p
type(psb_desc_type), intent(inout) :: desc_ac
integer, intent(out) :: info
end subroutine psb_zbldaggrmat
end interface
@ -105,7 +105,7 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
call psb_nullify_sp(ac)
if (.not.associated(p%iprcparm)) then
if (.not.allocated(p%iprcparm)) then
info = 2222
call psb_errpush(info,name)
goto 9999
@ -122,7 +122,7 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
& pre_smooth_,is_legal_ml_smooth_pos)
nullify(p%desc_data)
!!$ nullify(p%desc_data)
select case(p%iprcparm(f_type_))
case(f_ilu_n_)
call psb_check_def(p%iprcparm(ilu_fill_in_),'Level',0,is_legal_ml_lev)
@ -134,10 +134,6 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
& 1,is_legal_jac_sweeps)
nullify(p%d)
! Currently this is ignored by gen_aggrmap, but it could be
! changed in the future. Need to package nlaggr & mlia in a
! private data structure?
@ -150,22 +146,22 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
end if
if (debug) write(0,*) 'Out from genaggrmap',p%nlaggr
nullify(desc_p)
allocate(desc_p)
call psb_nullify_desc(desc_p)
call psb_bldaggrmat(a,desc_a,ac,p,desc_p,info)
nullify(desc_ac)
allocate(desc_ac)
call psb_nullify_desc(desc_ac)
call psb_bldaggrmat(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
info=4010
ch_err='psb_bld_aggrmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) 'Out from bldaggrmat',desc_p%matrix_data(:)
if (debug) write(0,*) 'Out from bldaggrmat',desc_ac%matrix_data(:)
call psb_baseprc_bld(ac,desc_p,p,info)
if (debug) write(0,*) 'Out from basaeprcbld',info
call psb_baseprc_bld(ac,desc_ac,p,info)
if (debug) write(0,*) 'Out from baseprcbld',info
if(info /= 0) then
info=4010
ch_err='psb_baseprc_bld'
@ -182,9 +178,9 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
! Hence a separate AC and a TRANSFER function at the end.
!
call psb_sp_transfer(ac,p%av(ac_),info)
call psb_cdfree(desc_p,info)
deallocate(desc_p)
p%base_a => p%av(ac_)
call psb_cdtransfer(desc_ac,p%desc_ac,info)
p%base_desc => p%desc_ac
call psb_erractionrestore(err_act)
return

2
src/prec/psb_zprc_aply.f90
Unescape Escape View File

@ -112,7 +112,7 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
end if
if (.not.(associated(prec%baseprecv))) then
if (.not.(allocated(prec%baseprecv))) then
write(0,*) 'Inconsistent preconditioner: neither SMTH nor BASE?'
end if
if (size(prec%baseprecv) >1) then

9
src/prec/psb_zprecbld.f90
Unescape Escape View File

@ -116,10 +116,10 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
iupd='F'
endif
if (.not.associated(p%baseprecv)) then
if (.not.allocated(p%baseprecv)) then
!! Error 1: should call precset
info=4010
ch_err='unassociated bpv'
ch_err='unallocated bpv'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
@ -128,7 +128,6 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
!
! ALso should define symbolic names for the preconditioners.
!
if (size(p%baseprecv) >= 1) then
call init_baseprc_av(p%baseprecv(1),info)
if (info /= 0) then
@ -190,11 +189,11 @@ contains
subroutine init_baseprc_av(p,info)
type(psb_zbaseprc_type), intent(inout) :: p
integer :: info
if (associated(p%av)) then
if (allocated(p%av)) then
! Have not decided what to do yet
end if
allocate(p%av(max_avsz),stat=info)
if (info /= 0) return
!!$ if (info /= 0) return
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))
end do

3
src/prec/psb_zprecfree.f90
Unescape Escape View File

@ -59,12 +59,11 @@ subroutine psb_zprecfree(p,info)
me=-1
if (associated(p%baseprecv)) then
if (allocated(p%baseprecv)) then
do i=1,size(p%baseprecv)
call psb_base_precfree(p%baseprecv(i),info)
end do
deallocate(p%baseprecv)
p%baseprecv => null()
end if
call psb_erractionrestore(err_act)
return

209
src/prec/psb_zprecset.f90
Unescape Escape View File

@ -34,7 +34,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zprecset(p,ptype,info,iv,rs,rv)
subroutine psb_zprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
use psb_serial_mod
use psb_descriptor_type
@ -45,148 +45,135 @@ subroutine psb_zprecset(p,ptype,info,iv,rs,rv)
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: iv(:)
integer, optional, intent(in) :: nlev,ilev
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
type(psb_zbaseprc_type), pointer :: bpv(:)=>null()
character(len=len(ptype)) :: typeup
integer :: isz, err
integer :: isz, err, nlev_, ilev_, i
info = 0
if (.not.associated(p%baseprecv)) then
allocate(p%baseprecv(1),stat=err)
call psb_nullify_baseprec(p%baseprecv(1))
if (present(ilev)) then
ilev_ = max(1, ilev)
else
ilev_ = 1
end if
if (present(nlev)) then
if (allocated(p%baseprecv)) then
write(0,*) 'Warning: NLEV is ignored when P is already allocated'
end if
nlev_ = max(1, nlev)
else
nlev_ = 1
end if
if (.not.allocated(p%baseprecv)) then
allocate(p%baseprecv(nlev_),stat=err)
else
nlev_ = size(p%baseprecv)
endif
if (.not.associated(p%baseprecv(1)%iprcparm)) then
allocate(p%baseprecv(1)%iprcparm(ifpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure',err
endif
end if
if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1
return
endif
call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
select case(toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(1)%iprcparm(p_type_) = noprec_
p%baseprecv(1)%iprcparm(f_type_) = f_none_
p%baseprecv(1)%iprcparm(restr_) = psb_none_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('DIAG','DIAGSC')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(1)%iprcparm(p_type_) = diagsc_
p%baseprecv(1)%iprcparm(f_type_) = f_none_
p%baseprecv(1)%iprcparm(restr_) = psb_none_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = diagsc_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('BJA','ILU')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(1)%iprcparm(p_type_) = bja_
p%baseprecv(1)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(1)%iprcparm(restr_) = psb_none_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 0
p%baseprecv(1)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bja_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('ASM','AS')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
! Defaults first
p%baseprecv(1)%iprcparm(p_type_) = asm_
p%baseprecv(1)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(1)%iprcparm(restr_) = psb_halo_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 1
p%baseprecv(1)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(1)%iprcparm(n_ovr_) = iv(1)
if (isz >= 2) p%baseprecv(1)%iprcparm(restr_) = iv(2)
if (isz >= 3) p%baseprecv(1)%iprcparm(prol_) = iv(3)
if (isz >= 4) p%baseprecv(1)%iprcparm(f_type_) = iv(4)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(n_ovr_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(restr_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(prol_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(4)
! Do not consider renum for the time being.
!!$ if (isz >= 5) p%baseprecv(1)%iprcparm(iren_) = iv(5)
!!$ if (isz >= 5) p%baseprecv(ilev_)%iprcparm(iren_) = iv(5)
end if
case ('ML', '2L','2LEV')
select case (size(p%baseprecv))
case(1)
! Reallocate
allocate(bpv(2),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure 2l:1',err
endif
bpv(1) = p%baseprecv(1)
call psb_nullify_baseprec(bpv(2))
deallocate(p%baseprecv)
p%baseprecv => bpv
nullify(bpv)
case(2)
! Do nothing
case default
! Error
end select
allocate(p%baseprecv(2)%iprcparm(ifpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure 2l:2',err
endif
allocate(p%baseprecv(2)%dprcparm(dfpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure 2l:3',err
endif
p%baseprecv(2)%iprcparm(:) = 0
p%baseprecv(2)%iprcparm(p_type_) = bja_
p%baseprecv(2)%iprcparm(restr_) = psb_none_
p%baseprecv(2)%iprcparm(prol_) = psb_none_
p%baseprecv(2)%iprcparm(iren_) = 0
p%baseprecv(2)%iprcparm(n_ovr_) = 0
p%baseprecv(2)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(2)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(2)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(2)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(2)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(2)%iprcparm(glb_smth_) = 1
p%baseprecv(2)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(2)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(2)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(2)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(2)%iprcparm(jac_sweeps_) = 1
case ('ML', '2L', '2LEV')
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bja_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(2)%iprcparm(ml_type_) = iv(1)
if (isz >= 2) p%baseprecv(2)%iprcparm(aggr_alg_) = iv(2)
if (isz >= 3) p%baseprecv(2)%iprcparm(coarse_mat_) = iv(3)
if (isz >= 4) p%baseprecv(2)%iprcparm(smth_pos_) = iv(4)
if (isz >= 5) p%baseprecv(2)%iprcparm(f_type_) = iv(5)
if (isz >= 6) p%baseprecv(2)%iprcparm(jac_sweeps_) = iv(6)
if (isz >= 7) p%baseprecv(2)%iprcparm(smth_kind_) = iv(7)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(ml_type_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(aggr_alg_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(coarse_mat_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(smth_pos_) = iv(4)
if (isz >= 5) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(5)
if (isz >= 6) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(6)
if (isz >= 7) p%baseprecv(ilev_)%iprcparm(smth_kind_) = iv(7)
end if
if (present(rs)) then
p%baseprecv(2)%iprcparm(om_choice_) = user_choice_
p%baseprecv(2)%dprcparm(smooth_omega_) = rs
p%baseprecv(ilev_)%iprcparm(om_choice_) = user_choice_
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = rs
end if

4
src/prec/psb_zsp_renum.f90
Unescape Escape View File

@ -58,8 +58,8 @@ subroutine psb_zsp_renum(a,desc_a,blck,p,atmp,info)
integer nztota, nztotb, nztmp, nzl, nnr, ir, mglob, mtype, n_row, &
& nrow_a,n_col, nhalo,lovr, ind, iind, pi,nr,ns,i,j,jj,k,kk
integer ::ictxt,np,me, err_act
integer, pointer :: itmp(:), itmp2(:)
complex(kind(1.d0)), pointer :: ztmp(:)
integer, allocatable :: itmp(:), itmp2(:)
complex(kind(1.d0)), allocatable :: ztmp(:)
real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,mpi_wtime, t7, t8
external mpi_wtime

4
src/psblas/psb_daxpby.f90
Unescape Escape View File

@ -135,7 +135,7 @@ subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
if ((in.ne.0)) then
if(desc_a%matrix_data(psb_n_row_).gt.0) then
call daxpby(desc_a%matrix_data(psb_n_col_),in,&
call daxpby(desc_a%matrix_data(psb_n_row_),in,&
& alpha,x(iix,jjx),size(x,1),beta,&
& y(iiy,jjy),size(y,1),info)
end if
@ -263,7 +263,7 @@ subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
end if
if(desc_a%matrix_data(psb_n_row_).gt.0) then
call daxpby(desc_a%matrix_data(psb_n_col_),ione,&
call daxpby(desc_a%matrix_data(psb_n_row_),ione,&
& alpha,x,size(x),beta,&
& y,size(y),info)
end if

101
src/psblas/psb_ddot.f90
Unescape Escape View File

@ -580,3 +580,104 @@ subroutine psb_dmdots(res, x, y, desc_a, info)
end if
return
end subroutine psb_dmdots
subroutine psb_ddot2v(res, x, y,w,z,desc_a, info)
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_penv_mod
implicit none
real(kind(1.d0)), intent(in) :: x(:), y(:),w(:), z(:)
real(kind(1.d0)), intent(out) :: res(:)
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
! locals
integer :: int_err(5), ictxt, np, me,&
& err_act, n, iix, jjx, ix, ijx, iy, ijy, iiy, jjy, i, m, j, k
real(kind(1.D0)) :: dot_local(2)
real(kind(1.d0)) :: ddot
character(len=20) :: name, ch_err
name='psb_ddot'
if(psb_get_errstatus().ne.0) return
info=0
call psb_erractionsave(err_act)
ictxt=desc_a%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
if (np == -ione) then
info = 2010
call psb_errpush(info,name)
goto 9999
endif
ix = ione
iy = ione
m = desc_a%matrix_data(psb_m_)
! check vector correctness
call psb_chkvect(m,ione,size(x,1),ix,ix,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ione,size(y,1),iy,iy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if ((iix.ne.ione).or.(iiy.ne.ione)) then
info=3040
call psb_errpush(info,name)
goto 9999
end if
if(m.ne.0) then
if(desc_a%matrix_data(psb_n_row_).gt.0) then
dot_local(1) = ddot(desc_a%matrix_data(psb_n_row_),&
& x,ione,y,ione)
dot_local(2) = ddot(desc_a%matrix_data(psb_n_row_),&
& w,ione,z,ione)
! adjust dot_local because overlapped elements are computed more than once
i=1
do while (desc_a%ovrlap_elem(i).ne.-ione)
dot_local(1) = dot_local(1) -&
& (desc_a%ovrlap_elem(i+1)-1)/desc_a%ovrlap_elem(i+1)*&
& x(desc_a%ovrlap_elem(i))*&
& y(desc_a%ovrlap_elem(i))
dot_local(2) = dot_local(2) -&
& (desc_a%ovrlap_elem(i+1)-1)/desc_a%ovrlap_elem(i+1)*&
& w(desc_a%ovrlap_elem(i))*&
& z(desc_a%ovrlap_elem(i))
i = i+2
end do
else
dot_local=0.d0
end if
else
dot_local=0.d0
end if
! compute global sum
call psb_sum(ictxt, dot_local)
res(1:2) = dot_local(1:2)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(ictxt)
return
end if
return
end subroutine psb_ddot2v

4
src/psblas/psb_dspsm.f90
Unescape Escape View File

@ -202,7 +202,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
end if
if (aliw) then
call psb_realloc(liwork,iwork,info)
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -493,7 +493,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
end if
if (aliw) then
call psb_realloc(liwork,iwork,info)
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'

2
src/psblas/psb_zspmm.f90
Unescape Escape View File

@ -190,7 +190,7 @@ subroutine psb_zspmm(alpha,a,x,beta,y,desc_a,info,&
end if
if (aliw) then
call psb_realloc(liwork,iwork,info)
allocate(iwork(liwork),stat=info)
if(info /= 0) then
info=4010
ch_err='psb_realloc'

4
src/psblas/psb_zspsm.f90
Unescape Escape View File

@ -206,7 +206,7 @@ subroutine psb_zspsm(alpha,a,x,beta,y,desc_a,info,&
end if
if (aliw) then
call psb_realloc(liwork,iwork,info)
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
@ -493,7 +493,7 @@ subroutine psb_zspsv(alpha,a,x,beta,y,desc_a,info,&
end if
if (aliw) then
call psb_realloc(liwork,iwork,info)
allocate(iwork(liwork),stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'

6
src/serial/f77/smmp.f
Unescape Escape View File

@ -22,7 +22,7 @@ c
* ib(*), jb(*), diagb,
* diagc,
* index(*)
integer, pointer :: ic(:),jc(:)
integer, allocatable :: ic(:),jc(:)
integer :: nze, info
integer, save :: iunit=11
c
@ -40,7 +40,11 @@ c$$$ write(iunit,*) 'Row:',i,' : ',jb(ib(i):ib(i+1)-1)
c$$$ enddo
if (size(ic) < n+1) then
write(0,*) 'Called realloc in SYMBMM '
call psb_realloc(n+1,ic,info)
if (info /=0) then
write(0,*) 'realloc failed in SYMBMM ',info
end if
endif
maxlmn = max(l,m,n)
do 10 i=1,maxlmn

10
src/serial/psb_dcsdp.f90
Unescape Escape View File

@ -195,11 +195,11 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (.not.associated(b%aspk).or.&
&.not.associated(b%ia1).or.&
&.not.associated(b%ia2).or.&
&.not.associated(b%pl).or.&
&.not.associated(b%pr)) then
if (.not.allocated(b%aspk).or.&
&.not.allocated(b%ia1).or.&
&.not.allocated(b%ia2).or.&
&.not.allocated(b%pl).or.&
&.not.allocated(b%pr)) then
call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info)
else if ((size(b%aspk) < aspk_size) .or.&
&(size(b%ia1) < ia1_size) .or.&

2
src/serial/psb_dcsrws.f90
Unescape Escape View File

@ -38,7 +38,7 @@ subroutine psb_dcsrws(rw,a,info,trans)
implicit none
type(psb_dspmat_type) :: a
real(kind(1.d0)), pointer :: rw(:)
real(kind(1.d0)), allocatable :: rw(:)
integer :: info
character, optional :: trans

7
src/serial/psb_dipcoo2csc.f90
Unescape Escape View File

@ -38,6 +38,7 @@ subroutine psb_dipcoo2csc(a,info,clshr)
use psb_serial_mod, only : psb_fixcoo
use psb_error_mod
use psb_string_mod
use psb_realloc_mod
implicit none
!....Parameters...
@ -45,7 +46,7 @@ subroutine psb_dipcoo2csc(a,info,clshr)
Integer, intent(out) :: info
logical, optional :: clshr
integer, pointer :: iaux(:), itemp(:)
integer, allocatable :: iaux(:), itemp(:)
!locals
logical :: clshr_
Integer :: nza, nr, i,j, idl,err_act,nc,icl
@ -76,8 +77,8 @@ subroutine psb_dipcoo2csc(a,info,clshr)
allocate(iaux(nc+1))
if(debug) write(0,*)'DIPCOO2CSC: out of fixcoo',nza,nc,size(a%ia2),size(iaux)
itemp => a%ia2
a%ia2 => iaux
call psb_transfer(a%ia2,itemp,info)
call psb_transfer(iaux,a%ia2,info)
!
! This routine can be used in two modes:

9
src/serial/psb_dipcoo2csr.f90
Unescape Escape View File

@ -35,6 +35,7 @@
subroutine psb_dipcoo2csr(a,info,rwshr)
use psb_spmat_type
use psb_const_mod
use psb_realloc_mod
use psb_serial_mod, only : psb_fixcoo
use psb_error_mod
use psb_string_mod
@ -45,7 +46,7 @@ subroutine psb_dipcoo2csr(a,info,rwshr)
Integer, intent(out) :: info
logical, optional :: rwshr
integer, pointer :: iaux(:), itemp(:)
integer, allocatable :: iaux(:), itemp(:)
!locals
logical :: rwshr_
Integer :: nza, nr, i,j,irw, idl,err_act
@ -81,9 +82,9 @@ subroutine psb_dipcoo2csr(a,info,rwshr)
if(debug) write(0,*)'DIPCOO2CSR: out of fixcoo',nza,nr,size(a%ia2),size(iaux)
itemp => a%ia1
a%ia1 => a%ia2
a%ia2 => iaux
call psb_transfer(a%ia1,itemp,info)
call psb_transfer(a%ia2,a%ia1,info)
call psb_transfer(iaux,a%ia2,info)
!
! This routine can be used in two modes:

17
src/serial/psb_dipcsr2coo.f90
Unescape Escape View File

@ -37,6 +37,7 @@ Subroutine psb_dipcsr2coo(a,info)
use psb_const_mod
use psb_error_mod
use psb_string_mod
use psb_realloc_mod
implicit none
!....Parameters...
@ -44,11 +45,11 @@ Subroutine psb_dipcsr2coo(a,info)
Integer, intent(out) :: info
!locals
Integer :: nza, nr
integer :: i,j,err_act
logical, parameter :: debug=.false.
integer, pointer :: iaux(:), itemp(:)
character(len=20) :: name
Integer :: nza, nr
integer :: i,j,err_act
logical, parameter :: debug=.false.
integer, allocatable :: iaux(:), itemp(:)
character(len=20) :: name, ch_err
name='psb_dipcsr2coo'
info = 0
@ -68,9 +69,9 @@ Subroutine psb_dipcsr2coo(a,info)
return
end if
!!$ write(0,*) 'ipcsr2coo ',a%m
itemp => a%ia2
a%ia2 => a%ia1
a%ia1 => iaux
call psb_transfer(a%ia2,itemp,info)
call psb_transfer(a%ia1,a%ia2,info)
call psb_transfer(iaux,a%ia1,info)
do i=1, nr
do j=itemp(i),itemp(i+1)-1

19
src/serial/psb_dneigh.f90
Unescape Escape View File

@ -43,12 +43,13 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev)
type(psb_dspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n, info ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer, optional :: lev ! level of neighbours to find
integer, pointer :: tmpn(:)=>null()
integer :: level, dim, i, j, k, n1, err_act, nn, nidx
character(len=20) :: name
integer, allocatable :: tmpn(:)
integer :: level, dim, i, j, k, r, c, brow,&
& elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx
character(len=20) :: name, ch_err
name='psb_dneigh'
info = 0
@ -110,7 +111,7 @@ contains
type(psb_dspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer :: dim, i, iidx
@ -143,7 +144,7 @@ contains
type(psb_dspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer :: dim, i, iidx, ip, nza
@ -207,10 +208,10 @@ contains
implicit none
type(psb_dspmat_type), intent(in) :: a ! the sparse matrix
type(psb_dspmat_type), intent(in),target :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer :: dim, i, iidx, ip, nza
integer, pointer :: ia1(:), ia2(:), ia3(:),&
@ -296,7 +297,7 @@ contains
type(psb_dspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
select case(toupper(a%fida(1:3)))

11
src/serial/psb_dnumbmm.f90
Unescape Escape View File

@ -77,6 +77,9 @@ subroutine psb_dnumbmm(a,b,c)
& c%ia2,c%ia1,0,c%aspk,temp)
else
call inner_numbmm(a,b,c,temp,info)
if (info /= 0) then
write(0,*) 'Error ',info,' from inner numbmm'
end if
end if
deallocate(temp)
return
@ -118,11 +121,16 @@ contains
ajj = aval(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
info = 2
return
else
temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
endif
@ -131,6 +139,8 @@ contains
do j = c%ia2(i),c%ia2(i+1)-1
if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
info = 3
return
else
c%aspk(j) = temp(c%ia1(j))
temp(c%ia1(j)) = dzero
@ -143,5 +153,4 @@ contains
end subroutine inner_numbmm
end subroutine psb_dnumbmm

2
src/serial/psb_dspinfo.f90
Unescape Escape View File

@ -46,7 +46,7 @@ subroutine psb_dspinfo(ireq,a,ires,info,iaux)
use psb_string_mod
implicit none
type(psb_dspmat_type), intent(in) :: a
type(psb_dspmat_type), intent(in), target :: a
integer, intent(in) :: ireq
integer, intent(out) :: ires, info
integer, intent(in), optional :: iaux

8
src/serial/psb_dsymbmm.f90
Unescape Escape View File

@ -52,7 +52,7 @@ subroutine psb_dsymbmm(a,b,c,info)
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, pointer :: ic(:),jc(:)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
interface psb_sp_getrow
@ -111,6 +111,8 @@ subroutine psb_dsymbmm(a,b,c,info)
else
call inner_symbmm(a,b,c,itemp,info)
endif
call psb_realloc(size(c%ia1),c%aspk,info)
c%pl(1) = 0
c%pr(1) = 0
c%m=a%m
@ -173,11 +175,15 @@ contains
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
else
if(index(ibcl(k)).eq.0) then
index(ibcl(k))=istart

12
src/serial/psb_dtransp.f90
Unescape Escape View File

@ -45,7 +45,7 @@ subroutine psb_dtransp(a,b,c,fmt)
character(len=5) :: fmt_
integer ::c_, info, nz
integer, pointer :: itmp(:)=>null()
integer, allocatable :: itmp(:)
type(psb_dspmat_type) :: tmp
if (present(c)) then
@ -60,7 +60,7 @@ subroutine psb_dtransp(a,b,c,fmt)
endif
if (.true.) then
if (associated(b%aspk)) call psb_sp_free(b,info)
if (allocated(b%aspk)) call psb_sp_free(b,info)
b%fida = 'COO'
b%descra = 'GUN'
call psb_csdp(a,b,info)
@ -70,7 +70,7 @@ subroutine psb_dtransp(a,b,c,fmt)
return
end if
else
if (associated(b%aspk)) call psb_sp_free(b,info)
if (allocated(b%aspk)) call psb_sp_free(b,info)
call psb_sp_clone(a,b,info)
if (b%fida=='CSR') then
@ -86,9 +86,9 @@ subroutine psb_dtransp(a,b,c,fmt)
!!$ write(0,*) 'TRANSP CHECKS:',a%m,a%k,&
!!$ &minval(b%ia1(1:nz)),maxval(b%ia1(1:nz)),&
!!$ &minval(b%ia2(1:nz)),maxval(b%ia2(1:nz))
itmp => b%ia1
b%ia1 => b%ia2
b%ia2 => itmp
call psb_transfer(b%ia1,itmp,info)
call psb_transfer(b%ia2,b%ia1,info)
call psb_transfer(itmp,b%ia2,info)
b%m = a%k
b%k = a%m

10
src/serial/psb_zcsdp.f90
Unescape Escape View File

@ -195,11 +195,11 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (.not.associated(b%aspk).or.&
&.not.associated(b%ia1).or.&
&.not.associated(b%ia2).or.&
&.not.associated(b%pl).or.&
&.not.associated(b%pr)) then
if (.not.allocated(b%aspk).or.&
&.not.allocated(b%ia1).or.&
&.not.allocated(b%ia2).or.&
&.not.allocated(b%pl).or.&
&.not.allocated(b%pr)) then
call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info)
else if ((size(b%aspk) < aspk_size) .or.&
&(size(b%ia1) < ia1_size) .or.&

2
src/serial/psb_zcsrws.f90
Unescape Escape View File

@ -38,7 +38,7 @@ subroutine psb_zcsrws(rw,a,info,trans)
implicit none
type(psb_zspmat_type) :: a
complex(kind(1.d0)), pointer :: rw(:)
complex(kind(1.d0)), allocatable :: rw(:)
integer :: info
character, optional :: trans

7
src/serial/psb_zipcoo2csc.f90
Unescape Escape View File

@ -38,6 +38,7 @@ subroutine psb_zipcoo2csc(a,info,clshr)
use psb_serial_mod, only : psb_fixcoo
use psb_error_mod
use psb_string_mod
use psb_realloc_mod
implicit none
!....Parameters...
@ -45,7 +46,7 @@ subroutine psb_zipcoo2csc(a,info,clshr)
Integer, intent(out) :: info
logical, optional :: clshr
integer, pointer :: iaux(:), itemp(:)
integer, allocatable :: iaux(:), itemp(:)
!locals
logical :: clshr_
Integer :: nza, i,j,irw, idl,err_act,nc,icl
@ -76,8 +77,8 @@ subroutine psb_zipcoo2csc(a,info,clshr)
allocate(iaux(nc+1))
if(debug) write(0,*)'DIPCOO2CSC: out of fixcoo',nza,nc,size(a%ia2),size(iaux)
itemp => a%ia2
a%ia2 => iaux
call psb_transfer(a%ia2,itemp,info)
call psb_transfer(iaux,a%ia2,info)
!
! This routine can be used in two modes:

9
src/serial/psb_zipcoo2csr.f90
Unescape Escape View File

@ -35,6 +35,7 @@
subroutine psb_zipcoo2csr(a,info,rwshr)
use psb_spmat_type
use psb_const_mod
use psb_realloc_mod
use psb_serial_mod, only : psb_fixcoo
use psb_error_mod
use psb_string_mod
@ -45,7 +46,7 @@ subroutine psb_zipcoo2csr(a,info,rwshr)
Integer, intent(out) :: info
logical, optional :: rwshr
integer, pointer :: iaux(:), itemp(:)
integer, allocatable :: iaux(:), itemp(:)
!locals
logical :: rwshr_
Integer :: nza, nr, i,j,irw, idl,err_act
@ -81,9 +82,9 @@ subroutine psb_zipcoo2csr(a,info,rwshr)
if(debug) write(0,*)'DIPCOO2CSR: out of fixcoo',nza,nr,size(a%ia2),size(iaux)
itemp => a%ia1
a%ia1 => a%ia2
a%ia2 => iaux
call psb_transfer(a%ia1,itemp,info)
call psb_transfer(a%ia2,a%ia1,info)
call psb_transfer(iaux,a%ia2,info)
!
! This routine can be used in two modes:

17
src/serial/psb_zipcsr2coo.f90
Unescape Escape View File

@ -37,6 +37,7 @@ Subroutine psb_zipcsr2coo(a,info)
use psb_const_mod
use psb_error_mod
use psb_string_mod
use psb_realloc_mod
implicit none
!....Parameters...
@ -44,11 +45,11 @@ Subroutine psb_zipcsr2coo(a,info)
Integer, intent(out) :: info
!locals
Integer :: nza, nr
integer :: i,j,err_act
logical, parameter :: debug=.false.
integer, pointer :: iaux(:), itemp(:)
character(len=20) :: name, ch_err
Integer :: nza, nr
integer :: i,j,err_act
logical, parameter :: debug=.false.
integer, allocatable :: iaux(:), itemp(:)
character(len=20) :: name, ch_err
name='psb_zipcsr2coo'
info = 0
@ -68,9 +69,9 @@ Subroutine psb_zipcsr2coo(a,info)
return
end if
!!$ write(0,*) 'ipcsr2coo ',a%m
itemp => a%ia2
a%ia2 => a%ia1
a%ia1 => iaux
call psb_transfer(a%ia2,itemp,info)
call psb_transfer(a%ia1,a%ia2,info)
call psb_transfer(iaux,a%ia1,info)
do i=1, nr
do j=itemp(i),itemp(i+1)-1

14
src/serial/psb_zneigh.f90
Unescape Escape View File

@ -43,10 +43,10 @@ subroutine psb_zneigh(a,idx,neigh,n,info,lev)
type(psb_zspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n, info ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer, optional :: lev ! level of neighbours to find
integer, pointer :: tmpn(:)=>null()
integer, allocatable :: tmpn(:)
integer :: level, dim, i, j, k, r, c, brow,&
& elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx
character(len=20) :: name, ch_err
@ -111,7 +111,7 @@ contains
type(psb_zspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer :: dim, i, iidx
@ -144,7 +144,7 @@ contains
type(psb_zspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer :: dim, i, iidx, ip, nza
@ -208,10 +208,10 @@ contains
implicit none
type(psb_zspmat_type), intent(in) :: a ! the sparse matrix
type(psb_zspmat_type), intent(in),target :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
integer :: dim, i, iidx, ip, nza
integer, pointer :: ia1(:), ia2(:), ia3(:),&
@ -297,7 +297,7 @@ contains
type(psb_zspmat_type), intent(in) :: a ! the sparse matrix
integer, intent(in) :: idx ! the index whose neighbours we want to find
integer, intent(out) :: n ! the number of neighbours and the info
integer, pointer :: neigh(:) ! the neighbours
integer, allocatable :: neigh(:) ! the neighbours
select case(toupper(a%fida(1:3)))

2
src/serial/psb_zspinfo.f90
Unescape Escape View File

@ -46,7 +46,7 @@ subroutine psb_zspinfo(ireq,a,ires,info,iaux)
use psb_string_mod
implicit none
type(psb_zspmat_type), intent(in) :: a
type(psb_zspmat_type), intent(in), target :: a
integer, intent(in) :: ireq
integer, intent(out) :: ires, info
integer, intent(in), optional :: iaux

5
src/serial/psb_zsymbmm.f90
Unescape Escape View File

@ -52,7 +52,7 @@ subroutine psb_zsymbmm(a,b,c,info)
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, pointer :: ic(:),jc(:)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
@ -100,7 +100,7 @@ subroutine psb_zsymbmm(a,b,c,info)
endif
nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info)
!!$ write(0,*) 'SYMBMM90 ',size(c%pl),size(c%pr)
!
! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code.
@ -112,6 +112,7 @@ subroutine psb_zsymbmm(a,b,c,info)
else
call inner_symbmm(a,b,c,itemp,info)
endif
call psb_realloc(size(c%ia1),c%aspk,info)
c%pl(1) = 0
c%pr(1) = 0
c%m=a%m

24
src/serial/psb_ztransc.f90
Unescape Escape View File

@ -36,6 +36,7 @@ subroutine psb_ztransc(a,b,c,fmt)
use psb_spmat_type
use psb_tools_mod
use psb_string_mod
use psb_realloc_mod
use psb_serial_mod, only : psb_ipcoo2csr, psb_ipcsr2coo, psb_fixcoo, psb_csdp
implicit none
@ -44,8 +45,8 @@ subroutine psb_ztransc(a,b,c,fmt)
character(len=*), optional :: fmt
character(len=5) :: fmt_
integer ::c_, info, nz, i
integer, pointer :: itmp(:)=>null()
integer :: c_, info, nz,i
integer, allocatable :: itmp(:)
type(psb_zspmat_type) :: tmp
if (present(c)) then
@ -59,10 +60,17 @@ subroutine psb_ztransc(a,b,c,fmt)
fmt_='CSR'
endif
if (.true.) then
b%fida = 'COO'
if (allocated(b%aspk)) call psb_sp_free(b,info)
b%fida = 'COO'
b%descra = 'GUN'
call psb_csdp(a,b,info)
!!$ write(0,*) 'Check from CSDP',b%m,b%k,b%fida,b%descra,b%infoa(psb_nnz_)
if (info /= 0) then
write(0,*) 'transp: info from CSDP ',info
return
end if
else
if (associated(b%aspk)) call psb_sp_free(b,info)
if (allocated(b%aspk)) call psb_sp_free(b,info)
call psb_sp_clone(a,b,info)
if (b%fida=='CSR') then
@ -73,9 +81,10 @@ subroutine psb_ztransc(a,b,c,fmt)
write(0,*) 'Unimplemented case in TRANSC '
endif
endif
itmp => b%ia1
b%ia1 => b%ia2
b%ia2 => itmp
call psb_transfer(b%ia1,itmp,info)
call psb_transfer(b%ia2,b%ia1,info)
call psb_transfer(itmp,b%ia2,info)
b%m = a%k
b%k = a%m
@ -84,7 +93,6 @@ subroutine psb_ztransc(a,b,c,fmt)
b%aspk(i) = conjg(b%aspk(i))
end do
!!$ write(0,*) 'Calling IPCOO2CSR from transc90 ',b%m,b%k
if (fmt_=='CSR') then
call psb_ipcoo2csr(b,info)
b%fida='CSR'

30
src/serial/psb_ztransp.f90
Unescape Escape View File

@ -36,6 +36,7 @@ subroutine psb_ztransp(a,b,c,fmt)
use psb_spmat_type
use psb_tools_mod
use psb_string_mod
use psb_realloc_mod
use psb_serial_mod, only : psb_ipcoo2csr, psb_ipcsr2coo, psb_fixcoo, psb_csdp
implicit none
@ -44,8 +45,8 @@ subroutine psb_ztransp(a,b,c,fmt)
character(len=*), optional :: fmt
character(len=5) :: fmt_
integer ::c_, info, nz
integer, pointer :: itmp(:)=>null()
integer :: c_, info, nz
integer, allocatable :: itmp(:)
type(psb_zspmat_type) :: tmp
if (present(c)) then
@ -58,11 +59,19 @@ subroutine psb_ztransp(a,b,c,fmt)
else
fmt_='CSR'
endif
if (.true.) then
b%fida = 'COO'
if (allocated(b%aspk)) call psb_sp_free(b,info)
b%fida = 'COO'
b%descra = 'GUN'
call psb_csdp(a,b,info)
!!$ write(0,*) 'Check from CSDP',b%m,b%k,b%fida,b%descra,b%infoa(psb_nnz_)
if (info /= 0) then
write(0,*) 'transp: info from CSDP ',info
return
end if
else
if (associated(b%aspk)) call psb_sp_free(b,info)
if (allocated(b%aspk)) call psb_sp_free(b,info)
call psb_sp_clone(a,b,info)
if (b%fida=='CSR') then
@ -73,17 +82,14 @@ subroutine psb_ztransp(a,b,c,fmt)
write(0,*) 'Unimplemented case in TRANSP '
endif
endif
!!$ nz = b%infoa(nnz_)
!!$ write(0,*) 'TRANSP CHECKS:',a%m,a%k,&
!!$ &minval(b%ia1(1:nz)),maxval(b%ia1(1:nz)),&
!!$ &minval(b%ia2(1:nz)),maxval(b%ia2(1:nz))
itmp => b%ia1
b%ia1 => b%ia2
b%ia2 => itmp
call psb_transfer(b%ia1,itmp,info)
call psb_transfer(b%ia2,b%ia1,info)
call psb_transfer(itmp,b%ia2,info)
b%m = a%k
b%k = a%m
!!$ write(0,*) 'Calling IPCOO2CSR from transp90 ',b%m,b%k
if (fmt_=='CSR') then
call psb_ipcoo2csr(b,info)
b%fida='CSR'

10
src/tools/psb_cdall.f90
Unescape Escape View File

@ -59,7 +59,7 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info)
Integer :: counter,i,j,np,me,loc_row,err,loc_col,nprocs,&
& l_ov_ix,l_ov_el,idx, err_act, itmpov, k, ns
integer :: int_err(5),exch(2)
integer, pointer :: prc_v(:), temp_ovrlap(:), ov_idx(:),ov_el(:)
integer, allocatable :: prc_v(:), temp_ovrlap(:), ov_idx(:),ov_el(:)
logical, parameter :: debug=.false.
character(len=20) :: name, char_err
@ -166,7 +166,6 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info)
endif
desc_a%glob_to_loc(i) = -(np+prc_v(1)+1)
j=1
!!$ do while ((j <= nprocs).and.(prc_v(j) /= me))
do
if (j > nprocs) exit
if (prc_v(j) == me) exit
@ -259,8 +258,8 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info)
l_ov_ix = l_ov_ix + 1
ov_idx(l_ov_ix) = -1
desc_a%ovrlap_index => ov_idx
desc_a%ovrlap_elem => ov_el
call psb_transfer(ov_idx,desc_a%ovrlap_index,info)
if (info == 0) call psb_transfer(ov_el,desc_a%ovrlap_elem,info)
deallocate(prc_v,temp_ovrlap,stat=info)
if (info /= no_err) then
info=4000
@ -269,7 +268,6 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info)
Goto 9999
endif
! estimate local cols number
!!$ loc_col=int((psb_colrow_+1.d0)*loc_row)+1
loc_col=min(2*loc_row,m)
allocate(desc_a%loc_to_glob(loc_col),&
&desc_a%lprm(1),stat=info)
@ -287,9 +285,7 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info)
desc_a%loc_to_glob(k) = i
endif
enddo
nullify(desc_a%bnd_elem,desc_a%halo_index)
!!$ if (debug) write(*,*) 'PSB_CDALL: Last bits in desc_a', loc_row,k
! set fields in desc_a%MATRIX_DATA....
desc_a%matrix_data(psb_n_row_) = loc_row
desc_a%matrix_data(psb_n_col_) = loc_row

8
src/tools/psb_cdalv.f90
Unescape Escape View File

@ -59,7 +59,7 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
& loc_col,nprocs,n,itmpov, k,&
& l_ov_ix,l_ov_el,idx, flag_, err_act
integer :: int_err(5),exch(2)
Integer, Pointer :: temp_ovrlap(:), ov_idx(:),ov_el(:)
Integer, allocatable :: temp_ovrlap(:), ov_idx(:),ov_el(:)
logical, parameter :: debug=.false.
character(len=20) :: name
@ -225,8 +225,8 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
l_ov_ix = l_ov_ix + 1
ov_idx(l_ov_ix) = -1
desc_a%ovrlap_index => ov_idx
desc_a%ovrlap_elem => ov_el
call psb_transfer(ov_idx,desc_a%ovrlap_index,info)
if (info == 0) call psb_transfer(ov_el,desc_a%ovrlap_elem,info)
deallocate(temp_ovrlap,stat=info)
if (info /= 0) then
info=4000
@ -254,9 +254,7 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
desc_a%loc_to_glob(k) = i
endif
enddo
nullify(desc_a%bnd_elem,desc_a%halo_index)
!!$ if (debug) write(*,*) 'PSB_CDALL: Last bits in desc_a', loc_row,k
! set fields in desc_a%MATRIX_DATA....
desc_a%matrix_data(psb_n_row_) = loc_row
desc_a%matrix_data(psb_n_col_) = loc_row

11
src/tools/psb_cdasb.f90
Unescape Escape View File

@ -52,8 +52,8 @@ subroutine psb_cdasb(desc_a,info)
!....Locals....
integer :: int_err(5), itemp(2)
integer,pointer :: ovrlap_index(:),halo_index(:)
integer :: i,err,np,me,&
integer,allocatable :: ovrlap_index(:),halo_index(:)
integer :: i,err,np,me,&
& lovrlap,lhalo,nhalo,novrlap,max_size,max_halo,n_col,ldesc_halo,&
& ldesc_ovrlap, dectype, err_act
integer :: ictxt,n_row
@ -104,11 +104,8 @@ subroutine psb_cdasb(desc_a,info)
endif
call psb_realloc(desc_a%matrix_data(psb_n_col_),desc_a%loc_to_glob,info)
ovrlap_index => desc_a%ovrlap_index
nullify(desc_a%ovrlap_index)
halo_index => desc_a%halo_index
nullify(desc_a%halo_index)
call psb_transfer(desc_a%ovrlap_index,ovrlap_index,info)
call psb_transfer(desc_a%halo_index,halo_index,info)
call psi_cnv_dsc(halo_index,ovrlap_index,desc_a,info)
if (info /= 0) then

8
src/tools/psb_cdcpy.f90
Unescape Escape View File

@ -63,14 +63,16 @@ subroutine psb_cdcpy(desc_in, desc_out, info)
call psb_erractionsave(err_act)
name = 'psb_cdcpy'
ictxt=desc_in%matrix_data(psb_ctxt_)
! check on blacs grid
call psb_info(ictxt, me, np)
if (debug) write(0,*) me,'Entered CDCPY'
if (np == -1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
endif
call psb_nullify_desc(desc_out)

39
src/tools/psb_cdfree.f90
Unescape Escape View File

@ -56,7 +56,7 @@ subroutine psb_cdfree(desc_a,info)
name = 'psb_cdfree'
if (.not.associated(desc_a%matrix_data)) then
if (.not.allocated(desc_a%matrix_data)) then
info=295
call psb_errpush(info,name)
return
@ -73,7 +73,7 @@ subroutine psb_cdfree(desc_a,info)
endif
!...deallocate desc_a....
if(.not.associated(desc_a%loc_to_glob)) then
if(.not.allocated(desc_a%loc_to_glob)) then
info=295
call psb_errpush(info,name)
goto 9999
@ -87,7 +87,7 @@ subroutine psb_cdfree(desc_a,info)
goto 9999
end if
if (.not.associated(desc_a%glob_to_loc)) then
if (.not.allocated(desc_a%glob_to_loc)) then
info=295
call psb_errpush(info,name)
goto 9999
@ -101,7 +101,7 @@ subroutine psb_cdfree(desc_a,info)
goto 9999
end if
if (.not.associated(desc_a%halo_index)) then
if (.not.allocated(desc_a%halo_index)) then
info=295
call psb_errpush(info,name)
goto 9999
@ -115,14 +115,13 @@ subroutine psb_cdfree(desc_a,info)
goto 9999
end if
!!$ if (.not.associated(desc_a%bnd_elem)) then
if (.not.allocated(desc_a%bnd_elem)) then
!!$ info=296
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!deallocate halo_index field
if (associated(desc_a%bnd_elem)) then
else
!deallocate halo_index field
deallocate(desc_a%bnd_elem,stat=info)
if (info /= 0) then
info=2054
@ -130,7 +129,8 @@ subroutine psb_cdfree(desc_a,info)
goto 9999
end if
end if
if (.not.associated(desc_a%ovrlap_index)) then
if (.not.allocated(desc_a%ovrlap_index)) then
info=295
call psb_errpush(info,name)
goto 9999
@ -160,7 +160,7 @@ subroutine psb_cdfree(desc_a,info)
goto 9999
end if
if (associated(desc_a%idx_space)) then
if (allocated(desc_a%idx_space)) then
deallocate(desc_a%idx_space,stat=info)
if (info /= 0) then
info=2056
@ -169,25 +169,6 @@ subroutine psb_cdfree(desc_a,info)
end if
end if
!!$ if (associated(desc_a%halo_pt)) then
!!$ deallocate(desc_a%halo_pt,stat=info)
!!$ if (info /= 0) then
!!$ info=2056
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!!$ end if
!!$
!!$ if (associated(desc_a%ovrlap_pt)) then
!!$ deallocate(desc_a%ovrlap_pt,stat=info)
!!$ if (info /= 0) then
!!$ info=2056
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!!$ end if
call psb_nullify_desc(desc_a)
call psb_erractionrestore(err_act)

3
src/tools/psb_cdins.f90
Unescape Escape View File

@ -98,8 +98,7 @@ subroutine psb_cdins(nz,ia,ja,desc_a,info)
goto 9999
end if
if (.not.associated(desc_a%halo_index)) then
if (.not.allocated(desc_a%halo_index)) then
allocate(desc_a%halo_index(relocsz))
desc_a%halo_index(:) = -1
endif

26
src/tools/psb_cdtransfer.f90
Unescape Escape View File

@ -50,7 +50,7 @@ subroutine psb_cdtransfer(desc_in, desc_out, info)
!....parameters...
type(psb_desc_type), intent(inout) :: desc_in
type(psb_desc_type), intent(out) :: desc_out
type(psb_desc_type), intent(out) :: desc_out
integer, intent(out) :: info
!locals
@ -74,19 +74,19 @@ subroutine psb_cdtransfer(desc_in, desc_out, info)
goto 9999
endif
call psb_nullify_desc(desc_out)
desc_out%matrix_data => desc_in%matrix_data
desc_out%halo_index => desc_in%halo_index
desc_out%bnd_elem => desc_in%bnd_elem
desc_out%ovrlap_elem => desc_in%ovrlap_elem
desc_out%ovrlap_index => desc_in%ovrlap_index
desc_out%loc_to_glob => desc_in%loc_to_glob
desc_out%glob_to_loc => desc_in%glob_to_loc
desc_out%lprm => desc_in%lprm
desc_out%idx_space => desc_in%idx_space
!!$ call psb_nullify_desc(desc_out)
!!$
call psb_transfer( desc_in%matrix_data , desc_out%matrix_data , info)
call psb_transfer( desc_in%halo_index , desc_out%halo_index , info)
call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info)
call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info)
call psb_transfer( desc_in%ovrlap_index, desc_out%ovrlap_index , info)
call psb_transfer( desc_in%loc_to_glob , desc_out%loc_to_glob , info)
call psb_transfer( desc_in%glob_to_loc , desc_out%glob_to_loc , info)
call psb_transfer( desc_in%lprm , desc_out%lprm , info)
call psb_transfer( desc_in%idx_space , desc_out%idx_space , info)
call psb_nullify_desc(desc_in)
call psb_erractionrestore(err_act)
return

4
src/tools/psb_dallc.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_dalloc(x, desc_a, info, n)
implicit none
!....parameters...
real(kind(1.d0)), pointer :: x(:,:)
real(kind(1.d0)), allocatable, intent(out) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n
@ -197,7 +197,7 @@ subroutine psb_dallocv(x, desc_a,info,n)
implicit none
!....parameters...
real(kind(1.d0)), pointer :: x(:)
real(kind(1.d0)), allocatable, intent(out) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer :: info
integer, optional, intent(in) :: n

11
src/tools/psb_dasb.f90
Unescape Escape View File

@ -48,7 +48,7 @@ subroutine psb_dasb(x, desc_a, info)
implicit none
type(psb_desc_type), intent(in) :: desc_a
real(kind(1.d0)), pointer :: x(:,:)
real(kind(1.d0)), allocatable, intent(inout) :: x(:,:)
integer, intent(out) :: info
! local variables
@ -63,7 +63,7 @@ subroutine psb_dasb(x, desc_a, info)
name='psb_dasb'
call psb_erractionsave(err_act)
if ((.not.associated(desc_a%matrix_data))) then
if ((.not.allocated(desc_a%matrix_data))) then
info=3110
call psb_errpush(info,name)
goto 9999
@ -97,7 +97,8 @@ subroutine psb_dasb(x, desc_a, info)
i1sz = size(x,dim=1)
i2sz = size(x,dim=2)
if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol
if (i1sz.lt.ncol) then
if (i1sz < ncol) then
call psb_realloc(ncol,i2sz,x,info)
if (info /= 0) then
info=2025
@ -178,7 +179,7 @@ subroutine psb_dasbv(x, desc_a, info)
implicit none
type(psb_desc_type), intent(in) :: desc_a
real(kind(1.d0)), pointer :: x(:)
real(kind(1.d0)), allocatable, intent(inout) :: x(:)
integer, intent(out) :: info
! local variables
@ -213,7 +214,7 @@ subroutine psb_dasbv(x, desc_a, info)
if (debug) write(*,*) name,' sizes: ',nrow,ncol
i1sz = size(x)
if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol
if (i1sz.lt.ncol) then
if (i1sz < ncol) then
call psb_realloc(ncol,x,info)
if (info /= 0) then
info=2025

6
src/tools/psb_dcdovr.f90
Unescape Escape View File

@ -170,9 +170,8 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
index_dim = size(desc_a%halo_index)
elem_dim = size(desc_a%halo_index)
nullify(desc_ov%ovrlap_index,desc_ov%halo_index,desc_ov%ovrlap_elem)
allocate(desc_ov%ovrlap_elem(novr*(Max(elem_dim,1)+3)),&
& desc_ov%matrix_data(psb_mdata_size_),STAT=INFO)
call psb_realloc(psb_mdata_size_,desc_ov%matrix_data,info)
if (info==0) call psb_realloc(novr*(Max(elem_dim,1)+3),desc_ov%ovrlap_elem,info)
if (info /= 0) then
info=4000
call psb_errpush(info,name)
@ -199,7 +198,6 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
Write(0,*)'Start cdovrbld',me,lworks,lworkr
call psb_barrier(ictxt)
endif
!
! The real work goes on in here....
!

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

Write Preview
Loading…
Cancel
Save