Written top-level Fortran95 version of SMMP.

src/serial/psb_dnumbmm.f90

@@ -31,6 +31,13 @@
 ! File:  psb_dnumbmm.f90 
 ! Subroutine: 
 ! Parameters:
+!
+!
+! Note: This subroutine performs the numerical product of two sparse matrices.
+!       It is modeled after the SMMP package by R. Bank and C. Douglas, but is 
+!       rewritten in Fortran 95 making use of our sparse matrix facilities.
+!
+!
 
 subroutine psb_dnumbmm(a,b,c)
   use psb_spmat_type
@@ -60,7 +67,11 @@ subroutine psb_dnumbmm(a,b,c)
     return
   endif
   call psb_realloc(size(c%ia1),c%aspk,info)
-  if (.true.) then 
+  !
+  ! Note: we still have to test about possible performance hits. 
+  !
+  !
+  if (.false.) then 
     call numbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
          & b%ia2,b%ia1,0,b%aspk,&
          & c%ia2,c%ia1,0,c%aspk,temp)
@@ -78,7 +89,56 @@ contains
     real(kind(1.d0))      :: temp(:)
     integer, allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
     real(kind(1.d0)), allocatable :: aval(:),bval(:)
-    integer  :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn
+    integer  :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn,minln
+    real(kind(1.d0))      :: ajj
+
+
+    n = a%m
+    m = a%k 
+    l = b%k 
+    maxlmn = max(l,m,n)
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & aval(maxlmn),bval(maxlmn), stat=info)
+    if (info /= 0) then 
+      return
+    endif
+
+    do i = 1,maxlmn
+      temp(i) = dzero
+    end do
+    minlm = min(l,m)
+    minln = min(l,n)
+    minmn = min(m,n)
+    do  i = 1,n
+
+      call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
+
+      do jj=1, nazr
+        j=iacl(jj)
+        ajj = aval(jj)
+        if ((j<1).or.(j>m)) then 
+          write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m
+        endif
+        call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
+          else
+            temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
+          endif
+        enddo
+      end do
+      do  j = c%ia2(i),c%ia2(i+1)-1
+        if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn))  then 
+          write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
+        else
+          c%aspk(j) = temp(c%ia1(j))
+          temp(c%ia1(j)) = dzero
+        endif
+      end do
+    end do
+
+
 
 
   end subroutine inner_numbmm

src/serial/psb_dsymbmm.f90

@@ -31,6 +31,12 @@
 ! File:  psb_dsymbmm.f90 
 ! Subroutine: 
 ! Parameters:
+!
+!
+! Note: This subroutine performs the symbolic product of two sparse matrices.
+!       It is modeled after the SMMP package by R. Bank and C. Douglas, but is 
+!       rewritten in Fortran 95 making use of our sparse matrix facilities.
+!
 
 subroutine psb_dsymbmm(a,b,c)
   use psb_spmat_type
@@ -73,6 +79,10 @@ subroutine psb_dsymbmm(a,b,c)
   nze = max(a%m+1,2*a%m)
   call psb_sp_reall(c,nze,info)
 !!$  write(0,*) 'SYMBMM90 ',size(c%pl),size(c%pr)
+  !
+  ! Note: we need to test whether there is a performance impact 
+  !       in not using the original Douglas & Bank code. 
+  !
   if (.false.) then 
     call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
          & b%ia2,b%ia1,0,&

src/serial/psb_znumbmm.f90

@@ -31,6 +31,13 @@
 ! File:  psb_znumbmm.f90 
 ! Subroutine: 
 ! Parameters:
+!
+!
+! Note: This subroutine performs the numerical product of two sparse matrices.
+!       It is modeled after the SMMP package by R. Bank and C. Douglas, but is 
+!       rewritten in Fortran 95 making use of our sparse matrix facilities.
+!
+!
 
 subroutine psb_znumbmm(a,b,c)
   use psb_spmat_type
@@ -40,12 +47,101 @@ subroutine psb_znumbmm(a,b,c)
   complex(kind(1.d0)), allocatable :: temp(:)
   integer                       :: info
   
