Support of xAy and traceproduct

README.md

@@ -36,6 +36,44 @@ See the brief [documentation](doc/documentation.md) for more details regarding t
 The brief documentation is available in the directory *doc* (see *mainpage.md*). An extended documentation can be generated with *Doxygen*.  
 
 
+## Usage examples  
+
+To __evaluate the quadratic form__ `x' A x` for a sparse matrix `A` and a dense vector `x`:  
+
+````  
+use modsparse, only: crssparse  
+type(crssparse) :: A  
+real(8) :: x(100), q  
+
+!... build A and set x ...  
+q = A%xAx(x)  
+````  
+
+To __evaluate the bilinear form__ `x' A y` for a sparse matrix `A` and two dense vectors `x` and `y` (`xAx` is the special case `x = y`):  
+
+````  
+use modsparse, only: crssparse  
+type(crssparse) :: A  
+real(8) :: x(100), y(100), q  
+
+!... build A and set x, y ...  
+q = A%xAy(x, y)  
+````  
+
+To __compute the trace of the product of a square sub-block of a sparse matrix with a dense matrix__, e.g. `trace(A(3:6, 3:6) * B)`:  
+
+````  
+use modsparse, only: crssparse  
+type(crssparse) :: A  
+real(8) :: B(4,4), t  
+
+!... build A and set B ...  
+t = A%traceproduct(3, 6, 3, 6, B)  
+````  
+
+*xAx*, *xAy*, and *traceproduct* work on matrices in COO and CRS format and run in O(nnz). `xAy` also works on non-square matrices (length of `x` = rows, length of `y` = columns). Indices in *traceproduct* are 1-based and inclusive.  
+
+
 ## Acknowledgements  
 
 The code in modspainv.f90 is based on code originally written by Karin Meyer

doc/documentation.md

@@ -202,6 +202,30 @@ call mat%multbyv(alpha, trans, v, val, y)
 where *alpha* and *val* are double-precision real values, *v* and *y* are vectors, and *trans* (= 'n' or 't') relates to the transposition of the matrix.   
 The method *multbyv* is based on the MKL Sparse BLAS library.  
 
+To __evaluate the quadratic form__ *x' A x* for a sparse (square) matrix *A* and a dense vector *x*, the method *xAx* must be used:  
+
+````  
+q = mat%xAx(x)  
+````  
+
+where *mat* is a `coosparse` or `crssparse` object and *x* has the same length as the matrix dimension. The result is computed as *dot(x, A*x)*, reusing the format-specific matrix-vector product, and therefore runs in O(nnz).
+
+To __evaluate the bilinear form__ *x' A y* for a sparse matrix *A* and two dense vectors *x* and *y* (different from each other), the method *xAy* must be used:
+
+````
+q = mat%xAy(x, y)
+````
+
+where *mat* is a `coosparse` or `crssparse` object, *x* has the length of the number of rows of *A*, and *y* has the length of the number of columns of *A* (both hold when *A* is square). `xAx` is the special case `xAx(x) = xAy(x, x)`. The result is computed as *dot(x, A*y)*, reusing the format-specific matrix-vector product, and therefore runs in O(nnz). Unlike `xAx`, `xAy` also works on non-square matrices.
+
+To __compute the trace of the product of a square sub-block of a sparse matrix with a dense matrix__, the method *traceproduct* must be used:
+
+````  
+t = mat%traceproduct(r1, r2, c1, c2, b)  
+````  
+
+where *(r1, r2, c1, c2)* is a 1-based, inclusive tuple defining the square sub-block *A(r1:r2, c1:c2)* (i.e. `r2-r1+1 == c2-c1+1`), and *b* is a dense square matrix of that block size. The returned value is `trace(A(r1:r2, c1:c2) * b)`. On CRS the cost is O(nnz) in the block; on COO it is O(nnz) overall; an informative error is raised on unsupported formats. Setting *b* to the identity gives the plain sub-block trace.
+
 To get the number of __non-zero elements__ of a sparse matrix, the method *nonzero* must be used, e.g.:  
 
 ````  

src/modsparse.f90

@@ -72,10 +72,16 @@ module modsparse
   procedure,public::issorted
   !> @brief Returns true if square matrix; else returns false
   procedure,public::issquare
-  !> @brief Multiplication with a vector or matrix
-  generic,public::mult=>multbyv,multbym
-  !> @brief Prints the sparse matrix to a file
-  procedure,public::printtofile=>printtofile_gen
+   !> @brief Multiplication with a vector or matrix
+   generic,public::mult=>multbyv,multbym
+    !> @brief Computes the quadratic form; e.g., rv = mat%xAx(x)
+    procedure,public::xAx=>xAx_gen
+    !> @brief Computes the bilinear form x'Ay; e.g., rv = mat%xAy(x,y)
+    procedure,public::xAy=>xAy_gen
+    !> @brief Computes the trace of a block product; e.g., rv = mat%traceproduct(r1,r2,c1,c2,b)
+    procedure,public::traceproduct=>traceproduct_gen
+   !> @brief Prints the sparse matrix to a file
+   procedure,public::printtofile=>printtofile_gen
   !> @brief Prints the sparse matrix in a rectangular/square format to the output mat\%unlog
   procedure,public::printsquaretofile=>printsquaretofile_gen
   !> @brief Prints the current status of a sparse matrix to the output mat\%unlog ; e.g., call mat%printstats()
