interface to, use and test 1-norm asum from BLAS

fortran-lang · Sep 19, 2024 · 93b37ff · 93b37ff
1 parent 37616ad
commit 93b37ff
Show file tree

Hide file tree

Showing 3 changed files with 67 additions and 3 deletions.
diff --git a/src/stdlib_linalg_blas.fypp b/src/stdlib_linalg_blas.fypp
@@ -9,6 +9,60 @@ module stdlib_linalg_blas
      implicit none(type,external)
      public
 
+          interface asum
+          !! ASUM takes the sum of the absolute values.
+#ifdef STDLIB_EXTERNAL_BLAS
+               pure real(dp) function dasum( n, x, incx )
+                    import sp,dp,qp,ilp,lk 
+                    implicit none(type,external) 
+                    integer(ilp), intent(in) :: incx,n
+                    real(dp), intent(in) :: x(*)
+               end function dasum
+#else 
+               module procedure stdlib_dasum
+#endif
+#ifdef STDLIB_EXTERNAL_BLAS
+               pure real(dp) function dzasum( n, x, incx )
+                    import sp,dp,qp,ilp,lk 
+                    implicit none(type,external) 
+                    integer(ilp), intent(in) :: incx,n
+                    complex(dp), intent(in) :: x(*)
+               end function dzasum
+#else 
+               module procedure stdlib_dzasum
+#endif
+#:for rk,rt,ri in REAL_KINDS_TYPES
+#:if not rk in ["sp","dp"]
+               module procedure stdlib_${ri}$asum
+#:endif
+#:endfor
+#:for rk,rt,ri in CMPLX_KINDS_TYPES
+#:if not rk in ["sp","dp"]
+               module procedure stdlib_${c2ri(ri)}$zasum
+#:endif
+#:endfor               
+#ifdef STDLIB_EXTERNAL_BLAS
+               pure real(sp) function sasum( n, x, incx )
+                    import sp,dp,qp,ilp,lk 
+                    implicit none(type,external) 
+                    integer(ilp), intent(in) :: incx,n
+                    real(sp), intent(in) :: x(*)
+               end function sasum
+#else 
+               module procedure stdlib_sasum
+#endif
+#ifdef STDLIB_EXTERNAL_BLAS
+               pure real(sp) function scasum( n, x, incx )
+                    import sp,dp,qp,ilp,lk 
+                    implicit none(type,external) 
+                    integer(ilp), intent(in) :: incx,n
+                    complex(sp), intent(in) :: x(*)
+               end function scasum
+#else 
+               module procedure stdlib_scasum
+#endif
+          end interface asum
+
           interface axpy
           !! AXPY constant times a vector plus a vector.
 #ifdef STDLIB_EXTERNAL_BLAS

diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp
@@ -8,7 +8,7 @@
 ! Vector norms
 submodule(stdlib_linalg) stdlib_linalg_norms
      use stdlib_linalg_constants
-     use stdlib_linalg_blas, only: nrm2
+     use stdlib_linalg_blas, only: nrm2,asum
      use stdlib_linalg_lapack, only: lange
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR     
@@ -171,7 +171,6 @@ submodule(stdlib_linalg) stdlib_linalg_norms
 
         integer(ilp) :: sze,norm_request
         real(${rk}$) :: rorder
-        ${rt}$, pointer :: a1d(:)
         intrinsic :: abs, sum, sqrt, maxval, minval, conjg
 
         sze = size(a,kind=ilp)
@@ -195,10 +194,14 @@ submodule(stdlib_linalg) stdlib_linalg_norms
 
         select case(norm_request)
             case(NORM_ONE)
+                #:if rank==1
+                nrm = asum(sze,a,incx=1_ilp)
+                #:else
                 nrm = sum( abs(a) )
+                #:endif
             case(NORM_TWO)            
                 #:if rank==1
-                nrm = nrm2(sze,a,incx=1)
+                nrm = nrm2(sze,a,incx=1_ilp)
                 #:elif rt.startswith('complex')
                 nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) )                
                 #:else

diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp
@@ -212,6 +212,13 @@ module test_linalg_norm
                if (allocated(error)) return                    
 
            end do 
+
+           ! Compare ND whole vector norm with unrolled vector norm
+           write(msg,"('Unrolled (1d) vs ${rank}$d order=',i0,' norm')") order
+           call check(error,abs(norm(a,order)-norm(b,order))<tol*max(1.0_${rk}$,norm(a,order)),&
+                      trim(msg))
+           if (allocated(error)) return              
+
 
         end do