16. Linear Calculations of Sparse Matrices
Note – This document was created using AI translation.
16.1 Overview
Matrices in which many of the elements are 0 are referred to as sparse matrices. Such matrices often appear as coefficients in systems of linear equations for numerical solutions to partial differential equations.
Representing sparse matrices as regular two-dimensional arrays consumes significant memory to store the zeros and requires considerable CPU time for operations involving zeros during computation. To address this, special storage formats for sparse matrices are used.
Linear calculations with sparse matrices not only involve different storage formats, but also require different algorithms due to the large scale of the problems applied (e.g., involving tens of thousands or more). As a result, conventional linear calculation programs cannot be directly used. Specifically, iterative methods (such as Krylov subspace methods) are often employed for systems of linear equations, while methods such as the Lanczos method or the Arnoldi method are frequently used for eigenvalues and eigenvectors.
The following sections will explain example problems in sparse matrix calculations, the storage formats for sparse matrices used in XLPack, basic calculation programs, sparse matrix checks, solutions for systems of linear equations, file input/output, and more.
16.2 Example Problem in Sparse Matrix Calculations
This section explains an example of solving Laplace’s equation numerically using the five-point finite difference approximation.
Consider solving the following partial differential equation in a square domain where x = 0 to x = 1, and y = 0 to y = 1:
-d2u/dx2 - d2u/dy2 = 0
with Dirichlet boundary conditions:
u(x, 0) = 0, u(0, y) = 0,
u(x, 1) = x, u(1, y) = y
Such partial differential equations are applied, for example, to describe the temperature distribution of a plate in thermal equilibrium.
To perform finite difference approximation, the region in both the x and y directions is divided into N equal parts. Setting the division width as h = 1/N, the value of u at the grid point x = ih, y = jh (i = 0 to N, j = 0 to N) is represented as u[i,j]. The derivatives are approximated using differences as follows:
d2u/dx2 = (u[i-1,j] - 2u[i,j] + u[i+1,j])/h2 d2u/dy2 = (u[i,j-1] - 2u[i,j] + u[i,j+1])/h2
From this, Laplace’s equation is approximated by the following formula:
-u[i,j-1] - u[i-1,j] + 4u[i,j] - u[i+1,j] - u[i,j+1] = 0
The boundary conditions (Dirichlet conditions) directly assign values to the grid points along the boundary as follows:
u[i,0] = 0 (i = 0 to N) u[0,j] = 0 (j = 0 to N) u[i,N] = ih (i = 0 to N) u[N,j] = jh (j = 0 to N)
Thus, for the square region, the values on the edges are known, leaving only the interior part to be computed. Here, let m = N – 1, and reassign symbols as u1 = u[1,1], u2 = u[2,1], …, um = u[m,1], um+1 = u[1,2], …, um*m = u[m,m].
For N = 7 (m = 6), the arrangement is illustrated as follows.
In this case, out of the 64 grid points, 28 points are known from the boundary conditions, while the values of u at the remaining 36 points are unknown. This leads to the following system of linear equations:
Au = b
where, A is a 36 x 36 square matrix, u is a vector consisting of u1, u2, …, u36, and b is the right-hand side vector.
Using the approximate formula of Laplace’s equation above for i = 1 to 6, j = 1 to 6, the matrix A and vector b can be determined from the 36 equations.
( 4 -1 -1 )
(-1 4 -1 -1 )
( -1 4 -1 -1 )
( -1 4 -1 -1 )
( -1 4 -1 -1 )
( -1 4 -1 )
(-1 4 -1 -1 )
( -1 -1 4 -1 -1 )
( -1 -1 4 -1 -1 )
( -1 -1 4 -1 -1 )
( -1 -1 4 -1 -1 )
( -1 -1 4 -1 )
A = ( . . . )
( -1 4 -1 )
( -1 -1 4 -1 )
( -1 -1 4 -1 )
( -1 -1 4 -1 )
( -1 -1 4 -1 )
( -1 -1 4 )
( u[1,0] + u[0,1] = 0 )
( u[2,0] = 0 )
( u[3,0] = 0 )
( u[4,0] = 0 )
( u[5,0] = 0 )
( u[6,0] + u[6,1] = h )
( u[0,2] = 0 )
b = ( . . . )
( u[5,7] = 5h )
( u[6,7] + u[7,6] = 12h )
In this way, A becomes a symmetric sparse matrix with at most 5 non-zero elements in each row.
While the example with N = 7 was explained, achieving a more precise solution requires increasing N. For instance, if N = 100, the system of equations will have approximately 10,000 unknowns. Furthermore, for three-dimensional problems, it would involve approximately 1,000,000 unknowns. In such cases, conventional linear calculation programs can no longer be used, and programs specifically designed for large-scale sparse matrices become necessary.
16.3 Storage Formats for Sparse Matrices
XLPack uses the following three storage formats:
(1) Compressed Sparse Row (CSR) format
(2) Compressed Sparse Column (CSC) format
(3) Coordinate (COO) format
Formats (1) and (2) are used for computation, while format (3) is used for input and output. The basic functions of XLPack utilize only formats (1) and (3).
As an example, consider the following 3 × 4 matrix.
( 1 0 -2 4 ) ( -1 2 0 -9 ) ( 0 5 0 -8 )
(1) Compressed Sparse Row (CSR) Format
In CSR format, only the nonzero elements are stored using one floating-point array, val, and two integer arrays, rowptr and colind. The val array contains the values of nonzero elements arranged in row-wise order. The rowptr array represents the position of the first element of each row within the val array. The colind array stores the column indices corresponding to the elements in val. The last element of rowptr is set to the final position in val plus one.
Row and column indices generally start from 1 in Fortran programs and 0 in C programs, so compatibility should be considered.
C language format:
val { 1, -2, 4, -1, 2, -9, 5, -8 }
rowptr { 0, 3, 6, 8 }
colind { 0, 2, 3, 0, 1, 3, 1, 3 }
Fortran format:
val { 1, -2, 4, -1, 2, -9, 5, -8 }
rowptr { 1, 4, 7, 9 }
colind { 1, 3, 4, 1, 2, 4, 2, 4 }
For symmetric matrices, only the diagonal elements and either the left or right part of the matrix may be stored to reduce memory usage. Here, such matrices are referred to as SSR format (not a commonly used term).
In the example above, column indices within each row are arranged in ascending order. While some programs may allow arbitrary ordering within rows, XLPack requires that the column indices be sorted in ascending order for computational efficiency.
(2) Compressed Sparse Column (CSC) Format
Since CSC format is not used in XLPack’s basic functions, its explanation is omitted. However, it stores nonzero elements in a column-wise order.
(3) Coordinate (COO) Format
In COO format, only the nonzero elements are stored using one floating-point array, val, and two integer arrays, rowind and colind. The val array holds the values of nonzero elements, while rowind and colind store the corresponding row and column indices.
Row and column indices generally start from 1 in Fortran programs and 0 in C programs, so compatibility should be considered.
C language format:
val { 1, -2, 4, -1, 2, -9, 5, -8 }
rowind { 0, 0, 0, 1, 1, 1, 2, 2 }
colind { 0, 2, 3, 0, 1, 3, 1, 3 }
Fortran format:
val { 1, -2, 4, -1, 2, -9, 5, -8 }
rowind { 1, 1, 1, 2, 2, 2, 3, 3 }
colind { 1, 3, 4, 1, 2, 4, 2, 4 }
In the example above, row and column indices are arranged in ascending order. However, some programs may allow arbitrary ordering, depending on implementation.
16.4 Basic Computational Programs
In XLPack (basic functions), the following fundamental computational programs are available for matrices in CSR format:
CsrDusmv y <- αAx + βy or y <- αA^Tx + βy CsrDussv Solves Ax = b or A^Tx = b (triangular matrix) CsrDusmm C <- αAB + βC or C <- αA^TB + βC CsrDussm Solves AX = B or A^TX = B (triangular matrix) SsrDusmv y <- αAx + βy (symmetric matrix) CsrTrans Transpose of a sparse matrix
Additionally, the following conversion programs are available:
CooCsr COO -> CSR CsrCoo CSR -> COO SsrCsr SSR (CSR symmetric matrix format) -> CSR (full symmetric matrix) CsrSsr CSR (full symmetric matrix) -> SSR (CSR symmetric matrix format) CsrDense CSR -> Dense matrix DenseCsr Dense matrix -> CSR
16.5 Sparse Matrix Checking
When using matrices in CSR format, special care must be taken to avoid errors related to pointers and indices. These values correspond to variable addresses, meaning that incorrect handling can lead to unintended memory reads or writes, causing access errors.
In VBA, pointers are generally not used, and array bounds are always checked, preventing crashes. However, errors in CSR matrix pointers or indices can still cause Excel to crash, so extra caution is required during program development.
To help reduce such errors, the sparse matrix checking program CsrCheck is provided. While it does not completely eliminate errors, it can be a useful tool. An example of its usage is shown below.
Sub Start_2()
Const M = 3, N = 4
Dim Val(M * N) As Double, Ptr(M) As Long, Ind(M * N) As Long
Dim Result(9) As Long, Info As Long
'-- Check (1)
Debug.Print "** Check 1 **"
Val(0) = 1: Val(1) = -2: Val(2) = 4: Val(3) = -1: Val(4) = 2
Val(5) = -9: Val(6) = 5: Val(7) = -8
Ptr(0) = 0: Ptr(1) = 3: Ptr(2) = 6: Ptr(3) = 8
Ind(0) = 0: Ind(1) = 2: Ind(2) = 3: Ind(3) = 0: Ind(4) = 1
Ind(5) = 3: Ind(6) = 1: Ind(7) = 3
Call CsrCheck(M, N, Val(), Ptr(), Ind(), Result(), Info)
Call PrintResult(Info, Result())
'-- Check (2)
Debug.Print "** Check 2 **"
Val(0) = 0: Val(1) = 4: Val(2) = -2: Val(3) = -1: Val(4) = 2
Val(5) = -9: Val(6) = 5: Val(7) = -8
Ptr(0) = 0: Ptr(1) = 3: Ptr(2) = 6: Ptr(3) = 8
Ind(0) = 0: Ind(1) = 3: Ind(2) = 2: Ind(3) = 0: Ind(4) = 1
Ind(5) = 3: Ind(6) = 1: Ind(7) = 5
Call CsrCheck(M, N, Val(), Ptr(), Ind(), Result(), Info)
Call PrintResult(Info, Result())
End Sub
Sub PrintResult(Info As Long, Result() As Long)
Debug.Print "Info =" + Str(Info)
Debug.Print " Sym =" + Str(Result(0))
Debug.Print " Nnz =" + Str(Result(1));
Debug.Print " Nnz(L) =" + Str(Result(2));
Debug.Print " Nnz(U) =" + Str(Result(3));
Debug.Print " Nnz(D) =" + Str(Result(4))
Debug.Print " Zero =" + Str(Result(5));
Debug.Print " Zero(D) =" + Str(Result(6))
Debug.Print " Empty rows =" + Str(Result(7));
Debug.Print " Unsorted rows =" + Str(Result(8))
Debug.Print " Invalid inds =" + Str(Result(9))
End SubThe output is as folows:
** Check 1 ** Info = 0 Sym = 0 Nnz = 8 Nnz(L) = 2 Nnz(U) = 4 Nnz(D) = 2 Zero = 0 Zero(D) = 0 Empty rows = 0 Unsorted rows = 0 Invalid inds = 0 ** Check 2 ** Info = 14 Sym = 0 Nnz = 8 Nnz(L) = 2 Nnz(U) = 4 Nnz(D) = 2 Zero = 1 Zero(D) = 1 Empty rows = 0 Unsorted rows = 1 Invalid inds = 1
Check 2 found one occurrence of an element with value 0 (harmless), one occurrence of rows not in ascending order (causing malfunctions in basic arithmetic programs and many other applications), and one occurrence of an invalid index (causing an access error).
16.6 Systems of Linear Equations
16.6.1 Solving Systems of Linear Equations
In sparse matrix linear calculations, iterative methods (such as Krylov subspace methods) are primarily used for solving systems of linear equations.
LU decomposition and Cholesky decomposition, commonly used for solving general systems of linear equations, allow coefficient matrices to have zero elements, but these elements may become nonzero as calculations progress. In contrast, iterative methods maintain zero elements in the coefficient matrix throughout the computation, making them well-suited for sparse matrix storage formats.
16.6.1.1 Conjugate Gradient (CG) Method
The conjugate gradient (CG) method is commonly used when the coefficient matrix A is symmetric and positive definite.
The solution to the equation Ax = b minimizes the following function:
f(x) = (1/2)(x, Ax) - (x, b)
In other words, solving this equation is equivalent to solving a nonlinear optimization problem with this as the objective function. The gradient descent method is applied (see Chapter 10).
Let pk be the direction vector and αk be the step size. The gradient descent method follows the following iterative formula:
xk+1 = xk + αkpk
To minimize f(xk+1), the step size and direction vector are determined as follows:
αk = (pk, rk)/(pk, Apk) rk+1 = rk - αkApk βk = -(rk+1, Apk)/(pk, Apk) pk+1 = rk+1 + βkpk
Here, rk represents the residual b – Axk.
Starting with an initial approximation x0, the iteration begins with p0 = r0 = b – Ax0.
In theory, this iterative process converges in at most n iterations. However, due to computational errors, it may not always do so in practice.
There is another interpretation of the CG method. Without going into detail, the Krylov subspace method generates a specialized space to approximate the solution of systems of linear equations. The CG method can be considered a variation of this approach.
16.6.1.2 Bi-Conjugate Gradient (BICG) Method
The conjugate gradient (CG) method cannot be used when the coefficient matrix A is a general (asymmetric) matrix. Therefore, various alternative methods have been proposed, with the bi-conjugate gradient (BICG) method being one of the most promising solutions.
Consider the equation Ax = b along with its dual equation ATx* = b*.
Using these, we define the following symmetric matrix and vectors:
A~ = ( A 0 )
( 0 AT )
x~ = ( x )
( x* )
b~ = ( b )
( b* )
H~ = ( 0 I )
( I 0 )
Next, we define the H~ inner product as follows:
<x~, y~> = (x~, H~y~) = (H~x~, y~)
Considering the equation A~x~ = b~, its solution minimizes the following function:
f(x~) = (1/2)<x~, A~x~> - <x~, b~>
By applying the computational procedures of the CG method to this problem, we derive the BICG method.
The BICG method can also be considered a variation of the Krylov subspace method.
16.6.2 Solving Systems of Linear Equations Using XLPack
In XLPack (basic functions), systems of linear equations can be solved using the conjugate gradient (CG) method for symmetric matrices (subroutine Cg1) and the bi-conjugate gradient (BICG) method for general (asymmetric) matrices (subroutine Bicg1).
The following example demonstrates solving the previously discussed Laplace equation using a five-point finite difference approximation. Since the coefficient matrix is symmetric, the CG method is used.
16.6.2.1 Generating the Coefficient Matrix and Right-Hand Side Vector
Following the problem setup, the coefficient matrix A and the right-hand side vector b are prepared accordingly.
Sub Mat5DF(M As Long, Val() As Double, Ptr() As Long, Ind() As Long)
Dim MM As Long, II As Long, I As Long, J As Long, K As Long
MM = M * M
K = 0
For II = 0 To MM - 1
I = Int(II / M)
J = II - I * M
Ptr(II) = K
If I > 0 Then
Ind(K) = II - M
Val(K) = -1
K = K + 1
End If
If J > 0 Then
Ind(K) = II - 1
Val(K) = -1
K = K + 1
End If
Ind(K) = II
Val(K) = 4
K = K + 1
If J < M - 1 Then
Ind(K) = II + 1
Val(K) = -1
K = K + 1
End If
If I < M - 1 Then
Ind(K) = II + M
Val(K) = -1
K = K + 1
End If
Next
Ptr(MM) = K
End Sub
Sub Rhs5DF(M As Long, B() As Double)
Dim H As Double, MM As Long, I As Long
H = 1# / (M + 1)
MM = M * M
For I = 0 To MM - 1
B(I) = 0
Next
For I = M - 1 To MM - 1 Step M
B(I) = H * (Int(I / M) + 1)
Next
For I = M * (M - 1) To MM - 1
B(I) = B(I) + H * (I - M * (M - 1) + 1)
Next
End Sub16.6.2.2 Computing the Solution Using the CG Method
Once A and b are prepared, the subroutine Cg1 is called to obtain the solution. The initial values are all set to zero. After computing the solution, it is output in an appropriate format.
Sub Start()
Const N = 7, M = N - 1, MM = M * M
Dim Val(5 * MM) As Double, Ptr(MM) As Long, Ind(5 * MM) As Long
Dim B(MM - 1) As Double, X(MM - 1) As Double
Dim Resid As Double, Iter As Long, Info As Long
Dim I As Long, J As Long
'-- Set coefficients
Call Mat5DF(M, Val(), Ptr(), Ind())
Call Rhs5DF(M, B())
'-- Set initial approximation
For I = 0 To MM - 1
X(I) = 0
Next
'-- Solve equations
Call Cg1(MM, Val(), Ptr(), Ind(), B(), X(), Info, Iter, Resid)
'-- Output result
If Info = 0 Then
Call Output(N, X())
Else
MsgBox ("Error in Cg1: Info =" + Str(Info))
End If
End SubWhen this program is executed, the following results are obtained. Here, the values of all lattice points, including boundary points, are displayed.
16.7 File Input and Output
One of the commonly used formats for storing sparse matrices as text files is the MM format (Matrix Market Exchange Format). This format is often used for distributing test matrices.
In XLPack (basic functionality), the following input and output programs are provided for MM format files:
MMRead Reads a Matrix Market format file MMReadInfo Reads matrix information from a Matrix Market format file MMWrite Writes to a Matrix Market format file
Appendix: Computational Example Using the Finite Element Method
XLPack includes programs related to partial differential equations as part of its sparse matrix computation features (currently in an experimental version and subject to change in the future). Here, we will use the finite element method program to solve the above example problem of the Laplace equation.
This example uses a simple mesh with triangular elements. By inputting the information of lattice points, finite elements, and boundary elements into Fem2p, it outputs a finite element matrix in CSR format along with a right-hand side vector. Since this is a symmetric matrix, solving it using the CG method yields a numerical solution to the partial differential equation.
Sub Start()
Const Nx = 7, Ny = 7, N = (Nx + 1) * (Ny + 1), Ne = 2 * Nx * Ny
Const Nb = 2 * (Nx + Ny), Nb2 = 0, LdKnc = 4, LdKs = 3
Dim X(N - 1) As Double, Y(N - 1) As Double
Dim P(N - 1) As Double, Q(N - 1) As Double, F(N - 1) As Double
Dim Knc(LdKnc - 1, Ne - 1) As Long, Ks(LdKs - 1, Nb - 1) As Long
Dim Lb(Nb - 1) As Long, Ib(Nb - 1) As Long, Bv(Nb - 1) As Double
Dim Ks2() As Long, Alpha() As Double, Beta() As Double
Dim Val(20 * N) As Double, Ptr(N) As Long, Ind(20 * N) As Long
Dim B(N - 1) As Double, U(N - 1) As Double
Dim Ux As Double, Err As Double
Dim Info As Long
Dim I As Long, J As Long
'-- Set mesh data
Call SetData(Nx, Ny, N, Ne, Nb, X(), Y(), Knc(), Ks(), Lb(), P(), Q(), F(), Ib(), Bv())
'-- Assemble FEM matrix
Call Fem2p(N, Ne, X(), Y(), Knc(), P(), Q(), F(), Nb, Ib(), Bv(), Nb2, Ks2(), Alpha(), Beta(), Val(), Ptr(), Ind(), B(), Info)
'-- Solve equation
Dim Iter As Long, Res As Double
Call Cg1(N, Val(), Ptr(), Ind(), B(), U(), Info, Iter, Res)
Debug.Print "Info =" + Str(Info) + ", Iter =" + Str(Iter) + ", Res =" + Str(Res)
'-- Display solution
For I = 0 To Nx
For J = 0 To Ny
Cells(4 + Nx - J, 1 + I) = U((Ny + 1) * J + I)
Next
Next
End Sub
Sub SetData(Nx As Long, Ny As Long, N As Long, Ne As Long, Nb As Long, X() As Double, Y() As Double, Knc() As Long, Ks() As Long, Lb() As Long, P() As Double, Q() As Double, F() As Double, Ib() As Long, Bv() As Double)
Dim Id() As Long
Dim I As Long, J As Long, K As Long
Call Mesh23(Nx, Ny, X(), Y(), Knc(), Ks(), Lb())
For I = 0 To N - 1
P(I) = 1
Q(I) = 0
F(I) = 0
Next
'-- Set boundary condition 1
ReDim Id(N - 1)
For I = 0 To Nb - 1
For J = 1 To 2
Id(Ks(J, I) - 1) = Lb(I)
Next
Next
K = 0
For I = 0 To N - 1
If Id(I) = 1 Or Id(I) = 4 Then
Bv(K) = 0
ElseIf Id(I) = 2 Then
Bv(K) = Y(I)
ElseIf Id(I) = 3 Then
Bv(K) = X(I)
Else
GoTo Continue
End If
Ib(K) = I + 1
K = K + 1
Continue:
Next
End SubBy executing this program, you can obtain the same results as those obtained using the difference method (within the margin of error).
As a function related to partial differential equations, the output program for VTK files, WriteVtkug, is provided. The final part of the above program is modified using this as below.
Sub Start2()
Const Nx = 7, Ny = 7, N = (Nx + 1) * (Ny + 1), Ne = 2 * Nx * Ny
Const Nb = 2 * (Nx + Ny), Nb2 = 0, LdKnc = 4, LdKs = 3
Dim X(N - 1) As Double, Y(N - 1) As Double
Dim P(N - 1) As Double, Q(N - 1) As Double, F(N - 1) As Double
Dim Knc(LdKnc - 1, Ne - 1) As Long, Ks(LdKs - 1, Nb - 1) As Long
Dim Lb(Nb - 1) As Long, Ib(Nb - 1) As Long, Bv(Nb - 1) As Double
Dim Ks2() As Long, Alpha() As Double, Beta() As Double
Dim Val(20 * N) As Double, Ptr(N) As Long, Ind(20 * N) As Long
Dim B(N - 1) As Double, U(N - 1) As Double
Dim Ux As Double, Err As Double
Dim Info As Long
Dim I As Long, J As Long
'-- Set mesh data
Call SetData(Nx, Ny, N, Ne, Nb, X(), Y(), Knc(), Ks(), Lb(), P(), Q(), F(), Ib(), Bv())
'-- Assemble FEM matrix
Call Fem2p(N, Ne, X(), Y(), Knc(), P(), Q(), F(), Nb, Ib(), Bv(), Nb2, Ks2(), Alpha(), Beta(), Val(), Ptr(), Ind(), B(), Info)
'-- Solve equation
Dim Iter As Long, Res As Double
Call Cg1(N, Val(), Ptr(), Ind(), B(), U(), Info, Iter, Res)
Debug.Print "Info =" + Str(Info) + ", Iter =" + Str(Iter) + ", Res =" + Str(Res)
'-- Output solution
Dim Z(N - 1) As Double
Call WriteVtkug("Test.vtk", N, X(), Y(), Z(), Ne, Knc(), U(), Info)
Debug.Print "Info =" + Str(Info)
End SubBy executing this program, you can obtain the Test.vtk file. You can input the vtk file into external programs to display the results. For example, when viewed in the ParaView program, it appears as follows:
In XLPack (basic functions), other input/output programs for files generated by the mesh creation software gmsh (ReadGmsh22, WriteGmsh22) are provided. Using these allows calculations in various shaped regions created with gmsh. Below is a simple example of calculation results.
