Merge branch 'unify_aggr_bld' into unify_aggr_bld_richardson to fix REPL instance

configure

@@ -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
+SHELL'
-am__quote'
 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>
+# <http://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/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; }
+am_make=${MAKE-make}
-cat > confinc.mk << 'END'
+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.
+_am_result=none
-echo '.include "confinc.mk" # ignored' > confmf.BSD
+# First try GNU make style include.
-# Other make implementations (GNU, Solaris 10, AIX) do it like this.
+echo "include confinc" > confmf
-echo 'include confinc.mk # ignored' > confmf.GNU
+# Ignore all kinds of additional output from 'make'.
-_am_result=no
+case `$am_make -s -f confmf 2> /dev/null` in #(
-for s in GNU BSD; do
+*the\ am__doit\ target*)
-  { echo "$as_me:$LINENO: ${MAKE-make} -f confmf.$s && cat confinc.out" >&5
+  am__include=include
-   (${MAKE-make} -f confmf.$s && cat confinc.out) >&5 2>&5
+  am__quote=
-   ac_status=$?
+  _am_result=GNU
-   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 ;; #(
-  *) :
-     ;;
 esac
-  if test "$am__include" != "#"; then
+# Now try BSD make style include.
-    _am_result="yes ($s style)"
+if test "$am__include" = "#"; then
-    break
+   echo '.include "confinc"' > confmf
-  fi
+   case `$am_make -s -f confmf 2> /dev/null` in #(
-done
+   *the\ am__doit\ target*)
-rm -f confinc.* confmf.*
+     am__include=.include
-{ $as_echo "$as_me:${as_lineno-$LINENO}: result: ${_am_result}" >&5
+     am__quote="\""
-$as_echo "${_am_result}" >&6; }
+     _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
+  case $CONFIG_FILES in
-  # requiring Autoconf 2.70 or later.
+  *\'*) eval set x "$CONFIG_FILES" ;;
-  case $CONFIG_FILES in #(
+  *)   set x $CONFIG_FILES ;;
-  *\'*) :
+  esac
-    eval set x "$CONFIG_FILES" ;; #(
-  *) :
-    set x $CONFIG_FILES ;; #(
-  *) :
-     ;;
-esac
   shift
-  # Used to flag and report bootstrapping failures.
+  for mf
-  am_rc=0
-  for am_mf
   do
     # Strip MF so we end up with the name of the file.
-    am_mf=`$as_echo "$am_mf" | sed -e 's/:.*$//'`
+    mf=`echo "$mf" | sed -e 's/:.*$//'`
-    # Check whether this is an Automake generated Makefile which includes
+    # Check whether this is an Automake generated Makefile or not.
-    # dependency-tracking related rules and includes.
+    # We used to match only the files named 'Makefile.in', but
-    # Grep'ing the whole file directly is not great: AIX grep has a line
+    # 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 \
+    if sed -n 's,^#.*generated by automake.*,X,p' "$mf" | grep X >/dev/null 2>&1; then
-      || continue
+      dirpart=`$as_dirname -- "$mf" ||
-    am_dirpart=`$as_dirname -- "$am_mf" ||
+$as_expr X"$mf" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
-$as_expr X"$am_mf" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
+	 X"$mf" : 'X\(//\)[^/]' \| \
-	 X"$am_mf" : 'X\(//\)[^/]' \| \
+	 X"$mf" : 'X\(//\)$' \| \
-	 X"$am_mf" : 'X\(//\)$' \| \
+	 X"$mf" : 'X\(/\)' \| . 2>/dev/null ||
-	 X"$am_mf" : 'X\(/\)' \| . 2>/dev/null ||
+$as_echo X"$mf" |
-$as_echo X"$am_mf" |
     sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{
 	    s//\1/
 	    q
@@ -10602,48 +10608,53 @@ $as_echo X"$am_mf" |
 	    q
 	  }
 	  s/.*/./; q'`
