Fixed bug in raw aggregation; this had already been discovered, but never

committed to svn.

mld2p4/psb_dbldaggrmat.f90

@@ -40,7 +40,7 @@ subroutine psb_dbldaggrmat(a,desc_a,ac,desc_ac,p,info)
   implicit none
 
   type(psb_dspmat_type), intent(in), target  :: a
-  type(psb_dspmat_type), intent(out), target :: ac
+  type(psb_dspmat_type), intent(inout), target :: ac
   type(psb_desc_type), intent(in)            :: desc_a
   type(psb_desc_type), intent(inout)         :: desc_ac
   type(psb_dbaseprc_type), intent(inout), target  :: p
@@ -100,9 +100,9 @@ contains
   subroutine raw_aggregate(info)
     use psb_base_mod
     use psb_prec_type
+    use mpi
     implicit none
 
-    include 'mpif.h'
     integer, intent(out)   :: info
     type(psb_dspmat_type)          :: b
     integer, pointer :: nzbr(:), idisp(:)
@@ -138,8 +138,8 @@ contains
       do i=1, nrow
         p%mlia(i) = p%mlia(i) + naggrm1
       end do
+      call psb_halo(p%mlia,desc_a,info)
     end if
-    call psb_halo(p%mlia,desc_a,info)
 
     if(info /= 0) then
       call psb_errpush(4010,name,a_err='psb_halo')
@@ -169,10 +169,16 @@ contains
 
     nzt = psb_sp_get_nnzeros(b)
 
+    j = 0 
     do i=1, nzt 
-      b%ia1(i) = p%mlia(b%ia1(i))
-      b%ia2(i) = p%mlia(b%ia2(i))
+      if ((1<=b%ia2(i)).and.(b%ia2(i)<=nrow)) then 
+        j = j + 1
+        b%aspk(j) = b%aspk(i)
+        b%ia1(j)  = p%mlia(b%ia1(i))
+        b%ia2(j)  = p%mlia(b%ia2(i))
+      end if
     enddo
+    b%infoa(psb_nnz_)=j
     call psb_fixcoo(b,info)
 
     nzt = psb_sp_get_nnzeros(b)
@@ -239,7 +245,6 @@ contains
         call psb_errpush(4010,name,a_err='psb_cdall')
         goto 9999
       end if
-      
       call psb_cdasb(desc_ac,info)
       if(info /= 0) then
         call psb_errpush(4010,name,a_err='psb_cdasb')

mld2p4/psb_zbldaggrmat.f90

@@ -41,7 +41,7 @@ subroutine psb_zbldaggrmat(a,desc_a,ac,desc_ac,p,info)
 
   type(psb_zspmat_type), intent(in), target  :: a
   type(psb_zbaseprc_type), intent(inout),target     :: p
-  type(psb_zspmat_type), intent(out), target :: ac
+  type(psb_zspmat_type), intent(inout), target :: ac
   type(psb_desc_type), intent(in)            :: desc_a
   type(psb_desc_type), intent(inout)         :: desc_ac
   integer, intent(out)                       :: info
@@ -99,11 +99,11 @@ contains
   subroutine raw_aggregate(info)
     use psb_base_mod
     use psb_prec_type
+    use mpi
     implicit none
 
-    include 'mpif.h'
     integer, intent(out)   :: info
-    type(psb_zspmat_type)          :: b, tmp
+    type(psb_zspmat_type)          :: b
     integer, pointer :: nzbr(:), idisp(:)
     integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
          & naggr, np, me, nzt,jl,nzl,nlr,&
@@ -137,8 +137,8 @@ contains
       do i=1, nrow
         p%mlia(i) = p%mlia(i) + naggrm1
       end do
+      call psb_halo(p%mlia,desc_a,info)
     end if
-    call psb_halo(p%mlia,desc_a,info)
 
     if(info /= 0) then
       call psb_errpush(4010,name,a_err='psb_halo')
@@ -168,10 +168,16 @@ contains
 
     nzt = psb_sp_get_nnzeros(b)
 
+    j = 0 
     do i=1, nzt 
-      b%ia1(i) = p%mlia(b%ia1(i))
-      b%ia2(i) = p%mlia(b%ia2(i))
+      if ((1<=b%ia2(i)).and.(b%ia2(i)<=nrow)) then 
+        j = j + 1
+        b%aspk(j) = b%aspk(i)
+        b%ia1(j)  = p%mlia(b%ia1(i))
+        b%ia2(j)  = p%mlia(b%ia2(i))
+      end if
     enddo
+    b%infoa(psb_nnz_)=j
     call psb_fixcoo(b,info)
 
     nzt = psb_sp_get_nnzeros(b)