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

Test merging of main/kinsol-stop

Browse Source
testmergemaint
sfilippone 1 week ago
parent 3a30c9495d
commit 1468c6ea2a
104 changed files with 2862 additions and 587 deletions
Show Stats Download Patch File Download Diff File

4
CMakeLists.txt
Unescape Escape View File

@ -428,7 +428,7 @@ if(METIS_FOUND)
include(CheckTypeSize)
message(STATUS "METIS PATH ${METIS_INCLUDES} and metis libraries ${METIS_LIBRARIES}")
# Make sure this path is correct
# set(METISINCFILE "metis.h") # Adjust this to your actual path
set(CMAKE_METIS_INCFILE "metis.h") # Adjust this to your actual path
# Specify the configuration file
# set(HEADER_TEMPLATE "${CMAKE_CURRENT_SOURCE_DIR}/util/psb_metis_int.h.in")
@ -509,7 +509,7 @@ if(METIS_FOUND)
# 1. Tell CMake where to find metis.h for the check
set(CMAKE_EXTRA_INCLUDE_FILES "${METIS_INCLUDES}/metis.h")
set(CMAKE_EXTRA_INCLUDE_FILES "${METIS_INCLUDES}/${CMAKE_METIS_INCFILE}")
# 2. Check the size of Metis's own type: real_t
# This replaces checking 'float' and 'double' separately

2
Makefile
Unescape Escape View File

@ -1,5 +1,5 @@
include Make.inc
all: dirs mods objs libd
@echo "====================================="
@echo "PSBLAS libraries Compilation Successful."

4
ReleaseNews
Unescape Escape View File

@ -1,4 +1,8 @@
WHAT'S NEW
Version 3.9.0-1
1. Fix licensing issues
2. Fix build and packaging
Version 3.9
1. PSBLAS3-EXT has been folded into the main library
2. Renamed GPU into CUDA.

3
base/CMakeLists.txt
Unescape Escape View File

@ -435,7 +435,7 @@ set(PSB_base_source_files
modules/comm/psi_s_comm_a_mod.f90
modules/comm/psb_i2_comm_a_mod.f90
modules/comm/psi_i2_comm_a_mod.f90
modules/comm/psi_i2_comm_v_mod.f90
modules/comm/psi_i2_comm_v_mod.f90
modules/comm/psi_m_comm_a_mod.f90
modules/comm/psi_l_comm_v_mod.f90
modules/comm/psb_comm_mod.f90
@ -460,6 +460,7 @@ set(PSB_base_source_files
modules/comm/psb_d_comm_mod.f90
modules/comm/psi_e_comm_a_mod.f90
modules/comm/psb_c_comm_a_mod.f90
modules/comm/psi_i2_comm_a_mod.f90
modules/comm/psb_linmap_mod.f90
modules/comm/psb_z_comm_a_mod.f90
modules/comm/psi_c_comm_a_mod.f90

2
base/modules/psi_i2_mod.F90
Unescape Escape View File

@ -31,7 +31,7 @@
!
module psi_i2_mod
use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_mpk_, psb_epk_, &
use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_mpk_, psb_epk_, &
& psb_lpk_, psb_i2pk_
use psi_m_comm_a_mod
use psi_e_comm_a_mod

2
base/modules/psi_i_mod.F90
Unescape Escape View File

@ -31,7 +31,7 @@
!
module psi_i_mod
use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_mpk_, psb_epk_, &
use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_mpk_, psb_epk_, &
& psb_lpk_, psb_i2pk_
use psi_m_comm_a_mod
use psi_e_comm_a_mod

2
base/modules/psi_l_mod.F90
Unescape Escape View File

@ -31,7 +31,7 @@
!
module psi_l_mod
use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_mpk_, psb_epk_, &
use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_mpk_, psb_epk_, &
& psb_lpk_, psb_i2pk_
use psi_m_comm_a_mod
use psi_e_comm_a_mod

10
base/modules/serial/psb_c_base_vect_mod.F90
Unescape Escape View File

@ -1288,21 +1288,19 @@ contains
end if
end function c_base_get_entry
subroutine c_base_set_entry(x, index, val)
subroutine c_base_set_entry(x, index, val)
implicit none
class(psb_c_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_spk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
x%v(index) = val
call x%set_host()
end if
end subroutine c_base_set_entry
!
! Overwrite with absolute value
!

2
base/modules/serial/psb_c_vect_mod.F90
Unescape Escape View File

@ -863,7 +863,7 @@ contains
if (allocated(x%v)) res = x%v%get_entry(index)
end function c_vect_get_entry
subroutine c_vect_set_entry(x,index,val)
subroutine c_vect_set_entry(x,index,val)
implicit none
class(psb_c_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index

10
base/modules/serial/psb_d_base_vect_mod.F90
Unescape Escape View File

@ -1295,21 +1295,19 @@ contains
end if
end function d_base_get_entry
subroutine d_base_set_entry(x, index, val)
subroutine d_base_set_entry(x, index, val)
implicit none
class(psb_d_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_dpk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
x%v(index) = val
call x%set_host()
end if
end subroutine d_base_set_entry
!
! Overwrite with absolute value
!

2
base/modules/serial/psb_d_vect_mod.F90
Unescape Escape View File

@ -870,7 +870,7 @@ contains
if (allocated(x%v)) res = x%v%get_entry(index)
end function d_vect_get_entry
subroutine d_vect_set_entry(x,index,val)
subroutine d_vect_set_entry(x,index,val)
implicit none
class(psb_d_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index

10
base/modules/serial/psb_s_base_vect_mod.F90
Unescape Escape View File

@ -1295,21 +1295,19 @@ contains
end if
end function s_base_get_entry
subroutine s_base_set_entry(x, index, val)
subroutine s_base_set_entry(x, index, val)
implicit none
class(psb_s_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_spk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
x%v(index) = val
call x%set_host()
end if
end subroutine s_base_set_entry
!
! Overwrite with absolute value
!

2
base/modules/serial/psb_s_vect_mod.F90
Unescape Escape View File

@ -870,7 +870,7 @@ contains
if (allocated(x%v)) res = x%v%get_entry(index)
end function s_vect_get_entry
subroutine s_vect_set_entry(x,index,val)
subroutine s_vect_set_entry(x,index,val)
implicit none
class(psb_s_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index

10
base/modules/serial/psb_z_base_vect_mod.F90
Unescape Escape View File

@ -1288,21 +1288,19 @@ contains
end if
end function z_base_get_entry
subroutine z_base_set_entry(x, index, val)
subroutine z_base_set_entry(x, index, val)
implicit none
class(psb_z_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_dpk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
x%v(index) = val
call x%set_host()
end if
end subroutine z_base_set_entry
!
! Overwrite with absolute value
!

2
base/modules/serial/psb_z_vect_mod.F90
Unescape Escape View File

@ -863,7 +863,7 @@ contains
if (allocated(x%v)) res = x%v%get_entry(index)
end function z_vect_get_entry
subroutine z_vect_set_entry(x,index,val)
subroutine z_vect_set_entry(x,index,val)
implicit none
class(psb_z_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index

11
base/modules/tools/psb_c_tools_mod.F90
Unescape Escape View File

@ -454,6 +454,17 @@ Module psb_c_tools_mod
end function
end interface
interface psb_setelem
subroutine psb_c_setelem(index,val,x,desc_a,info)
import
type(psb_c_vect_type), intent(inout) :: x
integer(psb_lpk_), intent(in) :: index
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
complex(psb_spk_) ::val
end subroutine psb_c_setelem
end interface
interface psb_remap
subroutine psb_c_remap(np_remap, desc_in, a_in, &
& ipd, isrc, nrsrc, naggr, desc_out, a_out, info)

11
base/modules/tools/psb_d_tools_mod.F90
Unescape Escape View File

@ -454,6 +454,17 @@ Module psb_d_tools_mod
end function
end interface
interface psb_setelem
subroutine psb_d_setelem(index,val,x,desc_a,info)
import
type(psb_d_vect_type), intent(inout) :: x
integer(psb_lpk_), intent(in) :: index
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_) ::val
end subroutine psb_d_setelem
end interface
interface psb_remap
subroutine psb_d_remap(np_remap, desc_in, a_in, &
& ipd, isrc, nrsrc, naggr, desc_out, a_out, info)

11
base/modules/tools/psb_s_tools_mod.F90
Unescape Escape View File

@ -454,6 +454,17 @@ Module psb_s_tools_mod
end function
end interface
interface psb_setelem
subroutine psb_s_setelem(index,val,x,desc_a,info)
import
type(psb_s_vect_type), intent(inout) :: x
integer(psb_lpk_), intent(in) :: index
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
real(psb_spk_) ::val
end subroutine psb_s_setelem
end interface
interface psb_remap
subroutine psb_s_remap(np_remap, desc_in, a_in, &
& ipd, isrc, nrsrc, naggr, desc_out, a_out, info)

11
base/modules/tools/psb_z_tools_mod.F90
Unescape Escape View File

@ -454,6 +454,17 @@ Module psb_z_tools_mod
end function
end interface
interface psb_setelem
subroutine psb_z_setelem(index,val,x,desc_a,info)
import
type(psb_z_vect_type), intent(inout) :: x
integer(psb_lpk_), intent(in) :: index
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
complex(psb_dpk_) ::val
end subroutine psb_z_setelem
end interface
interface psb_remap
subroutine psb_z_remap(np_remap, desc_in, a_in, &
& ipd, isrc, nrsrc, naggr, desc_out, a_out, info)

3
base/tools/Makefile
Unescape Escape View File

@ -28,7 +28,8 @@ FOBJS = psb_cdall.o psb_cdals.o psb_cdalv.o psb_cd_inloc.o psb_cdins.o psb_cdprt
psb_s_map.o psb_d_map.o psb_c_map.o psb_z_map.o \
psb_s_par_csr_spspmm.o psb_d_par_csr_spspmm.o psb_c_par_csr_spspmm.o psb_z_par_csr_spspmm.o \
psb_s_glob_transpose.o psb_d_glob_transpose.o psb_c_glob_transpose.o psb_z_glob_transpose.o \
psb_cgetelem.o psb_dgetelem.o psb_sgetelem.o psb_zgetelem.o
psb_cgetelem.o psb_dgetelem.o psb_sgetelem.o psb_zgetelem.o \
psb_csetelem.o psb_dsetelem.o psb_ssetelem.o psb_zsetelem.o
MPFOBJS = psb_icdasb.o psb_ssphalo.o psb_dsphalo.o psb_csphalo.o psb_zsphalo.o \
psb_dcdbldext.o psb_zcdbldext.o psb_scdbldext.o psb_ccdbldext.o \

112
base/tools/psb_csetelem.f90
Unescape Escape View File

@ -0,0 +1,112 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
! Function: psb_c_setelem
! Set entries into a dense vector. Note: the row indices in index
! are assumed to be in global numbering and are converted on the fly.
! Row indices not belonging to the current process have to be in the halo,
! othewise failure is ensured.
!
! Arguments:
! x - type(psb_c_vect_type) The source vector
! desc_a - type(psb_desc_type). The communication descriptor.
! index - integer. Row index of x of the value to extract
! iam - integer. Index of the process requesting the value
! info - integer. return code
subroutine psb_c_setelem(index,val,x,desc_a,info)
use psb_base_mod, psb_protect_name => psb_c_setelem
use psi_mod
implicit none
type(psb_c_vect_type), intent(inout) :: x
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_lpk_), intent(in) :: index
integer(psb_ipk_), intent(out) :: info
complex(psb_spk_) :: val
!locals
integer(psb_ipk_) :: localindex(1)
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me, err_act
integer(psb_lpk_) :: gindex(1)
integer(psb_lpk_), allocatable :: myidx(:),mylocal(:)
character(len=20) :: name
logical, parameter :: debug = .false.
info = 0
gindex(1) = index
if (psb_errstatus_fatal()) return
info=psb_success_
call psb_erractionsave(err_act)
name = 'psb_c_getelem'
if (.not.desc_a%is_ok()) then
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
goto 9999
end if
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
call desc_a%indxmap%g2l(gindex,localindex,info,owned=.false.)
if(debug.and.(localindex(1) < 1)) then
write(*,*)"Process ",me," owns ",desc_a%get_local_rows()," rows"," Global index is ",gindex,"Local index is ",localindex
myidx = desc_a%get_global_indices(owned=.false.)
mylocal = desc_a%get_global_indices(owned=.true.)
write(*,*)"My (local+halo) indexes are: ",myidx
write(*,*)"My (local) indexes are: ",mylocal
end if
if ( localindex(1) < 1) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err="Index not in the HALO")
goto 9999
else
call x%set_entry(localindex(1),val)
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
end subroutine psb_c_setelem

112
base/tools/psb_dsetelem.f90
Unescape Escape View File

@ -0,0 +1,112 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
! Function: psb_d_setelem
! Set entries into a dense vector. Note: the row indices in index
! are assumed to be in global numbering and are converted on the fly.
! Row indices not belonging to the current process have to be in the halo,
! othewise failure is ensured.
!
! Arguments:
! x - type(psb_d_vect_type) The source vector
! desc_a - type(psb_desc_type). The communication descriptor.
! index - integer. Row index of x of the value to extract
! iam - integer. Index of the process requesting the value
! info - integer. return code
subroutine psb_d_setelem(index,val,x,desc_a,info)
use psb_base_mod, psb_protect_name => psb_d_setelem
use psi_mod
implicit none
type(psb_d_vect_type), intent(inout) :: x
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_lpk_), intent(in) :: index
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_) :: val
!locals
integer(psb_ipk_) :: localindex(1)
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me, err_act
integer(psb_lpk_) :: gindex(1)
integer(psb_lpk_), allocatable :: myidx(:),mylocal(:)
character(len=20) :: name
logical, parameter :: debug = .false.
info = 0
gindex(1) = index
if (psb_errstatus_fatal()) return
info=psb_success_
call psb_erractionsave(err_act)
name = 'psb_d_getelem'
if (.not.desc_a%is_ok()) then
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
goto 9999
end if
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
call desc_a%indxmap%g2l(gindex,localindex,info,owned=.false.)
if(debug.and.(localindex(1) < 1)) then
write(*,*)"Process ",me," owns ",desc_a%get_local_rows()," rows"," Global index is ",gindex,"Local index is ",localindex
myidx = desc_a%get_global_indices(owned=.false.)
mylocal = desc_a%get_global_indices(owned=.true.)
write(*,*)"My (local+halo) indexes are: ",myidx
write(*,*)"My (local) indexes are: ",mylocal
end if
if ( localindex(1) < 1) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err="Index not in the HALO")
goto 9999
else
call x%set_entry(localindex(1),val)
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
end subroutine psb_d_setelem

112
base/tools/psb_ssetelem.f90
Unescape Escape View File

@ -0,0 +1,112 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
! Function: psb_s_setelem
! Set entries into a dense vector. Note: the row indices in index
! are assumed to be in global numbering and are converted on the fly.
! Row indices not belonging to the current process have to be in the halo,
! othewise failure is ensured.
!
! Arguments:
! x - type(psb_s_vect_type) The source vector
! desc_a - type(psb_desc_type). The communication descriptor.
! index - integer. Row index of x of the value to extract
! iam - integer. Index of the process requesting the value
! info - integer. return code
subroutine psb_s_setelem(index,val,x,desc_a,info)
use psb_base_mod, psb_protect_name => psb_s_setelem
use psi_mod
implicit none
type(psb_s_vect_type), intent(inout) :: x
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_lpk_), intent(in) :: index
integer(psb_ipk_), intent(out) :: info
real(psb_spk_) :: val
!locals
integer(psb_ipk_) :: localindex(1)
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me, err_act
integer(psb_lpk_) :: gindex(1)
integer(psb_lpk_), allocatable :: myidx(:),mylocal(:)
character(len=20) :: name
logical, parameter :: debug = .false.
info = 0
gindex(1) = index
if (psb_errstatus_fatal()) return
info=psb_success_
call psb_erractionsave(err_act)
name = 'psb_s_getelem'
if (.not.desc_a%is_ok()) then
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
goto 9999
end if
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
call desc_a%indxmap%g2l(gindex,localindex,info,owned=.false.)
if(debug.and.(localindex(1) < 1)) then
write(*,*)"Process ",me," owns ",desc_a%get_local_rows()," rows"," Global index is ",gindex,"Local index is ",localindex
myidx = desc_a%get_global_indices(owned=.false.)
mylocal = desc_a%get_global_indices(owned=.true.)
write(*,*)"My (local+halo) indexes are: ",myidx
write(*,*)"My (local) indexes are: ",mylocal
end if
if ( localindex(1) < 1) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err="Index not in the HALO")
goto 9999
else
call x%set_entry(localindex(1),val)
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
end subroutine psb_s_setelem

112
base/tools/psb_zsetelem.f90
Unescape Escape View File

@ -0,0 +1,112 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
! Function: psb_z_setelem
! Set entries into a dense vector. Note: the row indices in index
! are assumed to be in global numbering and are converted on the fly.
! Row indices not belonging to the current process have to be in the halo,
! othewise failure is ensured.
!
! Arguments:
! x - type(psb_z_vect_type) The source vector
! desc_a - type(psb_desc_type). The communication descriptor.
! index - integer. Row index of x of the value to extract
! iam - integer. Index of the process requesting the value
! info - integer. return code
subroutine psb_z_setelem(index,val,x,desc_a,info)
use psb_base_mod, psb_protect_name => psb_z_setelem
use psi_mod
implicit none
type(psb_z_vect_type), intent(inout) :: x
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_lpk_), intent(in) :: index
integer(psb_ipk_), intent(out) :: info
complex(psb_dpk_) :: val
!locals
integer(psb_ipk_) :: localindex(1)
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me, err_act
integer(psb_lpk_) :: gindex(1)
integer(psb_lpk_), allocatable :: myidx(:),mylocal(:)
character(len=20) :: name
logical, parameter :: debug = .false.
info = 0
gindex(1) = index
if (psb_errstatus_fatal()) return
info=psb_success_
call psb_erractionsave(err_act)
name = 'psb_z_getelem'
if (.not.desc_a%is_ok()) then
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
goto 9999
end if
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
call desc_a%indxmap%g2l(gindex,localindex,info,owned=.false.)
if(debug.and.(localindex(1) < 1)) then
write(*,*)"Process ",me," owns ",desc_a%get_local_rows()," rows"," Global index is ",gindex,"Local index is ",localindex
myidx = desc_a%get_global_indices(owned=.false.)
mylocal = desc_a%get_global_indices(owned=.true.)
write(*,*)"My (local+halo) indexes are: ",myidx
write(*,*)"My (local) indexes are: ",mylocal
end if
if ( localindex(1) < 1) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err="Index not in the HALO")
goto 9999
else
call x%set_entry(localindex(1),val)
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
end subroutine psb_z_setelem

