diff --git a/README.md b/README.md
index fda14b8..58e98a5 100644
--- a/README.md
+++ b/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
diff --git a/doc/documentation.md b/doc/documentation.md
index 702acc4..882fc7f 100644
--- a/doc/documentation.md
+++ b/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.:  
 
 ````  
diff --git a/src/modsparse.f90 b/src/modsparse.f90
index 34ec13f..1ba2649 100644
--- a/src/modsparse.f90
+++ b/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
diff --git a/src/modsparse_gen.f90 b/src/modsparse_gen.f90
index 422b1a8..1b45221 100644
--- a/src/modsparse_gen.f90
+++ b/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
diff --git a/test/modtest_coo.f90 b/test/modtest_coo.f90
index de806c1..30a64a6 100644
--- a/test/modtest_coo.f90
+++ b/test/modtest_coo.f90
@@ -5,9 +5,10 @@ module modtest_coo
  use, intrinsic :: iso_fortran_env, only: int64, output_unit, wp => real64
 #endif
  use testdrive, only: new_unittest, unittest_type, error_type, check
- use modsparse, only: coosparse
- use modtest_common, only: tol_wp, verbose, ia, ja, a, aspsd&
-                       , addval => addval_coo, getmat, matcheck, printmat
+  use modsparse, only: coosparse
+  use modrandom, only: setseed, rand_stdnormal
+  use modtest_common, only: tol_wp, verbose, ia, ja, a, aspsd&
+                        , addval => addval_coo, getmat, matcheck, printmat
  implicit none
  private
 
@@ -60,6 +61,9 @@ subroutine collect_coo(testsuite)
     , new_unittest("coo multbymat_y_y_n", test_multbymat_y_y_n) &
     , new_unittest("coo multbymat_y_y_y", test_multbymat_y_y_y) &
     , new_unittest("coo multbymat_sym", test_multbymat_sym) &
+    , new_unittest("coo xAx", test_xAx_coo) &
+   , new_unittest("coo xAy", test_xAy_coo) &
+    , new_unittest("coo traceproduct", test_traceproduct_coo) &
     , new_unittest("coo nonzero", test_nonzero) &
     , new_unittest("coo nonzero_sym", test_nonzero_sym) &
     , new_unittest("coo scale", test_scale) &
@@ -1164,9 +1168,176 @@ subroutine test_submatrix_upper_full(error)
  call check(error, coosub%getdim(2), colend - colstart + 1, 'submatrix_dim2')
  if(allocated(error))return
 