-  allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
+  interface psb_sp_getrow
+    subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
+      use psb_spmat_type
+      type(psb_zspmat_type), intent(in) :: a
+      integer, intent(in)       :: irw
+      integer, intent(out)      :: nz
+      integer, intent(inout)    :: ia(:), ja(:)
+      complex(kind(1.d0)),  intent(inout)    :: val(:)
+      integer, intent(in), target, optional :: iren(:)
+      integer, intent(in), optional :: lrw
+      integer, intent(out)  :: info
+    end subroutine psb_zspgetrow
+  end interface
+
 
+  allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
+  if (info /= 0) then
+    return
+  endif
   call psb_realloc(size(c%ia1),c%aspk,info)
+  !
+  ! Note: we still have to test about possible performance hits. 
+  !
+  !
+  if (.false.) then 
     call znumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
          & b%ia2,b%ia1,0,b%aspk,&
          & c%ia2,c%ia1,0,c%aspk,temp)
+  else
+    call inner_numbmm(a,b,c,temp,info)
+  end if
   deallocate(temp) 
   return
+
+contains 
+
+  subroutine inner_numbmm(a,b,c,temp,info)
+    type(psb_zspmat_type) :: a,b,c
+    integer               :: info
+    complex(kind(1.d0))      :: temp(:)
+    integer, allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    complex(kind(1.d0)), allocatable :: aval(:),bval(:)
+    integer  :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn,minln
+    complex(kind(1.d0))      :: ajj
+
+
+    n = a%m
+    m = a%k 
+    l = b%k 
+    maxlmn = max(l,m,n)
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & aval(maxlmn),bval(maxlmn), stat=info)
+    if (info /= 0) then 
+      return
+    endif
+
+    do i = 1,maxlmn
+      temp(i) = dzero
+    end do
+    minlm = min(l,m)
+    minln = min(l,n)
+    minmn = min(m,n)
+    do  i = 1,n
+
+      call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
+
+      do jj=1, nazr
+        j=iacl(jj)
+        ajj = aval(jj)
+        if ((j<1).or.(j>m)) then 
+          write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m
+        endif
+        call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
+          else
+            temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
+          endif
+        enddo
+      end do
+      do  j = c%ia2(i),c%ia2(i+1)-1
+        if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn))  then 
+          write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
+        else
+          c%aspk(j) = temp(c%ia1(j))
+          temp(c%ia1(j)) = dzero
+        endif
+      end do
+    end do
+
+
+
+
+  end subroutine inner_numbmm
+
+
 end subroutine psb_znumbmm

src/serial/psb_zsymbmm.f90

@@ -31,6 +31,12 @@
 ! File:  psb_zsymbmm.f90 
 ! Subroutine: 
 ! Parameters:
+!
+!
+! Note: This subroutine performs the symbolic product of two sparse matrices.
+!       It is modeled after the SMMP package by R. Bank and C. Douglas, but is 
+!       rewritten in Fortran 95 making use of our sparse matrix facilities.
+!
 
 subroutine psb_zsymbmm(a,b,c)
   use psb_spmat_type
@@ -73,6 +79,10 @@ subroutine psb_zsymbmm(a,b,c)
   nze = max(a%m+1,2*a%m)
   call psb_sp_reall(c,nze,info)
 !!$  write(0,*) 'SYMBMM90 ',size(c%pl),size(c%pr)
+  !
+  ! Note: we need to test whether there is a performance impact 
+  !       in not using the original Douglas & Bank code. 
+  !
   if (.false.) then 
     call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
          & b%ia2,b%ia1,0,&
@@ -88,6 +98,7 @@ subroutine psb_zsymbmm(a,b,c)
   c%descra='GUN'
   deallocate(itemp) 
   return
+
 contains
   subroutine inner_symbmm(a,b,c,index,info)
     type(psb_zspmat_type) :: a,b,c
@@ -166,4 +177,5 @@ contains
       end do main
 
   end subroutine inner_symbmm
+
 end subroutine psb_zsymbmm