3
cbind/base/psb_c_base.c
Unescape Escape View File

@ -58,4 +58,7 @@ char *psb_c_pop_errmsg()
return(tmp);
}
void psb_c_print_pointer(void *p){
fprintf(stderr,"psb_c_print_pointer %p\n",p);
}
// Convertire il comunicatore fortran in comunicatore c

5
cbind/base/psb_c_base.h
Unescape Escape View File

@ -30,6 +30,10 @@ extern "C" {
psb_i_t *ctxt;
} psb_c_ctxt;
typedef struct PSB_C_OBJTYPE {
void *item;
} psb_c_objtype;
void psb_c_check_error(psb_c_ctxt cctxt);
@ -42,6 +46,7 @@ extern "C" {
void psb_c_seterraction_print();
void psb_c_seterraction_abort();
void psb_c_print_pointer(void *p);
/* Environment routines */
void psb_c_init(psb_c_ctxt *cctxt);
void psb_c_init_from_fint(psb_c_ctxt *cctxt, psb_i_t f_comm);

4
cbind/base/psb_c_cbase.h
Unescape Escape View File

@ -40,6 +40,8 @@ psb_i_t psb_c_cgefree(psb_c_cvector *xh, psb_c_descriptor *cdh);
psb_i_t psb_c_cgereinit(psb_c_cvector *xh, psb_c_descriptor *cdh, bool clear);
psb_c_t psb_c_cgetelem(psb_c_cvector *xh,psb_l_t index,psb_c_descriptor *cd);
psb_c_t psb_c_cmatgetelem(psb_c_cspmat *ah,psb_l_t rowindex,psb_l_t colindex,psb_c_descriptor *cdh);
psb_i_t psb_c_csetelem(psb_l_t index, psb_c_t val,
psb_c_cvector *xh, psb_c_descriptor *cd);
/* sparse matrices*/
psb_c_cspmat* psb_c_new_cspmat();
@ -65,6 +67,8 @@ psb_i_t psb_c_cspasb_opt(psb_c_cspmat *mh, psb_c_descriptor *cdh,
psb_i_t psb_c_csprn(psb_c_cspmat *mh, psb_c_descriptor *cdh, _Bool clear);
psb_i_t psb_c_cmat_name_print(psb_c_cspmat *mh, char *name);
psb_i_t psb_c_cvect_set_scal(psb_c_cvector *xh, psb_c_t val);
psb_i_t psb_c_cvect_set_scal_bound(psb_c_cvector *xh, psb_c_t val,
psb_i_t ifirst, psb_i_t ilast);
psb_i_t psb_c_cvect_set_vect(psb_c_cvector *xh, psb_c_t *val, psb_i_t n);
psb_c_t psb_c_cvect_get_entry(psb_c_cvector *xh, psb_i_t index);
psb_i_t psb_c_cvect_set_entry(psb_c_cvector *xh, psb_i_t index, psb_c_t val);

5
cbind/base/psb_c_dbase.h
Unescape Escape View File

@ -23,7 +23,6 @@ psb_i_t psb_c_dvect_f_get_cpy(psb_d_t *v, psb_c_dvector *xh);
psb_i_t psb_c_dvect_zero(psb_c_dvector *xh);
psb_d_t *psb_c_dvect_f_get_pnt( psb_c_dvector *xh);
psb_i_t psb_c_dvect_clone(psb_c_dvector *xh,psb_c_dvector *yh);
psb_i_t psb_c_dgeall(psb_c_dvector *xh, psb_c_descriptor *cdh);
psb_i_t psb_c_dgeall_remote(psb_c_dvector *xh, psb_c_descriptor *cdh);
psb_i_t psb_c_dgeall_remote_options(psb_c_dvector *xh, psb_c_descriptor *cdh,
@ -40,7 +39,9 @@ psb_i_t psb_c_dgefree(psb_c_dvector *xh, psb_c_descriptor *cdh);
psb_i_t psb_c_dgereinit(psb_c_dvector *xh, psb_c_descriptor *cdh, bool clear);
psb_d_t psb_c_dgetelem(psb_c_dvector *xh,psb_l_t index,psb_c_descriptor *cd);
psb_d_t psb_c_dmatgetelem(psb_c_dspmat *ah,psb_l_t rowindex,psb_l_t colindex,psb_c_descriptor *cdh);
psb_d_t psb_c_dgetelem(psb_c_dvector *xh, psb_l_t index,psb_c_descriptor *cd);
psb_i_t psb_c_dsetelem(psb_l_t index, psb_d_t val,
psb_c_dvector *xh, psb_c_descriptor *cd);
/* sparse matrices*/
psb_c_dspmat* psb_c_new_dspmat();

6
cbind/base/psb_c_sbase.h
Unescape Escape View File

@ -40,7 +40,8 @@ psb_i_t psb_c_sgefree(psb_c_svector *xh, psb_c_descriptor *cdh);
psb_i_t psb_c_sgereinit(psb_c_svector *xh, psb_c_descriptor *cdh, bool clear);
psb_s_t psb_c_sgetelem(psb_c_svector *xh,psb_l_t index,psb_c_descriptor *cd);
psb_s_t psb_c_smatgetelem(psb_c_sspmat *ah,psb_l_t rowindex,psb_l_t colindex,psb_c_descriptor *cdh);
psb_i_t psb_c_ssetelem(psb_l_t index, psb_s_t val,
psb_c_svector *xh, psb_c_descriptor *cd);
/* sparse matrices*/
psb_c_sspmat* psb_c_new_sspmat();
@ -69,7 +70,8 @@ psb_i_t psb_c_svect_set_scal(psb_c_svector *xh, psb_s_t val);
psb_i_t psb_c_svect_set_vect(psb_c_svector *xh, psb_s_t *val, psb_i_t n);
psb_s_t psb_c_svect_get_entry(psb_c_svector *xh, psb_i_t index);
psb_i_t psb_c_svect_set_entry(psb_c_svector *xh, psb_i_t index, psb_s_t val);
psb_i_t psb_c_svect_set_scal_bound(psb_c_svector *xh, psb_s_t val,
psb_i_t ifirst, psb_i_t ilast);
/* psblas computational routines */
psb_s_t psb_c_sgedot(psb_c_svector *xh, psb_c_svector *yh, psb_c_descriptor *cdh);

259
cbind/base/psb_c_serial_cbind_mod.F90
Unescape Escape View File

@ -176,6 +176,30 @@ contains
end function psb_c_cvect_set_scal
function psb_c_cvect_set_scal_bound(x,val,ifirst,ilast) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_cvector) :: x
type(psb_c_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: ifirst, ilast
complex(c_float_complex) :: val
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
call xp%set(val,first=ifirst,last=ilast)
info = 0
end function psb_c_cvect_set_scal_bound
function psb_c_cvect_set_vect(x,val,n) bind(c) result(info)
use psb_base_mod
implicit none
@ -203,7 +227,6 @@ contains
function psb_c_cvect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_cvector) :: x
type(psb_c_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
@ -222,20 +245,20 @@ contains
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_cvect_set_entry
function psb_c_cvect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_cvector) :: x
type(psb_c_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
complex(c_float_complex) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
@ -270,4 +293,232 @@ contains
end function psb_c_cvect_clone
function psb_c_cnnz(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = 0
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = psb_nnz(ap,descp,info)
end function psb_c_cnnz
function psb_c_cis_matupd(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_upd()
end function
function psb_c_cis_matasb(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_asb()
end function
function psb_c_cis_matbld(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_bld()
end function
function psb_c_cset_matupd(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_upd()
res = psb_success_
end function
function psb_c_cset_matasb(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
res = -1;
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_asb()
res = psb_success_
end function
function psb_c_cset_matbld(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_cspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_bld()
res = psb_success_
end function
function psb_c_ccopy_mat(ah,bh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_cspmat) :: ah,bh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap,bp
integer(psb_c_ipk_) :: info
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
if (c_associated(bh%item)) then
call c_f_pointer(bh%item,bp)
else
return
end if
call ap%clone(bp,info)
res = info
end function
end module psb_c_serial_cbind_mod

36
cbind/base/psb_c_tools_cbind_mod.F90
Unescape Escape View File

@ -688,7 +688,6 @@ contains
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
res = psb_getelem(ap,rowindex,colindex,descp,info)
@ -696,8 +695,43 @@ contains
res = psb_getelem(ap,rowindex+(1-ixb),colindex+(1-ixb),descp,info)
end if
res=info
return
end function psb_c_cmatgetelem
function psb_c_csetelem(index,val,xh,cdh) bind(c) result(res)
implicit none
type(psb_c_cvector) :: xh
integer(psb_c_lpk_), value :: index
type(psb_c_descriptor) :: cdh
complex(c_float_complex), value :: val
integer(psb_c_ipk_) :: res
type(psb_c_vect_type), pointer :: xp
type(psb_desc_type), pointer :: descp
integer(psb_c_ipk_) :: info, ixb
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(xh%item)) then
call c_f_pointer(xh%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_setelem(index,val,xp,descp,info)
else
call psb_setelem(index+(1-ixb),val,xp,descp,info)
end if
res = info
return
end function psb_c_csetelem
end module psb_c_tools_cbind_mod

5
cbind/base/psb_c_zbase.h
Unescape Escape View File

@ -40,7 +40,8 @@ psb_i_t psb_c_zgefree(psb_c_zvector *xh, psb_c_descriptor *cdh);
psb_i_t psb_c_zgereinit(psb_c_zvector *xh, psb_c_descriptor *cdh, bool clear);
psb_z_t psb_c_zgetelem(psb_c_zvector *xh,psb_l_t index,psb_c_descriptor *cd);
psb_z_t psb_c_zmatgetelem(psb_c_zspmat *ah,psb_l_t rowindex,psb_l_t colindex,psb_c_descriptor *cdh);
psb_i_t psb_c_zsetelem(psb_l_t index, psb_z_t val,
psb_c_zvector *xh, psb_c_descriptor *cd);
/* sparse matrices*/
psb_c_zspmat* psb_c_new_zspmat();
@ -68,6 +69,8 @@ psb_i_t psb_c_zsprn(psb_c_zspmat *mh, psb_c_descriptor *cdh, _Bool clear);
psb_i_t psb_c_zmat_name_print(psb_c_zspmat *mh, char *name);
psb_i_t psb_c_zvect_set_scal(psb_c_zvector *xh, psb_z_t val);
psb_i_t psb_c_zvect_set_vect(psb_c_zvector *xh, psb_z_t *val, psb_i_t n);
psb_i_t psb_c_zvect_set_scal_bound(psb_c_zvector *xh, psb_z_t val,
psb_i_t ifirst, psb_i_t ilast);
psb_z_t psb_c_zvect_get_entry(psb_c_zvector *xh, psb_i_t index);
psb_i_t psb_c_zvect_set_entry(psb_c_zvector *xh, psb_i_t index, psb_z_t val);

259
cbind/base/psb_d_serial_cbind_mod.F90
Unescape Escape View File

@ -176,6 +176,30 @@ contains
end function psb_c_dvect_set_scal
function psb_c_dvect_set_scal_bound(x,val,ifirst,ilast) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_dvector) :: x
type(psb_d_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: ifirst, ilast
real(c_double) :: val
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
call xp%set(val,first=ifirst,last=ilast)
info = 0
end function psb_c_dvect_set_scal_bound
function psb_c_dvect_set_vect(x,val,n) bind(c) result(info)
use psb_base_mod
implicit none
@ -203,7 +227,6 @@ contains
function psb_c_dvect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_dvector) :: x
type(psb_d_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
@ -222,20 +245,20 @@ contains
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_dvect_set_entry
function psb_c_dvect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_dvector) :: x
type(psb_d_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
real(c_double) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
@ -270,4 +293,232 @@ contains
end function psb_c_dvect_clone
function psb_c_dnnz(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = 0
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = psb_nnz(ap,descp,info)
end function psb_c_dnnz
function psb_c_dis_matupd(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_upd()
end function
function psb_c_dis_matasb(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_asb()
end function
function psb_c_dis_matbld(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_bld()
end function
function psb_c_dset_matupd(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_upd()
res = psb_success_
end function
function psb_c_dset_matasb(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
res = -1;
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_asb()
res = psb_success_
end function
function psb_c_dset_matbld(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_dspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_bld()
res = psb_success_
end function
function psb_c_dcopy_mat(ah,bh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_dspmat) :: ah,bh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap,bp
integer(psb_c_ipk_) :: info
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
if (c_associated(bh%item)) then
call c_f_pointer(bh%item,bp)
else
return
end if
call ap%clone(bp,info)
res = info
end function
end module psb_d_serial_cbind_mod

36
cbind/base/psb_d_tools_cbind_mod.F90
Unescape Escape View File

@ -698,7 +698,6 @@ contains
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
res = psb_getelem(ap,rowindex,colindex,descp,info)
@ -706,8 +705,43 @@ contains
res = psb_getelem(ap,rowindex+(1-ixb),colindex+(1-ixb),descp,info)
end if
res=info
return
end function psb_c_dmatgetelem
function psb_c_dsetelem(index,val,xh,cdh) bind(c) result(res)
implicit none
type(psb_c_dvector) :: xh
integer(psb_c_lpk_), value :: index
type(psb_c_descriptor) :: cdh
real(c_double), value :: val
integer(psb_c_ipk_) :: res
type(psb_d_vect_type), pointer :: xp
type(psb_desc_type), pointer :: descp
integer(psb_c_ipk_) :: info, ixb
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(xh%item)) then
call c_f_pointer(xh%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_setelem(index,val,xp,descp,info)
else
call psb_setelem(index+(1-ixb),val,xp,descp,info)
end if
res = info
return
end function psb_c_dsetelem
end module psb_d_tools_cbind_mod

11
cbind/base/psb_objhandle_mod.F90
Unescape Escape View File

@ -42,4 +42,15 @@ module psb_objhandle_mod
type(c_ptr) :: item = c_null_ptr
end type psb_c_zspmat
interface
subroutine psb_c_print_pointer(p) bind(c,name='psb_c_print_pointer')
use iso_c_binding
type(c_ptr), value :: p
end subroutine psb_c_print_pointer
end interface
contains
function psb_c_get_new_object() result(res)
type(psb_c_object_type) :: res
res%item = c_null_ptr
end function psb_c_get_new_object
end module psb_objhandle_mod

259
cbind/base/psb_s_serial_cbind_mod.F90
Unescape Escape View File

@ -176,6 +176,30 @@ contains
end function psb_c_svect_set_scal
function psb_c_svect_set_scal_bound(x,val,ifirst,ilast) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_svector) :: x
type(psb_s_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: ifirst, ilast
real(c_float) :: val
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
call xp%set(val,first=ifirst,last=ilast)
info = 0
end function psb_c_svect_set_scal_bound
function psb_c_svect_set_vect(x,val,n) bind(c) result(info)
use psb_base_mod
implicit none
@ -203,7 +227,6 @@ contains
function psb_c_svect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_svector) :: x
type(psb_s_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
@ -222,20 +245,20 @@ contains
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_svect_set_entry
function psb_c_svect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_svector) :: x
type(psb_s_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
real(c_float) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
@ -270,4 +293,232 @@ contains
end function psb_c_svect_clone
function psb_c_snnz(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = 0
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = psb_nnz(ap,descp,info)
end function psb_c_snnz
function psb_c_sis_matupd(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_upd()
end function
function psb_c_sis_matasb(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_asb()
end function
function psb_c_sis_matbld(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_bld()
end function
function psb_c_sset_matupd(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_upd()
res = psb_success_
end function
function psb_c_sset_matasb(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
res = -1;
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_asb()
res = psb_success_
end function
function psb_c_sset_matbld(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_sspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_bld()
res = psb_success_
end function
function psb_c_scopy_mat(ah,bh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_sspmat) :: ah,bh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap,bp
integer(psb_c_ipk_) :: info
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
if (c_associated(bh%item)) then
call c_f_pointer(bh%item,bp)
else
return
end if
call ap%clone(bp,info)
res = info
end function
end module psb_s_serial_cbind_mod

36
cbind/base/psb_s_tools_cbind_mod.F90
Unescape Escape View File

@ -698,7 +698,6 @@ contains
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
res = psb_getelem(ap,rowindex,colindex,descp,info)
@ -706,8 +705,43 @@ contains
res = psb_getelem(ap,rowindex+(1-ixb),colindex+(1-ixb),descp,info)
end if
res=info
return
end function psb_c_smatgetelem
function psb_c_ssetelem(index,val,xh,cdh) bind(c) result(res)
implicit none
type(psb_c_svector) :: xh
integer(psb_c_lpk_), value :: index
type(psb_c_descriptor) :: cdh
real(c_float), value :: val
integer(psb_c_ipk_) :: res
type(psb_s_vect_type), pointer :: xp
type(psb_desc_type), pointer :: descp
integer(psb_c_ipk_) :: info, ixb
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(xh%item)) then
call c_f_pointer(xh%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_setelem(index,val,xp,descp,info)
else
call psb_setelem(index+(1-ixb),val,xp,descp,info)
end if
res = info
return
end function psb_c_ssetelem
end module psb_s_tools_cbind_mod

259
cbind/base/psb_z_serial_cbind_mod.F90
Unescape Escape View File

@ -176,6 +176,30 @@ contains
end function psb_c_zvect_set_scal
function psb_c_zvect_set_scal_bound(x,val,ifirst,ilast) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_zvector) :: x
type(psb_z_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: ifirst, ilast
complex(c_double_complex) :: val
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
call xp%set(val,first=ifirst,last=ilast)
info = 0
end function psb_c_zvect_set_scal_bound
function psb_c_zvect_set_vect(x,val,n) bind(c) result(info)
use psb_base_mod
implicit none
@ -203,7 +227,6 @@ contains
function psb_c_zvect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_zvector) :: x
type(psb_z_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
@ -222,20 +245,20 @@ contains
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_zvect_set_entry
function psb_c_zvect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_zvector) :: x
type(psb_z_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
complex(c_double_complex) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
@ -270,4 +293,232 @@ contains
end function psb_c_zvect_clone
function psb_c_znnz(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = 0
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = psb_nnz(ap,descp,info)
end function psb_c_znnz
function psb_c_zis_matupd(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_upd()
end function
function psb_c_zis_matasb(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_asb()
end function
function psb_c_zis_matbld(ah,cdh) bind(c) result(res)
implicit none
logical(c_bool) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
integer(psb_c_ipk_) :: info
res = .false.
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
res = ap%is_bld()
end function
function psb_c_zset_matupd(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_upd()
res = psb_success_
end function
function psb_c_zset_matasb(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
res = -1;
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_asb()
res = psb_success_
end function
function psb_c_zset_matbld(ah,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_zspmat) :: ah
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
call ap%set_bld()
res = psb_success_
end function
function psb_c_zcopy_mat(ah,bh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
type(psb_c_zspmat) :: ah,bh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap,bp
integer(psb_c_ipk_) :: info
res = -1
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(ah%item)) then
call c_f_pointer(ah%item,ap)
else
return
end if
if (c_associated(bh%item)) then
call c_f_pointer(bh%item,bp)
else
return
end if
call ap%clone(bp,info)
res = info
end function
end module psb_z_serial_cbind_mod

36
cbind/base/psb_z_tools_cbind_mod.F90
Unescape Escape View File

@ -688,7 +688,6 @@ contains
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
res = psb_getelem(ap,rowindex,colindex,descp,info)
@ -696,8 +695,43 @@ contains
res = psb_getelem(ap,rowindex+(1-ixb),colindex+(1-ixb),descp,info)
end if
res=info
return
end function psb_c_zmatgetelem
function psb_c_zsetelem(index,val,xh,cdh) bind(c) result(res)
implicit none
type(psb_c_zvector) :: xh
integer(psb_c_lpk_), value :: index
type(psb_c_descriptor) :: cdh
complex(c_double_complex), value :: val
integer(psb_c_ipk_) :: res
type(psb_z_vect_type), pointer :: xp
type(psb_desc_type), pointer :: descp
integer(psb_c_ipk_) :: info, ixb
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
if (c_associated(xh%item)) then
call c_f_pointer(xh%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_setelem(index,val,xp,descp,info)
else
call psb_setelem(index+(1-ixb),val,xp,descp,info)
end if
res = info
return
end function psb_c_zsetelem
end module psb_z_tools_cbind_mod

24
cbind/linsolve/psb_base_linsolve_cbind_mod.f90
Unescape Escape View File

@ -6,6 +6,7 @@ module psb_base_linsolve_cbind_mod
type, bind(c) :: solveroptions
integer(psb_c_ipk_) :: iter, itmax, itrace, irst, istop
real(c_double) :: eps, err
type(psb_c_object_type) :: s1, s2
end type solveroptions
contains
@ -20,24 +21,27 @@ contains
options%istop = 2
options%irst = 10
options%eps = 1.d-6
options%s1 = psb_c_get_new_object()
options%s2 = psb_c_get_new_object()
res = 0
end function psb_c_DefaultSolverOptions
function psb_c_PrintSolverOptions(options)&
& bind(c,name='psb_c_PrintSolverOptions') result(res)
implicit none
type(solveroptions) :: options
type(solveroptions), value :: options
integer(psb_c_ipk_) :: res
write(*,*) 'PSBLAS C Interface Solver Options '
write(*,*) ' Maximum number of iterations :', options%itmax
write(*,*) ' Tracing :', options%itrace
write(*,*) ' Stopping Criterion :', options%istop
write(*,*) ' Restart :', options%irst
write(*,*) ' EPS (tolerance) :', options%eps
write(0,*) 'PSBLAS C Interface Solver Options '
write(0,*) ' Maximum number of iterations :', options%itmax
write(0,*) ' Tracing :', options%itrace
write(0,*) ' Stopping Criterion :', options%istop
write(0,*) ' Restart :', options%irst
write(0,*) ' EPS (tolerance) :', options%eps
write(0,*) ' S1 scaling :', c_associated(options%s1%item)
write(0,*) ' S2 scaling :', c_associated(options%s2%item)
res = 0
end function psb_c_PrintSolverOptions
end module psb_base_linsolve_cbind_mod

51
cbind/linsolve/psb_clinsolve_cbind_mod.f90
Unescape Escape View File

@ -23,14 +23,16 @@ contains
res= psb_c_ckrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_ckrylov
function psb_c_ckrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istop,s1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
use psb_linsolve_mod
use psb_objhandle_mod
@ -46,12 +48,14 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
type(psb_cprec_type), pointer :: precp
type(psb_c_vect_type), pointer :: xp, bp
type(psb_c_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
real(psb_spk_) :: feps,ferr
@ -81,6 +85,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call psb_stringc2f(methd,fmethd)
@ -89,14 +103,33 @@ contains
fitrace = itrace
first = irst
fistop = istop
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
end function psb_c_ckrylov_opt

51
cbind/linsolve/psb_dlinsolve_cbind_mod.f90
Unescape Escape View File

@ -23,14 +23,16 @@ contains
res= psb_c_dkrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_dkrylov
function psb_c_dkrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istop,s1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
use psb_linsolve_mod
use psb_objhandle_mod
@ -46,12 +48,14 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
type(psb_dprec_type), pointer :: precp
type(psb_d_vect_type), pointer :: xp, bp
type(psb_d_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
real(psb_dpk_) :: feps,ferr
@ -81,6 +85,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call psb_stringc2f(methd,fmethd)
@ -89,14 +103,33 @@ contains
fitrace = itrace
first = irst
fistop = istop
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
end function psb_c_dkrylov_opt

5
cbind/linsolve/psb_linsolve_cbind.h
Unescape Escape View File

@ -22,11 +22,12 @@ typedef struct psb_c_solveroptions {
int istop; /* Stopping criterion: 1:backward error 2: ||r||_2/||b||_2 */
double eps; /* Stopping tolerance */
double err; /* Convergence indicator on exit */
void *s1;
void *s2;
} psb_c_SolverOptions;
int psb_c_DefaultSolverOptions(psb_c_SolverOptions *opt);
int psb_c_PrintSolverOptions(psb_c_SolverOptions *opt);
int psb_c_PrintSolverOptions(psb_c_SolverOptions opt);
int psb_c_skrylov(const char *method, psb_c_sspmat *ah, psb_c_sprec *ph,
psb_c_svector *bh, psb_c_svector *xh,
psb_c_descriptor *cdh, psb_c_SolverOptions *opt);

51
cbind/linsolve/psb_slinsolve_cbind_mod.f90
Unescape Escape View File

@ -23,14 +23,16 @@ contains
res= psb_c_skrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_skrylov
function psb_c_skrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istop,s1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
use psb_linsolve_mod
use psb_objhandle_mod
@ -46,12 +48,14 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
type(psb_sprec_type), pointer :: precp
type(psb_s_vect_type), pointer :: xp, bp
type(psb_s_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
real(psb_spk_) :: feps,ferr
@ -81,6 +85,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call psb_stringc2f(methd,fmethd)
@ -89,14 +103,33 @@ contains
fitrace = itrace
first = irst
fistop = istop
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
end function psb_c_skrylov_opt

51
cbind/linsolve/psb_zlinsolve_cbind_mod.f90
Unescape Escape View File

@ -23,14 +23,16 @@ contains
res= psb_c_zkrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_zkrylov
function psb_c_zkrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istop,s1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
use psb_linsolve_mod
use psb_objhandle_mod
@ -46,12 +48,14 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
type(psb_zprec_type), pointer :: precp
type(psb_z_vect_type), pointer :: xp, bp
type(psb_z_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
real(psb_dpk_) :: feps,ferr
@ -81,6 +85,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call psb_stringc2f(methd,fmethd)
@ -89,14 +103,33 @@ contains
fitrace = itrace
first = irst
fistop = istop
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
end function psb_c_zkrylov_opt

4
cuda/License-spgpu
Unescape Escape View File

@ -24,3 +24,7 @@
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
<<<<<<< HEAD
=======
>>>>>>> kinsol-stop

5
linsolve/Makefile
Unescape Escape View File

@ -9,6 +9,8 @@ MODDIR=../modules
MODOBJS= psb_base_linsolve_conv_mod.o \
psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o \
psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o \
psb_s_linsolve_mod.o psb_c_linsolve_mod.o \
psb_d_linsolve_mod.o psb_z_linsolve_mod.o \
psb_linsolve_mod.o
OBJS=$(MODOBJS)
@ -37,6 +39,9 @@ impld: $(OBJS)
psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o: psb_base_linsolve_conv_mod.o
psb_linsolve_conv_mod.o: psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o
psb_linsolve_mod.o: psb_s_linsolve_mod.o psb_c_linsolve_mod.o psb_d_linsolve_mod.o psb_z_linsolve_mod.o
$(F90OBJS): $(MODOBJS)
$(OBJS): $(MODDIR)/$(PRECMODNAME)$(.mod) $(MODDIR)/$(BASEMODNAME)$(.mod)

10
linsolve/impl/psb_cbicg.f90
Unescape Escape View File

@ -95,7 +95,7 @@
!
subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -111,6 +111,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_spk_), optional, intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! !$ local data
complex(psb_spk_), allocatable, target :: aux(:)
type(psb_c_vect_type), allocatable, target :: wwrk(:)
@ -236,7 +237,8 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -262,7 +264,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
rho = czero
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -316,7 +318,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-alpha,q,cone,r,desc_a,info)
call psb_geaxpby(-alpha,qt,cone,rt,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_ccg.F90
Unescape Escape View File

@ -96,7 +96,7 @@
!
!
subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -112,6 +112,8 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = Local data
complex(psb_spk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:)
integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@ -245,7 +247,8 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
rho = czero
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
&desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -289,7 +292,7 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(alpha,p,cone,x,desc_a,info)
call psb_geaxpby(-alpha,q,cone,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_ccgs.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! estimate of) residual.
!
Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_spk_), allocatable, target :: aux(:)
type(psb_c_vect_type), allocatable, target :: wwrk(:)
@ -223,7 +224,8 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -246,7 +248,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -320,7 +322,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_ccgstab.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! where r is the (preconditioned, recursive
! estimate of) residual.
!
Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = Local data
complex(psb_spk_), allocatable, target :: aux(:),wwrk(:,:)
type(psb_c_vect_type) :: q, r, p, v, s, t, z, f
@ -235,7 +236,8 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
End If
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -252,7 +254,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(cone,r,czero,q,desc_a,info)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) then
info=psb_err_from_subroutine_
@ -372,7 +374,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(omega,z,cone,x,desc_a,info)
call psb_geaxpby(cone,s,czero,r,desc_a,info)
call psb_geaxpby(-omega,t,cone,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')

10
linsolve/impl/psb_ccgstabl.f90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -120,6 +120,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_spk_), allocatable, target :: aux(:), gamma(:),&
& gamma1(:), gamma2(:), taum(:,:), sigma(:)
@ -267,7 +268,8 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
rt0 => wwrk(10)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -305,7 +307,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',b%get_nrows()
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -409,7 +411,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-gamma1(j),rh(j),cone,rh(0),desc_a,info)
enddo
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_cfcg.F90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -120,6 +120,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = Local data
type(psb_c_vect_type) :: v, w, d , q, r
complex(psb_spk_) :: alpha, beta, delta, gamma, theta
@ -227,7 +228,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
& a,x,b,eps,desc_a,stopdat,info)
& a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
itx = 0
restart: do
@ -246,7 +247,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart'
exit restart
end if
@ -302,7 +303,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
itx = itx + 1
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration'
exit restart
end if

10
linsolve/impl/psb_cgcr.f90
Unescape Escape View File

@ -106,7 +106,7 @@
!
subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace, irst, istop)
& itmax,iter,err,itrace, irst, istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -124,6 +124,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_spk_), allocatable :: alpha(:), h(:,:)
type(psb_c_vect_type), allocatable :: z(:), c(:), c_scale(:)
@ -253,7 +254,8 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
restart: do
if (itx>= itmax_) exit restart
h = czero
@ -276,7 +278,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
nrst = nrst + 1
@ -307,7 +309,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(cone, r, czero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), cone, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (j >= irst) exit iteration

30
linsolve/impl/psb_ckrylov.f90
Unescape Escape View File

@ -80,7 +80,7 @@
! estimate of) residual
!
Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod,only : psb_cprec_type
@ -97,11 +97,12 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_c_vect_type), intent(inout), optional :: s1, s2
abstract interface
subroutine psb_ckryl_vect(a,prec,b,x,eps,&
& desc_a,info,itmax,iter,err,itrace,istop)
& desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
import :: psb_ipk_, psb_spk_, psb_desc_type, &
& psb_cspmat_type, psb_cprec_type, psb_c_vect_type
type(psb_cspmat_type), intent(in) :: a
@ -114,9 +115,10 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_spk_), optional, intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_ckryl_vect
Subroutine psb_ckryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,irst,istop)
&itmax,iter,err, itrace,irst,istop,s1,s2)
import :: psb_ipk_, psb_spk_, psb_desc_type, &
& psb_cspmat_type, psb_cprec_type, psb_c_vect_type
Type(psb_cspmat_type), Intent(in) :: a
@ -129,9 +131,10 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_ckryl_rest_vect
Subroutine psb_ckryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,istop,cond)
&itmax,iter,err, itrace,istop,cond,s1,s2)
import :: psb_ipk_, psb_spk_, psb_desc_type, &
& psb_cspmat_type, psb_cprec_type, psb_c_vect_type
Type(psb_cspmat_type), Intent(in) :: a
@ -144,6 +147,7 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err, cond
type(psb_c_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_ckryl_cond_vect
end interface
@ -180,37 +184,37 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
select case(psb_toupper(method))
case('CG')
call psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('FCG')
call psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('GCR')
call psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('CGS')
call psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICG')
call psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICGSTAB')
call psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('RGMRES','GMRES')
call psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case('MINRES','PMINRES')
call psb_cminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case default
if (me == 0) write(psb_err_unit,*) trim(name),&
& ': Warning: Unknown method ',method,&
& ', defaulting to BiCGSTAB'
call psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
end select
if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)

52
linsolve/impl/psb_crgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! 3: Same as 2 but with X and B scaled
! by s1 and s2
! where r is the (preconditioned, recursive
! estimate of) residual.
! irst - integer(optional) Input: restart parameter
!
subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
@ -124,6 +126,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_c_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_spk_), allocatable :: aux(:)
complex(psb_spk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:)
@ -268,9 +271,20 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(cone,s1,b,czero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(cone,s1,b,czero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-cone,a,x,cone,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_errpush(info,name)
goto 9999
end if
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
rs(1) = psb_genrm2(v(1),desc_a,info)
rs(2:) = czero
if (info /= psb_success_) then
@ -378,8 +394,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(cone,x,czero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(cone,b,czero,w1,desc_a,info)
call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(cone,xt,czero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
complex(psb_spk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_c_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_cprec_type) :: prec
type(psb_c_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),cone,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(cone,w1,cone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_crgmres_vect

11
linsolve/impl/psb_crichardson.f90
Unescape Escape View File

@ -120,8 +120,17 @@ Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else

10
linsolve/impl/psb_dbicg.f90
Unescape Escape View File

@ -95,7 +95,7 @@
!
subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -111,6 +111,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_dpk_), optional, intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! !$ local data
real(psb_dpk_), allocatable, target :: aux(:)
type(psb_d_vect_type), allocatable, target :: wwrk(:)
@ -236,7 +237,8 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -262,7 +264,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
rho = dzero
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -316,7 +318,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-alpha,q,done,r,desc_a,info)
call psb_geaxpby(-alpha,qt,done,rt,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_dcg.F90
Unescape Escape View File

@ -96,7 +96,7 @@
!
!
subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -112,6 +112,8 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = Local data
real(psb_dpk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:)
integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@ -253,7 +255,8 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
rho = dzero
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
&desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -306,7 +309,7 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(alpha,p,done,x,desc_a,info)
call psb_geaxpby(-alpha,q,done,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_dcgs.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! estimate of) residual.
!
Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_dpk_), allocatable, target :: aux(:)
type(psb_d_vect_type), allocatable, target :: wwrk(:)
@ -223,7 +224,8 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -246,7 +248,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -320,7 +322,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_dcgstab.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! where r is the (preconditioned, recursive
! estimate of) residual.
!
Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = Local data
real(psb_dpk_), allocatable, target :: aux(:),wwrk(:,:)
type(psb_d_vect_type) :: q, r, p, v, s, t, z, f
@ -235,7 +236,8 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
End If
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -252,7 +254,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(done,r,dzero,q,desc_a,info)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) then
info=psb_err_from_subroutine_
@ -372,7 +374,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(omega,z,done,x,desc_a,info)
call psb_geaxpby(done,s,dzero,r,desc_a,info)
call psb_geaxpby(-omega,t,done,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')

10
linsolve/impl/psb_dcgstabl.f90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -120,6 +120,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_dpk_), allocatable, target :: aux(:), gamma(:),&
& gamma1(:), gamma2(:), taum(:,:), sigma(:)
@ -267,7 +268,8 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
rt0 => wwrk(10)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -305,7 +307,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',b%get_nrows()
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -409,7 +411,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-gamma1(j),rh(j),done,rh(0),desc_a,info)
enddo
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_dfcg.F90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -120,6 +120,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = Local data
type(psb_d_vect_type) :: v, w, d , q, r
real(psb_dpk_) :: alpha, beta, delta, gamma, theta
@ -227,7 +228,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
& a,x,b,eps,desc_a,stopdat,info)
& a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
itx = 0
restart: do
@ -246,7 +247,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart'
exit restart
end if
@ -302,7 +303,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
itx = itx + 1
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration'
exit restart
end if

10
linsolve/impl/psb_dgcr.f90
Unescape Escape View File

@ -106,7 +106,7 @@
!
subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace, irst, istop)
& itmax,iter,err,itrace, irst, istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -124,6 +124,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_dpk_), allocatable :: alpha(:), h(:,:)
type(psb_d_vect_type), allocatable :: z(:), c(:), c_scale(:)
@ -253,7 +254,8 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
restart: do
if (itx>= itmax_) exit restart
h = dzero
@ -276,7 +278,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
nrst = nrst + 1
@ -307,7 +309,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(done, r, dzero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), done, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (j >= irst) exit iteration

30
linsolve/impl/psb_dkrylov.f90
Unescape Escape View File

@ -80,7 +80,7 @@
! estimate of) residual
!
Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod,only : psb_dprec_type
@ -97,11 +97,12 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_d_vect_type), intent(inout), optional :: s1, s2
abstract interface
subroutine psb_dkryl_vect(a,prec,b,x,eps,&
& desc_a,info,itmax,iter,err,itrace,istop)
& desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_dspmat_type, psb_dprec_type, psb_d_vect_type
type(psb_dspmat_type), intent(in) :: a
@ -114,9 +115,10 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_dpk_), optional, intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_dkryl_vect
Subroutine psb_dkryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,irst,istop)
&itmax,iter,err, itrace,irst,istop,s1,s2)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_dspmat_type, psb_dprec_type, psb_d_vect_type
Type(psb_dspmat_type), Intent(in) :: a
@ -129,9 +131,10 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_dkryl_rest_vect
Subroutine psb_dkryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,istop,cond)
&itmax,iter,err, itrace,istop,cond,s1,s2)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_dspmat_type, psb_dprec_type, psb_d_vect_type
Type(psb_dspmat_type), Intent(in) :: a
@ -144,6 +147,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err, cond
type(psb_d_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_dkryl_cond_vect
end interface
@ -180,37 +184,37 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
select case(psb_toupper(method))
case('CG')
call psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('FCG')
call psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('GCR')
call psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('CGS')
call psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICG')
call psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICGSTAB')
call psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('RGMRES','GMRES')
call psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case('MINRES','PMINRES')
call psb_dminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case default
if (me == 0) write(psb_err_unit,*) trim(name),&
& ': Warning: Unknown method ',method,&
& ', defaulting to BiCGSTAB'
call psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
end select
if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)

52
linsolve/impl/psb_drgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! 3: Same as 2 but with X and B scaled
! by s1 and s2
! where r is the (preconditioned, recursive
! estimate of) residual.
! irst - integer(optional) Input: restart parameter
!
subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
@ -124,6 +126,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_d_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_dpk_), allocatable :: aux(:)
real(psb_dpk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:)
@ -268,9 +271,20 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(done,s1,b,dzero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(done,s1,b,dzero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-done,a,x,done,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_errpush(info,name)
goto 9999
end if
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
rs(1) = psb_genrm2(v(1),desc_a,info)
rs(2:) = dzero
if (info /= psb_success_) then
@ -378,8 +394,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(done,x,dzero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(done,b,dzero,w1,desc_a,info)
call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(done,xt,dzero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
real(psb_dpk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_d_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_dprec_type) :: prec
type(psb_d_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),done,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(done,w1,done,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_drgmres_vect

11
linsolve/impl/psb_drichardson.f90
Unescape Escape View File

@ -120,8 +120,17 @@ Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else

10
linsolve/impl/psb_sbicg.f90
Unescape Escape View File

@ -95,7 +95,7 @@
!
subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -111,6 +111,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_spk_), optional, intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! !$ local data
real(psb_spk_), allocatable, target :: aux(:)
type(psb_s_vect_type), allocatable, target :: wwrk(:)
@ -236,7 +237,8 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -262,7 +264,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
rho = szero
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -316,7 +318,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-alpha,q,sone,r,desc_a,info)
call psb_geaxpby(-alpha,qt,sone,rt,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_scg.F90
Unescape Escape View File

@ -96,7 +96,7 @@
!
!
subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -112,6 +112,8 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = Local data
real(psb_spk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:)
integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@ -253,7 +255,8 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
rho = szero
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
&desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -306,7 +309,7 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(alpha,p,sone,x,desc_a,info)
call psb_geaxpby(-alpha,q,sone,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_scgs.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! estimate of) residual.
!
Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_spk_), allocatable, target :: aux(:)
type(psb_s_vect_type), allocatable, target :: wwrk(:)
@ -223,7 +224,8 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -246,7 +248,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -320,7 +322,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_scgstab.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! where r is the (preconditioned, recursive
! estimate of) residual.
!
Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = Local data
real(psb_spk_), allocatable, target :: aux(:),wwrk(:,:)
type(psb_s_vect_type) :: q, r, p, v, s, t, z, f
@ -235,7 +236,8 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
End If
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -252,7 +254,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(sone,r,szero,q,desc_a,info)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) then
info=psb_err_from_subroutine_
@ -372,7 +374,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(omega,z,sone,x,desc_a,info)
call psb_geaxpby(sone,s,szero,r,desc_a,info)
call psb_geaxpby(-omega,t,sone,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')

10
linsolve/impl/psb_scgstabl.f90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -120,6 +120,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_spk_), allocatable, target :: aux(:), gamma(:),&
& gamma1(:), gamma2(:), taum(:,:), sigma(:)
@ -267,7 +268,8 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
rt0 => wwrk(10)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -305,7 +307,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',b%get_nrows()
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -409,7 +411,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-gamma1(j),rh(j),sone,rh(0),desc_a,info)
enddo
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_sfcg.F90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -120,6 +120,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = Local data
type(psb_s_vect_type) :: v, w, d , q, r
real(psb_spk_) :: alpha, beta, delta, gamma, theta
@ -227,7 +228,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
& a,x,b,eps,desc_a,stopdat,info)
& a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
itx = 0
restart: do
@ -246,7 +247,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart'
exit restart
end if
@ -302,7 +303,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
itx = itx + 1
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration'
exit restart
end if

10
linsolve/impl/psb_sgcr.f90
Unescape Escape View File

@ -106,7 +106,7 @@
!
subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace, irst, istop)
& itmax,iter,err,itrace, irst, istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -124,6 +124,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_spk_), Optional, Intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_spk_), allocatable :: alpha(:), h(:,:)
type(psb_s_vect_type), allocatable :: z(:), c(:), c_scale(:)
@ -253,7 +254,8 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
restart: do
if (itx>= itmax_) exit restart
h = szero
@ -276,7 +278,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
nrst = nrst + 1
@ -307,7 +309,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(sone, r, szero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), sone, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (j >= irst) exit iteration

30
linsolve/impl/psb_skrylov.f90
Unescape Escape View File

@ -80,7 +80,7 @@
! estimate of) residual
!
Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod,only : psb_sprec_type
@ -97,11 +97,12 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_s_vect_type), intent(inout), optional :: s1, s2
abstract interface
subroutine psb_skryl_vect(a,prec,b,x,eps,&
& desc_a,info,itmax,iter,err,itrace,istop)
& desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
import :: psb_ipk_, psb_spk_, psb_desc_type, &
& psb_sspmat_type, psb_sprec_type, psb_s_vect_type
type(psb_sspmat_type), intent(in) :: a
@ -114,9 +115,10 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_spk_), optional, intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_skryl_vect
Subroutine psb_skryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,irst,istop)
&itmax,iter,err, itrace,irst,istop,s1,s2)
import :: psb_ipk_, psb_spk_, psb_desc_type, &
& psb_sspmat_type, psb_sprec_type, psb_s_vect_type
Type(psb_sspmat_type), Intent(in) :: a
@ -129,9 +131,10 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_skryl_rest_vect
Subroutine psb_skryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,istop,cond)
&itmax,iter,err, itrace,istop,cond,s1,s2)
import :: psb_ipk_, psb_spk_, psb_desc_type, &
& psb_sspmat_type, psb_sprec_type, psb_s_vect_type
Type(psb_sspmat_type), Intent(in) :: a
@ -144,6 +147,7 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err, cond
type(psb_s_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_skryl_cond_vect
end interface
@ -180,37 +184,37 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
select case(psb_toupper(method))
case('CG')
call psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('FCG')
call psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('GCR')
call psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('CGS')
call psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICG')
call psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICGSTAB')
call psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('RGMRES','GMRES')
call psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case('MINRES','PMINRES')
call psb_sminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case default
if (me == 0) write(psb_err_unit,*) trim(name),&
& ': Warning: Unknown method ',method,&
& ', defaulting to BiCGSTAB'
call psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
end select
if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)

