Fix AMG cbind interface

amgprec/amg_c_prec_type.f90

@@ -636,6 +636,11 @@ contains
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free(info)
+        if (psb_errstatus_fatal()) then
+          info=psb_err_internal_error_
+          call psb_errpush(info,name)
+          goto 9999
+        end if
       end do
       deallocate(prec%precv,stat=info)
     end if

amgprec/amg_d_prec_type.f90

@@ -636,6 +636,11 @@ contains
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free(info)
+        if (psb_errstatus_fatal()) then
+          info=psb_err_internal_error_
+          call psb_errpush(info,name)
+          goto 9999
+        end if
       end do
       deallocate(prec%precv,stat=info)
     end if

amgprec/amg_s_prec_type.f90

@@ -636,6 +636,11 @@ contains
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free(info)
+        if (psb_errstatus_fatal()) then
+          info=psb_err_internal_error_
+          call psb_errpush(info,name)
+          goto 9999
+        end if
       end do
       deallocate(prec%precv,stat=info)
     end if

amgprec/amg_z_prec_type.f90

@@ -636,6 +636,11 @@ contains
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free(info)
+        if (psb_errstatus_fatal()) then
+          info=psb_err_internal_error_
+          call psb_errpush(info,name)
+          goto 9999
+        end if
       end do
       deallocate(prec%precv,stat=info)
     end if

amgprec/impl/amg_c_hierarchy_bld.F90

@@ -366,14 +366,10 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
 
     if (i>2) then
       if (sizeratio < mnaggratio) then
-        if (sizeratio > 1) then
-          newsz = i
-        else
-          !
-          ! We are not gaining 
-          !
-          newsz = i-1
-        end if
+        !
+        ! We are not gaining 
+        !
+        newsz = i-1
       end if
 
       if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then 

amgprec/impl/amg_cfile_prec_descr.f90

@@ -88,6 +88,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix)
   logical             :: is_symgs
   character(len=20), parameter :: name='amg_file_prec_descr'
   integer(psb_ipk_)  :: iout_, root_, verbosity_
+  integer(psb_lpk_)  :: gl_nrows,gl_ncols,gl_nzeros
   character(1024)    :: prefix_
 
   info = psb_success_
@@ -122,6 +123,12 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix)
     if (root_ == -1) root_ = me
 
     if (verbosity_ >=0) then 
+      gl_nrows = prec%precv(1)%base_a%get_nrows()
+      gl_ncols = prec%precv(1)%base_a%get_ncols()
+      gl_nzeros = prec%precv(1)%base_a%get_nzeros()
+      call psb_sum(ctxt,gl_nrows)
+      call psb_sum(ctxt,gl_ncols)
+      call psb_sum(ctxt,gl_nzeros)
       !
       ! The preconditioner description is printed by processor psb_root_.
       ! This agrees with the fact that all the parameters defining the
@@ -141,7 +148,11 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix)
 
         write(iout_,*) 
         write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
-
+        write(iout_,*) 
+        write(iout_,*) trim(prefix_),' Base matrix      : ',&
+             & gl_nrows, gl_ncols, gl_nzeros
+        write(iout_,*)
+        
         if (nlev == 1) then
           !
           ! Here we have a gigantic kludge just to handle Symmetrized Gauss-Seidel.

amgprec/impl/amg_d_hierarchy_bld.F90

@@ -366,14 +366,10 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
 
     if (i>2) then
       if (sizeratio < mnaggratio) then
-        if (sizeratio > 1) then
-          newsz = i
-        else
-          !
-          ! We are not gaining 
-          !
-          newsz = i-1
-        end if
+        !
+        ! We are not gaining 
+        !
+        newsz = i-1
       end if
 
       if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then 
@@ -434,7 +430,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
            & ilaggr,nlaggr,op_prol,info)
       if (do_timings) call psb_toc(idx_matasb)      
     end if
-
     if (info /= psb_success_) then 
       write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i
       call psb_errpush(psb_err_internal_error_,name,&

amgprec/impl/amg_dfile_prec_descr.f90

@@ -88,6 +88,7 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix)
   logical             :: is_symgs
   character(len=20), parameter :: name='amg_file_prec_descr'
   integer(psb_ipk_)  :: iout_, root_, verbosity_
+  integer(psb_lpk_)  :: gl_nrows,gl_ncols,gl_nzeros
   character(1024)    :: prefix_
 
   info = psb_success_
@@ -122,6 +123,12 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix)
     if (root_ == -1) root_ = me
 
     if (verbosity_ >=0) then 
+      gl_nrows = prec%precv(1)%base_a%get_nrows()
+      gl_ncols = prec%precv(1)%base_a%get_ncols()
+      gl_nzeros = prec%precv(1)%base_a%get_nzeros()
+      call psb_sum(ctxt,gl_nrows)
+      call psb_sum(ctxt,gl_ncols)
+      call psb_sum(ctxt,gl_nzeros)
       !
       ! The preconditioner description is printed by processor psb_root_.
       ! This agrees with the fact that all the parameters defining the
@@ -141,7 +148,11 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix)
 
         write(iout_,*) 
         write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