- call check(error, all(getmat(coosub) == mat(rowstart:rowend, colstart:colend)) &
+  call check(error, all(getmat(coosub) == mat(rowstart:rowend, colstart:colend)) &
             , 'submatrix_upper_full')
 
 end subroutine
 
+!QUADRATIC FORM
+subroutine test_xAx_coo(error)
+  type(error_type), allocatable, intent(out) :: error
+
+  integer :: i
+  integer, parameter :: n = 6
+
+  type(coosparse) :: coo
+  real(wp) :: mat(n,n), x(n), expected, quad
+  real(wp) :: tol
+
+  !Known SPD matrix (aspsd + a(2,2)=202 => SPD)
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+
+  mat = getmat(coo)
+
+  !fixed x
+  x = [1._wp, 2._wp, 3._wp, 4._wp, 5._wp, 6._wp]
+  quad = coo%xAx(x)
+  expected = dot_product(x, matmul(mat, x))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(quad - expected) < tol, more = 'coo xAx fixed')
+  if(allocated(error))return
+
+  !random x
+  call setseed(42)
+  do i = 1, n
+    x(i) = rand_stdnormal()
+  enddo
+  quad = coo%xAx(x)
+  expected = dot_product(x, matmul(mat, x))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(quad - expected) < tol, more = 'coo xAx random')
+  if(allocated(error))return
+
+  !empty matrix -> 0
+  coo = coosparse(4, unlog = sparse_unit)
+  call check(error, abs(coo%xAx([1._wp, 2._wp, -1._wp, 3._wp])) < tol_wp&
+           , more = 'coo xAx empty')
+
+end subroutine test_xAx_coo
+
+!BILINEAR FORM
+subroutine test_xAy_coo(error)
+  type(error_type), allocatable, intent(out) :: error
+
+  integer :: i
+  integer, parameter :: n = 6
+
+  type(coosparse) :: coo
+  real(wp) :: mat(n,n), x(n), y(n), expected, bilin
+  real(wp) :: tol
+
+  !Known SPD matrix (aspsd + a(2,2)=202 => SPD)
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+
+  mat = getmat(coo)
+
+  !fixed x, y (x != y)
+  x = [1._wp, 2._wp, 3._wp, 4._wp, 5._wp, 6._wp]
+  y = [2._wp, 1._wp, 4._wp, 3._wp, 6._wp, 5._wp]
+  bilin = coo%xAy(x, y)
+  expected = dot_product(x, matmul(mat, y))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(bilin - expected) < tol, more = 'coo xAy fixed')
+  if(allocated(error))return
+
+  !identity: xAx(x) == xAy(x, x)
+  call check(error, abs(coo%xAx(x) - coo%xAy(x, x)) < tol_wp&
+           , more = 'coo xAx identity xAy(x,x)')
+  if(allocated(error))return
+
+  !random x, y
+  call setseed(42)
+  do i = 1, n
+    x(i) = rand_stdnormal()
+    y(i) = rand_stdnormal()
+  enddo
+  bilin = coo%xAy(x, y)
+  expected = dot_product(x, matmul(mat, y))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(bilin - expected) < tol, more = 'coo xAy random')
+  if(allocated(error))return
+
+  !empty matrix -> 0
+  coo = coosparse(4, unlog = sparse_unit)
+  call check(error, abs(coo%xAy([1._wp, 2._wp, -1._wp, 3._wp], [3._wp, -1._wp, 2._wp, 1._wp])) < tol_wp&
+           , more = 'coo xAy empty')
+
+end subroutine test_xAy_coo
+
+!TRACE OF A BLOCK PRODUCT
+subroutine test_traceproduct_coo(error)
+  type(error_type), allocatable, intent(out) :: error
+
+  integer :: id
+  integer, parameter :: n = 6
+
+  type(coosparse) :: coo
+  real(wp) :: mat(n,n), b(n,n), expected, result
+  real(wp) :: tol
+
+  !Known SPD matrix (aspsd + a(2,2)=202)
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+
+  mat = getmat(coo)
+
+  !Case 1: whole-matrix block, b = identity => expected = trace(mat)
+  b = 0._wp
+  do id = 1, n
+    b(id, id) = 1._wp
+  enddo
+  expected = sum([(mat(id, id), id = 1, n)])
+  result = coo%traceproduct(1, n, 1, n, b)
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'coo traceproduct whole I')
+  if(allocated(error))return
+
+  !Case 2: whole-matrix block, deterministic dense b
+  b = 0._wp
+  do id = 0, n*n - 1
+    b(mod(id, n) + 1, id/n + 1) = real(id, kind=wp)
+  enddo
+  expected = sum(mat * transpose(b))
+  result = coo%traceproduct(1, n, 1, n, b)
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'coo traceproduct whole dense')
+  if(allocated(error))return
+
+  !Case 3: 2x2 block at (1:2, 1:2)
+  b(1, 1) = 1._wp; b(1, 2) = 2._wp; b(2, 1) = 3._wp; b(2, 2) = 4._wp
+  expected = mat(1,1)*b(1,1) + mat(1,2)*b(2,1) &
+           + mat(2,1)*b(1,2) + mat(2,2)*b(2,2)
+  result = coo%traceproduct(1, 2, 1, 2, b(1:2, 1:2))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'coo traceproduct 2x2 diag')
+  if(allocated(error))return
+
+  !Case 4: 1x1 off-diagonal at (1,3): A(1,3)=13
+  expected = mat(1,3) * 7._wp
+  result = coo%traceproduct(1, 1, 3, 3, reshape([7._wp], [1,1]))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'coo traceproduct 1x1 off')
+  if(allocated(error))return
+
+  !Case 5: empty matrix -> 0
+  coo = coosparse(4, unlog = sparse_unit)
+  b(1:4, 1:4) = 0._wp
+  do id = 1, 4
+    b(id, id) = 1._wp
+  enddo
+  call check(error, abs(coo%traceproduct(1, 4, 1, 4, b(1:4, 1:4))) < tol_wp&
+           , more = 'coo traceproduct empty')
+
+end subroutine test_traceproduct_coo
+
+
+
 end module modtest_coo
diff --git a/test/modtest_crs.f90 b/test/modtest_crs.f90
index bcaf236..d76b1f5 100644
--- a/test/modtest_crs.f90
+++ b/test/modtest_crs.f90
@@ -5,8 +5,9 @@ module modtest_crs
  use, intrinsic :: iso_fortran_env, only: int64, wp => real64
 #endif
  use testdrive, only: new_unittest, unittest_type, error_type, check