@@ -156,16 +162,36 @@ subroutine printsquare_gen(sparse,output)
   module elemental subroutine destroy_gen_gen(sparse)
    class(gen_sparse),intent(inout)::sparse
   end subroutine
-  !**CONJUGATE GRADIENT
-  module subroutine cg_gen(sparse,x,y,maxiter,tol)
-   !sparse*x=y
-   class(gen_sparse),intent(in)::sparse
-   integer(kind=int32),intent(inout),optional::maxiter
-   real(kind=wp),intent(inout)::x(:)
-   real(kind=wp),intent(in)::y(:)
-   real(kind=wp),intent(inout),optional::tol
-  end subroutine
-  !**GET ELEMENTS
+   !**CONJUGATE GRADIENT
+   module subroutine cg_gen(sparse,x,y,maxiter,tol)
+    !sparse*x=y
+    class(gen_sparse),intent(in)::sparse
+    integer(kind=int32),intent(inout),optional::maxiter
+    real(kind=wp),intent(inout)::x(:)
+    real(kind=wp),intent(in)::y(:)
+    real(kind=wp),intent(inout),optional::tol
+   end subroutine
+   !**QUADRATIC FORM
+   module function xAx_gen(sparse,x) result(rv)
+    class(gen_sparse),intent(in)::sparse
+    real(kind=wp),intent(in)::x(:)
+    real(kind=wp)::rv
+   end function
+   !**BILINEAR FORM
+   module function xAy_gen(sparse,x,y) result(rv)
+    class(gen_sparse),intent(in)::sparse
+    real(kind=wp),intent(in)::x(:)
+    real(kind=wp),intent(in)::y(:)
+    real(kind=wp)::rv
+   end function
+   !**TRACE OF A BLOCK PRODUCT
+   module function traceproduct_gen(sparse,r1,r2,c1,c2,b) result(rv)
+    class(gen_sparse),intent(in)::sparse
+    integer(kind=int32),intent(in)::r1,r2,c1,c2
+    real(kind=wp),intent(in)::b(:,:)
+    real(kind=wp)::rv
+   end function
+   !**GET ELEMENTS
   pure module function getdim_gen(sparse,dim1) result(dimget)
    class(gen_sparse),intent(in)::sparse
    integer(kind=int32),intent(in)::dim1

src/modsparse_gen.f90

@@ -75,11 +75,163 @@ module subroutine cg_gen(sparse,x,y,maxiter,tol)
   rsold = rsnew
  enddo
 
- if(present(maxiter)) maxiter = i
- if(present(tol)) tol = sqrt(rsnew) / ynorm
+  if(present(maxiter)) maxiter = i
+  if(present(tol)) tol = sqrt(rsnew) / ynorm
 
 end subroutine
 
