recover changes (fortran-lang#1117)

Author: jalvesz

doc/specs/stdlib_sparse.md

@@ -338,7 +338,7 @@ If the `diagonal` array has not been previously allocated, the `diag` subroutine
 
 ### Syntax
 
-`call ` [[stdlib_sparse_conversion(module):csr2sellc(interface)]] `(csr,ell[,num_nz_rows])`
+`call ` [[stdlib_sparse_conversion(module):csr2ell(interface)]] `(csr,ell[,num_nz_rows])`
 
 ### Arguments
 
@@ -361,4 +361,39 @@ If the `diagonal` array has not been previously allocated, the `diag` subroutine
 ### Example
 ```fortran
 {!example/linalg/example_sparse_spmv.f90!}
-```
+```
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+## Operator overloading (`+`, `-`, `*`, `/`) {#operators}
+
+### Status
+
+Experimental
+
+### Description
+
+The definition of all standard arithmetic operators have been overloaded to be applicable for the matrix types defined by `stdlib_sparse`. The operators have been overloaded to support the following tuple combinations of the left-hand-side of the operation: same type and kind matrix-matrix, matrix-scalar and scalar-matrix.
+
+### Syntax
+
+- Matrix-matrix operators :
+
+`C = A + B`
+
+`C = A - B`
+
+`C = A * B`
+
+`C = A / B`
+
+- Matrix scalar operators : 
+
+`B = A + alpha` or `B = alpha + A` 
+
+`B = A - alpha` or `B = alpha - A`
+
+`B = A * alpha` or `B = alpha * A`
+
+`B = A / alpha` or `B = alpha / A`
+
+*Note*: scalar addition and subtraction operators perform element-wise operations only on the stored (non-zero) values, not on the full mathematical matrix. Meaning, the sparsity pattern is preserved.

src/sparse/CMakeLists.txt

@@ -2,6 +2,7 @@ set(sparse_fppFiles
     stdlib_sparse_constants.fypp
     stdlib_sparse_conversion.fypp
     stdlib_sparse_kinds.fypp
+    stdlib_sparse_operators.fypp
     stdlib_sparse_spmv.fypp
     )
 

src/sparse/stdlib_sparse_conversion.fypp

@@ -161,6 +161,7 @@ module stdlib_sparse_conversion
         #:for k1, t1, s1 in (KINDS_TYPES)
         module procedure :: coo_from_ijv_${s1}$
         module procedure :: csr_from_ijv_${s1}$
+        module procedure :: csc_from_ijv_${s1}$
         module procedure :: ell_from_ijv_${s1}$
         module procedure :: sellc_from_ijv_${s1}$
         #:endfor
@@ -258,8 +259,9 @@ contains
         temp(1,1:COO%nnz) = COO%index(2,1:COO%nnz)
         temp(2,1:COO%nnz) = COO%index(1,1:COO%nnz)
         allocate(data, source = COO%data )
+
         nnz = COO%nnz
-        call sort_coo_unique_${s1}$( temp, data, nnz, COO%nrows, COO%ncols )
+        call sort_coo_unique_${s1}$( temp, data, nnz, COO%ncols, COO%nrows )
 
         if( allocated(CSC%row) ) then
             CSC%row(1:COO%nnz)  = temp(2,1:COO%nnz)
@@ -740,6 +742,40 @@ contains
     end subroutine
     #:endfor
 
+    #:for k1, t1, s1 in (KINDS_TYPES)
+    subroutine csc_from_ijv_${s1}$(CSC,row,col,data,nrows,ncols)
+        type(CSC_${s1}$_type), intent(inout) :: CSC
+        integer(ilp), intent(in) :: row(:)
+        integer(ilp), intent(in) :: col(:)
+        ${t1}$, intent(in), optional :: data(:)
+        integer(ilp), intent(in), optional :: nrows
+        integer(ilp), intent(in), optional :: ncols
+
+        integer(ilp) :: nrows_, ncols_
+        !---------------------------------------------------------
+        if(present(nrows)) then
+            nrows_ = nrows
+        else 
+            nrows_ = maxval(row)
+        end if
+        if(present(ncols)) then
+            ncols_ = ncols
+        else 
+            ncols_ = maxval(col)
+        end if
+        !---------------------------------------------------------
+        block
+            type(COO_${s1}$_type) :: COO
+            if(present(data)) then
+                call from_ijv(COO,row,col,data=data,nrows=nrows_,ncols=ncols_)
+            else 
+                call from_ijv(COO,row,col,nrows=nrows_,ncols=ncols_)
+            end if
+            call coo2csc(COO,CSC)
+        end block
+    end subroutine
+    #:endfor
+
     #:for k1, t1, s1 in (KINDS_TYPES)
     subroutine ell_from_ijv_${s1}$(ELL,row,col,data,nrows,ncols,num_nz_rows)
         type(ELL_${s1}$_type), intent(inout) :: ELL

src/sparse/stdlib_sparse_kinds.fypp

@@ -3,6 +3,10 @@
 #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
 #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
 #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
+#:set OP_NAMES = ["add","sub","mul","div"]
+#:set OP_SYMBOLS = ["+","-","*","/"]
+#:set OPERATORS = list(zip(OP_NAMES, OP_SYMBOLS))
+
 !! The `stdlib_sparse_kinds` module provides derived type definitions for different sparse matrices
 !!
 ! This code was modified from https://github.com/jalvesz/FSPARSE by its author: Alves Jose
@@ -13,6 +17,8 @@ module stdlib_sparse_kinds
     private
     public :: sparse_full, sparse_lower, sparse_upper
     public :: sparse_op_none, sparse_op_transpose, sparse_op_hermitian
+    public :: operator(+), operator(-), operator(*), operator(/)
+
     !! version: experimental
     !!
     !! Base sparse type holding the meta data related to the storage capacity of a matrix.
@@ -128,6 +134,34 @@ module stdlib_sparse_kinds
     end type
     #:endfor
 
+    #:for op, sym in OPERATORS
+    !! Overload the `${sym}$` operator for sparse matrices
+    !! [Specifications](../page/specs/stdlib_sparse.html#operators)    
+    interface operator(${sym}$)
+        #:for matrix in SPARSE_KINDS
+        #:for k, t, s in (KINDS_TYPES)
+        pure module function sparse_${op}$_${matrix}$_${s}$(a, b) result(c)
+        type(${matrix}$_${s}$_type), intent(in) :: a, b
+        type(${matrix}$_${s}$_type) :: c
+        end function
+
+        pure module function sparse_${op}$_${matrix}$_scalar_${s}$(a, b) result(c)
+        ${t}$, intent(in) :: a
+        type(${matrix}$_${s}$_type), intent(in) :: b
+        type(${matrix}$_${s}$_type) :: c
+        end function
+
+        pure module function sparse_${op}$_scalar_${matrix}$_${s}$(a, b) result(c)
+        type(${matrix}$_${s}$_type), intent(in) :: a
+        ${t}$, intent(in) :: b
+        type(${matrix}$_${s}$_type) :: c
+        end function
+        #:endfor
+        #:endfor
+    end interface
+
+    #:endfor
+
 contains
 
     !! (re)Allocate matrix memory for the COO type

src/sparse/stdlib_sparse_operators.fypp

@@ -0,0 +1,44 @@
+#:include "common.fypp"
+#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
+#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
+#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
+#:set OP_NAMES = ["add","sub","mul","div"]
+#:set OP_SYMBOLS = ["+","-","*","/"]
+#:set OPERATORS = list(zip(OP_NAMES, OP_SYMBOLS))
+
+submodule(stdlib_sparse_kinds) stdlib_sparse_operators
+  implicit none
+
+  contains
+
+#:for op, sym in OPERATORS
+#:for matrix in SPARSE_KINDS
+#:for k, t, s in (KINDS_TYPES)
+    pure module function sparse_${op}$_${matrix}$_${s}$(a, b) result(c)
+      type(${matrix}$_${s}$_type), intent(in) :: a, b
+      type(${matrix}$_${s}$_type) :: c
+      c = a
+      c%data = c%data ${sym}$ b%data
+    end function
+
+    pure module function sparse_${op}$_${matrix}$_scalar_${s}$(a, b) result(c)
+      ${t}$, intent(in) :: a
+      type(${matrix}$_${s}$_type), intent(in) :: b
+      type(${matrix}$_${s}$_type) :: c
+      c = b
+      c%data = a ${sym}$ c%data
+    end function
+
+    pure module function sparse_${op}$_scalar_${matrix}$_${s}$(a, b) result(c)
+      type(${matrix}$_${s}$_type), intent(in) :: a
+      ${t}$, intent(in) :: b
+      type(${matrix}$_${s}$_type) :: c
+      c = a
+      c%data = c%data ${sym}$ b
+    end function
+
+#:endfor
+#:endfor
+#:endfor
+
+end submodule stdlib_sparse_operators

test/linalg/test_linalg_sparse.fypp

@@ -1,6 +1,9 @@
 #:include "common.fypp"
 #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
 #:set KINDS_TYPES = R_KINDS_TYPES
+#:set OP_NAMES = ["add","sub","mul","div"]
+#:set OP_SYMBOLS = ["+","-","*","/"]
+#:set OPERATORS = list(zip(OP_NAMES, OP_SYMBOLS))
 module test_sparse_spmv
     use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
     use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64
@@ -25,7 +28,8 @@ contains
             new_unittest('sellc', test_sellc),  &
             new_unittest('symmetries', test_symmetries), &
             new_unittest('diagonal', test_diagonal), &
-            new_unittest('add_get_values', test_add_get_values) &
+            new_unittest('add_get_values', test_add_get_values),  &
+            new_unittest('sparse_operators', test_sparse_operators) &
         ]
     end subroutine
 
@@ -332,6 +336,7 @@ contains
             diagonal = 0.0
             call coo2csc( COO, CSC )
             call diag( CSC , diagonal )
+            print *, 'diagonal csc:', diagonal
             call check(error, all(diagonal == [1,2,3,4]) )
             if (allocated(error)) return
         end block
@@ -383,6 +388,105 @@ contains
         #:endfor
     end subroutine
 
+    subroutine test_sparse_operators(error)
+        !> Error handling
+        type(error_type), allocatable, intent(out) :: error
+        #:for matrix in ["COO","CSR","CSC"] # ELL and SELLC are supported but since they have explicit zeros in data array, the tests are less straightforward
+        #:for k, t, s in (KINDS_TYPES)
+        block
+            integer, parameter :: wp = ${k}$
+            real(wp), parameter :: tol = 1000*epsilon(0._wp)
+            real(wp) :: scalar
+            integer :: row(10), col(10)
+            real(wp) :: data(10)
+            type(${matrix}$_${s}$_type) :: a, b, c
+            ${t}$:: err
+
+            data(:)   = real([9,-3,4,7,8,-1,8,4,5,6],kind=wp)
+            col(:)    = [1,5,1,2,2,3,4,1,3,4]
+            row(:)    = [1,1,2,2,3,3,3,4,4,4]
+            
+            call from_ijv(a, row, col, data)
+            call from_ijv(b, row, col, data)
+            scalar = 2.0_wp
+            
+            ! Test sparse + sparse
+            c = a + b
+            err = sum( abs( c%data - 2.0_wp*a%data ) ) / size(data)
+            call check(error, err <= tol, "error in sparse + sparse" )
+            if (allocated(error)) return
+
+            ! Test scalar + sparse
+            c = scalar + a
+            err = sum( abs( c%data - (scalar + a%data) ) ) / size(data)
+            call check(error, err <= tol, "error in scalar + sparse" )
+            if (allocated(error)) return
+
+            ! Test sparse + scalar
+            c = a + scalar
+            err = sum( abs( c%data - (a%data + scalar) ) ) / size(data)
+            call check(error, err <= tol, "error in sparse + scalar" )
+            if (allocated(error)) return
+
+            ! Test sparse * sparse
+            c = a * b
+            err = sum( abs( c%data - a%data*b%data ) ) / size(data)
+            call check(error, err <= tol, "error in sparse * sparse" )
+            if (allocated(error)) return
+
+            ! Test scalar * sparse
+            c = scalar * a
+            err = sum( abs( c%data - (scalar*a%data) ) ) / size(data)
+            call check(error, err <= tol, "error in scalar * sparse" )
+            if (allocated(error)) return
+
+            ! Test sparse * scalar
+            c = a * scalar
+            err = sum( abs( c%data - (a%data*scalar) ) ) / size(data)
+            call check(error, err <= tol, "error in sparse * scalar" )
+            if (allocated(error)) return
+
+            ! Test sparse - sparse
+            c = a - b
+            err = sum( abs( c%data ) ) / size(data)
+            call check(error, err <= tol, "error in sparse - sparse" )
+            if (allocated(error)) return
+
+            ! Test scalar - sparse
+            c = scalar - a
+            err = sum( abs( c%data - (scalar - a%data) ) ) / size(data)
+            call check(error, err <= tol, "error in scalar - sparse" )
+            if (allocated(error)) return
+
+            ! Test sparse - scalar
+            c = a - scalar
+            err = sum( abs( c%data - (a%data - scalar) ) ) / size(data)
+            call check(error, err <= tol, "error in sparse - scalar" )
+            if (allocated(error)) return
+
+            ! Test sparse / sparse
+            c = a / b
+            err = sum( abs( c%data - 1.0_wp ) ) / size(data)
+            call check(error, err <= tol, "error in sparse / sparse" )
+            if (allocated(error)) return
+
+            ! Test scalar / sparse
+            c = scalar / a
+            err = sum( abs( c%data - (scalar / a%data) ) ) / size(data)
+            call check(error, err <= tol, "error in scalar / sparse" )
+            if (allocated(error)) return
+
+            ! Test sparse / scalar
+            c = a / scalar
+            err = sum( abs( c%data - (a%data / scalar) ) ) / size(data)
+            call check(error, err <= tol, "error in sparse / scalar" )
+            if (allocated(error)) return
+
+        end block
+        #:endfor
+        #:endfor
+    end subroutine
+
 end module