mirror of
https://github.com/sfilippone/amg4psblas.git
synced 2026-10-07 15:15:07 +00:00
Compare commits
9
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Cirdans-Home
|
971363ccf3 | ||
|
Cirdans-Home
|
55b6b63b72 | ||
|
Cirdans-Home
|
6bfd227a80 | ||
 Salvatore Filippone
|
be4e7b4837 | ||
|
Salvatore Filippone
|
83c14397eb | ||
|
Salvatore Filippone
|
f67b61d23f | ||
|
Salvatore Filippone
|
a72690f9c3 | ||
|
Salvatore Filippone
|
5c2eef76ab | ||
|
Salvatore Filippone
|
c9e1719098 |
@@ -666,6 +666,7 @@ am__nodep
|
||||
AMDEPBACKSLASH
|
||||
AMDEP_FALSE
|
||||
AMDEP_TRUE
|
||||
am__quote
|
||||
am__include
|
||||
DEPDIR
|
||||
ac_ct_CC
|
||||
@@ -724,6 +725,7 @@ infodir
|
||||
docdir
|
||||
oldincludedir
|
||||
includedir
|
||||
runstatedir
|
||||
localstatedir
|
||||
sharedstatedir
|
||||
sysconfdir
|
||||
@@ -742,8 +744,7 @@ PACKAGE_VERSION
|
||||
PACKAGE_TARNAME
|
||||
PACKAGE_NAME
|
||||
PATH_SEPARATOR
|
||||
SHELL
|
||||
am__quote'
|
||||
SHELL'
|
||||
ac_subst_files=''
|
||||
ac_user_opts='
|
||||
enable_option_checking
|
||||
@@ -835,6 +836,7 @@ datadir='${datarootdir}'
|
||||
sysconfdir='${prefix}/etc'
|
||||
sharedstatedir='${prefix}/com'
|
||||
localstatedir='${prefix}/var'
|
||||
runstatedir='${localstatedir}/run'
|
||||
includedir='${prefix}/include'
|
||||
oldincludedir='/usr/include'
|
||||
docdir='${datarootdir}/doc/${PACKAGE_TARNAME}'
|
||||
@@ -1087,6 +1089,15 @@ do
|
||||
| -silent | --silent | --silen | --sile | --sil)
|
||||
silent=yes ;;
|
||||
|
||||
-runstatedir | --runstatedir | --runstatedi | --runstated \
|
||||
| --runstate | --runstat | --runsta | --runst | --runs \
|
||||
| --run | --ru | --r)
|
||||
ac_prev=runstatedir ;;
|
||||
-runstatedir=* | --runstatedir=* | --runstatedi=* | --runstated=* \
|
||||
| --runstate=* | --runstat=* | --runsta=* | --runst=* | --runs=* \
|
||||
| --run=* | --ru=* | --r=*)
|
||||
runstatedir=$ac_optarg ;;
|
||||
|
||||
-sbindir | --sbindir | --sbindi | --sbind | --sbin | --sbi | --sb)
|
||||
ac_prev=sbindir ;;
|
||||
-sbindir=* | --sbindir=* | --sbindi=* | --sbind=* | --sbin=* \
|
||||
@@ -1224,7 +1235,7 @@ fi
|
||||
for ac_var in exec_prefix prefix bindir sbindir libexecdir datarootdir \
|
||||
datadir sysconfdir sharedstatedir localstatedir includedir \
|
||||
oldincludedir docdir infodir htmldir dvidir pdfdir psdir \
|
||||
libdir localedir mandir
|
||||
libdir localedir mandir runstatedir
|
||||
do
|
||||
eval ac_val=\$$ac_var
|
||||
# Remove trailing slashes.
|
||||
@@ -1377,6 +1388,7 @@ Fine tuning of the installation directories:
|
||||
--sysconfdir=DIR read-only single-machine data [PREFIX/etc]
|
||||
--sharedstatedir=DIR modifiable architecture-independent data [PREFIX/com]
|
||||
--localstatedir=DIR modifiable single-machine data [PREFIX/var]
|
||||
--runstatedir=DIR modifiable per-process data [LOCALSTATEDIR/run]
|
||||
--libdir=DIR object code libraries [EPREFIX/lib]
|
||||
--includedir=DIR C header files [PREFIX/include]
|
||||
--oldincludedir=DIR C header files for non-gcc [/usr/include]
|
||||
@@ -2665,7 +2677,7 @@ fi
|
||||
|
||||
{ $as_echo "$as_me:${as_lineno-$LINENO}: Loaded $pac_cv_status_file $FC $MPIFC $BLACS_LIBS" >&5
|
||||
$as_echo "$as_me: Loaded $pac_cv_status_file $FC $MPIFC $BLACS_LIBS" >&6;}
|
||||
am__api_version='1.16'
|
||||
am__api_version='1.15'
|
||||
|
||||
ac_aux_dir=
|
||||
for ac_dir in "$srcdir" "$srcdir/.." "$srcdir/../.."; do
|
||||
@@ -3210,8 +3222,8 @@ MAKEINFO=${MAKEINFO-"${am_missing_run}makeinfo"}
|
||||
|
||||
# For better backward compatibility. To be removed once Automake 1.9.x
|
||||
# dies out for good. For more background, see:
|
||||
# <https://lists.gnu.org/archive/html/automake/2012-07/msg00001.html>
|
||||
# <https://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
|
||||
# <http://lists.gnu.org/archive/html/automake/2012-07/msg00001.html>
|
||||
# <http://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
|
||||
mkdir_p='$(MKDIR_P)'
|
||||
|
||||
# We need awk for the "check" target (and possibly the TAP driver). The
|
||||
@@ -3262,7 +3274,7 @@ END
|
||||
Aborting the configuration process, to ensure you take notice of the issue.
|
||||
|
||||
You can download and install GNU coreutils to get an 'rm' implementation
|
||||
that behaves properly: <https://www.gnu.org/software/coreutils/>.
|
||||
that behaves properly: <http://www.gnu.org/software/coreutils/>.
|
||||
|
||||
If you want to complete the configuration process using your problematic
|
||||
'rm' anyway, export the environment variable ACCEPT_INFERIOR_RM_PROGRAM
|
||||
@@ -4161,45 +4173,45 @@ DEPDIR="${am__leading_dot}deps"
|
||||
|
||||
ac_config_commands="$ac_config_commands depfiles"
|
||||
|
||||
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking whether ${MAKE-make} supports the include directive" >&5
|
||||
$as_echo_n "checking whether ${MAKE-make} supports the include directive... " >&6; }
|
||||
cat > confinc.mk << 'END'
|
||||
|
||||
am_make=${MAKE-make}
|
||||
cat > confinc << 'END'
|
||||
am__doit:
|
||||
@echo this is the am__doit target >confinc.out
|
||||
@echo this is the am__doit target
|
||||
.PHONY: am__doit
|
||||
END
|
||||
# If we don't find an include directive, just comment out the code.
|
||||
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for style of include used by $am_make" >&5
|
||||
$as_echo_n "checking for style of include used by $am_make... " >&6; }
|
||||
am__include="#"
|
||||
am__quote=
|
||||
# BSD make does it like this.
|
||||
echo '.include "confinc.mk" # ignored' > confmf.BSD
|
||||
# Other make implementations (GNU, Solaris 10, AIX) do it like this.
|
||||
echo 'include confinc.mk # ignored' > confmf.GNU
|
||||
_am_result=no
|
||||
for s in GNU BSD; do
|
||||
{ echo "$as_me:$LINENO: ${MAKE-make} -f confmf.$s && cat confinc.out" >&5
|
||||
(${MAKE-make} -f confmf.$s && cat confinc.out) >&5 2>&5
|
||||
ac_status=$?
|
||||
echo "$as_me:$LINENO: \$? = $ac_status" >&5
|
||||
(exit $ac_status); }
|
||||
case $?:`cat confinc.out 2>/dev/null` in #(
|
||||
'0:this is the am__doit target') :
|
||||
case $s in #(
|
||||
BSD) :
|
||||
am__include='.include' am__quote='"' ;; #(
|
||||
*) :
|
||||
am__include='include' am__quote='' ;;
|
||||
esac ;; #(
|
||||
*) :
|
||||
;;
|
||||
_am_result=none
|
||||
# First try GNU make style include.
|
||||
echo "include confinc" > confmf
|
||||
# Ignore all kinds of additional output from 'make'.
|
||||
case `$am_make -s -f confmf 2> /dev/null` in #(
|
||||
*the\ am__doit\ target*)
|
||||
am__include=include
|
||||
am__quote=
|
||||
_am_result=GNU
|
||||
;;
|
||||
esac
|
||||
if test "$am__include" != "#"; then
|
||||
_am_result="yes ($s style)"
|
||||
break
|
||||
fi
|
||||
done
|
||||
rm -f confinc.* confmf.*
|
||||
{ $as_echo "$as_me:${as_lineno-$LINENO}: result: ${_am_result}" >&5
|
||||
$as_echo "${_am_result}" >&6; }
|
||||
# Now try BSD make style include.
|
||||
if test "$am__include" = "#"; then
|
||||
echo '.include "confinc"' > confmf
|
||||
case `$am_make -s -f confmf 2> /dev/null` in #(
|
||||
*the\ am__doit\ target*)
|
||||
am__include=.include
|
||||
am__quote="\""
|
||||
_am_result=BSD
|
||||
;;
|
||||
esac
|
||||
fi
|
||||
|
||||
|
||||
{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $_am_result" >&5
|
||||
$as_echo "$_am_result" >&6; }
|
||||
rm -f confinc confmf
|
||||
|
||||
# Check whether --enable-dependency-tracking was given.
|
||||
if test "${enable_dependency_tracking+set}" = set; then :
|
||||
@@ -10111,7 +10123,7 @@ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
|
||||
#
|
||||
# INIT-COMMANDS
|
||||
#
|
||||
AMDEP_TRUE="$AMDEP_TRUE" MAKE="${MAKE-make}"
|
||||
AMDEP_TRUE="$AMDEP_TRUE" ac_aux_dir="$ac_aux_dir"
|
||||
|
||||
_ACEOF
|
||||
|
||||
@@ -10556,35 +10568,29 @@ $as_echo "$as_me: executing $ac_file commands" >&6;}
|
||||
# Older Autoconf quotes --file arguments for eval, but not when files
|
||||
# are listed without --file. Let's play safe and only enable the eval
|
||||
# if we detect the quoting.
|
||||
# TODO: see whether this extra hack can be removed once we start
|
||||
# requiring Autoconf 2.70 or later.
|
||||
case $CONFIG_FILES in #(
|
||||
*\'*) :
|
||||
eval set x "$CONFIG_FILES" ;; #(
|
||||
*) :
|
||||
set x $CONFIG_FILES ;; #(
|
||||
*) :
|
||||
;;
|
||||
esac
|
||||
case $CONFIG_FILES in
|
||||
*\'*) eval set x "$CONFIG_FILES" ;;
|
||||
*) set x $CONFIG_FILES ;;
|
||||
esac
|
||||
shift
|
||||
# Used to flag and report bootstrapping failures.
|
||||
am_rc=0
|
||||
for am_mf
|
||||
for mf
|
||||
do
|
||||
# Strip MF so we end up with the name of the file.
|
||||
am_mf=`$as_echo "$am_mf" | sed -e 's/:.*$//'`
|
||||
# Check whether this is an Automake generated Makefile which includes
|
||||
# dependency-tracking related rules and includes.
|
||||
# Grep'ing the whole file directly is not great: AIX grep has a line
|
||||
mf=`echo "$mf" | sed -e 's/:.*$//'`
|
||||
# Check whether this is an Automake generated Makefile or not.
|
||||
# We used to match only the files named 'Makefile.in', but
|
||||
# some people rename them; so instead we look at the file content.
|
||||
# Grep'ing the first line is not enough: some people post-process
|
||||
# each Makefile.in and add a new line on top of each file to say so.
|
||||
# Grep'ing the whole file is not good either: AIX grep has a line
|
||||
# limit of 2048, but all sed's we know have understand at least 4000.
|
||||
sed -n 's,^am--depfiles:.*,X,p' "$am_mf" | grep X >/dev/null 2>&1 \
|
||||
|| continue
|
||||
am_dirpart=`$as_dirname -- "$am_mf" ||
|
||||
$as_expr X"$am_mf" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
|
||||
X"$am_mf" : 'X\(//\)[^/]' \| \
|
||||
X"$am_mf" : 'X\(//\)$' \| \
|
||||
X"$am_mf" : 'X\(/\)' \| . 2>/dev/null ||
|
||||
$as_echo X"$am_mf" |
|
||||
if sed -n 's,^#.*generated by automake.*,X,p' "$mf" | grep X >/dev/null 2>&1; then
|
||||
dirpart=`$as_dirname -- "$mf" ||
|
||||
$as_expr X"$mf" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
|
||||
X"$mf" : 'X\(//\)[^/]' \| \
|
||||
X"$mf" : 'X\(//\)$' \| \
|
||||
X"$mf" : 'X\(/\)' \| . 2>/dev/null ||
|
||||
$as_echo X"$mf" |
|
||||
sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{
|
||||
s//\1/
|
||||
q
|
||||
@@ -10602,48 +10608,53 @@ $as_echo X"$am_mf" |
|
||||
q
|
||||
}
|
||||
s/.*/./; q'`
|
||||
am_filepart=`$as_basename -- "$am_mf" ||
|
||||
$as_expr X/"$am_mf" : '.*/\([^/][^/]*\)/*$' \| \
|
||||
X"$am_mf" : 'X\(//\)$' \| \
|
||||
X"$am_mf" : 'X\(/\)' \| . 2>/dev/null ||
|
||||
$as_echo X/"$am_mf" |
|
||||
sed '/^.*\/\([^/][^/]*\)\/*$/{
|
||||
else
|
||||
continue
|
||||
fi
|
||||
# Extract the definition of DEPDIR, am__include, and am__quote
|
||||
# from the Makefile without running 'make'.
|
||||
DEPDIR=`sed -n 's/^DEPDIR = //p' < "$mf"`
|
||||
test -z "$DEPDIR" && continue
|
||||
am__include=`sed -n 's/^am__include = //p' < "$mf"`
|
||||
test -z "$am__include" && continue
|
||||
am__quote=`sed -n 's/^am__quote = //p' < "$mf"`
|
||||
# Find all dependency output files, they are included files with
|
||||
# $(DEPDIR) in their names. We invoke sed twice because it is the
|
||||
# simplest approach to changing $(DEPDIR) to its actual value in the
|
||||
# expansion.
|
||||
for file in `sed -n "
|
||||
s/^$am__include $am__quote\(.*(DEPDIR).*\)$am__quote"'$/\1/p' <"$mf" | \
|
||||
sed -e 's/\$(DEPDIR)/'"$DEPDIR"'/g'`; do
|
||||
# Make sure the directory exists.
|
||||
test -f "$dirpart/$file" && continue
|
||||
fdir=`$as_dirname -- "$file" ||
|
||||
$as_expr X"$file" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
|
||||
X"$file" : 'X\(//\)[^/]' \| \
|
||||
X"$file" : 'X\(//\)$' \| \
|
||||
X"$file" : 'X\(/\)' \| . 2>/dev/null ||
|
||||
$as_echo X"$file" |
|
||||
sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{
|
||||
s//\1/
|
||||
q
|
||||
}
|
||||
/^X\/\(\/\/\)$/{
|
||||
/^X\(\/\/\)[^/].*/{
|
||||
s//\1/
|
||||
q
|
||||
}
|
||||
/^X\/\(\/\).*/{
|
||||
/^X\(\/\/\)$/{
|
||||
s//\1/
|
||||
q
|
||||
}
|
||||
/^X\(\/\).*/{
|
||||
s//\1/
|
||||
q
|
||||
}
|
||||
s/.*/./; q'`
|
||||
{ echo "$as_me:$LINENO: cd "$am_dirpart" \
|
||||
&& sed -e '/# am--include-marker/d' "$am_filepart" \
|
||||
| $MAKE -f - am--depfiles" >&5
|
||||
(cd "$am_dirpart" \
|
||||
&& sed -e '/# am--include-marker/d' "$am_filepart" \
|
||||
| $MAKE -f - am--depfiles) >&5 2>&5
|
||||
ac_status=$?
|
||||
echo "$as_me:$LINENO: \$? = $ac_status" >&5
|
||||
(exit $ac_status); } || am_rc=$?
|
||||
as_dir=$dirpart/$fdir; as_fn_mkdir_p
|
||||
# echo "creating $dirpart/$file"
|
||||
echo '# dummy' > "$dirpart/$file"
|
||||
done
|
||||
done
|
||||
if test $am_rc -ne 0; then
|
||||
{ { $as_echo "$as_me:${as_lineno-$LINENO}: error: in \`$ac_pwd':" >&5
|
||||
$as_echo "$as_me: error: in \`$ac_pwd':" >&2;}
|
||||
as_fn_error $? "Something went wrong bootstrapping makefile fragments
|
||||
for automatic dependency tracking. Try re-running configure with the
|
||||
'--disable-dependency-tracking' option to at least be able to build
|
||||
the package (albeit without support for automatic dependency tracking).
|
||||
See \`config.log' for more details" "$LINENO" 5; }
|
||||
fi
|
||||
{ am_dirpart=; unset am_dirpart;}
|
||||
{ am_filepart=; unset am_filepart;}
|
||||
{ am_mf=; unset am_mf;}
|
||||
{ am_rc=; unset am_rc;}
|
||||
rm -f conftest-deps.mk
|
||||
}
|
||||
;;
|
||||
|
||||
|
||||
@@ -105,116 +105,218 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,y,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
endif
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,89 +107,140 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,y,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -197,47 +248,125 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_c_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
sm%pa => a
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_c_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,6 +61,7 @@ subroutine mld_c_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,146 +102,63 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -250,6 +167,12 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_c_vect_type),intent(inout) :: y
|
||||
complex(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
complex(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_c_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,190 +108,90 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if (sweeps > 0) then
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_c_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
end associate
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,12 +71,11 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_c_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_c_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,12 +72,11 @@ subroutine mld_c_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_c_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -105,116 +105,218 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,y,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
endif
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,89 +107,140 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,y,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -197,47 +248,125 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_d_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
sm%pa => a
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_d_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,6 +61,7 @@ subroutine mld_d_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,146 +102,63 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -250,6 +167,12 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_d_vect_type),intent(inout) :: y
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_d_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,190 +108,90 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if (sweeps > 0) then
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_d_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
end associate
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,12 +71,11 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_d_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_d_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,12 +72,11 @@ subroutine mld_d_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_d_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -105,116 +105,218 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,y,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
endif
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,89 +107,140 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,y,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -197,47 +248,125 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_s_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
sm%pa => a
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_s_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,6 +61,7 @@ subroutine mld_s_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,146 +102,63 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -250,6 +167,12 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_s_vect_type),intent(inout) :: y
|
||||
real(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
real(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_s_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,190 +108,90 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if (sweeps > 0) then
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_s_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
end associate
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,12 +71,11 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_s_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_s_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,12 +72,11 @@ subroutine mld_s_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_s_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -105,116 +105,218 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,y,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
endif
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,89 +107,140 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,y,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -197,47 +248,125 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_z_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
sm%pa => a
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_z_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,6 +61,7 @@ subroutine mld_z_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,146 +102,63 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -250,6 +167,12 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_z_vect_type),intent(inout) :: y
|
||||
complex(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
complex(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_z_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,190 +108,90 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if (sweeps > 0) then
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_z_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
end associate
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,12 +71,11 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_z_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,6 +77,7 @@ subroutine mld_z_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,12 +72,11 @@ subroutine mld_z_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
if( sm%checkres ) sm%pa => a
|
||||
sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_z_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -65,6 +65,7 @@ module mld_c_as_smoother
|
||||
! parent type.
|
||||
! class(mld_c_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_cspmat_type), pointer :: pa => null()
|
||||
type(psb_cspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
@@ -65,6 +65,7 @@ module mld_d_as_smoother
|
||||
! parent type.
|
||||
! class(mld_d_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_dspmat_type), pointer :: pa => null()
|
||||
type(psb_dspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
@@ -65,6 +65,7 @@ module mld_s_as_smoother
|
||||
! parent type.
|
||||
! class(mld_s_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_sspmat_type), pointer :: pa => null()
|
||||
type(psb_sspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
@@ -65,6 +65,7 @@ module mld_z_as_smoother
|
||||
! parent type.
|
||||
! class(mld_z_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_zspmat_type), pointer :: pa => null()
|
||||
type(psb_zspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
Reference in New Issue
Block a user