-    am_filepart=`$as_basename -- "$am_mf" ||
+    else
-$as_expr X/"$am_mf" : '.*/\([^/][^/]*\)/*$' \| \
+      continue
-	 X"$am_mf" : 'X\(//\)$' \| \
+    fi
-	 X"$am_mf" : 'X\(/\)' \| . 2>/dev/null ||
+    # Extract the definition of DEPDIR, am__include, and am__quote
-$as_echo X/"$am_mf" |
+    # from the Makefile without running 'make'.
-    sed '/^.*\/\([^/][^/]*\)\/*$/{
+    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" \
+      as_dir=$dirpart/$fdir; as_fn_mkdir_p
-      && sed -e '/# am--include-marker/d' "$am_filepart" \
+      # echo "creating $dirpart/$file"
-        | $MAKE -f - am--depfiles" >&5
+      echo '# dummy' > "$dirpart/$file"
-   (cd "$am_dirpart" \
+    done
-      && 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=$?
   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
 }
  ;;
 

mlprec/impl/smoother/mld_c_as_smoother_apply.f90

@@ -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 
+  if (.false.) then 
-    !
+    if (sweeps > 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) 
+      ! 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,&
+      !  Unroll  the first iteration and fold it inside SELECT CASE
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !  this will save one SPMM when INIT=Z, and will be
-      goto 9999
+      !  significant when sweeps=1 (a common case)
-    endif
+      !
+      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 
+      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,&
-    ! Apply multiple sweeps of an AS solver
+             & work=aux,trans=trans_)
-    ! to compute an approximate solution of a linear system.
+        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')             
-    !
-    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('U')
-    !  Unroll  the first iteration and fold it inside SELECT CASE
+        if (.not.present(initu)) then
-    !  this will save one SPMM when INIT=Z, and will be
+          call psb_errpush(psb_err_internal_error_,name,&
-    !  significant when sweeps=1 (a common case)
+               & a_err='missing initu to smoother_apply')
-    !
+          goto 9999
-    call psb_geaxpby(cone,x,czero,tx,desc_data,info)    
+        end if
-    if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
+        call psb_geaxpby(cone,x,czero,ww,desc_data,info)    
-    if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%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 default
-    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')
+             & a_err='wrong  init to smoother_apply')
         goto 9999
-      end if
+      end select
-      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)
+
+      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')
+           & i_err=(/itwo,sweeps,izero,izero,izero/))
-      goto 9999      
+      goto 9999
-    end if
+
-    
+    endif
-    !
-    ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
-    !
-    call psb_geaxpby(alpha,ty,beta,y,desc_data,info) 
-    
-    
   else
-    
+    if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then 
-    info = psb_err_iarg_neg_
+      !
-    call psb_errpush(info,name,&
+      ! Shortcut: in this case there is nothing else to be done. 
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      !
-    goto 9999
+      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

mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90

@@ -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 
+  if (.true.) then
-    !
+    
-    ! Shortcut: in this case there is nothing else to be done. 
+    if (sweeps > 0) then 
-    !
+      !
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) 
+      !
-
+      ! Apply multiple sweeps of an AS solver
-    if (info /= psb_success_) then
+      ! to compute an approximate solution of a linear system.
-      call psb_errpush(psb_err_internal_error_,name,&
+      !
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !
-      goto 9999
+      if (size(wv) < 3) then 
-    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 
         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
+      end if
-      if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
+
+      !
+      ! 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
+    else if (sweeps >= 0) then 
-        !
+      !
-        ! 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
+      ! Apply multiple sweeps of an AS solver
-        ! and Y(j) is the approximate solution at sweep j.
+      ! to compute an approximate solution of a linear system.
-        !
+      !
-        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,&
+      if (size(wv) < 3) then 
-             & work=aux,trans=trans_)
+        call psb_errpush(psb_err_internal_error_,name,&
-
+             & a_err='invalid wv size in smoother_apply')
-        if (info /= psb_success_) exit
+        goto 9999
-
-        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)
+      ! 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) 
+      do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
-    end associate
+           & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
+           & (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
 
-  else
+      if (do_realloc_wv) then
-    
+        call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
-    info = psb_err_iarg_neg_
+        call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-    call psb_errpush(info,name,&
+        call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      end if
-    goto 9999
 
-  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)

mlprec/impl/smoother/mld_c_as_smoother_bld.f90

@@ -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_

mlprec/impl/smoother/mld_c_as_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_c_as_smoother_free.f90

@@ -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

mlprec/impl/smoother/mld_c_jac_smoother_apply.f90

@@ -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 (sweeps > 0) then 
-    !  if .not.sv%is_iterative, there's no need to pass init
+    call psb_geasb(tx,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
+    call psb_geasb(ty,desc_data,info)
-    
+    select case (init_)
-    if (info /= psb_success_) then
+    case('Z') 
-      call psb_errpush(psb_err_internal_error_,&
+
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      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
+    end select
-    
+
-  else if (sweeps >= 0) then 
+    do i=1, sweeps-1
-    if (associated(sm%pa)) then
       !
-      ! 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_geaxpby(cone,x,czero,tx,desc_data,info)
-      call psb_geasb(ty,desc_data,info)
+      call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
-      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)
 
-      !
+      if (info /= psb_success_) exit
-      !  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') 
 
-        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')
+      if (info /= psb_success_) exit
-        call psb_geaxpby(cone,x,czero,tx,desc_data,info)
+    end do
-        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') 
 
-      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
+    if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
-        call psb_errpush(psb_err_internal_error_,name,&
-             & a_err='wrong  init to smoother_apply')
-        goto 9999
-      end select
 
-      do i=1, sweeps-1
+    if (info /= psb_success_) then 
-        !
+      info=psb_err_internal_error_
-        ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
+      call psb_errpush(info,name,& 
-        ! block diagonal part and the remaining part of the local matrix
+           & a_err='subsolve with Jacobi sweeps > 1')
-        ! and Y(j) is the approximate solution at sweep j.
+      goto 9999      
-        !
-        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
     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_

mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90

@@ -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
+  if (sweeps > 0) 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 ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
+    if(sm%checkres) then
-    !  if .not.sv%is_iterative, there's no need to pass init
+      call psb_geall(r,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
+      call psb_geasb(r,desc_data,info)
+      resdenum = psb_genrm2(x,desc_data,info)
+    end if
 
-    if (info /= psb_success_) then
+    associate(tx => wv(1), ty => wv(2))
-      call psb_errpush(psb_err_internal_error_,&
+      select case (init_)
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      case('Z')
-      goto 9999
-    endif
 
-  else if (sweeps >= 0) then
+        call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
-    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')
+      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')
+      case('U')
-          call psb_geaxpby(cone,x,czero,tx,desc_data,info)
+        if (.not.present(initu)) then
-          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
           call psb_errpush(psb_err_internal_error_,name,&
-               & a_err='wrong  init to smoother_apply')
+               & a_err='missing initu 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')
           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
+      if (info /= psb_success_) then
-      !
+        info=psb_err_internal_error_
-      !
-      ! 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_
         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))
 
-        !
+    end associate
-        !  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')
 
-          call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
+  else  if (sweeps == 0) then
-
+    !
-        case('Y')
+    ! 0 sweeps of smoother is the identity operator
-          call psb_geaxpby(cone,x,czero,tx,desc_data,info)
+    !
-          call psb_geaxpby(cone,y,czero,ty,desc_data,info)
+    call psb_geaxpby(alpha,x,beta,y,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
 

mlprec/impl/smoother/mld_c_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_c_jac_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_c_l1_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_d_as_smoother_apply.f90

@@ -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 
+  if (.false.) then 
-    !
+    if (sweeps > 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) 
+      ! 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,&
+      !  Unroll  the first iteration and fold it inside SELECT CASE
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !  this will save one SPMM when INIT=Z, and will be
-      goto 9999
+      !  significant when sweeps=1 (a common case)
-    endif
+      !
+      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 
+      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,&
-    ! Apply multiple sweeps of an AS solver
+             & work=aux,trans=trans_)
-    ! to compute an approximate solution of a linear system.
+        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')             
-    !
-    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('U')
-    !  Unroll  the first iteration and fold it inside SELECT CASE
+        if (.not.present(initu)) then
-    !  this will save one SPMM when INIT=Z, and will be
+          call psb_errpush(psb_err_internal_error_,name,&
-    !  significant when sweeps=1 (a common case)
+               & a_err='missing initu to smoother_apply')
-    !
+          goto 9999
-    call psb_geaxpby(done,x,dzero,tx,desc_data,info)    
+        end if
-    if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
+        call psb_geaxpby(done,x,dzero,ww,desc_data,info)    
-    if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%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 default
-    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')
+             & a_err='wrong  init to smoother_apply')
         goto 9999
-      end if
+      end select
-      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)
+
+      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')
+           & i_err=(/itwo,sweeps,izero,izero,izero/))
-      goto 9999      
+      goto 9999
-    end if
+
-    
+    endif
-    !
-    ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
-    !
-    call psb_geaxpby(alpha,ty,beta,y,desc_data,info) 
-    
-    
   else
-    
+    if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then 
-    info = psb_err_iarg_neg_
+      !
-    call psb_errpush(info,name,&
+      ! Shortcut: in this case there is nothing else to be done. 
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      !
-    goto 9999
+      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

mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90

@@ -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 
+  if (.true.) then
-    !
+    
-    ! Shortcut: in this case there is nothing else to be done. 
+    if (sweeps > 0) then 
-    !
+      !
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) 
+      !
-
+      ! Apply multiple sweeps of an AS solver
-    if (info /= psb_success_) then
+      ! to compute an approximate solution of a linear system.
-      call psb_errpush(psb_err_internal_error_,name,&
+      !
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !
-      goto 9999
+      if (size(wv) < 3) then 
-    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 
         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
+      end if
-      if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
+
+      !
+      ! 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
+    else if (sweeps >= 0) then 
-        !
+      !
-        ! 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
+      ! Apply multiple sweeps of an AS solver
-        ! and Y(j) is the approximate solution at sweep j.
+      ! to compute an approximate solution of a linear system.
-        !
+      !
-        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,&
+      if (size(wv) < 3) then 
-             & work=aux,trans=trans_)
+        call psb_errpush(psb_err_internal_error_,name,&
-
+             & a_err='invalid wv size in smoother_apply')
-        if (info /= psb_success_) exit
+        goto 9999
-
-        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)
+      ! 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) 
+      do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
-    end associate
+           & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
+           & (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
 
-  else
+      if (do_realloc_wv) then
-    
+        call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
-    info = psb_err_iarg_neg_
+        call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-    call psb_errpush(info,name,&
+        call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      end if
-    goto 9999
 
-  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)

mlprec/impl/smoother/mld_d_as_smoother_bld.f90

@@ -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_

mlprec/impl/smoother/mld_d_as_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_d_as_smoother_free.f90

@@ -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

mlprec/impl/smoother/mld_d_jac_smoother_apply.f90

@@ -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 (sweeps > 0) then 
-    !  if .not.sv%is_iterative, there's no need to pass init
+    call psb_geasb(tx,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
+    call psb_geasb(ty,desc_data,info)
-    
+    select case (init_)
-    if (info /= psb_success_) then
+    case('Z') 
-      call psb_errpush(psb_err_internal_error_,&
+
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      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
+    end select
-    
+
-  else if (sweeps >= 0) then 
+    do i=1, sweeps-1
-    if (associated(sm%pa)) then
       !
-      ! 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_geaxpby(done,x,dzero,tx,desc_data,info)
-      call psb_geasb(ty,desc_data,info)
+      call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
-      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)
 
-      !
+      if (info /= psb_success_) exit
-      !  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') 
 
-        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')
+      if (info /= psb_success_) exit
-        call psb_geaxpby(done,x,dzero,tx,desc_data,info)
+    end do
-        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') 
 
-      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
+    if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
-        call psb_errpush(psb_err_internal_error_,name,&
-             & a_err='wrong  init to smoother_apply')
-        goto 9999
-      end select
 
-      do i=1, sweeps-1
+    if (info /= psb_success_) then 
-        !
+      info=psb_err_internal_error_
-        ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
+      call psb_errpush(info,name,& 
-        ! block diagonal part and the remaining part of the local matrix
+           & a_err='subsolve with Jacobi sweeps > 1')
-        ! and Y(j) is the approximate solution at sweep j.
+      goto 9999      
-        !
-        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
     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_

mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90

@@ -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
+  if (sweeps > 0) 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 ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
+    if(sm%checkres) then
-    !  if .not.sv%is_iterative, there's no need to pass init
+      call psb_geall(r,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
+      call psb_geasb(r,desc_data,info)
+      resdenum = psb_genrm2(x,desc_data,info)
+    end if
 
-    if (info /= psb_success_) then
+    associate(tx => wv(1), ty => wv(2))
-      call psb_errpush(psb_err_internal_error_,&
+      select case (init_)
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      case('Z')
-      goto 9999
-    endif
 
-  else if (sweeps >= 0) then
+        call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
-    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')
+      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')
+      case('U')
-          call psb_geaxpby(done,x,dzero,tx,desc_data,info)
+        if (.not.present(initu)) then
-          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
           call psb_errpush(psb_err_internal_error_,name,&
-               & a_err='wrong  init to smoother_apply')
+               & a_err='missing initu 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')
           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
+      if (info /= psb_success_) then
-      !
+        info=psb_err_internal_error_
-      !
-      ! 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_
         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))
 
-        !
+    end associate
-        !  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')
 
-          call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
+  else  if (sweeps == 0) then
-
+    !
-        case('Y')
+    ! 0 sweeps of smoother is the identity operator
-          call psb_geaxpby(done,x,dzero,tx,desc_data,info)
+    !
-          call psb_geaxpby(done,y,dzero,ty,desc_data,info)
+    call psb_geaxpby(alpha,x,beta,y,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
 

mlprec/impl/smoother/mld_d_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_d_jac_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_d_l1_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_s_as_smoother_apply.f90

@@ -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 
+  if (.false.) then 
-    !
+    if (sweeps > 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) 
+      ! 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,&
+      !  Unroll  the first iteration and fold it inside SELECT CASE
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !  this will save one SPMM when INIT=Z, and will be
-      goto 9999
+      !  significant when sweeps=1 (a common case)
-    endif
+      !
+      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 
+      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,&
-    ! Apply multiple sweeps of an AS solver
+             & work=aux,trans=trans_)
-    ! to compute an approximate solution of a linear system.
+        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')             
-    !
-    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('U')
-    !  Unroll  the first iteration and fold it inside SELECT CASE
+        if (.not.present(initu)) then
-    !  this will save one SPMM when INIT=Z, and will be
+          call psb_errpush(psb_err_internal_error_,name,&
-    !  significant when sweeps=1 (a common case)
+               & a_err='missing initu to smoother_apply')
-    !
+          goto 9999
-    call psb_geaxpby(sone,x,szero,tx,desc_data,info)    
+        end if
-    if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
+        call psb_geaxpby(sone,x,szero,ww,desc_data,info)    
-    if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%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 default
-    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')
+             & a_err='wrong  init to smoother_apply')
         goto 9999
-      end if
+      end select
-      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)
+
+      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')
+           & i_err=(/itwo,sweeps,izero,izero,izero/))
-      goto 9999      
+      goto 9999
-    end if
+
-    
+    endif
-    !
-    ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
-    !
-    call psb_geaxpby(alpha,ty,beta,y,desc_data,info) 
-    
-    
   else
-    
+    if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then 
-    info = psb_err_iarg_neg_
+      !
-    call psb_errpush(info,name,&
+      ! Shortcut: in this case there is nothing else to be done. 
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      !
-    goto 9999
+      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

mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90

@@ -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 
+  if (.true.) then
-    !
+    
-    ! Shortcut: in this case there is nothing else to be done. 
+    if (sweeps > 0) then 
-    !
+      !
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) 
+      !
-
+      ! Apply multiple sweeps of an AS solver
-    if (info /= psb_success_) then
+      ! to compute an approximate solution of a linear system.
-      call psb_errpush(psb_err_internal_error_,name,&
+      !
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !
-      goto 9999
+      if (size(wv) < 3) then 
-    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 
         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
+      end if
-      if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
+
+      !
+      ! 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
+    else if (sweeps >= 0) then 
-        !
+      !
-        ! 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
+      ! Apply multiple sweeps of an AS solver
-        ! and Y(j) is the approximate solution at sweep j.
+      ! to compute an approximate solution of a linear system.
-        !
+      !
-        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,&
+      if (size(wv) < 3) then 
-             & work=aux,trans=trans_)
+        call psb_errpush(psb_err_internal_error_,name,&
-
+             & a_err='invalid wv size in smoother_apply')
-        if (info /= psb_success_) exit
+        goto 9999
-
-        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)
+      ! 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) 
+      do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
-    end associate
+           & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
+           & (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
 
-  else
+      if (do_realloc_wv) then
-    
+        call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
-    info = psb_err_iarg_neg_
+        call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-    call psb_errpush(info,name,&
+        call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      end if
-    goto 9999
 
-  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)

mlprec/impl/smoother/mld_s_as_smoother_bld.f90

@@ -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_

mlprec/impl/smoother/mld_s_as_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_s_as_smoother_free.f90

@@ -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

mlprec/impl/smoother/mld_s_jac_smoother_apply.f90

@@ -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 (sweeps > 0) then 
-    !  if .not.sv%is_iterative, there's no need to pass init
+    call psb_geasb(tx,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
+    call psb_geasb(ty,desc_data,info)
-    
+    select case (init_)
-    if (info /= psb_success_) then
+    case('Z') 
-      call psb_errpush(psb_err_internal_error_,&
+
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      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
+    end select
-    
+
-  else if (sweeps >= 0) then 
+    do i=1, sweeps-1
-    if (associated(sm%pa)) then
       !
-      ! 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_geaxpby(sone,x,szero,tx,desc_data,info)
-      call psb_geasb(ty,desc_data,info)
+      call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
-      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)
 
-      !
+      if (info /= psb_success_) exit
-      !  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') 
 
-        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')
+      if (info /= psb_success_) exit
-        call psb_geaxpby(sone,x,szero,tx,desc_data,info)
+    end do
-        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') 
 
-      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
+    if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
-        call psb_errpush(psb_err_internal_error_,name,&
-             & a_err='wrong  init to smoother_apply')
-        goto 9999
-      end select
 
-      do i=1, sweeps-1
+    if (info /= psb_success_) then 
-        !
+      info=psb_err_internal_error_
-        ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
+      call psb_errpush(info,name,& 
-        ! block diagonal part and the remaining part of the local matrix
+           & a_err='subsolve with Jacobi sweeps > 1')
-        ! and Y(j) is the approximate solution at sweep j.
+      goto 9999      
-        !
-        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
     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_

mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90

@@ -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
+  if (sweeps > 0) 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 ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
+    if(sm%checkres) then
-    !  if .not.sv%is_iterative, there's no need to pass init
+      call psb_geall(r,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
+      call psb_geasb(r,desc_data,info)
+      resdenum = psb_genrm2(x,desc_data,info)
+    end if
 
-    if (info /= psb_success_) then
+    associate(tx => wv(1), ty => wv(2))
-      call psb_errpush(psb_err_internal_error_,&
+      select case (init_)
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      case('Z')
-      goto 9999
-    endif
 
-  else if (sweeps >= 0) then
+        call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
-    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')
+      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')
+      case('U')
-          call psb_geaxpby(sone,x,szero,tx,desc_data,info)
+        if (.not.present(initu)) then
-          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
           call psb_errpush(psb_err_internal_error_,name,&
-               & a_err='wrong  init to smoother_apply')
+               & a_err='missing initu 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')
           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
+      if (info /= psb_success_) then
-      !
+        info=psb_err_internal_error_
-      !
-      ! 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_
         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))
 
-        !
+    end associate
-        !  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')
 
-          call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
+  else  if (sweeps == 0) then
-
+    !
-        case('Y')
+    ! 0 sweeps of smoother is the identity operator
-          call psb_geaxpby(sone,x,szero,tx,desc_data,info)
+    !
-          call psb_geaxpby(sone,y,szero,ty,desc_data,info)
+    call psb_geaxpby(alpha,x,beta,y,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
 

mlprec/impl/smoother/mld_s_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_s_jac_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_s_l1_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_z_as_smoother_apply.f90

@@ -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 
+  if (.false.) then 
-    !
+    if (sweeps > 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) 
+      ! 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,&
+      !  Unroll  the first iteration and fold it inside SELECT CASE
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !  this will save one SPMM when INIT=Z, and will be
-      goto 9999
+      !  significant when sweeps=1 (a common case)
-    endif
+      !
+      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 
+      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,&
-    ! Apply multiple sweeps of an AS solver
+             & work=aux,trans=trans_)
-    ! to compute an approximate solution of a linear system.
+        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')             
-    !
-    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('U')
-    !  Unroll  the first iteration and fold it inside SELECT CASE
+        if (.not.present(initu)) then
-    !  this will save one SPMM when INIT=Z, and will be
+          call psb_errpush(psb_err_internal_error_,name,&
-    !  significant when sweeps=1 (a common case)
+               & a_err='missing initu to smoother_apply')
-    !
+          goto 9999
-    call psb_geaxpby(zone,x,zzero,tx,desc_data,info)    
+        end if
-    if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
+        call psb_geaxpby(zone,x,zzero,ww,desc_data,info)    
-    if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%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 default
-    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')
+             & a_err='wrong  init to smoother_apply')
         goto 9999
-      end if
+      end select
-      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)
+
+      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')
+           & i_err=(/itwo,sweeps,izero,izero,izero/))
-      goto 9999      
+      goto 9999
-    end if
+
-    
+    endif
-    !
-    ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
-    !
-    call psb_geaxpby(alpha,ty,beta,y,desc_data,info) 
-    
-    
   else
-    
+    if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then 
-    info = psb_err_iarg_neg_
+      !
-    call psb_errpush(info,name,&
+      ! Shortcut: in this case there is nothing else to be done. 
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      !
-    goto 9999
+      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

mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90

@@ -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 
+  if (.true.) then
-    !
+    
-    ! Shortcut: in this case there is nothing else to be done. 
+    if (sweeps > 0) then 
-    !
+      !
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) 
+      !
-
+      ! Apply multiple sweeps of an AS solver
-    if (info /= psb_success_) then
+      ! to compute an approximate solution of a linear system.
-      call psb_errpush(psb_err_internal_error_,name,&
+      !
-           & a_err='Error in sub_aply Jacobi Sweeps = 1')
+      !
-      goto 9999
+      if (size(wv) < 3) then 
-    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 
         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
+      end if
-      if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
+
+      !
+      ! 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
+    else if (sweeps >= 0) then 
-        !
+      !
-        ! 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
+      ! Apply multiple sweeps of an AS solver
-        ! and Y(j) is the approximate solution at sweep j.
+      ! to compute an approximate solution of a linear system.
-        !
+      !
-        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,&
+      if (size(wv) < 3) then 
-             & work=aux,trans=trans_)
+        call psb_errpush(psb_err_internal_error_,name,&
-
+             & a_err='invalid wv size in smoother_apply')
-        if (info /= psb_success_) exit
+        goto 9999
-
-        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)
+      ! 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) 
+      do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
-    end associate
+           & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
+           & (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
 
-  else
+      if (do_realloc_wv) then
-    
+        call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
-    info = psb_err_iarg_neg_
+        call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-    call psb_errpush(info,name,&
+        call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
-         & i_err=(/itwo,sweeps,izero,izero,izero/))
+      end if
-    goto 9999
 
-  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)

mlprec/impl/smoother/mld_z_as_smoother_bld.f90

@@ -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_

mlprec/impl/smoother/mld_z_as_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_z_as_smoother_free.f90

@@ -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

mlprec/impl/smoother/mld_z_jac_smoother_apply.f90

@@ -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 (sweeps > 0) then 
-    !  if .not.sv%is_iterative, there's no need to pass init
+    call psb_geasb(tx,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
+    call psb_geasb(ty,desc_data,info)
-    
+    select case (init_)
-    if (info /= psb_success_) then
+    case('Z') 
-      call psb_errpush(psb_err_internal_error_,&
+
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      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
+    end select
-    
+
-  else if (sweeps >= 0) then 
+    do i=1, sweeps-1
-    if (associated(sm%pa)) then
       !
-      ! 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_geaxpby(zone,x,zzero,tx,desc_data,info)
-      call psb_geasb(ty,desc_data,info)
+      call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
-      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)
 
-      !
+      if (info /= psb_success_) exit
-      !  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') 
 
-        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')
+      if (info /= psb_success_) exit
-        call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
+    end do
-        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') 
 
-      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
+    if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
-        call psb_errpush(psb_err_internal_error_,name,&
-             & a_err='wrong  init to smoother_apply')
-        goto 9999
-      end select
 
-      do i=1, sweeps-1
+    if (info /= psb_success_) then 
-        !
+      info=psb_err_internal_error_
-        ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
+      call psb_errpush(info,name,& 
-        ! block diagonal part and the remaining part of the local matrix
+           & a_err='subsolve with Jacobi sweeps > 1')
-        ! and Y(j) is the approximate solution at sweep j.
+      goto 9999      
-        !
-        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
     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_

mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90

@@ -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
+  if (sweeps > 0) 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 ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
+    if(sm%checkres) then
-    !  if .not.sv%is_iterative, there's no need to pass init
+      call psb_geall(r,desc_data,info)
-    call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
+      call psb_geasb(r,desc_data,info)
+      resdenum = psb_genrm2(x,desc_data,info)
+    end if
 
-    if (info /= psb_success_) then
+    associate(tx => wv(1), ty => wv(2))
-      call psb_errpush(psb_err_internal_error_,&
+      select case (init_)
-           & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
+      case('Z')
-      goto 9999
-    endif
 
-  else if (sweeps >= 0) then
+        call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
-    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')
+      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')
+      case('U')
-          call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
+        if (.not.present(initu)) then
-          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
           call psb_errpush(psb_err_internal_error_,name,&
-               & a_err='wrong  init to smoother_apply')
+               & a_err='missing initu 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')
           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
+      if (info /= psb_success_) then
-      !
+        info=psb_err_internal_error_
-      !
-      ! 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_
         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))
 
-        !
+    end associate
-        !  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')
 
-          call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
+  else  if (sweeps == 0) then
-
+    !
-        case('Y')
+    ! 0 sweeps of smoother is the identity operator
-          call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
+    !
-          call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
+    call psb_geaxpby(alpha,x,beta,y,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
 

mlprec/impl/smoother/mld_z_jac_smoother_bld.f90

@@ -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)

mlprec/impl/smoother/mld_z_jac_smoother_clone.f90

@@ -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_

mlprec/impl/smoother/mld_z_l1_jac_smoother_bld.f90

@@ -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)

mlprec/mld_c_as_smoother.f90

@@ -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

mlprec/mld_d_as_smoother.f90

@@ -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

mlprec/mld_s_as_smoother.f90

@@ -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

mlprec/mld_z_as_smoother.f90

@@ -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