-
+        write(iout_,*) 
+        write(iout_,*) trim(prefix_),' Base matrix      : ',&
+             & gl_nrows, gl_ncols, gl_nzeros
+        write(iout_,*)
+        
         if (nlev == 1) then
           !
           ! Here we have a gigantic kludge just to handle Symmetrized Gauss-Seidel.

amgprec/impl/amg_s_hierarchy_bld.F90

@@ -366,14 +366,10 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
 
     if (i>2) then
       if (sizeratio < mnaggratio) then
-        if (sizeratio > 1) then
-          newsz = i
-        else
-          !
-          ! We are not gaining 
-          !
-          newsz = i-1
-        end if
+        !
+        ! We are not gaining 
+        !
+        newsz = i-1
       end if
 
       if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then 

amgprec/impl/amg_sfile_prec_descr.f90

@@ -88,6 +88,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix)
   logical             :: is_symgs
   character(len=20), parameter :: name='amg_file_prec_descr'
   integer(psb_ipk_)  :: iout_, root_, verbosity_
+  integer(psb_lpk_)  :: gl_nrows,gl_ncols,gl_nzeros
   character(1024)    :: prefix_
 
   info = psb_success_
@@ -122,6 +123,12 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix)
     if (root_ == -1) root_ = me
 
     if (verbosity_ >=0) then 
+      gl_nrows = prec%precv(1)%base_a%get_nrows()
+      gl_ncols = prec%precv(1)%base_a%get_ncols()
+      gl_nzeros = prec%precv(1)%base_a%get_nzeros()
+      call psb_sum(ctxt,gl_nrows)
+      call psb_sum(ctxt,gl_ncols)
+      call psb_sum(ctxt,gl_nzeros)
       !
       ! The preconditioner description is printed by processor psb_root_.
       ! This agrees with the fact that all the parameters defining the
@@ -141,7 +148,11 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix)
 
         write(iout_,*) 
         write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
-
+        write(iout_,*) 
+        write(iout_,*) trim(prefix_),' Base matrix      : ',&
+             & gl_nrows, gl_ncols, gl_nzeros
+        write(iout_,*)
+        
         if (nlev == 1) then
           !
           ! Here we have a gigantic kludge just to handle Symmetrized Gauss-Seidel.

amgprec/impl/amg_z_hierarchy_bld.F90

@@ -366,14 +366,10 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
 
     if (i>2) then
       if (sizeratio < mnaggratio) then
-        if (sizeratio > 1) then
-          newsz = i
-        else
-          !
-          ! We are not gaining 
-          !
-          newsz = i-1
-        end if
+        !
+        ! We are not gaining 
+        !
+        newsz = i-1
       end if
 
       if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then 

amgprec/impl/amg_zfile_prec_descr.f90

@@ -88,6 +88,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix)
   logical             :: is_symgs
   character(len=20), parameter :: name='amg_file_prec_descr'
   integer(psb_ipk_)  :: iout_, root_, verbosity_
+  integer(psb_lpk_)  :: gl_nrows,gl_ncols,gl_nzeros
   character(1024)    :: prefix_
 
   info = psb_success_
@@ -122,6 +123,12 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix)
     if (root_ == -1) root_ = me
 
     if (verbosity_ >=0) then 
+      gl_nrows = prec%precv(1)%base_a%get_nrows()
+      gl_ncols = prec%precv(1)%base_a%get_ncols()
+      gl_nzeros = prec%precv(1)%base_a%get_nzeros()
+      call psb_sum(ctxt,gl_nrows)
+      call psb_sum(ctxt,gl_ncols)
+      call psb_sum(ctxt,gl_nzeros)
       !
       ! The preconditioner description is printed by processor psb_root_.
       ! This agrees with the fact that all the parameters defining the
@@ -141,7 +148,11 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix)
 
         write(iout_,*) 
         write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
-
+        write(iout_,*) 
+        write(iout_,*) trim(prefix_),' Base matrix      : ',&
+             & gl_nrows, gl_ncols, gl_nzeros
+        write(iout_,*)
+        
         if (nlev == 1) then
           !
           ! Here we have a gigantic kludge just to handle Symmetrized Gauss-Seidel.

amgprec/impl/solver/amg_c_diag_solver_bld.f90

@@ -91,6 +91,10 @@ subroutine amg_c_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
@@ -172,6 +176,10 @@ subroutine amg_c_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_c_jac_solver_bld.f90

@@ -100,6 +100,10 @@ subroutine amg_c_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_c_l1_jac_solver_bld.f90

@@ -103,6 +103,10 @@ subroutine amg_c_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_c_mumps_solver_bld.F90

@@ -80,24 +80,19 @@ subroutine c_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   call psb_info(ctxt, iam, np)      
   if (sv%ipar(1) == amg_local_solver_ ) then
     call psb_init(ctxt1,np=1,basectxt=ctxt,ids=(/iam/))
-    icomm = psb_get_mpi_comm(ctxt1)
-    allocate(sv%local_ctxt,stat=info)
-    sv%local_ctxt = ctxt1
-    !write(*,*)iam,'mumps_bld: local +++++>',icomm,sv%local_ctxt
-    call psb_info(ctxt1, me, np)
-    npr  = np
   else if (sv%ipar(1) == amg_global_solver_ ) then
