psblas3

config/pac.m4
 configure.ac
 configure
 util/Makefile
 util/psb_amd_interface.c
 util/psb_d_renum_impl.F90

First version of amd interface. To be debugged.

config/pac.m4

@@ -1583,6 +1583,7 @@ fi
 
 
 if test "x$pac_amd_header_ok" == "xyes" ; then 
+      psblas_cv_amd_includes="$AMD_INCLUDES"
       AMD_LIBS="$psblas_cv_amd $AMD_LIBDIR"
       LIBS="$AMD_LIBS -lm $LIBS";
       AC_MSG_CHECKING([for amd_order in $AMD_LIBS])

configure

@@ -9152,6 +9152,7 @@ fi
 
 
 if test "x$pac_amd_header_ok" == "xyes" ; then
+      psblas_cv_amd_includes="$AMD_INCLUDES"
       AMD_LIBS="$psblas_cv_amd $AMD_LIBDIR"
       LIBS="$AMD_LIBS -lm $LIBS";
       { $as_echo "$as_me:${as_lineno-$LINENO}: checking for amd_order in $AMD_LIBS" >&5
@@ -9282,6 +9283,7 @@ CPPFLAGS="$SAVE_CPPFLAGS";
 
 if test "x$psblas_cv_have_amd" == "xyes" ; then
    FDEFINES="$mld_cv_define_prepend-DHAVE_AMD $FDEFINES"
+   CDEFINES="-DHAVE_AMD_ $psblas_cv_amd_includes $CDEFINES"
 fi
 
 

configure.ac

@@ -684,6 +684,7 @@ fi
 PAC_CHECK_AMD
 if test "x$psblas_cv_have_amd" == "xyes" ; then 
    FDEFINES="$mld_cv_define_prepend-DHAVE_AMD $FDEFINES"
+   CDEFINES="-DHAVE_AMD_ $psblas_cv_amd_includes $CDEFINES"
 fi
 
 

util/Makefile

@@ -17,7 +17,8 @@ IMPLOBJS= psb_s_hbio_impl.o psb_d_hbio_impl.o \
 	psb_s_renum_impl.o psb_d_renum_impl.o \
 	psb_c_renum_impl.o psb_z_renum_impl.o 
 MODOBJS=psb_util_mod.o $(BASEOBJS)
-OBJS=$(MODOBJS) $(IMPLOBJS)
+COBJS=psb_amd_interface.o
+OBJS=$(MODOBJS) $(IMPLOBJS) $(COBJS)
 LIBMOD=psb_util_mod$(.mod)
 LOCAL_MODS=$(MODOBJS:.o=$(.mod))
 LIBNAME=$(UTILLIBNAME)

util/psb_amd_interface.c

@@ -0,0 +1,45 @@
+/* 
+              Parallel Sparse BLAS  version 3.0
+    (C) Copyright 2006, 2007, 2008, 2009, 2010
+                       Salvatore Filippone    University of Rome Tor Vergata
+                       Alfredo Buttari        CNRS-IRIT, Toulouse
+ 
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions
+  are met:
+    1. Redistributions of source code must retain the above copyright
+       notice, this list of conditions and the following disclaimer.
+    2. Redistributions in binary form must reproduce the above copyright
+       notice, this list of conditions, and the following disclaimer in the
+       documentation and/or other materials provided with the distribution.
+    3. The name of the PSBLAS group or the names of its contributors may
+       not be used to endorse or promote products derived from this
+       software without specific written permission.
+ 
+  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+  TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+  POSSIBILITY OF SUCH DAMAGE.
+ 
+*/
+#ifdef HAVE_AMD_
+#include "amd.h"
+#endif
+
+int psb_amd_interface(int n, int Ap[], int Ai[], int P[])
+{ 
+  int i;
+#ifdef HAVE_AMD_
+  i= amd_order(n,Ap,Ai, P,(double *)NULL, (double *)NULL);
+  if ((i==AMD_OK)||(AMD_OK_BUT_JUMBLED)) return(0);
+#endif
+  return(-1);
+}
+

util/psb_d_renum_impl.F90

@@ -180,20 +180,33 @@ contains
 
 
   subroutine psb_mat_renum_amd(a,info,operm)
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    use iso_c_binding
+#endif
     use psb_base_mod    
     implicit none 
     type(psb_dspmat_type), intent(inout) :: a
     integer, intent(out) :: info
     integer, allocatable, optional, intent(out) :: operm(:)
-
+       
     ! 
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    interface 
+      function psb_amd_interface(n,ap,ai,p)&
+           & result(res) bind(c,name='psb_amd_interface')
+        use iso_c_binding
+        integer(c_int) :: res, n
+        integer(c_int) :: ap(*), ai(*), p(*)
+      end function psb_amd_interface
+    end interface
+#endif
+
+    type(psb_d_csc_sparse_mat)  :: acsc
     class(psb_d_base_sparse_mat), allocatable :: aa
-    type(psb_d_csr_sparse_mat)  :: acsr
     type(psb_d_coo_sparse_mat)  :: acoo
     
     integer :: err_act
     character(len=20)           :: name
-    integer, allocatable :: ndstk(:,:), iold(:), ndeg(:), perm(:) 
     integer :: i, j, k, ideg, nr, ibw, ipf, idpth
 
     info = psb_success_
@@ -204,67 +217,48 @@ contains
 
     info = psb_success_
     nr = a%get_nrows()
-    allocate(operm(nr))
+    allocate(perm(nr))
-    do i=1, nr
+
-      operm(i) = i
+    call a%mold(aa)
+    call a%mv_to(aa)
+    call aa%mv_to_fmt(acsc,info)
+
+    acsc%ia(:)  = acsc%ia(:) - 1
+    acsc%icp(:) = acsc%icp(:) - 1
+    info = psb_amd_interface(nr,acsc%icp,acsc%ia,perm)
+    if (info /= psb_success_) then 
+      info = psb_err_from_subroutine_
+      call psb_errpush(info,name) 
+      goto 9999
+    end if
+    perm(:)     = perm(:) + 1
+    acsc%ia(:)  = acsc%ia(:) - 1
+    acsc%icp(:) = acsc%icp(:) - 1
+
+    call acsc%mv_to_coo(acoo,info)
+    do i=1, acoo%get_nzeros()
+      acoo%ia(i) = perm(acoo%ia(i))
+      acoo%ja(i) = perm(acoo%ja(i))
     end do
-!!$    call a%mold(aa)
+    call acoo%fix(info) 
-!!$    call a%mv_to(aa)
+
-!!$    call aa%mv_to_fmt(acsr,info)
+    ! Get back to where we started from
-!!$    ! Insert call to gps_reduce
+    call aa%mv_from_coo(acoo,info)
-!!$    nr   = acsr%get_nrows()
+    call a%mv_from(aa)
-!!$    ideg = 0
+    if (present(operm)) then 
-!!$    do i=1, nr
+      call psb_realloc(nr,operm,info)
-!!$      ideg = max(ideg,acsr%irp(i+1)-acsr%irp(i))
+      if (info /= psb_success_) then 
-!!$    end do
+        info = psb_err_alloc_dealloc_
-!!$    allocate(ndstk(nr,ideg), iold(nr), perm(nr+1), ndeg(nr),stat=info)
+        call psb_errpush(info,name) 
-!!$    if (info /= 0) then
+        goto 9999
-!!$      info = psb_err_alloc_dealloc_
+      end if
-!!$      call psb_errpush(info, name) 
+      operm(1:nr) = perm(1:nr)
-!!$      goto 9999
+    end if
-!!$    end if
+    
-!!$    do i=1, nr
+    deallocate(aa,perm)
-!!$      iold(i) = i 
+    
-!!$      ndstk(i,:) = 0
-!!$      k  = 0
-!!$      do j=acsr%irp(i),acsr%irp(i+1)-1
-!!$        k = k + 1
-!!$        ndstk(i,k) = acsr%ja(j)
-!!$      end do
-!!$    end do
-!!$    perm = 0
-!!$
-!!$    call psb_gps_reduce(ndstk,nr,ideg,iold,perm,ndeg,ibw,ipf,idpth)
-!!$
-!!$    if (.not.psb_isaperm(nr,perm)) then 
-!!$      write(0,*) 'Something wrong: bad perm from gps_reduce'
-!!$      info = psb_err_from_subroutine_
-!!$      call psb_errpush(info,name) 
-!!$      goto 9999
-!!$    end if
-!!$    ! Move to coordinate to apply renumbering
-!!$    call acsr%mv_to_coo(acoo,info)
-!!$    do i=1, acoo%get_nzeros()
-!!$      acoo%ia(i) = perm(acoo%ia(i))
-!!$      acoo%ja(i) = perm(acoo%ja(i))
-!!$    end do
-!!$    call acoo%fix(info) 
-!!$
-!!$    ! Get back to where we started from
-!!$    call aa%mv_from_coo(acoo,info)
-!!$    call a%mv_from(aa)
-!!$    if (present(operm)) then 
-!!$      call psb_realloc(nr,operm,info)
-!!$      if (info /= psb_success_) then 
-!!$        info = psb_err_alloc_dealloc_
-!!$        call psb_errpush(info,name) 
-!!$        goto 9999
-!!$      end if
-!!$      operm(1:nr) = perm(1:nr)
-!!$    end if
-!!$
-!!$    deallocate(aa)
 #else 
+
     info = psb_err_missing_aux_lib_
     call psb_errpush(info, name) 
     goto 9999