- use modsparse, only: coosparse, crssparse, assignment(=)
- use modtest_common, only: tol_wp, verbose, ia, ja, a, aspsd, addval => addval_coo&
+  use modsparse, only: coosparse, crssparse, assignment(=)
+  use modrandom, only: setseed, rand_stdnormal
+  use modtest_common, only: tol_wp, verbose, ia, ja, a, aspsd, addval => addval_coo&
                        , getmat, matcheck, printmat&
                        , iaspsdf, jaspsdf, aspsdf&
                        , iaspsdf1, jaspsdf1, aspsdf1
@@ -38,6 +39,9 @@ subroutine collect_crs(testsuite)
     , new_unittest("crs test external", test_external) &
     , new_unittest("crs test harville", test_harville) &
     , new_unittest("crs cg", test_cg) &
+    , new_unittest("crs xAx", test_xAx_crs) &
+    , new_unittest("crs xAy", test_xAy_crs) &
+    , new_unittest("crs traceproduct", test_traceproduct_crs) &
     , new_unittest("crs chol", test_chol) &
     , new_unittest("crs diag vect", test_diag_vect) &
     , new_unittest("crs diag vect lupper", test_diag_vect_lupper) &
@@ -479,10 +483,202 @@ subroutine test_cg(error)
  mat_l = reshape([(merge(1._wp, 0._wp, i/nrow.eq.mod(i,nrow)), i = 0, nrow**2 - 1)], [nrow, nrow])
  where(getmat(crs) <= tol_wp) mat_l = 0._wp
 
- call check(error, all(abs(matmul(getmat(crs), xmat) - mat_l) < tol_wp), 'crs cg')
+  call check(error, all(abs(matmul(getmat(crs), xmat) - mat_l) < tol_wp), 'crs cg')
 
 end subroutine
 