-    icomm = psb_get_mpi_comm(ctxt)
-    !write(*,*)iam,'mumps_bld: global +++++>',icomm,ctxt
-    call psb_info(ctxt, iam, np)
-    me = iam 
-    npr  = np
+    call psb_init(ctxt1,basectxt=ctxt)
   else
     info = psb_err_internal_error_
     call psb_errpush(info,name,& 
          & a_err='Invalid local/global solver in MUMPS')
     goto 9999
   end if
+  !allocate(sv%local_ctxt,stat=info)
+  icomm = ctxt1%get_mpic()
+  sv%local_ctxt = ctxt1
+  !write(*,*)iam,'mumps_bld: local +++++>',icomm!,sv%local_ctxt%ctxt
+  call psb_info(ctxt1, me, npr)
   npc  = 1
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),' start'
@@ -115,8 +110,8 @@ subroutine c_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     end if
   end if
 
+  call mpi_comm_dup(icomm,sv%id%comm,info)
 
-  sv%id%comm    =  icomm
   sv%id%job     = -1
   sv%id%par     =  1
   if (sv%ipar(3) == 2) then

amgprec/impl/solver/amg_d_diag_solver_bld.f90

@@ -91,6 +91,10 @@ subroutine amg_d_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
@@ -172,6 +176,10 @@ subroutine amg_d_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_d_jac_solver_bld.f90

@@ -100,6 +100,10 @@ subroutine amg_d_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_d_l1_jac_solver_bld.f90

@@ -103,6 +103,10 @@ subroutine amg_d_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_d_mumps_solver_bld.F90

@@ -80,24 +80,19 @@ subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   call psb_info(ctxt, iam, np)      
   if (sv%ipar(1) == amg_local_solver_ ) then
     call psb_init(ctxt1,np=1,basectxt=ctxt,ids=(/iam/))
-    icomm = psb_get_mpi_comm(ctxt1)
-    allocate(sv%local_ctxt,stat=info)
-    sv%local_ctxt = ctxt1
-    !write(*,*)iam,'mumps_bld: local +++++>',icomm,sv%local_ctxt
-    call psb_info(ctxt1, me, np)
-    npr  = np
   else if (sv%ipar(1) == amg_global_solver_ ) then
-    icomm = psb_get_mpi_comm(ctxt)
-    !write(*,*)iam,'mumps_bld: global +++++>',icomm,ctxt
-    call psb_info(ctxt, iam, np)
-    me = iam 
-    npr  = np
+    call psb_init(ctxt1,basectxt=ctxt)
   else
     info = psb_err_internal_error_
     call psb_errpush(info,name,& 
          & a_err='Invalid local/global solver in MUMPS')
     goto 9999
   end if
+  !allocate(sv%local_ctxt,stat=info)
+  icomm = ctxt1%get_mpic()
+  sv%local_ctxt = ctxt1
+  !write(*,*)iam,'mumps_bld: local +++++>',icomm!,sv%local_ctxt%ctxt
+  call psb_info(ctxt1, me, npr)
   npc  = 1
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),' start'
@@ -115,8 +110,8 @@ subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     end if
   end if
 
+  call mpi_comm_dup(icomm,sv%id%comm,info)
 
-  sv%id%comm    =  icomm
   sv%id%job     = -1
   sv%id%par     =  1
   if (sv%ipar(3) == 2) then

amgprec/impl/solver/amg_s_diag_solver_bld.f90

@@ -91,6 +91,10 @@ subroutine amg_s_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
@@ -172,6 +176,10 @@ subroutine amg_s_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_s_jac_solver_bld.f90

@@ -100,6 +100,10 @@ subroutine amg_s_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_s_l1_jac_solver_bld.f90

@@ -103,6 +103,10 @@ subroutine amg_s_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_s_mumps_solver_bld.F90

@@ -80,24 +80,19 @@ subroutine s_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   call psb_info(ctxt, iam, np)      
   if (sv%ipar(1) == amg_local_solver_ ) then
     call psb_init(ctxt1,np=1,basectxt=ctxt,ids=(/iam/))
-    icomm = psb_get_mpi_comm(ctxt1)
-    allocate(sv%local_ctxt,stat=info)
-    sv%local_ctxt = ctxt1
-    !write(*,*)iam,'mumps_bld: local +++++>',icomm,sv%local_ctxt
-    call psb_info(ctxt1, me, np)
-    npr  = np
   else if (sv%ipar(1) == amg_global_solver_ ) then
-    icomm = psb_get_mpi_comm(ctxt)
-    !write(*,*)iam,'mumps_bld: global +++++>',icomm,ctxt
-    call psb_info(ctxt, iam, np)
-    me = iam 
-    npr  = np
+    call psb_init(ctxt1,basectxt=ctxt)
   else
     info = psb_err_internal_error_
     call psb_errpush(info,name,& 
          & a_err='Invalid local/global solver in MUMPS')
     goto 9999
   end if