52
linsolve/impl/psb_srgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! 3: Same as 2 but with X and B scaled
! by s1 and s2
! where r is the (preconditioned, recursive
! estimate of) residual.
! irst - integer(optional) Input: restart parameter
!
subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
@ -124,6 +126,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
type(psb_s_vect_type), intent(inout), optional :: s1, s2
! = local data
real(psb_spk_), allocatable :: aux(:)
real(psb_spk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:)
@ -268,9 +271,20 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(sone,s1,b,szero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(sone,s1,b,szero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-sone,a,x,sone,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_errpush(info,name)
goto 9999
end if
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
rs(1) = psb_genrm2(v(1),desc_a,info)
rs(2:) = szero
if (info /= psb_success_) then
@ -378,8 +394,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(sone,x,szero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(sone,b,szero,w1,desc_a,info)
call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(sone,xt,szero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
real(psb_spk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_s_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_sprec_type) :: prec
type(psb_s_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),sone,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(sone,w1,sone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_srgmres_vect

11
linsolve/impl/psb_srichardson.f90
Unescape Escape View File

@ -120,8 +120,17 @@ Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else

10
linsolve/impl/psb_zbicg.f90
Unescape Escape View File

@ -95,7 +95,7 @@
!
subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -111,6 +111,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace, istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_dpk_), optional, intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! !$ local data
complex(psb_dpk_), allocatable, target :: aux(:)
type(psb_z_vect_type), allocatable, target :: wwrk(:)
@ -236,7 +237,8 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -262,7 +264,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
rho = zzero
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -316,7 +318,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-alpha,q,zone,r,desc_a,info)
call psb_geaxpby(-alpha,qt,zone,rt,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_zcg.F90
Unescape Escape View File

@ -96,7 +96,7 @@
!
!
subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -112,6 +112,8 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = Local data
complex(psb_dpk_), allocatable, target :: aux(:),td(:),tu(:),eig(:),ewrk(:)
integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@ -245,7 +247,8 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
rho = zzero
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
&desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -289,7 +292,7 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(alpha,p,zone,x,desc_a,info)
call psb_geaxpby(-alpha,q,zone,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_zcgs.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! estimate of) residual.
!
Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
& itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_dpk_), allocatable, target :: aux(:)
type(psb_z_vect_type), allocatable, target :: wwrk(:)
@ -223,7 +224,8 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -246,7 +248,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -320,7 +322,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

10
linsolve/impl/psb_zcgstab.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! where r is the (preconditioned, recursive
! estimate of) residual.
!
Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -109,6 +109,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = Local data
complex(psb_dpk_), allocatable, target :: aux(:),wwrk(:,:)
type(psb_z_vect_type) :: q, r, p, v, s, t, z, f
@ -235,7 +236,8 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
End If
itx = 0
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -252,7 +254,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(zone,r,zzero,q,desc_a,info)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) then
info=psb_err_from_subroutine_
@ -372,7 +374,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
call psb_geaxpby(omega,z,zone,x,desc_a,info)
call psb_geaxpby(zone,s,zzero,r,desc_a,info)
call psb_geaxpby(-omega,t,zone,r,desc_a,info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (psb_errstatus_fatal()) Then
call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')

10
linsolve/impl/psb_zcgstabl.f90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -120,6 +120,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_dpk_), allocatable, target :: aux(:), gamma(:),&
& gamma1(:), gamma2(:), taum(:,:), sigma(:)
@ -267,7 +268,8 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
rt0 => wwrk(10)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -305,7 +307,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',b%get_nrows()
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999
@ -409,7 +411,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(-gamma1(j),rh(j),zone,rh(0),desc_a,info)
enddo
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
if (info /= psb_success_) Then
call psb_errpush(psb_err_from_subroutine_non_,name)
goto 9999

9
linsolve/impl/psb_zfcg.F90
Unescape Escape View File

@ -104,7 +104,7 @@
!
!
subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
& itmax,iter,err,itrace,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -120,6 +120,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = Local data
type(psb_z_vect_type) :: v, w, d , q, r
complex(psb_dpk_) :: alpha, beta, delta, gamma, theta
@ -227,7 +228,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
& scratch=.true.,mold=x%v)
call psb_init_conv(methdname,istop_,itrace_,itmax_,&
& a,x,b,eps,desc_a,stopdat,info)
& a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
itx = 0
restart: do
@ -246,7 +247,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from restart'
exit restart
end if
@ -302,7 +303,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
itx = itx + 1
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
if (debug.and.(me==0)) write(0,*) name,' Exit on convergence from iteration'
exit restart
end if

10
linsolve/impl/psb_zgcr.f90
Unescape Escape View File

@ -106,7 +106,7 @@
!
subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace, irst, istop)
& itmax,iter,err,itrace, irst, istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -124,6 +124,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
real(psb_dpk_), Optional, Intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_dpk_), allocatable :: alpha(:), h(:,:)
type(psb_z_vect_type), allocatable :: z(:), c(:), c_scale(:)
@ -253,7 +254,8 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
itx = 0
nrst = -1
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
& desc_a,stopdat,info,s1=s1,s2=s2)
restart: do
if (itx>= itmax_) exit restart
h = zzero
@ -276,7 +278,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
nrst = nrst + 1
@ -307,7 +309,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(zone, r, zzero, r, desc_a, info)
call psb_geaxpby(-alpha(j), c_scale(j), zone, r, desc_a, info)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
if (j >= irst) exit iteration