+!**BILINEAR FORM
+module function xAy_gen(sparse,x,y) result(rv)
+ !x' A y (A need not be square; size(x)=rows, size(y)=cols)
+ !Accumulates in-place over the stored non-zero elements (no dense A*y, O(nnz)).
+ !OpenMP parallel reduction over the independent entries.
+ class(gen_sparse),intent(in)::sparse
+ real(kind=wp),intent(in)::x(:)
+ real(kind=wp),intent(in)::y(:)
+ real(kind=wp)::rv
+
+ integer(kind=int32)::i,c
+ integer(kind=int64)::ic
+ logical::lsym
+
+ if(size(x).ne.sparse%getdim(1))then
+  write(sparse%unlog,'(a)')' ERROR (xAy): the length of x does not match the number of rows of the matrix'
+  error stop
+ endif
+ if(size(y).ne.sparse%getdim(2))then
+  write(sparse%unlog,'(a)')' ERROR (xAy): the length of y does not match the number of columns of the matrix'
+  error stop
+ endif
+
+ rv=0._wp
+ !A stored element (i,k) contributes A(i,k) x(i) y(k); a symmetric matrix also
+ !contributes its mirror A(k,i)=A(i,k), giving the extra A(i,k) x(k) y(i).
+ select type(sparse)
+   type is(coosparse)
+    lsym=sparse%lsymmetric .and. sparse%lupperstorage
+    !$omp parallel do reduction(+:rv)
+    do ic=1_int64,sparse%nel
+     if(sparse%ij(1,ic).eq.0)cycle
+     i=sparse%ij(1,ic)
+     c=sparse%ij(2,ic)
+     rv=rv+sparse%a(ic)*x(i)*y(c)
+     if(lsym.and.i.ne.c) rv=rv+sparse%a(ic)*x(c)*y(i)
+    enddo
+    !$omp end parallel do
+   type is(crssparse)
+    lsym=sparse%lsymmetric .and. sparse%lupperstorage
+    !$omp parallel do reduction(+:rv)
+    do i=1,sparse%dim1
+     do c=sparse%ia(i),sparse%ia(i+1)-1
+      rv=rv+sparse%a(c)*(x(i)*y(sparse%ja(c))&
+            +merge(x(sparse%ja(c))*y(i),0._wp,lsym.and.i.ne.sparse%ja(c)))
+     enddo
+    enddo
+    !$omp end parallel do
+  class default
+   write(sparse%unlog,'(a)')' ERROR (xAy): unsupported format'
+   call sparse%printstats
+   error stop
+  end select
+
+end function
+
+!**QUADRATIC FORM
+module function xAx_gen(sparse,x) result(rv)
+ !x' A x (A must be square); special case of xAy with x=y
+ class(gen_sparse),intent(in)::sparse
+ real(kind=wp),intent(in)::x(:)
+ real(kind=wp)::rv
+
+ if(.not.sparse%issquare())then
+  write(sparse%unlog,'(a)')' ERROR (xAx): the matrix must be square'
+  error stop
+ endif
+
+ rv=xAy_gen(sparse,x,x)
+
+end function
+
+!**TRACE OF A BLOCK PRODUCT
+module function traceproduct_gen(sparse,r1,r2,c1,c2,b) result(rv)
+ !trace( A(r1:r2, c1:c2) * b ) where A(r1:r2, c1:c2) is square
+ class(gen_sparse),intent(in)::sparse
+ integer(kind=int32),intent(in)::r1,r2,c1,c2
+ real(kind=wp),intent(in)::b(:,:)
+ real(kind=wp)::rv
+
+ integer(kind=int32)::i,j,k,nrow,ncol,kblk
+ integer(kind=int64)::ic
+ integer(kind=int32)::rmin,rmax
+ logical::lsym
+
+ rv=0._wp
+ nrow=sparse%getdim(1)
+ ncol=sparse%getdim(2)
+ kblk=r2-r1+1
+
+ if(r1.lt.1 .or. r2.gt.nrow .or. r1.gt.r2)then
+  write(sparse%unlog,'(a)')' ERROR (traceproduct): the block row range is out of bounds'
+  error stop
+ endif
+ if(c1.lt.1 .or. c2.gt.ncol .or. c1.gt.c2)then
+  write(sparse%unlog,'(a)')' ERROR (traceproduct): the block column range is out of bounds'
+  error stop
+ endif
+ if(r2-r1.ne.c2-c1)then
+  write(sparse%unlog,'(a)')' ERROR (traceproduct): the block is not square'
+  error stop
+ endif
+ if(size(b,1).ne.kblk .or. size(b,2).ne.kblk)then
+  write(sparse%unlog,'(a)')' ERROR (traceproduct): b does not have the block size'
+  error stop
+ endif
+
+  select type(sparse)
+   type is(coosparse)
+    lsym=sparse%lsymmetric .and. sparse%lupperstorage
+    !$omp parallel do reduction(+:rv)
+    do ic=1_int64,sparse%nel
+     if(sparse%ij(1,ic).eq.0)cycle
+     j=sparse%ij(1,ic)
+     k=sparse%ij(2,ic)
+     !stored element A(j,k), position (k-c1+1, j-r1+1) of b if within the block
+     if(j.ge.r1 .and. j.le.r2 .and. k.ge.c1 .and. k.le.c2)then
+      rv=rv+sparse%a(ic)*b(k-c1+1,j-r1+1)
+     endif
+     !symmetric mirror A(k,j), if stored as upper triangle
+     if(lsym .and. j.ne.k .and. k.ge.r1 .and. k.le.r2 .and. j.ge.c1 .and. j.le.c2)then
+      rv=rv+sparse%a(ic)*b(j-c1+1,k-r1+1)
+     endif
+    enddo
+    !$omp end parallel do
+   type is(crssparse)
+    lsym=sparse%lsymmetric .and. sparse%lupperstorage
+    rmin=min(r1,c1)
+    rmax=max(r2,c2)
+    !$omp parallel do reduction(+:rv)
+    do i=rmin,rmax
+     do j=sparse%ia(i),sparse%ia(i+1)-1
+      k=sparse%ja(j)
+      !stored element A(i,k), position (k-c1+1, i-r1+1) of b if within the block
+      if(i.ge.r1 .and. i.le.r2 .and. k.ge.c1 .and. k.le.c2)then
+       rv=rv+sparse%a(j)*b(k-c1+1,i-r1+1)
+      endif
+      !symmetric mirror A(k,i), if stored as upper triangle
+      if(lsym .and. i.ne.k .and. k.ge.r1 .and. k.le.r2 .and. i.ge.c1 .and. i.le.c2)then
+       rv=rv+sparse%a(j)*b(i-c1+1,k-r1+1)
+      endif
+     enddo
+    enddo
+    !$omp end parallel do
+   class default
+    write(sparse%unlog,'(a)')' ERROR (traceproduct): unsupported format'
+    call sparse%printstats
+    error stop
+  end select
+
+end function
+
 !**GET ELEMENTS
 pure module function getdim_gen(sparse,dim1) result(dimget)
  class(gen_sparse),intent(in)::sparse