+  !allocate(sv%local_ctxt,stat=info)
+  icomm = ctxt1%get_mpic()
+  sv%local_ctxt = ctxt1
+  !write(*,*)iam,'mumps_bld: local +++++>',icomm!,sv%local_ctxt%ctxt
+  call psb_info(ctxt1, me, npr)
   npc  = 1
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),' start'
@@ -115,8 +110,8 @@ subroutine s_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     end if
   end if
 
+  call mpi_comm_dup(icomm,sv%id%comm,info)
 
-  sv%id%comm    =  icomm
   sv%id%job     = -1
   sv%id%par     =  1
   if (sv%ipar(3) == 2) then

amgprec/impl/solver/amg_z_diag_solver_bld.f90

@@ -91,6 +91,10 @@ subroutine amg_z_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
@@ -172,6 +176,10 @@ subroutine amg_z_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_z_jac_solver_bld.f90

@@ -100,6 +100,10 @@ subroutine amg_z_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_z_l1_jac_solver_bld.f90

@@ -103,6 +103,10 @@ subroutine amg_z_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
   allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)

amgprec/impl/solver/amg_z_mumps_solver_bld.F90

@@ -80,24 +80,19 @@ subroutine z_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   call psb_info(ctxt, iam, np)      
   if (sv%ipar(1) == amg_local_solver_ ) then
     call psb_init(ctxt1,np=1,basectxt=ctxt,ids=(/iam/))
-    icomm = psb_get_mpi_comm(ctxt1)
-    allocate(sv%local_ctxt,stat=info)
-    sv%local_ctxt = ctxt1
-    !write(*,*)iam,'mumps_bld: local +++++>',icomm,sv%local_ctxt
-    call psb_info(ctxt1, me, np)
-    npr  = np
   else if (sv%ipar(1) == amg_global_solver_ ) then
-    icomm = psb_get_mpi_comm(ctxt)
-    !write(*,*)iam,'mumps_bld: global +++++>',icomm,ctxt
-    call psb_info(ctxt, iam, np)
-    me = iam 
-    npr  = np
+    call psb_init(ctxt1,basectxt=ctxt)
   else
     info = psb_err_internal_error_
     call psb_errpush(info,name,& 
          & a_err='Invalid local/global solver in MUMPS')
     goto 9999
   end if
+  !allocate(sv%local_ctxt,stat=info)
+  icomm = ctxt1%get_mpic()
+  sv%local_ctxt = ctxt1
+  !write(*,*)iam,'mumps_bld: local +++++>',icomm!,sv%local_ctxt%ctxt
+  call psb_info(ctxt1, me, npr)
   npc  = 1
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),' start'
@@ -115,8 +110,8 @@ subroutine z_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     end if
   end if
 
+  call mpi_comm_dup(icomm,sv%id%comm,info)
 
-  sv%id%comm    =  icomm
   sv%id%job     = -1
   sv%id%par     =  1
   if (sv%ipar(3) == 2) then

cbind/amgprec/amg_c_dprec.h

@@ -36,8 +36,9 @@ extern "C" {
   psb_i_t amg_c_dallocate_wrk(amg_c_dprec *ph, const char *chfmt);
 
   psb_i_t amg_c_dkrylov(const char *method, psb_c_dspmat *ah, amg_c_dprec *ph, 
-		  psb_c_dvector *bh, psb_c_dvector *xh,
-		  psb_c_descriptor *cdh, psb_c_SolverOptions *opt);
+			psb_c_dvector *bh, psb_c_dvector *xh,
+			psb_c_descriptor *cdh, psb_c_dvector *s1,
+			psb_c_dvector *s2, psb_c_SolverOptions *opt);
 
 
 #ifdef __cplusplus

cbind/amgprec/amg_c_zprec.h

@@ -40,7 +40,8 @@ extern "C"
 
   psb_i_t amg_c_zkrylov(const char *method, psb_c_zspmat *ah, amg_c_zprec *ph,
                         psb_c_zvector *bh, psb_c_zvector *xh,
-                        psb_c_descriptor *cdh, psb_c_SolverOptions *opt);
+                        psb_c_descriptor *cdh, , psb_c_zvector *s1,
+			psb_c_zvector *s2, psb_c_SolverOptions *opt);
 
 #ifdef __cplusplus
 }

cbind/amgprec/amg_dprec_cbind_mod.F90

@@ -3,6 +3,7 @@ module amg_dprec_cbind_mod
   use iso_c_binding
   use amg_prec_mod
   use psb_base_cbind_mod
+  use psb_dlinsolve_cbind_mod
 
   type, bind(c) :: amg_c_dprec
     type(c_ptr) :: item = c_null_ptr