30
linsolve/impl/psb_zkrylov.f90
Unescape Escape View File

@ -80,7 +80,7 @@
! estimate of) residual
!
Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod
use psb_prec_mod,only : psb_zprec_type
@ -97,11 +97,12 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_z_vect_type), intent(inout), optional :: s1, s2
abstract interface
subroutine psb_zkryl_vect(a,prec,b,x,eps,&
& desc_a,info,itmax,iter,err,itrace,istop)
& desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_zspmat_type, psb_zprec_type, psb_z_vect_type
type(psb_zspmat_type), intent(in) :: a
@ -114,9 +115,10 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), optional, intent(in) :: itmax, itrace,istop
integer(psb_ipk_), optional, intent(out) :: iter
real(psb_dpk_), optional, intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_zkryl_vect
Subroutine psb_zkryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,irst,istop)
&itmax,iter,err, itrace,irst,istop,s1,s2)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_zspmat_type, psb_zprec_type, psb_z_vect_type
Type(psb_zspmat_type), Intent(in) :: a
@ -129,9 +131,10 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_zkryl_rest_vect
Subroutine psb_zkryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err, itrace,istop,cond)
&itmax,iter,err, itrace,istop,cond,s1,s2)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_zspmat_type, psb_zprec_type, psb_z_vect_type
Type(psb_zspmat_type), Intent(in) :: a
@ -144,6 +147,7 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err, cond
type(psb_z_vect_type), intent(inout), optional :: s1, s2
end subroutine psb_zkryl_cond_vect
end interface
@ -180,37 +184,37 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
select case(psb_toupper(method))
case('CG')
call psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('FCG')
call psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
&itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
case('GCR')
call psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('CGS')
call psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICG')
call psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('BICGSTAB')
call psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
case('RGMRES','GMRES')
call psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case('MINRES','PMINRES')
call psb_zminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
case default
if (me == 0) write(psb_err_unit,*) trim(name),&
& ': Warning: Unknown method ',method,&
& ', defaulting to BiCGSTAB'
call psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
&itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
end select
if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)