+!QUADRATIC FORM
+subroutine test_xAx_crs(error)
+  type(error_type), allocatable, intent(out) :: error
+
+  integer :: i
+  integer, parameter :: n = 6
+
+  type(coosparse) :: coo
+  type(crssparse) :: crs
+  real(wp) :: mat(n,n), x(n), expected, quad
+  real(wp) :: tol
+
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+  crs = coo
+
+  mat = getmat(crs)
+
+  !fixed x
+  x = [1._wp, 2._wp, 3._wp, 4._wp, 5._wp, 6._wp]
+  quad = crs%xAx(x)
+  expected = dot_product(x, matmul(mat, x))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(quad - expected) < tol, more = 'crs xAx fixed')
+  if(allocated(error))return
+
+  !random x
+  call setseed(42)
+  do i = 1, n
+    x(i) = rand_stdnormal()
+  enddo
+  quad = crs%xAx(x)
+  expected = dot_product(x, matmul(mat, x))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(quad - expected) < tol, more = 'crs xAx random')
+  if(allocated(error))return
+
+  !empty matrix -> 0
+  coo = coosparse(4, unlog = sparse_unit)
+  crs = coo
+  call check(error, abs(crs%xAx([1._wp, 2._wp, -1._wp, 3._wp])) < tol_wp&
+           , more = 'crs xAx empty')
+
+end subroutine test_xAx_crs
+
+!BILINEAR FORM
+subroutine test_xAy_crs(error)
+  type(error_type), allocatable, intent(out) :: error
+
+  integer :: i
+  integer, parameter :: n = 6
+
+  type(coosparse) :: coo
+  type(crssparse) :: crs
+  real(wp) :: mat(n,n), x(n), y(n), expected, bilin
+  real(wp) :: tol
+
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+  crs = coo
+
+  mat = getmat(crs)
+
+  !fixed x, y (x != y)
+  x = [1._wp, 2._wp, 3._wp, 4._wp, 5._wp, 6._wp]
+  y = [2._wp, 1._wp, 4._wp, 3._wp, 6._wp, 5._wp]
+  bilin = crs%xAy(x, y)
+  expected = dot_product(x, matmul(mat, y))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(bilin - expected) < tol, more = 'crs xAy fixed')
+  if(allocated(error))return
+
+  !identity: xAx(x) == xAy(x, x)
+  call check(error, abs(crs%xAx(x) - crs%xAy(x, x)) < tol_wp&
+           , more = 'crs xAx identity xAy(x,x)')
+  if(allocated(error))return
+
+  !random x, y
+  call setseed(42)
+  do i = 1, n
+    x(i) = rand_stdnormal()
+    y(i) = rand_stdnormal()
+  enddo
+  bilin = crs%xAy(x, y)
+  expected = dot_product(x, matmul(mat, y))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(bilin - expected) < tol, more = 'crs xAy random')
+  if(allocated(error))return
+
+  !CRS result equals COO result for the same input (before reassigning the objects)
+  !different accumulation order -> allow a small relative tolerance
+  call check(error, abs(crs%xAy(x, y) - coo%xAy(x, y)) < tol_wp * max(abs(crs%xAy(x, y)), 1._wp)&
+           , more = 'crs-xAy close to coo')
+  if(allocated(error))return
+
+  !empty matrix -> 0
+  coo = coosparse(4, unlog = sparse_unit)
+  crs = coo
+  call check(error, abs(crs%xAy([1._wp, 2._wp, -1._wp, 3._wp], [3._wp, -1._wp, 2._wp, 1._wp])) < tol_wp&
+           , more = 'crs xAy empty')
+
+end subroutine test_xAy_crs
+
+!TRACE OF A BLOCK PRODUCT
+subroutine test_traceproduct_crs(error)
+  type(error_type), allocatable, intent(out) :: error
+
+  integer :: id
+  integer, parameter :: n = 6
+
+  type(coosparse) :: coo
+  type(crssparse) :: crs
+  real(wp) :: mat(n,n), b(n,n), expected, result
+  real(wp) :: tol
+
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+  crs = coo
+
+  mat = getmat(crs)
+
+  !Case 1: whole-matrix block, b = identity => expected = trace(mat)
+  b = 0._wp
+  do id = 1, n
+    b(id, id) = 1._wp
+  enddo
+  expected = sum([(mat(id, id), id = 1, n)])
+  result = crs%traceproduct(1, n, 1, n, b)
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'crs traceproduct whole I')
+  if(allocated(error))return
+
+  !Case 2: whole-matrix block, deterministic dense b
+  b = 0._wp
+  do id = 0, n*n - 1
+    b(mod(id, n) + 1, id/n + 1) = real(id, kind=wp)
+  enddo
+  expected = sum(mat * transpose(b))
+  result = crs%traceproduct(1, n, 1, n, b)
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'crs traceproduct whole dense')
+  if(allocated(error))return
+
+  !Case 3: 2x2 block at (1:2, 1:2)
+  b(1, 1) = 1._wp; b(1, 2) = 2._wp; b(2, 1) = 3._wp; b(2, 2) = 4._wp
+  expected = mat(1,1)*b(1,1) + mat(1,2)*b(2,1) &
+           + mat(2,1)*b(1,2) + mat(2,2)*b(2,2)
+  result = crs%traceproduct(1, 2, 1, 2, b(1:2, 1:2))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'crs traceproduct 2x2 diag')
+  if(allocated(error))return
+
+  !Case 4: 1x1 off-diagonal at (1,3): A(1,3)=13
+  expected = mat(1,3) * 7._wp
+  result = crs%traceproduct(1, 1, 3, 3, reshape([7._wp], [1,1]))
+  tol = 1.e-3_wp * max(abs(expected), 1._wp)
+  call check(error, abs(result - expected) < tol, more = 'crs traceproduct 1x1 off')
+  if(allocated(error))return
+
+  !Case 5: empty matrix -> 0
+  coo = coosparse(4, unlog = sparse_unit)
+  crs = coo
+  b(1:4, 1:4) = 0._wp
+  do id = 1, 4
+    b(id, id) = 1._wp
+  enddo
+  call check(error, abs(crs%traceproduct(1, 4, 1, 4, b(1:4, 1:4))) < tol_wp&
+           , more = 'crs traceproduct empty')
+  if(allocated(error))return
+
+  !Case 6: CRS result equals COO result for the same input
+  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
+  call coo%setsymmetric()
+  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
+  call coo%add(2, 2, 202._wp)
+  crs = coo
+  b(1:4, 1:4) = 0._wp
+  do id = 1, 4
+    b(id, id) = 1._wp
+  enddo
+  call check(error, crs%traceproduct(1, 3, 1, 3, b(1:3, 1:3)) == &
+                  coo%traceproduct(1, 3, 1, 3, b(1:3, 1:3)) &
+           , more = 'crs-traceproduct equal to coo')
+
+end subroutine test_traceproduct_crs
+
 !CHOLESKY FACTOR
 subroutine test_chol(error)
  type(error_type), allocatable, intent(out) :: error