@@ -172,15 +173,17 @@ contains
   end function amg_c_dprecbld
 
   function  amg_c_dhierarchy_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_dprec_type), pointer  :: precp
     type(psb_dspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret, act
+
 
     res = -1
 
@@ -204,11 +207,16 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
 
     return
   end function amg_c_dhierarchy_build
 
   function  amg_c_dsmoothers_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
@@ -217,7 +225,7 @@ contains
     type(psb_dspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
-    integer(psb_ipk_)     :: iret
+    integer(psb_ipk_)     :: iret, act
 
     res = -1
 
@@ -241,8 +249,10 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
-
-    return
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
   end function amg_c_dsmoothers_build
 
   function  amg_c_dsmoothers_build_opt(ah,cdh,ph,afmt,cdfmt) bind(c) result(res)
@@ -360,7 +370,7 @@ contains
   end function amg_c_dsmoothers_build_opt
 
   function  amg_c_dkrylov(methd,&
-       & ah,ph,bh,xh,cdh,options) bind(c) result(res)
+       & ah,ph,bh,xh,cdh,s1,s2,options) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
     use psb_linsolve_mod
@@ -368,20 +378,21 @@ contains
     use psb_dlinsolve_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1,s2
     character(c_char)       :: methd(*)
     type(solveroptions)     :: options
 
     res= amg_c_dkrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh,  &
          & itmax=options%itmax, iter=options%iter,&
          & itrace=options%itrace, istop=options%istop,&
-         & irst=options%irst, err=options%err)
+         & irst=options%irst, err=options%err, s1=s1,s2=s2)
 
   end function amg_c_dkrylov
 
 
   function  amg_c_dkrylov_opt(methd,&
-       & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
+       & ah,ph,bh,xh,eps,cdh,itmax,iter,err,&
+       & itrace,irst,istop,s1,s2) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
     use psb_linsolve_mod
@@ -389,16 +400,17 @@ contains
     use psb_prec_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1, s2
     integer(psb_c_ipk_), value :: itmax,itrace,irst,istop
     real(c_double), value :: eps
     integer(psb_c_ipk_)        :: iter
     real(c_double)        :: err
     character(c_char)       :: methd(*)
+
     type(psb_desc_type), pointer   :: descp
     type(psb_dspmat_type), pointer :: ap
     type(amg_dprec_type), pointer  :: precp
-    type(psb_d_vect_type), pointer :: xp, bp
+    type(psb_d_vect_type), pointer :: xp, bp, s1p, s2p
 
     integer(psb_ipk_)     :: iret,fitmax,fitrace,first,fistop,fiter
     character(len=20)     :: fmethd
@@ -431,6 +443,16 @@ contains
       return
     end if
 
+    if (c_associated(s1%item)) then
+      call c_f_pointer(s1%item,s1p)
+    else
+      nullify(s1p)
+    end if
+    if (c_associated(s2%item)) then
+      call c_f_pointer(s2%item,s2p)
+    else
+      nullify(s2p)
+    end if
 
     call psb_stringc2f(methd,fmethd)
     feps    = eps
@@ -439,10 +461,28 @@ contains
     first   = irst
     fistop  = istop
 
-    call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
-         & descp, iret,&
-         & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
-         & irst=first, err=ferr)
+    if (associated(s1p).and.associated(s2p)) then
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr,s1=s1p,s2=s2p)
+    else if  (associated(s1p)) then
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr,s1=s1p)
+    else  if (associated(s2p)) then
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr,s2=s2p)
+    else
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr)
+    end if
+    
     iter = fiter
     err  = ferr
     res = min(iret,0)
@@ -497,7 +537,7 @@ contains
     res = AMGC_ERR_FILTER(info)
     AMGC_ERR_HANDLE(res)
     return
-end function amg_c_dprecapply
+  end function amg_c_dprecapply
 
   function amg_c_dprecapply_opt(ph,bc,xc,cdh,ctrans) bind(c,name="amg_c_dprecapply_opt") result(res)
     use psb_base_mod
@@ -553,7 +593,7 @@ end function amg_c_dprecapply
     res = AMGC_ERR_FILTER(info)
     AMGC_ERR_HANDLE(res)
     return
-end function amg_c_dprecapply_opt
+  end function amg_c_dprecapply_opt
 
   function  amg_c_dprecfree(ph) bind(c) result(res)
     implicit none

cbind/amgprec/amg_zprec_cbind_mod.F90

@@ -140,11 +140,11 @@ contains
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_zprec_type), pointer  :: precp
     type(psb_zspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret
 
     res = -1
 
@@ -173,15 +173,17 @@ contains
   end function amg_c_zprecbld
 
   function  amg_c_zhierarchy_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_zprec_type), pointer  :: precp
     type(psb_zspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret, act
+
 
     res = -1
 
@@ -205,20 +207,25 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
-
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
+      
     return
   end function amg_c_zhierarchy_build
 
   function  amg_c_zsmoothers_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_zprec_type), pointer  :: precp
     type(psb_zspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret, act
 
     res = -1
 
@@ -242,11 +249,15 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
 
     return
   end function amg_c_zsmoothers_build
 
-  function  amg_c_zsmoothers_build_format(ah,cdh,ph,afmt,cdfmt) bind(c) result(res)
+  function  amg_c_zsmoothers_build_opt(ah,cdh,ph,afmt,cdfmt) bind(c) result(res)
 #if defined (PSB_HAVE_CUDA)
     use psb_cuda_mod
 #endif
@@ -352,34 +363,40 @@ contains
     AMGC_ERR_HANDLE(res)
 
     return
-  end function amg_c_zsmoothers_build_format
+  end function amg_c_zsmoothers_build_opt
 
   function  amg_c_zkrylov(methd,&
-       & ah,ph,bh,xh,cdh,options) bind(c) result(res)
+       & ah,ph,bh,xh,cdh,s1,s2,options) bind(c) result(res)
+    use psb_base_mod
+    use psb_prec_mod
     use psb_linsolve_mod
     use psb_prec_cbind_mod
     use psb_zlinsolve_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1,s2
     character(c_char)       :: methd(*)
     type(solveroptions)     :: options
 
     res= amg_c_zkrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh,  &
          & itmax=options%itmax, iter=options%iter,&
          & itrace=options%itrace, istop=options%istop,&
