parallelism; however, since we implement a batched routine for independent operations, we will use the serial block Jacobi algorithm for individual matrices and rely on the parallelism of the batch processing. The parallel version, where multiple blocks are processed simultaneously, can still be used when the batch size is very small, but we will focus on the serial version. In this section we will discuss the implementation details for two global memory block Jacobi algorithms that differ only in the way block columns are orthogonalized and compare their performance with parallel streamed calls to the cuSOLVER 8 [@nvidia-cusolver] library routines.

Gram Matrix Block Jacobi SVD
----------------------------

The block Jacobi algorithm is very similar to the vector Algorithm \[alg:jacobi\], orthogonalizing pairs of blocks columns instead of vectors. The first method of orthogonalizing pairs of block columns is based on the SVD of their Gram matrix. During the $p$-th sweep, each pair of $m \times k$ block columns $A^{(p)}_i$ and $A^{(p)}_j$ is orthogonalized by forming a $2k \times 2k$ Gram matrix $G^{(p)}_{ij} = {[A^{(p)}_i A^{(p)}_j]}^T [A^{(p)}_i A^{(p)}_j] = {A^{(p)}_{ij}}^T A^{(p)}_{ij}$ and generating a block rotation matrix $U^{(p)}_{ij}$, computed as the left singular vectors of $G^{(p)}_{ij}$ (or equivalently its eigenvectors, since it is 