52
linsolve/impl/psb_zrgmres.f90
Unescape Escape View File

@ -102,13 +102,15 @@
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! 3: Same as 2 but with X and B scaled
! by s1 and s2
! where r is the (preconditioned, recursive
! estimate of) residual.
! irst - integer(optional) Input: restart parameter
!
subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop)
& itmax,iter,err,itrace,irst,istop,s1,s2)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
@ -124,6 +126,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
type(psb_z_vect_type), intent(inout), optional :: s1, s2
! = local data
complex(psb_dpk_), allocatable :: aux(:)
complex(psb_dpk_), allocatable :: c(:), s(:), h(:,:), rs(:), rst(:)
@ -268,9 +271,20 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(zone,s1,b,zzero,v(1),desc_a,info)
bni = psb_geamax(v(1),desc_a,info)
else
bni = psb_geamax(b,desc_a,info)
end if
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
if (present(s1)) then
call psb_gemlt(zone,s1,b,zzero,v(1),desc_a,info)
bn2 = psb_genrm2(v(1),desc_a,info)
else
bn2 = psb_genrm2(b,desc_a,info)
end if
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
@ -282,6 +296,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
end if
call psb_spmm(-zone,a,x,zone,v(1),desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
@ -323,7 +338,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_errpush(info,name)
goto 9999
end if
if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
rs(1) = psb_genrm2(v(1),desc_a,info)
rs(2:) = zzero
if (info /= psb_success_) then
@ -378,8 +394,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
@ -391,10 +413,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
rst = rs
call psb_geaxpby(zone,x,zzero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
call psb_geaxpby(zone,b,zzero,w1,desc_a,info)
call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
rni = psb_geamax(w1,desc_a,info)
xni = psb_geamax(xt,desc_a,info)
errnum = rni
@ -432,7 +455,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
! = x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
!
end select
@ -452,7 +476,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) !
end select
if (itx >= itmax_) then
@ -523,11 +547,12 @@ contains
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
complex(psb_dpk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_z_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_zprec_type) :: prec
type(psb_z_vect_type), intent(inout), optional :: s2
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
@ -539,12 +564,13 @@ contains
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
call w%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
call psb_geaxpby(rs(k),v(k),zone,w,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
if (present(s2)) call psb_gediv(s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
call psb_geaxpby(zone,w1,zone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_zrgmres_vect

11
linsolve/impl/psb_zrichardson.f90
Unescape Escape View File

@ -120,8 +120,17 @@ Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else

68
linsolve/psb_c_linsolve_conv_mod.f90
Unescape Escape View File

@ -83,15 +83,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (3)
case (psb_istop_rn2_abs_)
! Do nothing
case (psb_istop_rrn2_)
call psb_geall(r,desc_a,info)
call psb_geaxpby(cone,b,czero,r,desc_a,info)
call psb_spmm(-cone,a,x,cone,r,desc_a,info)
@ -108,8 +111,8 @@ contains
end if
stopdat%values(psb_ik_eps_) = eps
stopdat%values(psb_ik_errnum_) = dzero
stopdat%values(psb_ik_errden_) = done
stopdat%values(psb_ik_errnum_) = szero
stopdat%values(psb_ik_errden_) = sone
if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
& call log_header(methdname)
@ -123,7 +126,6 @@ contains
end subroutine psb_c_init_conv
function psb_c_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
use psb_base_mod
implicit none
@ -149,19 +151,26 @@ contains
res = .false.
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) =&
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = sone
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@ -201,8 +210,8 @@ contains
end function psb_c_check_conv
subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,&
& a,x,b,eps,desc_a,stopdat,info,s1,s2)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -213,6 +222,7 @@ contains
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer(psb_ipk_), intent(out) :: info
type(psb_c_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -236,15 +246,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (psb_istop_rn2_abs_)
! Do nothing
case (3)
case (psb_istop_rrn2_)
call psb_geasb(r,desc_a,info,scratch=.true.)
call psb_geaxpby(cone,b,czero,r,desc_a,info)
call psb_spmm(-cone,a,x,cone,r,desc_a,info)
@ -261,8 +274,8 @@ contains
end if
stopdat%values(psb_ik_eps_) = eps
stopdat%values(psb_ik_errnum_) = dzero
stopdat%values(psb_ik_errden_) = done
stopdat%values(psb_ik_errnum_) = szero
stopdat%values(psb_ik_errden_) = sone
if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
& call log_header(methdname)
@ -276,7 +289,8 @@ contains
end subroutine psb_c_init_conv_vect
function psb_c_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
function psb_c_check_conv_vect(methdname,it,x,r,&
& desc_a,stopdat,info,s1,s2) result(res)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -286,6 +300,7 @@ contains
type(psb_itconv_type) :: stopdat
logical :: res
integer(psb_ipk_), intent(out) :: info
type(psb_c_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -303,19 +318,26 @@ contains
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) = &
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = sone
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)

81
linsolve/psb_c_linsolve_mod.f90
Unescape Escape View File

@ -0,0 +1,81 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psb_linsolve_mod.f90
! Interfaces for linear solvers.
!
Module psb_c_linsolve_mod
use psb_const_mod
public
interface psb_krylov
Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, &
& psb_spk_, psb_c_vect_type
use psb_prec_mod, only : psb_cprec_type
character(len=*) :: method
Type(psb_cspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_c_vect_type), optional :: s1,s2
end Subroutine psb_ckrylov_vect
end interface
interface psb_richardson
Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, &
& psb_spk_, psb_c_vect_type
use psb_prec_mod, only : psb_cprec_type
Type(psb_cspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
end Subroutine psb_crichardson_vect
end interface
end module psb_c_linsolve_mod

60
linsolve/psb_d_linsolve_conv_mod.f90
Unescape Escape View File

@ -83,15 +83,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (3)
case (psb_istop_rn2_abs_)
! Do nothing
case (psb_istop_rrn2_)
call psb_geall(r,desc_a,info)
call psb_geaxpby(done,b,dzero,r,desc_a,info)
call psb_spmm(-done,a,x,done,r,desc_a,info)
@ -123,7 +126,6 @@ contains
end subroutine psb_d_init_conv
function psb_d_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
use psb_base_mod
implicit none
@ -149,19 +151,26 @@ contains
res = .false.
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) =&
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = done
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@ -201,8 +210,8 @@ contains
end function psb_d_check_conv
subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,&
& a,x,b,eps,desc_a,stopdat,info,s1,s2)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -213,6 +222,7 @@ contains
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer(psb_ipk_), intent(out) :: info
type(psb_d_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -236,15 +246,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (psb_istop_rn2_abs_)
! Do nothing
case (3)
case (psb_istop_rrn2_)
call psb_geasb(r,desc_a,info,scratch=.true.)
call psb_geaxpby(done,b,dzero,r,desc_a,info)
call psb_spmm(-done,a,x,done,r,desc_a,info)
@ -276,7 +289,8 @@ contains
end subroutine psb_d_init_conv_vect
function psb_d_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
function psb_d_check_conv_vect(methdname,it,x,r,&
& desc_a,stopdat,info,s1,s2) result(res)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -286,6 +300,7 @@ contains
type(psb_itconv_type) :: stopdat
logical :: res
integer(psb_ipk_), intent(out) :: info
type(psb_d_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -303,19 +318,26 @@ contains
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) = &
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = done
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)

81
linsolve/psb_d_linsolve_mod.f90
Unescape Escape View File

@ -0,0 +1,81 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psb_linsolve_mod.f90
! Interfaces for linear solvers.
!
Module psb_d_linsolve_mod
use psb_const_mod
public
interface psb_krylov
Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, &
& psb_dpk_, psb_d_vect_type
use psb_prec_mod, only : psb_dprec_type
character(len=*) :: method
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_d_vect_type), optional :: s1,s2
end Subroutine psb_dkrylov_vect
end interface
interface psb_richardson
Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, &
& psb_dpk_, psb_d_vect_type
use psb_prec_mod, only : psb_dprec_type
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
end Subroutine psb_drichardson_vect
end interface
end module psb_d_linsolve_mod