-         & irst=options%irst, err=options%err)
+         & irst=options%irst, err=options%err, s1=s1,s2=s2)
 
   end function amg_c_zkrylov
 
 
   function  amg_c_zkrylov_opt(methd,&
-       & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
+       & ah,ph,bh,xh,eps,cdh,itmax,iter,err,&
+       & itrace,irst,istop,s1,s2) bind(c) result(res)
+    use psb_base_mod
+    use psb_prec_mod
     use psb_linsolve_mod
     use psb_objhandle_mod
+    use psb_prec_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1, s2
     integer(psb_c_ipk_), value :: itmax,itrace,irst,istop
     real(c_double), value :: eps
     integer(psb_c_ipk_)        :: iter
@@ -388,7 +405,7 @@ contains
     type(psb_desc_type), pointer   :: descp
     type(psb_zspmat_type), pointer :: ap
     type(amg_zprec_type), pointer  :: precp
-    type(psb_z_vect_type), pointer :: xp, bp
+    type(psb_z_vect_type), pointer :: xp, bp, s1p, s2p
 
     integer(psb_ipk_)     :: iret,fitmax,fitrace,first,fistop,fiter
     character(len=20)     :: fmethd
@@ -420,6 +437,16 @@ contains
     else
       return
     end if
+    if (c_associated(s1%item)) then
+      call c_f_pointer(s1%item,s1p)
+    else
+      nullify(s1p)
+    end if
+    if (c_associated(s2%item)) then
+      call c_f_pointer(s2%item,s2p)
+    else
+      nullify(s2p)
+    end if
 
 
     call psb_stringc2f(methd,fmethd)
@@ -429,10 +456,28 @@ contains
     first   = irst
     fistop  = istop
 
-    call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
-         & descp, iret,&
-         & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
-         & irst=first, err=ferr)
+    if (associated(s1p).and.associated(s2p)) then
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr,s1=s1p,s2=s2p)
+    else if  (associated(s1p)) then
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr,s1=s1p)
+    else  if (associated(s2p)) then
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr,s2=s2p)
+    else
+      call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
+           & descp, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr)
+    end if
+
     iter = fiter
     err  = ferr
     res = min(iret,0)

configure

@@ -3451,7 +3451,7 @@ fi
 
 { printf "%s\n" "$as_me:${as_lineno-$LINENO}: Loaded $pac_cv_status_file $FC $MPIFC $BLACS_LIBS" >&5
 printf "%s\n" "$as_me: Loaded $pac_cv_status_file $FC $MPIFC $BLACS_LIBS" >&6;}
-am__api_version='1.17'
+am__api_version='1.18'
 
 
 
@@ -3721,10 +3721,14 @@ am_lf='
 '
 case `pwd` in
   *[\\\"\#\$\&\'\`$am_lf]*)
+    { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5
+printf "%s\n" "no" >&6; }
     as_fn_error $? "unsafe absolute working directory name" "$LINENO" 5;;
 esac
 case $srcdir in
   *[\\\"\#\$\&\'\`$am_lf\ \	]*)
+    { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5
+printf "%s\n" "no" >&6; }
     as_fn_error $? "unsafe srcdir value: '$srcdir'" "$LINENO" 5;;
 esac
 
@@ -4189,9 +4193,133 @@ AMTAR='$${TAR-tar}'
 
 
 # We'll loop over all known methods to create a tar archive until one works.
-_am_tools='gnutar  pax cpio none'
+_am_tools='gnutar plaintar pax cpio none'
 
-am__tar='$${TAR-tar} chof - "$$tardir"' am__untar='$${TAR-tar} xf -'
+# The POSIX 1988 'ustar' format is defined with fixed-size fields.
+      # There is notably a 21 bits limit for the UID and the GID.  In fact,
+      # the 'pax' utility can hang on bigger UID/GID (see automake bug#8343
+      # and bug#13588).
+      am_max_uid=2097151 # 2^21 - 1
+      am_max_gid=$am_max_uid
+      # The $UID and $GID variables are not portable, so we need to resort
+      # to the POSIX-mandated id(1) utility.  Errors in the 'id' calls
+      # below are definitely unexpected, so allow the users to see them
+      # (that is, avoid stderr redirection).
+      am_uid=`id -u || echo unknown`
+      am_gid=`id -g || echo unknown`
+      { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking whether UID '$am_uid' is supported by ustar format" >&5
+printf %s "checking whether UID '$am_uid' is supported by ustar format... " >&6; }
+      if test x$am_uid = xunknown; then
+        { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: ancient id detected; assuming current UID is ok, but dist-ustar might not work" >&5
+printf "%s\n" "$as_me: WARNING: ancient id detected; assuming current UID is ok, but dist-ustar might not work" >&2;}
+      elif test $am_uid -le $am_max_uid; then
+        { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5
+printf "%s\n" "yes" >&6; }
+      else
+        { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5
+printf "%s\n" "no" >&6; }
+        _am_tools=none
+      fi
+      { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking whether GID '$am_gid' is supported by ustar format" >&5
+printf %s "checking whether GID '$am_gid' is supported by ustar format... " >&6; }
+      if test x$gm_gid = xunknown; then
+        { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: ancient id detected; assuming current GID is ok, but dist-ustar might not work" >&5
+printf "%s\n" "$as_me: WARNING: ancient id detected; assuming current GID is ok, but dist-ustar might not work" >&2;}
+      elif test $am_gid -le $am_max_gid; then
+        { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5
+printf "%s\n" "yes" >&6; }
+      else
+        { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5
+printf "%s\n" "no" >&6; }
+        _am_tools=none
+      fi
+
+  { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking how to create a ustar tar archive" >&5
+printf %s "checking how to create a ustar tar archive... " >&6; }
+
+  # Go ahead even if we have the value already cached.  We do so because we
+  # need to set the values for the 'am__tar' and 'am__untar' variables.
+  _am_tools=${am_cv_prog_tar_ustar-$_am_tools}
+
+  for _am_tool in $_am_tools; do
+    case $_am_tool in
+    gnutar)
+      for _am_tar in tar gnutar gtar; do
+        { echo "$as_me:$LINENO: $_am_tar --version" >&5
+   ($_am_tar --version) >&5 2>&5
+   ac_status=$?
+   echo "$as_me:$LINENO: \$? = $ac_status" >&5
+   (exit $ac_status); } && break
+      done
+      am__tar="$_am_tar --format=ustar -chf - "'"$$tardir"'
+      am__tar_="$_am_tar --format=ustar -chf - "'"$tardir"'
+      am__untar="$_am_tar -xf -"
+      ;;
+    plaintar)
+      # Must skip GNU tar: if it does not support --format= it doesn't create
+      # ustar tarball either.
+      (tar --version) >/dev/null 2>&1 && continue
+      am__tar='tar chf - "$$tardir"'
+      am__tar_='tar chf - "$tardir"'
+      am__untar='tar xf -'
+      ;;
+    pax)
+      am__tar='pax -L -x ustar -w "$$tardir"'
+      am__tar_='pax -L -x ustar -w "$tardir"'
+      am__untar='pax -r'
+      ;;
+    cpio)
+      am__tar='find "$$tardir" -print | cpio -o -H ustar -L'
+      am__tar_='find "$tardir" -print | cpio -o -H ustar -L'
+      am__untar='cpio -i -H ustar -d'
+      ;;
+    none)
+      am__tar=false
+      am__tar_=false
+      am__untar=false
+      ;;
+    esac
+
+    # If the value was cached, stop now.  We just wanted to have am__tar
+    # and am__untar set.
+    test -n "${am_cv_prog_tar_ustar}" && break
+
+    # tar/untar a dummy directory, and stop if the command works.
+    rm -rf conftest.dir
+    mkdir conftest.dir
+    echo GrepMe > conftest.dir/file
+    { echo "$as_me:$LINENO: tardir=conftest.dir && eval $am__tar_ >conftest.tar" >&5
+   (tardir=conftest.dir && eval $am__tar_ >conftest.tar) >&5 2>&5
+   ac_status=$?
+   echo "$as_me:$LINENO: \$? = $ac_status" >&5
+   (exit $ac_status); }
+    rm -rf conftest.dir
+    if test -s conftest.tar; then
+      { echo "$as_me:$LINENO: $am__untar <conftest.tar" >&5
+   ($am__untar <conftest.tar) >&5 2>&5
+   ac_status=$?
+   echo "$as_me:$LINENO: \$? = $ac_status" >&5
+   (exit $ac_status); }
+      { echo "$as_me:$LINENO: cat conftest.dir/file" >&5
+   (cat conftest.dir/file) >&5 2>&5
+   ac_status=$?
+   echo "$as_me:$LINENO: \$? = $ac_status" >&5
+   (exit $ac_status); }
+      grep GrepMe conftest.dir/file >/dev/null 2>&1 && break
+    fi
+  done
+  rm -rf conftest.dir
+
+  if test ${am_cv_prog_tar_ustar+y}
+then :
+  printf %s "(cached) " >&6
+else case e in #(
+  e) am_cv_prog_tar_ustar=$_am_tool ;;
+esac
+fi
+
+  { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: $am_cv_prog_tar_ustar" >&5
+printf "%s\n" "$am_cv_prog_tar_ustar" >&6; }
 
 
 
@@ -5210,7 +5338,10 @@ _ACEOF
       break
     fi
   done
-  rm -f core conftest*
+  # aligned with autoconf, so not including core; see bug#72225.
+  rm -f -r a.out a.exe b.out conftest.$ac_ext conftest.$ac_objext \
+    conftest.dSYM conftest1.$ac_ext conftest1.$ac_objext conftest1.dSYM \
+    conftest2.$ac_ext conftest2.$ac_objext conftest2.dSYM
   unset am_i ;;
 esac
 fi