178
linsolve/psb_linsolve_mod.f90
Unescape Escape View File

@ -36,179 +36,9 @@
Module psb_linsolve_mod
use psb_const_mod
public
interface psb_krylov
Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, &
& psb_spk_, psb_s_vect_type
use psb_prec_mod, only : psb_sprec_type
character(len=*) :: method
Type(psb_sspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
end Subroutine psb_skrylov_vect
Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, &
& psb_spk_, psb_c_vect_type
use psb_prec_mod, only : psb_cprec_type
character(len=*) :: method
Type(psb_cspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
end Subroutine psb_ckrylov_vect
Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, &
& psb_dpk_, psb_d_vect_type
use psb_prec_mod, only : psb_dprec_type
character(len=*) :: method
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
end Subroutine psb_dkrylov_vect
Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, &
& psb_dpk_, psb_z_vect_type
use psb_prec_mod, only : psb_zprec_type
character(len=*) :: method
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
end Subroutine psb_zkrylov_vect
end interface
interface psb_richardson
Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, &
& psb_spk_, psb_s_vect_type
use psb_prec_mod, only : psb_sprec_type
Type(psb_sspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
end Subroutine psb_srichardson_vect
Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_cspmat_type, &
& psb_spk_, psb_c_vect_type
use psb_prec_mod, only : psb_cprec_type
Type(psb_cspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
end Subroutine psb_crichardson_vect
Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, &
& psb_dpk_, psb_d_vect_type
use psb_prec_mod, only : psb_dprec_type
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
end Subroutine psb_drichardson_vect
Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, &
& psb_dpk_, psb_z_vect_type
use psb_prec_mod, only : psb_zprec_type
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
end Subroutine psb_zrichardson_vect
end interface
use psb_s_linsolve_mod
use psb_d_linsolve_mod
use psb_c_linsolve_mod
use psb_z_linsolve_mod
end module psb_linsolve_mod

68
linsolve/psb_s_linsolve_conv_mod.f90
Unescape Escape View File

@ -83,15 +83,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (3)
case (psb_istop_rn2_abs_)
! Do nothing
case (psb_istop_rrn2_)
call psb_geall(r,desc_a,info)
call psb_geaxpby(sone,b,szero,r,desc_a,info)
call psb_spmm(-sone,a,x,sone,r,desc_a,info)
@ -108,8 +111,8 @@ contains
end if
stopdat%values(psb_ik_eps_) = eps
stopdat%values(psb_ik_errnum_) = dzero
stopdat%values(psb_ik_errden_) = done
stopdat%values(psb_ik_errnum_) = szero
stopdat%values(psb_ik_errden_) = sone
if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
& call log_header(methdname)
@ -123,7 +126,6 @@ contains
end subroutine psb_s_init_conv
function psb_s_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
use psb_base_mod
implicit none
@ -149,19 +151,26 @@ contains
res = .false.
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) =&
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = sone
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@ -201,8 +210,8 @@ contains
end function psb_s_check_conv
subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,&
& a,x,b,eps,desc_a,stopdat,info,s1,s2)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -213,6 +222,7 @@ contains
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer(psb_ipk_), intent(out) :: info
type(psb_s_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -236,15 +246,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (psb_istop_rn2_abs_)
! Do nothing
case (3)
case (psb_istop_rrn2_)
call psb_geasb(r,desc_a,info,scratch=.true.)
call psb_geaxpby(sone,b,szero,r,desc_a,info)
call psb_spmm(-sone,a,x,sone,r,desc_a,info)
@ -261,8 +274,8 @@ contains
end if
stopdat%values(psb_ik_eps_) = eps
stopdat%values(psb_ik_errnum_) = dzero
stopdat%values(psb_ik_errden_) = done
stopdat%values(psb_ik_errnum_) = szero
stopdat%values(psb_ik_errden_) = sone
if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
& call log_header(methdname)
@ -276,7 +289,8 @@ contains
end subroutine psb_s_init_conv_vect
function psb_s_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
function psb_s_check_conv_vect(methdname,it,x,r,&
& desc_a,stopdat,info,s1,s2) result(res)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -286,6 +300,7 @@ contains
type(psb_itconv_type) :: stopdat
logical :: res
integer(psb_ipk_), intent(out) :: info
type(psb_s_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -303,19 +318,26 @@ contains
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) = &
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = sone
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)

81
linsolve/psb_s_linsolve_mod.f90
Unescape Escape View File

@ -0,0 +1,81 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psb_linsolve_mod.f90
! Interfaces for linear solvers.
!
Module psb_s_linsolve_mod
use psb_const_mod
public
interface psb_krylov
Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, &
& psb_spk_, psb_s_vect_type
use psb_prec_mod, only : psb_sprec_type
character(len=*) :: method
Type(psb_sspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err,cond
type(psb_s_vect_type), optional :: s1,s2
end Subroutine psb_skrylov_vect
end interface
interface psb_richardson
Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_sspmat_type, &
& psb_spk_, psb_s_vect_type
use psb_prec_mod, only : psb_sprec_type
Type(psb_sspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err
end Subroutine psb_srichardson_vect
end interface
end module psb_s_linsolve_mod

60
linsolve/psb_z_linsolve_conv_mod.f90
Unescape Escape View File

@ -83,15 +83,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (3)
case (psb_istop_rn2_abs_)
! Do nothing
case (psb_istop_rrn2_)
call psb_geall(r,desc_a,info)
call psb_geaxpby(zone,b,zzero,r,desc_a,info)
call psb_spmm(-zone,a,x,zone,r,desc_a,info)
@ -123,7 +126,6 @@ contains
end subroutine psb_z_init_conv
function psb_z_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
use psb_base_mod
implicit none
@ -149,19 +151,26 @@ contains
res = .false.
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) =&
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = done
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@ -201,8 +210,8 @@ contains
end function psb_z_check_conv
subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,&
& a,x,b,eps,desc_a,stopdat,info,s1,s2)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -213,6 +222,7 @@ contains
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer(psb_ipk_), intent(out) :: info
type(psb_z_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -236,15 +246,18 @@ contains
stopdat%controls(psb_ik_itmax_) = itmax
select case(stopdat%controls(psb_ik_stopc_))
case (1)
case (psb_istop_ani_)
stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
if (info == psb_success_)&
& stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
case (2)
case (psb_istop_bn2_)
stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
case (psb_istop_rn2_abs_)
! Do nothing
case (3)
case (psb_istop_rrn2_)
call psb_geasb(r,desc_a,info,scratch=.true.)
call psb_geaxpby(zone,b,zzero,r,desc_a,info)
call psb_spmm(-zone,a,x,zone,r,desc_a,info)
@ -276,7 +289,8 @@ contains
end subroutine psb_z_init_conv_vect
function psb_z_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
function psb_z_check_conv_vect(methdname,it,x,r,&
& desc_a,stopdat,info,s1,s2) result(res)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
@ -286,6 +300,7 @@ contains
type(psb_itconv_type) :: stopdat
logical :: res
integer(psb_ipk_), intent(out) :: info
type(psb_z_vect_type), optional :: s1, s2
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me, np, err_act
@ -303,19 +318,26 @@ contains
select case(stopdat%controls(psb_ik_stopc_))
case(1)
case(psb_istop_ani_)
stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
if (info == psb_success_) &
& stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
stopdat%values(psb_ik_errden_) = &
& (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
& +stopdat%values(psb_ik_bni_))
case(2)
case(psb_istop_bn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
case(3)
case (psb_istop_rn2_abs_)
stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_abs_)
stopdat%values(psb_ik_errden_) = done
case(psb_istop_rrn2_)
stopdat%values(psb_ik_rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)

81
linsolve/psb_z_linsolve_mod.f90
Unescape Escape View File

@ -0,0 +1,81 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psb_linsolve_mod.f90
! Interfaces for linear solvers.
!
Module psb_z_linsolve_mod
use psb_const_mod
public
interface psb_krylov
Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,irst,istop,cond,s1,s2)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, &
& psb_dpk_, psb_z_vect_type
use psb_prec_mod, only : psb_zprec_type
character(len=*) :: method
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err,cond
type(psb_z_vect_type), optional :: s1,s2
end Subroutine psb_zkrylov_vect
end interface
interface psb_richardson
Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop)
use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_zspmat_type, &
& psb_dpk_, psb_z_vect_type
use psb_prec_mod, only : psb_zprec_type
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err
end Subroutine psb_zrichardson_vect
end interface
end module psb_z_linsolve_mod

4
prec/impl/psb_c_bjacprec_impl.f90
Unescape Escape View File

@ -644,6 +644,10 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold)
enddo
deallocate(prec%av,stat=info)
endif
if (allocated (prec%dv)) then
call prec%dv%free(info)
if (info == 0) deallocate(prec%dv)
end if
end if
if (.not.allocated(prec%av)) then
allocate(prec%av(psb_max_avsz),stat=info)

2
prec/impl/psb_c_diagprec_impl.f90
Unescape Escape View File

@ -270,7 +270,7 @@ subroutine psb_c_diag_precbld(a,desc_a,prec,info,amold,vmold,imold)
prec%d(i) = cone
end do
allocate(prec%dv,stat=info)
if (.not.allocated(prec%dv)) allocate(prec%dv,stat=info)
if (info == 0) then
if (present(vmold)) then
allocate(prec%dv%v,mold=vmold,stat=info)

4
prec/impl/psb_d_bjacprec_impl.f90
Unescape Escape View File

@ -644,6 +644,10 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold)
enddo
deallocate(prec%av,stat=info)
endif
if (allocated (prec%dv)) then
call prec%dv%free(info)
if (info == 0) deallocate(prec%dv)
end if
end if
if (.not.allocated(prec%av)) then
allocate(prec%av(psb_max_avsz),stat=info)

2
prec/impl/psb_d_diagprec_impl.f90
Unescape Escape View File

@ -270,7 +270,7 @@ subroutine psb_d_diag_precbld(a,desc_a,prec,info,amold,vmold,imold)
prec%d(i) = done
end do
allocate(prec%dv,stat=info)
if (.not.allocated(prec%dv)) allocate(prec%dv,stat=info)
if (info == 0) then
if (present(vmold)) then
allocate(prec%dv%v,mold=vmold,stat=info)

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

Write Preview
Loading…
Cancel
Save