@@ -10651,6 +10782,12 @@ if test "x$amg4psblas_cv_have_mumps" == "xyes" ; then
       { printf "%s\n" "$as_me:${as_lineno-$LINENO}: PSBLAS defines PSB_LPK_ as $pac_cv_psblas_lpk. MUMPS interfacing will fail when called in global mode on very large matrices. " >&5
 printf "%s\n" "$as_me: PSBLAS defines PSB_LPK_ as $pac_cv_psblas_lpk. MUMPS interfacing will fail when called in global mode on very large matrices. " >&6;}
    fi
+   MUMPS_LIBS="-lsmumps -ldmumps -lcmumps -lzmumps -lmumps_common -lpord"
+   if test "x$amg4psblas_cv_mumpslibdir" != "x" ; then
+      { printf "%s\n" "$as_me:${as_lineno-$LINENO}: MUMPSLIBDIR  $amg4psblas_cv_mumpslibdir .." >&5
+printf "%s\n" "$as_me: MUMPSLIBDIR  $amg4psblas_cv_mumpslibdir .." >&6;}
+      MUMPS_LIBS="${MUMPS_LIBS} -L$amg4psblas_cv_mumpslibdir"
+   fi
    if test "x$pac_mumps_fmods_ok" == "xyes" ; then
       FDEFINES="$amg_cv_define_prepend-DAMG_HAVE_MUMPS $amg_cv_define_prepend-DAMG_HAVE_MUMPS_MODULES $MUMPS_MODULES $FDEFINES"
        MUMPS_FLAGS="-DAMG_HAVE_MUMPS $MUMPS_MODULES"

configure.ac

@@ -763,7 +763,11 @@ dnl fi
 if test "x$amg4psblas_cv_have_mumps" == "xyes" ; then
    if test "x$pac_cv_psblas_lpk" == "x8" ; then
       AC_MSG_NOTICE([PSBLAS defines PSB_LPK_ as $pac_cv_psblas_lpk. MUMPS interfacing will fail when called in global mode on very large matrices. ])
-   fi	
+   fi
+   MUMPS_LIBS="-lsmumps -ldmumps -lcmumps -lzmumps -lmumps_common -lpord"
+   if test "x$amg4psblas_cv_mumpslibdir" != "x" ; then
+      MUMPS_LIBS="${MUMPS_LIBS} -L$amg4psblas_cv_mumpslibdir"
+   fi 		   
    if test "x$pac_mumps_fmods_ok" == "xyes" ; then 
       FDEFINES="$amg_cv_define_prepend-DAMG_HAVE_MUMPS $amg_cv_define_prepend-DAMG_HAVE_MUMPS_MODULES $MUMPS_MODULES $FDEFINES"
        MUMPS_FLAGS="-DAMG_HAVE_MUMPS $MUMPS_MODULES"

samples/advanced/pdegen/amg_d_genpde_mod.F90

@@ -176,7 +176,7 @@ contains
     else
       partition_ = 3
     end if
-    deltah   = done/(idim+2)
+    deltah   = done/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_dpk_* deltah
 
@@ -412,9 +412,9 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
           ! x, y, z coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
-          z = (iz-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
+          z = (iz)*deltah
           zt(k) = f_(x,y,z)
           ! internal point: build discretization
           !
@@ -643,7 +643,7 @@ contains
       f_ => d_null_func_2d
     end if
 
-    deltah   = done/(idim+2)
+    deltah   = done/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_dpk_* deltah
 
@@ -875,8 +875,8 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,glob_row,idim,idim)
           ! x, y coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
 
           zt(k) = f_(x,y)
           ! internal point: build discretization

samples/advanced/pdegen/amg_s_genpde_mod.F90

@@ -176,7 +176,7 @@ contains
     else
       partition_ = 3
     end if
-    deltah   = sone/(idim+2)
+    deltah   = sone/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_spk_* deltah
 
@@ -412,9 +412,9 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
           ! x, y, z coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
-          z = (iz-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
+          z = (iz)*deltah
           zt(k) = f_(x,y,z)
           ! internal point: build discretization
           !
@@ -643,7 +643,7 @@ contains
       f_ => s_null_func_2d
     end if
 
-    deltah   = sone/(idim+2)
+    deltah   = sone/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_spk_* deltah
 
@@ -875,8 +875,8 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,glob_row,idim,idim)
           ! x, y coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
 
           zt(k) = f_(x,y)
           ! internal point: build discretization

samples/simple/pdegen/amg_dpde_mod.f90

@@ -369,9 +369,9 @@ contains
         ! compute gridpoint coordinates
         call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
         ! x, y, z coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
-        z = (iz-1)*deltah
+        x = (ix)*deltah
+        y = (iy)*deltah
+        z = (iz)*deltah
         zt(k) = f_(x,y,z)
         ! internal point: build discretization
         !   

samples/simple/pdegen/amg_spde_mod.f90

@@ -369,9 +369,9 @@ contains
         ! compute gridpoint coordinates
         call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
         ! x, y, z coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
-        z = (iz-1)*deltah
+        x = (ix)*deltah
+        y = (iy)*deltah
+        z = (iz)*deltah
         zt(k) = f_(x,y,z)
         ! internal point: build discretization
         !