linsolve: add linsolve package
diff --git a/linsolve/bicg.go b/linsolve/bicg.go
new file mode 100644
index 0000000..042f75e
--- /dev/null
+++ b/linsolve/bicg.go
@@ -0,0 +1,128 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "math"
+
+ "gonum.org/v1/gonum/mat"
+)
+
+// BiCG implements the Bi-Conjugate Gradient method with preconditioning for
+// solving systems of linear equations
+// A * x = b,
+// where A is a nonsymmetric, nonsingular matrix. It uses limited memory storage
+// but the convergence may be irregular and the method may break down. BiCG
+// requires a multiplication with A and Aᵀ at each iteration. BiCGStab is a
+// related method that does not use Aᵀ.
+//
+// References:
+// - Barrett, R. et al. (1994). Section 2.3.5 BiConjugate Gradient (BiCG).
+// In Templates for the Solution of Linear Systems: Building Blocks
+// for Iterative Methods (2nd ed.) (pp. 19-20). Philadelphia, PA: SIAM.
+// Retrieved from http://www.netlib.org/templates/templates.pdf
+type BiCG struct {
+ x mat.VecDense
+ r, rt mat.VecDense
+ p, pt mat.VecDense
+ z, zt mat.VecDense
+
+ rho, rhoPrev float64
+
+ resume int
+}
+
+// Init initializes the data for a linear solve. See the Method interface for more details.
+func (b *BiCG) Init(x, residual *mat.VecDense) {
+ dim := x.Len()
+ if residual.Len() != dim {
+ panic("bicg: vector length mismatch")
+ }
+
+ b.x.CloneVec(x)
+ b.r.CloneVec(residual)
+ b.rt.CloneVec(&b.r)
+
+ b.p.Reset()
+ b.p.ReuseAsVec(dim)
+ b.pt.Reset()
+ b.pt.ReuseAsVec(dim)
+ b.z.Reset()
+ b.z.ReuseAsVec(dim)
+ b.zt.Reset()
+ b.zt.ReuseAsVec(dim)
+
+ b.rhoPrev = 1
+
+ b.resume = 1
+}
+
+// Iterate performs an iteration of the linear solve. See the Method interface for more details.
+//
+// BiCG will command the following operations:
+// MulVec
+// MulVec|Trans
+// PreconSolve
+// PreconSolve|Trans
+// CheckResidualNorm
+// MajorIteration
+// NoOperation
+func (b *BiCG) Iterate(ctx *Context) (Operation, error) {
+ switch b.resume {
+ case 1:
+ // Solve M^{-1} * r_{i-1}.
+ ctx.Src.CopyVec(&b.r)
+ b.resume = 2
+ return PreconSolve, nil
+ case 2:
+ b.z.CopyVec(ctx.Dst)
+ // Solve M^{-T} * rt_{i-1}.
+ ctx.Src.CopyVec(&b.rt)
+ b.resume = 3
+ return PreconSolve | Trans, nil
+ case 3:
+ b.zt.CopyVec(ctx.Dst)
+ b.rho = mat.Dot(&b.z, &b.rt)
+ if math.Abs(b.rho) < breakdownTol {
+ b.resume = 0
+ return NoOperation, &BreakdownError{math.Abs(b.rho), breakdownTol}
+ }
+ beta := b.rho / b.rhoPrev
+ b.p.AddScaledVec(&b.z, beta, &b.p)
+ b.pt.AddScaledVec(&b.zt, beta, &b.pt)
+ // Compute A * p.
+ ctx.Src.CopyVec(&b.p)
+ b.resume = 4
+ return MulVec, nil
+ case 4:
+ // z is overwritten and reused.
+ b.z.CopyVec(ctx.Dst)
+ // Compute Aᵀ * pt.
+ ctx.Src.CopyVec(&b.pt)
+ b.resume = 5
+ return MulVec | Trans, nil
+ case 5:
+ // zt is overwritten and reused.
+ b.zt.CopyVec(ctx.Dst)
+ alpha := b.rho / mat.Dot(&b.pt, &b.z)
+ ctx.X.AddScaledVec(ctx.X, alpha, &b.p)
+ b.rt.AddScaledVec(&b.rt, -alpha, &b.zt)
+ b.r.AddScaledVec(&b.r, -alpha, &b.z)
+ ctx.ResidualNorm = mat.Norm(&b.r, 2)
+ b.resume = 6
+ return CheckResidualNorm, nil
+ case 6:
+ if ctx.Converged {
+ b.resume = 0
+ return MajorIteration, nil
+ }
+ b.rhoPrev = b.rho
+ b.resume = 1
+ return MajorIteration, nil
+
+ default:
+ panic("bicg: Init not called")
+ }
+}
diff --git a/linsolve/bicgstab.go b/linsolve/bicgstab.go
new file mode 100644
index 0000000..e9046b4
--- /dev/null
+++ b/linsolve/bicgstab.go
@@ -0,0 +1,153 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "math"
+
+ "gonum.org/v1/gonum/mat"
+)
+
+// BiCGStab implements the BiConjugate Gradient Stabilized method with
+// preconditioning for solving systems of linear equations
+// A * x = b,
+// where A is a nonsymmetric, nonsingular matrix. The method is a variant of
+// BiCG but offers a smoother convergence and does not require multiplication
+// with Aᵀ.
+//
+// References:
+// - Barrett, R. et al. (1994). Section 2.3.8 BiConjugate Gradient Stabilized (Bi-CGSTAB).
+// In Templates for the Solution of Linear Systems: Building Blocks
+// for Iterative Methods (2nd ed.) (pp. 24-25). Philadelphia, PA: SIAM.
+// Retrieved from http://www.netlib.org/templates/templates.pdf
+type BiCGStab struct {
+ x mat.VecDense
+ r, rt mat.VecDense
+ p mat.VecDense
+ phat mat.VecDense
+ shat mat.VecDense
+ t mat.VecDense
+ v mat.VecDense
+
+ rho, rhoPrev float64
+ alpha float64
+ omega float64
+
+ resume int
+}
+
+// Init initializes the data for a linear solve. See the Method interface for more details.
+func (b *BiCGStab) Init(x, residual *mat.VecDense) {
+ dim := x.Len()
+ if residual.Len() != dim {
+ panic("bicgstab: vector length mismatch")
+ }
+
+ b.x.CloneVec(x)
+ b.r.CloneVec(residual)
+ b.rt.CloneVec(&b.r)
+
+ b.p.Reset()
+ b.p.ReuseAsVec(dim)
+ b.phat.Reset()
+ b.phat.ReuseAsVec(dim)
+ b.shat.Reset()
+ b.shat.ReuseAsVec(dim)
+ b.t.Reset()
+ b.t.ReuseAsVec(dim)
+ b.v.Reset()
+ b.v.ReuseAsVec(dim)
+
+ b.rhoPrev = 1
+ b.alpha = 0
+ b.omega = 1
+
+ b.resume = 1
+}
+
+// Iterate performs an iteration of the linear solve. See the Method interface for more details.
+//
+// BiCGStab will command the following operations:
+// MulVec
+// PreconSolve
+// CheckResidualNorm
+// MajorIteration
+// NoOperation
+func (b *BiCGStab) Iterate(ctx *Context) (Operation, error) {
+ switch b.resume {
+ case 1:
+ b.rho = mat.Dot(&b.rt, &b.r)
+ if math.Abs(b.rho) < breakdownTol {
+ b.resume = 0
+ return NoOperation, &BreakdownError{math.Abs(b.rho), breakdownTol}
+ }
+ // p_i = r_{i-1} + beta*(p_{i-1} - omega * v_{i-1})
+ beta := (b.rho / b.rhoPrev) * (b.alpha / b.omega)
+ b.p.AddScaledVec(&b.p, -b.omega, &b.v)
+ b.p.AddScaledVec(&b.r, beta, &b.p)
+ // Solve M^{-1} * p_i.
+ ctx.Src.CopyVec(&b.p)
+ b.resume = 2
+ return PreconSolve, nil
+ case 2:
+ b.phat.CopyVec(ctx.Dst)
+ // Compute A * \hat{p}_i.
+ ctx.Src.CopyVec(&b.phat)
+ b.resume = 3
+ return MulVec, nil
+ case 3:
+ b.v.CopyVec(ctx.Dst)
+ rtv := mat.Dot(&b.rt, &b.v)
+ if rtv == 0 {
+ b.resume = 0
+ return NoOperation, &BreakdownError{}
+ }
+ b.alpha = b.rho / rtv
+ // Form the residual and X so that we can check for tolerance early.
+ ctx.X.AddScaledVec(ctx.X, b.alpha, &b.phat)
+ b.r.AddScaledVec(&b.r, -b.alpha, &b.v)
+ ctx.ResidualNorm = mat.Norm(&b.r, 2)
+ b.resume = 4
+ return CheckResidualNorm, nil
+ case 4:
+ if ctx.Converged {
+ b.resume = 0
+ return MajorIteration, nil
+ }
+ // Solve M^{-1} * r_i.
+ ctx.Src.CopyVec(&b.r)
+ b.resume = 5
+ return PreconSolve, nil
+ case 5:
+ b.shat.CopyVec(ctx.Dst)
+ // Compute A * \hat{s}_i.
+ ctx.Src.CopyVec(&b.shat)
+ b.resume = 6
+ return MulVec, nil
+ case 6:
+ b.t.CopyVec(ctx.Dst)
+ b.omega = mat.Dot(&b.t, &b.r) / mat.Dot(&b.t, &b.t)
+ ctx.X.AddScaledVec(ctx.X, b.omega, &b.shat)
+ b.r.AddScaledVec(&b.r, -b.omega, &b.t)
+ ctx.ResidualNorm = mat.Norm(&b.r, 2)
+ b.resume = 7
+ return CheckResidualNorm, nil
+ case 7:
+ if ctx.Converged {
+ b.resume = 0
+ return MajorIteration, nil
+ }
+ if math.Abs(b.omega) < breakdownTol {
+ b.resume = 0
+ return NoOperation, &BreakdownError{math.Abs(b.omega), breakdownTol}
+ }
+ b.rhoPrev = b.rho
+ b.resume = 1
+ return MajorIteration, nil
+
+ default:
+ panic("bicgstab: Init not called")
+ }
+}
diff --git a/linsolve/cg.go b/linsolve/cg.go
new file mode 100644
index 0000000..97e13f6
--- /dev/null
+++ b/linsolve/cg.go
@@ -0,0 +1,99 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "gonum.org/v1/gonum/mat"
+)
+
+// CG implements the Conjugate Gradient iterative method with preconditioning
+// for solving systems of linear equations
+// A * x = b,
+// where A is a symmetric positive definite matrix. It requires minimal memory
+// storage and is a good choice for symmetric positive definite problems.
+//
+// References:
+// - Barrett, Richard et al. (1994). Section 2.3.1 Conjugate Gradient Method (CG).
+// In Templates for the Solution of Linear Systems: Building Blocks for
+// Iterative Methods (2nd ed.) (pp. 12-15). Philadelphia, PA: SIAM.
+// Retrieved from http://www.netlib.org/templates/templates.pdf
+// - Hestenes, M., and Stiefel, E. (1952). Methods of conjugate gradients for
+// solving linear systems. Journal of Research of the National Bureau of
+// Standards, 49(6), 409. doi:10.6028/jres.049.044
+// - Málek, J. and Strakoš, Z. (2015). Preconditioning and the Conjugate Gradient
+// Method in the Context of Solving PDEs. Philadelphia, PA: SIAM.
+type CG struct {
+ x mat.VecDense
+ r mat.VecDense
+ p mat.VecDense
+
+ rho, rhoPrev float64
+
+ resume int
+}
+
+// Init initializes the data for a linear solve. See the Method interface for more details.
+func (cg *CG) Init(x, residual *mat.VecDense) {
+ dim := x.Len()
+ if residual.Len() != dim {
+ panic("cg: vector length mismatch")
+ }
+
+ cg.x.CloneVec(x)
+ cg.r.CloneVec(residual)
+
+ cg.p.Reset()
+ cg.p.ReuseAsVec(dim)
+
+ cg.rhoPrev = 1
+
+ cg.resume = 1
+}
+
+// Iterate performs an iteration of the linear solve. See the Method interface for more details.
+//
+// CG will command the following operations:
+// MulVec
+// PreconSolve
+// CheckResidualNorm
+// MajorIteration
+func (cg *CG) Iterate(ctx *Context) (Operation, error) {
+ switch cg.resume {
+ case 1:
+ ctx.Src.CopyVec(&cg.r)
+ cg.resume = 2
+ // Compute z_{i-1} = M^{-1} * r_{i-1}.
+ return PreconSolve, nil
+ case 2:
+ z := ctx.Dst
+ cg.rho = mat.Dot(&cg.r, z) // ρ_{i-1} = r_{i-1} · z_{i-1}
+ beta := cg.rho / cg.rhoPrev // β_{i-1} = ρ_{i-1} / ρ_{i-2}
+ cg.p.AddScaledVec(z, beta, &cg.p) // p_i = z_{i-1} + β p_{i-1}
+ ctx.Src.CopyVec(&cg.p)
+ cg.resume = 3
+ // Compute A * p_i.
+ return MulVec, nil
+ case 3:
+ ap := ctx.Dst
+ alpha := cg.rho / mat.Dot(&cg.p, ap) // α_i = ρ_{i-1} / (p_i · A p_i)
+ cg.x.AddScaledVec(&cg.x, alpha, &cg.p) // x_i = x_{i-1} + α p_i
+ cg.r.AddScaledVec(&cg.r, -alpha, ap) // r_i = r_{i-1} - α A p_i
+ ctx.ResidualNorm = mat.Norm(&cg.r, 2)
+ cg.resume = 4
+ return CheckResidualNorm, nil
+ case 4:
+ ctx.X.CopyVec(&cg.x)
+ if ctx.Converged {
+ cg.resume = 0
+ return MajorIteration, nil
+ }
+ cg.rhoPrev = cg.rho
+ cg.resume = 1
+ return MajorIteration, nil
+
+ default:
+ panic("cg: Init not called")
+ }
+}
diff --git a/linsolve/doc.go b/linsolve/doc.go
new file mode 100644
index 0000000..d4b3364
--- /dev/null
+++ b/linsolve/doc.go
@@ -0,0 +1,118 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+/*
+Package linsolve provides iterative methods for solving linear systems.
+
+Background
+
+A system of linear equations can be written as
+
+ A * x = b,
+
+where A is a given n×n non-singular matrix, b is a given n-vector (the
+right-hand side), and x is an unknown n-vector.
+
+Direct methods such as the LU or QR decomposition compute (in the absence of
+roundoff errors) the exact solution after a finite number of steps. For a
+general matrix A they require O(n^3) arithmetic operations which becomes
+infeasible for large n due to excessive memory and time cost.
+
+Iterative methods, in contrast, generally do not compute the exact solution x.
+Starting from an initial estimate x_0, they instead compute a sequence x_i of
+increasingly accurate approximations to x. This iterative process is stopped
+when the estimated difference between x_i and the true x becomes smaller than a
+prescribed threshold. The number of iterations thus depends on the value of the
+threshold. If the desired threshold is very small, then the iterative methods
+may take as long or longer than a direct method. However, for many problems a
+decent approximation can be found in a small number of steps.
+
+The iterative methods implemented in this package do not access the elements of
+A directly, they instead ask for the result of matrix-vector products with A.
+For a general n×n matrix this requires O(n^2) operations, but can be much
+cheaper depending on the structure of the matrix (sparse, banded,
+block-stuctured, etc.). Such structure often arises in practical applications.
+An iterative method can thus be significantly cheaper than a direct method, by
+using a small number of iterations and taking advantage of matrix structure.
+
+Iterative methods are most often useful in the following situations:
+
+ - The system matrix A is sparse, blocked or has other special structure,
+ - The problem size is sufficiently large that a dense factorization of A is
+ costly in terms of compute time and/or memory storage,
+ - Computing the product of A (or A^T, if necessary) with a vector can be done
+ efficiently,
+ - An approximate solution is all that is required.
+
+Using linsolve
+
+The two most important elements of the API are the MulVecToer interface and the
+Iterative function.
+
+MulVecToer interface
+
+The MulVecToer interface represents the system matrix A. This abstracts the
+details of any particular matrix storage, and allows the user to exploit the
+properties of their particular matrix. Matrix types provided by gonum/mat and
+github.com/james-bowman/sparse packages implement this interface.
+
+Note that methods in this package have only limited means for checking whether
+the provided MulVecToer represents a matrix that satisfies all assumptions made
+by the chosen Method, for example if the matrix is actually symmetric positive
+definite.
+
+Iterative function
+
+The Iterative function is the entry point to the functionality provided by this
+package. It takes as parameters the matrix A (via the MulVecToer interface as
+discussed above), the right-hand side vector b, the iterative method and
+settings that control the iterative process and provide a way for reusing
+memory.
+
+Choosing an iterative method
+
+The choice of an iterative method is typically guided by the properties of the
+matrix A including symmetry, definiteness, sparsity, conditioning, and block
+structure. In general, performance on symmetric matrices is well understood (see
+the references below), with the conjugate gradient method being a good starting
+point. Non-symmetric matrices are much more difficult to assess, where any
+suggestion of a 'best' method is usually accompanied by a recommendation to use
+trial-and-error.
+
+Preconditioning
+
+Preconditioning is a family of techniques that attempt to transform the linear
+system into one that has the same solution but more favorable eigenspectrum. The
+transformation matrix is called a preconditioner. A good preconditioner will
+reduce the number of iterations needed to find a good approximate solution
+(hopefully enough to overcome the cost of applying the preconditioning!), and in
+some cases preconditioning is necessary to get any kind of convergence. In
+linsolve a preconditioner is specified by Settings.PreconSolve.
+
+Implementing Method interface
+
+This package allows external implementations of iterative solvers by means of
+the Method interface. It uses a reverse-communication style of API to
+"outsource" operations such as matrix-vector multiplication, preconditioner
+solve or convergence checks to the caller. The caller performs the commanded
+operation and passes the result again to Method. The matrix A and the right-hand
+side b are not directly available to Methods which encourages their cleaner
+implementation. See the documentation for Method, Operation, and Context for
+more information.
+
+References
+
+Further details about computational practice and mathematical theory of
+iterative methods can be found in the following references:
+
+ - Barrett, Richard et al. (1994). Templates for the Solution of Linear Systems:
+ Building Blocks for Iterative Methods (2nd ed.). Philadelphia, PA: SIAM.
+ Retrieved from http://www.netlib.org/templates/templates.pdf
+ - Saad, Yousef (2003). Iterative methods for sparse linear systems (2nd ed.).
+ Philadelphia, PA: SIAM. Retrieved from
+ http://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf
+ - Greenbaum, A. (1997). Iterative methods for solving linear systems.
+ Philadelphia, PA: SIAM.
+*/
+package linsolve
diff --git a/linsolve/gmres.go b/linsolve/gmres.go
new file mode 100644
index 0000000..d1a174b
--- /dev/null
+++ b/linsolve/gmres.go
@@ -0,0 +1,271 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "math"
+
+ "gonum.org/v1/gonum/blas"
+ "gonum.org/v1/gonum/blas/blas64"
+ "gonum.org/v1/gonum/mat"
+)
+
+// GMRES implements the Generalized Minimum Residual method with the modified
+// Gram-Schmidt orthogonalization for solving systems of linear equations
+// A * x = b,
+// where A is a nonsymmetric, nonsingular matrix. It may find a solution in
+// fewer iterations and with fewer matrix-vector products compared to BiCG or
+// BiCGStab at the price of much large memory storage. This implementation uses
+// restarts to limit the memory requirements. GMRES does not need the
+// multiplication with Aᵀ.
+//
+// References:
+// - Barrett, R. et al. (1994). Section 2.3.4 Generalized Minimal Residual
+// (GMRES). In Templates for the Solution of Linear Systems: Building Blocks
+// for Iterative Methods (2nd ed.) (pp. 17-19). Philadelphia, PA: SIAM.
+// Retrieved from http://www.netlib.org/templates/templates.pdf
+// - Saad, Y., and Schultz, M. (1986). GMRES: A generalized minimal residual
+// algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat.
+// Comput., 7(3), 856. doi:10.6028/jres.049.044
+// Retrieved from https://web.stanford.edu/class/cme324/saad-schultz.pdf
+type GMRES struct {
+ // Restart is the restart parameter which limits the computation and
+ // storage costs. It must hold that
+ // 1 <= Restart <= n
+ // where n is the dimension of the problem. If Restart is 0, n will be
+ // used instead. This guarantees convergence of GMRES and increases
+ // robustness. Many specific problems however, particularly for large
+ // n, will benefit in efficiency by setting Restart to
+ // a problem-dependent value less than n.
+ Restart int
+
+ // m is the used value of Restart.
+ m int
+ // v is an n×(m+1) matrix V whose columns form an orthonormal basis of the
+ // Krylov subspace.
+ v mat.Dense
+ // h is an (m+1)×m upper Hessenberg matrix H.
+ h mat.Dense
+ // givs holds Givens rotations that are used to reduce H to upper triangular
+ // form.
+ givs []givens
+
+ x mat.VecDense
+ y mat.VecDense
+ s mat.VecDense
+
+ k int // Loop variable for inner iterations.
+ resume int
+}
+
+// Init initializes the data for a linear solve. See the Method interface for more details.
+func (g *GMRES) Init(x, residual *mat.VecDense) {
+ dim := x.Len()
+ if residual.Len() != dim {
+ panic("gmres: vector length mismatch")
+ }
+
+ g.m = g.Restart
+ if g.m == 0 {
+ g.m = dim
+ }
+ if g.m <= 0 || dim < g.m {
+ panic("gmres: invalid value of Restart")
+ }
+
+ g.v.Reset()
+ g.v.ReuseAs(dim, g.m+1)
+ // Store the residual in the first column of V.
+ g.vcol(0).CopyVec(residual)
+
+ g.h.Reset()
+ g.h.ReuseAs(g.m+1, g.m)
+
+ if cap(g.givs) < g.m {
+ g.givs = make([]givens, g.m)
+ } else {
+ g.givs = g.givs[:g.m]
+ for i := range g.givs {
+ g.givs[i].c = 0
+ g.givs[i].s = 0
+ }
+ }
+
+ g.x.CloneVec(x)
+ g.y.Reset()
+ g.y.ReuseAsVec(g.m + 1)
+ g.s.Reset()
+ g.s.ReuseAsVec(g.m + 1)
+
+ g.resume = 1
+}
+
+// Iterate performs an iteration of the linear solve. See the Method interface for more details.
+//
+// GMRES will command the following operations:
+// MulVec
+// PreconSolve
+// CheckResidualNorm
+// MajorIteration
+// NoOperation
+func (g *GMRES) Iterate(ctx *Context) (Operation, error) {
+ switch g.resume {
+ case 1:
+ // The initial residual is in the first column of V.
+ ctx.Src.CopyVec(g.vcol(0))
+ g.resume = 2
+ // Solve M^{-1} * r_0.
+ return PreconSolve, nil
+ case 2:
+ // v_0 = M^{-1} * r_0
+ v0 := g.vcol(0)
+ v0.CopyVec(ctx.Dst)
+ // Normalize v_0.
+ norm := mat.Norm(v0, 2)
+ v0.ScaleVec(1/norm, v0)
+ // Initialize s to the elementary vector e_1 scaled by norm.
+ g.s.Zero()
+ g.s.SetVec(0, norm)
+
+ // Begin the inner for-loop for k going from 0 to m-1.
+ g.k = 0
+ fallthrough
+ case 3:
+ ctx.Src.CopyVec(g.vcol(g.k))
+ g.resume = 4
+ // Compute A * v_k.
+ return MulVec, nil
+ case 4:
+ ctx.Src.CopyVec(ctx.Dst)
+ g.resume = 5
+ // Solve M^{-1} * (A * v_k).
+ return PreconSolve, nil
+ case 5:
+ // v_{k+1} = M^{-1} * (A * v_k)
+ vk1 := g.vcol(g.k + 1)
+ vk1.CopyVec(ctx.Dst)
+ // Construct the k-th column of the upper Hessenberg matrix H
+ // using the modified Gram-Schmidt process to make v_{k+1}
+ // orthonormal to the first k+1 columns of V.
+ g.modifiedGS(g.k, &g.h, &g.v, vk1)
+ // Reduce H back to upper triangular form and update the vector s.
+ g.qr(g.k, g.givs, &g.h, &g.s)
+ // Check the approximate residual norm.
+ ctx.ResidualNorm = math.Abs(g.s.AtVec(g.k + 1))
+ g.resume = 6
+ return CheckResidualNorm, nil
+ case 6:
+ g.k++
+ if g.k < g.m && !ctx.Converged {
+ // Continue the inner for-loop.
+ g.resume = 3
+ return NoOperation, nil
+ }
+ // Either restarting or converged, we have to update the solution.
+ // Solve the upper triangular system H*y=s.
+ g.solveLeastSquares(g.k, &g.y, &g.h, &g.s)
+ // Compute x as a linear combination of columns of V.
+ g.updateSolution(g.k, &g.x, &g.v, &g.y)
+ ctx.X.CopyVec(&g.x)
+ if ctx.Converged {
+ g.resume = 0
+ return MajorIteration, nil
+ }
+ // We are restarting, so we have to also compute the residual.
+ g.resume = 7
+ return ComputeResidual, nil
+ case 7:
+ // Store the residual again in the first column of V.
+ g.vcol(0).CopyVec(ctx.Dst)
+ g.resume = 1
+ return MajorIteration, nil
+
+ default:
+ panic("gmres: Init not called")
+ }
+}
+
+// modifiedGS orthonormalizes the vector w with respect to the first k+1 columns
+// of V using the modified Gram-Schmidt algorithm, and stores the computed
+// coefficients in the k-th column of H.
+func (g *GMRES) modifiedGS(k int, h, v *mat.Dense, w *mat.VecDense) {
+ hk := h.ColView(k).(*mat.VecDense)
+ for j := 0; j <= k; j++ {
+ vj := v.ColView(j).(*mat.VecDense)
+ hkj := mat.Dot(vj, w)
+ hk.SetVec(j, hkj) // H[j,k] = v_j · w
+ w.AddScaledVec(w, -hkj, vj) // w -= H[j,k] * v_j
+ }
+ norm := mat.Norm(w, 2)
+ hk.SetVec(k+1, norm) // H[k+1,k] = |w|
+ w.ScaleVec(1/norm, w) // Normalize w.
+}
+
+// qr applies previous Givens rotations to the k-th column of H, computes the
+// next Givens rotation to zero out H[k+1,k] and applies it also to the vector s.
+func (g *GMRES) qr(k int, givs []givens, h *mat.Dense, s *mat.VecDense) {
+ // Apply previous Givens rotations to the k-th column of H.
+ hk := h.ColView(k).(*mat.VecDense)
+ for i, giv := range givs[:k] {
+ hki, hki1 := giv.apply(hk.AtVec(i), hk.AtVec(i+1))
+ hk.SetVec(i, hki)
+ hk.SetVec(i+1, hki1)
+ }
+
+ // Compute the k-th Givens rotation that zeros H[k+1,k].
+ givs[k].c, givs[k].s, _, _ = blas64.Rotg(hk.AtVec(k), hk.AtVec(k+1))
+
+ // Apply the k-th Givens rotation to (H[k,k], H[k+1,k]).
+ hkk, _ := givs[k].apply(hk.AtVec(k), hk.AtVec(k+1))
+ hk.SetVec(k, hkk)
+ // We don't call hk.SetVec(k+1, 0) because the element won't be accessed.
+
+ // Apply the k-th Givens rotation to (s[k], s[k+1]).
+ sk, sk1 := givs[k].apply(s.AtVec(k), s.AtVec(k+1))
+ s.SetVec(k, sk)
+ s.SetVec(k+1, sk1)
+}
+
+// solveLeastSquares solves the upper triangular linear system
+// H * y = s
+func (g *GMRES) solveLeastSquares(k int, y *mat.VecDense, h *mat.Dense, s *mat.VecDense) {
+ // Copy the first k elements of s into y.
+ y.CopyVec(s.SliceVec(0, k))
+ // Convert H into an upper triangular matrix.
+ hraw := h.RawMatrix()
+ htri := blas64.Triangular{
+ Uplo: blas.Upper,
+ Diag: blas.NonUnit,
+ N: k,
+ Data: hraw.Data,
+ Stride: hraw.Stride,
+ }
+ // Solve the upper triangular system storing the solution into y.
+ blas64.Trsv(blas.NoTrans, htri, y.RawVector())
+}
+
+// updateSolution updates the current solution vector x with a linear
+// combination of the first k columns of V:
+// x = x + V * y = x + \sum y_j * v_j
+func (g *GMRES) updateSolution(k int, x *mat.VecDense, v *mat.Dense, y *mat.VecDense) {
+ for j := 0; j < k; j++ {
+ vj := v.ColView(j)
+ x.AddScaledVec(x, y.AtVec(j), vj)
+ }
+}
+
+// vcol returns a view of the j-th column of the matrix V.
+func (g *GMRES) vcol(j int) *mat.VecDense {
+ return g.v.ColView(j).(*mat.VecDense)
+}
+
+// givens is a Givens rotation.
+type givens struct {
+ c, s float64
+}
+
+func (giv *givens) apply(x, y float64) (float64, float64) {
+ return giv.c*x + giv.s*y, giv.c*y - giv.s*x
+}
diff --git a/linsolve/internal/mmarket/reader.go b/linsolve/internal/mmarket/reader.go
new file mode 100644
index 0000000..eaf92d4
--- /dev/null
+++ b/linsolve/internal/mmarket/reader.go
@@ -0,0 +1,105 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package mmarket provides a type to read Matrix Market format files.
+package mmarket
+
+import (
+ "bufio"
+ "errors"
+ "fmt"
+ "io"
+ "strings"
+
+ "gonum.org/v1/exp/linsolve/internal/triplet"
+)
+
+var (
+ errBadFormat = errors.New("mmarket: bad file format")
+ errUnsupported = errors.New("mmarket: matrix type not supported")
+)
+
+type Reader struct {
+ s *bufio.Scanner
+}
+
+func NewReader(r io.Reader) *Reader {
+ return &Reader{
+ s: bufio.NewScanner(r),
+ }
+}
+
+// Read a real matrix in coordinate format and return its triplet representation.
+func (r *Reader) Read() (*triplet.Matrix, error) {
+ r.s.Scan()
+ if err := r.s.Err(); err != nil {
+ return nil, err
+ }
+ header := strings.Fields(r.s.Text())
+ if header[0] != "%%MatrixMarket" {
+ return nil, errBadFormat
+ }
+ if header[1] != "matrix" {
+ return nil, errBadFormat
+ }
+ if header[2] != "coordinate" {
+ return nil, errBadFormat
+ }
+ if header[3] != "real" {
+ return nil, errUnsupported
+ }
+ sym := header[4] == "symmetric"
+
+ var nr, nc, nnz int
+ for r.s.Scan() {
+ line := r.s.Text()
+ if line[0] == '%' {
+ continue
+ }
+ n, err := fmt.Sscan(line, &nr, &nc, &nnz)
+ if err != nil {
+ return nil, err
+ }
+ if n != 3 {
+ return nil, errBadFormat
+ }
+ break
+ }
+ if err := r.s.Err(); err != nil {
+ return nil, err
+ }
+
+ if sym && nr != nc {
+ return nil, errBadFormat
+ }
+
+ m := triplet.NewMatrix(nr, nc)
+ for i := 0; i < nnz; i++ {
+ if !r.s.Scan() {
+ return nil, errBadFormat
+ }
+ var (
+ i, j int
+ v float64
+ )
+ n, err := fmt.Sscan(r.s.Text(), &i, &j, &v)
+ if err != nil {
+ return nil, err
+ }
+ if n != 3 {
+ return nil, errBadFormat
+ }
+ if i < 1 || nr < i {
+ return nil, errBadFormat
+ }
+ if j < 1 || nc < j {
+ return nil, errBadFormat
+ }
+ m.Append(i-1, j-1, v)
+ if sym && i != j {
+ m.Append(j-1, i-1, v)
+ }
+ }
+ return m, nil
+}
diff --git a/linsolve/internal/triplet/triplet.go b/linsolve/internal/triplet/triplet.go
new file mode 100644
index 0000000..e7718d4
--- /dev/null
+++ b/linsolve/internal/triplet/triplet.go
@@ -0,0 +1,74 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package triplet provides triplet representation for sparse matrices.
+package triplet
+
+import "gonum.org/v1/gonum/mat"
+
+type triplet struct {
+ i, j int
+ v float64
+}
+
+type Matrix struct {
+ r, c int
+ data []triplet
+}
+
+func NewMatrix(r, c int) *Matrix {
+ if r <= 0 || c <= 0 {
+ panic("triplet: invalid shape")
+ }
+ return &Matrix{
+ r: r,
+ c: c,
+ }
+}
+
+func (m *Matrix) Dims() (r, c int) {
+ return m.r, m.c
+}
+
+// Append appends a non-zero element to the list of matrix elements without
+// checking whether it already exists.
+func (m *Matrix) Append(i, j int, v float64) {
+ if i < 0 || m.r <= i {
+ panic("triplet: row index out of range")
+ }
+ if j < 0 || m.c <= j {
+ panic("triplet: column index out of range")
+ }
+ if v == 0 {
+ return
+ }
+ m.data = append(m.data, triplet{i, j, v})
+}
+
+func (m *Matrix) MulVecTo(dst *mat.VecDense, trans bool, x mat.Vector) {
+ dst.Zero()
+ if trans {
+ if m.c != dst.Len() || m.r != x.Len() {
+ panic("triplet: dimension mismatch")
+ }
+ for _, aij := range m.data {
+ dst.SetVec(aij.j, dst.AtVec(aij.j)+aij.v*x.AtVec(aij.i))
+ }
+ return
+ }
+ if m.c != x.Len() || m.r != dst.Len() {
+ panic("triplet: dimension mismatch")
+ }
+ for _, aij := range m.data {
+ dst.SetVec(aij.i, dst.AtVec(aij.i)+aij.v*x.AtVec(aij.j))
+ }
+}
+
+func (m *Matrix) DenseCopy() *mat.Dense {
+ d := mat.NewDense(m.r, m.c, nil)
+ for _, aij := range m.data {
+ d.Set(aij.i, aij.j, aij.v)
+ }
+ return d
+}
diff --git a/linsolve/iterative.go b/linsolve/iterative.go
new file mode 100644
index 0000000..c49f1c7
--- /dev/null
+++ b/linsolve/iterative.go
@@ -0,0 +1,266 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "errors"
+
+ "gonum.org/v1/gonum/mat"
+)
+
+const defaultTolerance = 1e-8
+
+// ErrIterationLimit is returned when a maximum number of iterations were done
+// without converging to a solution.
+var ErrIterationLimit = errors.New("linsolve: iteration limit reached")
+
+// MulVecToer represents a square matrix A by means of a matrix-vector
+// multiplication.
+type MulVecToer interface {
+ // MulVecTo computes A*x or Aᵀ*x and stores the result into dst.
+ MulVecTo(dst *mat.VecDense, trans bool, x mat.Vector)
+}
+
+// Settings holds settings for solving a linear system.
+type Settings struct {
+ // InitX holds the initial guess. If it is nil or empty, the zero vector
+ // will be used, otherwise its length must be equal to the dimension of
+ // the system.
+ InitX *mat.VecDense
+
+ // Dst, if not nil, will be used for storing the approximate solution,
+ // otherwise a new vector will be allocated. In both cases the vector will
+ // also be returned in Result. If Dst is not empty, its length must be equal
+ // to the dimension of the system.
+ Dst *mat.VecDense
+
+ // Tolerance specifies error tolerance for the final (approximate)
+ // solution produced by the iterative method. The iteration will be
+ // stopped when
+ // |r_i| < Tolerance * |b|
+ // where r_i is the residual at i-th iteration.
+ //
+ // If Tolerance is zero, a default value of 1e-8 will be used, otherwise
+ // it must be positive and less than 1.
+ Tolerance float64
+
+ // MaxIterations is the limit on the number of iterations. If it is
+ // zero, a default value of twice the dimension of the system will be
+ // used.
+ MaxIterations int
+
+ // PreconSolve describes a preconditioner solve that stores into dst the
+ // solution of the system
+ // M * dst = rhs, or
+ // Mᵀ * dst = rhs,
+ // where M is the preconditioning matrix. If PreconSolve is nil, no
+ // preconditioning will be used (M is the identity).
+ PreconSolve func(dst *mat.VecDense, trans bool, rhs mat.Vector) error
+
+ // Work context can be provided to reduce memory allocation when solving
+ // multiple linear systems. If Work is not nil, its fields must be
+ // either empty or their length must be equal to the dimension of the
+ // system.
+ Work *Context
+}
+
+// defaultSettings fills zero fields of s with default values.
+func defaultSettings(s *Settings, dim int) {
+ if s.InitX != nil && s.InitX.Len() == 0 {
+ s.InitX.ReuseAsVec(dim)
+ }
+ if s.Dst == nil {
+ s.Dst = mat.NewVecDense(dim, nil)
+ } else if s.Dst.Len() == 0 {
+ s.Dst.ReuseAsVec(dim)
+ }
+ if s.Tolerance == 0 {
+ s.Tolerance = defaultTolerance
+ }
+ if s.MaxIterations == 0 {
+ s.MaxIterations = 4 * dim
+ }
+ if s.PreconSolve == nil {
+ s.PreconSolve = NoPreconditioner
+ }
+ if s.Work == nil {
+ s.Work = NewContext(dim)
+ } else {
+ if s.Work.X.Len() == 0 {
+ s.Work.X.ReuseAsVec(dim)
+ }
+ if s.Work.Src.Len() == 0 {
+ s.Work.Src.ReuseAsVec(dim)
+ }
+ if s.Work.Dst.Len() == 0 {
+ s.Work.Dst.ReuseAsVec(dim)
+ }
+ }
+}
+
+func checkSettings(s *Settings, dim int) {
+ if s.InitX != nil && s.InitX.Len() != dim {
+ panic("linsolve: mismatched length of initial guess")
+ }
+ if s.Dst.Len() != dim {
+ panic("linsolve: mismatched destination length")
+ }
+ if s.Tolerance <= 0 || 1 <= s.Tolerance {
+ panic("linsolve: invalid tolerance")
+ }
+ if s.MaxIterations <= 0 {
+ panic("linsolve: negative iteration limit")
+ }
+ if s.Work.X.Len() != dim || s.Work.Src.Len() != dim || s.Work.Dst.Len() != dim {
+ panic("linsolve: mismatched work context length")
+ }
+}
+
+// Result holds the result of an iterative solve.
+type Result struct {
+ // X is the approximate solution.
+ X *mat.VecDense
+
+ // ResidualNorm is an approximation to the norm of the final residual.
+ ResidualNorm float64
+
+ // Stats holds statistics about the iterative solve.
+ Stats Stats
+}
+
+// Stats holds statistics about an iterative solve.
+type Stats struct {
+ // Iterations is the number of iterations performed by Method.
+ Iterations int
+
+ // MulVec is the number of MulVec operations commanded by Method.
+ MulVec int
+
+ // PreconSolve is the number of PreconSolve operations commanded by Method.
+ PreconSolve int
+}
+
+// Iterative finds an approximate solution of the system of n linear equations
+// A*x = b,
+// where A is a nonsingular square matrix of order n and b is the right-hand
+// side vector, using an iterative method m. If m is nil, default GMRES will be
+// used.
+//
+// settings provide means for adjusting parameters of the iterative process. See
+// the Settings documentation for more information. Iterative will not modify
+// the fields of Settings. If settings is nil, default settings will be used.
+//
+// Note that the default choices of Method and Settings were chosen to provide
+// accuracy and robustness, rather than speed. There are many algorithms for
+// iterative linear solutions that have different tradeoffs, and can exploit
+// special structure in the A matrix. Similarly, in many cases the number of
+// iterations can be significantly reduced by using an appropriate
+// preconditioner or increasing the error tolerance. Combined, these choices can
+// significantly reduce computation time. Thus, while Iterative has supplied
+// defaults, users are strongly encouraged to adjust these defaults for their
+// problem.
+func Iterative(a MulVecToer, b *mat.VecDense, m Method, settings *Settings) (*Result, error) {
+ n := b.Len()
+
+ var s Settings
+ if settings != nil {
+ s = *settings
+ }
+ defaultSettings(&s, n)
+ checkSettings(&s, n)
+
+ var stats Stats
+ ctx := s.Work
+ rInit := mat.NewVecDense(n, nil)
+ if s.InitX != nil {
+ // Initial x is provided.
+ ctx.X.CloneVec(s.InitX)
+ computeResidual(rInit, a, b, ctx.X, &stats)
+ } else {
+ // Initial x is the zero vector.
+ ctx.X.Zero()
+ // Residual b-A*x is then equal to b.
+ rInit.CopyVec(b)
+ }
+
+ if m == nil {
+ m = &GMRES{}
+ }
+
+ var err error
+ ctx.ResidualNorm = mat.Norm(rInit, 2)
+ if ctx.ResidualNorm >= s.Tolerance {
+ err = iterate(a, b, rInit, s, m, &stats)
+ } else {
+ s.Dst.CopyVec(ctx.X)
+ }
+
+ return &Result{
+ X: s.Dst,
+ ResidualNorm: ctx.ResidualNorm,
+ Stats: stats,
+ }, err
+}
+
+func iterate(a MulVecToer, b, initRes *mat.VecDense, settings Settings, method Method, stats *Stats) error {
+ bNorm := mat.Norm(b, 2)
+ if bNorm == 0 {
+ bNorm = 1
+ }
+
+ ctx := settings.Work
+ settings.Dst.CopyVec(ctx.X)
+
+ method.Init(ctx.X, initRes)
+ for {
+ op, err := method.Iterate(ctx)
+ if err != nil {
+ return err
+ }
+ switch op {
+ case NoOperation:
+ case MulVec, MulVec | Trans:
+ stats.MulVec++
+ a.MulVecTo(ctx.Dst, op&Trans == Trans, ctx.Src)
+ case PreconSolve, PreconSolve | Trans:
+ stats.PreconSolve++
+ err = settings.PreconSolve(ctx.Dst, op&Trans == Trans, ctx.Src)
+ if err != nil {
+ return err
+ }
+ case CheckResidualNorm:
+ ctx.Converged = ctx.ResidualNorm < settings.Tolerance*bNorm
+ case ComputeResidual:
+ computeResidual(ctx.Dst, a, b, ctx.X, stats)
+ case MajorIteration:
+ stats.Iterations++
+ if ctx.Converged {
+ settings.Dst.CopyVec(ctx.X)
+ return nil
+ }
+ if stats.Iterations == settings.MaxIterations {
+ settings.Dst.CopyVec(ctx.X)
+ return ErrIterationLimit
+ }
+ default:
+ panic("linsolve: invalid operation")
+ }
+ }
+}
+
+// NoPreconditioner implements the identity preconditioner.
+func NoPreconditioner(dst *mat.VecDense, trans bool, rhs mat.Vector) error {
+ if dst.Len() != rhs.Len() {
+ panic("linsolve: mismatched vector length")
+ }
+ dst.CloneVec(rhs)
+ return nil
+}
+
+func computeResidual(dst *mat.VecDense, a MulVecToer, b, x *mat.VecDense, stats *Stats) {
+ stats.MulVec++
+ a.MulVecTo(dst, false, x)
+ dst.AddScaledVec(b, -1, dst)
+}
diff --git a/linsolve/iterative_example_test.go b/linsolve/iterative_example_test.go
new file mode 100644
index 0000000..70b8d54
--- /dev/null
+++ b/linsolve/iterative_example_test.go
@@ -0,0 +1,81 @@
+// Copyright ©2019 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve_test
+
+import (
+ "fmt"
+ "math"
+
+ "gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/linsolve"
+ "gonum.org/v1/gonum/mat"
+)
+
+// System represents a linear system with a symmetric band matrix
+// A*x = b
+type System struct {
+ A *mat.SymBandDense
+ B *mat.VecDense
+}
+
+// L2Projection returns a linear system whose solution is the L2 projection of f
+// into the space of piecewise linear functions defined on the given grid.
+//
+// References:
+// - M. Larson, F. Bengzon, The Finite Element Method: Theory,
+// Implementations, and Applications. Springer (2013), Section 1.3, also
+// available at:
+// http://www.springer.com/cda/content/document/cda_downloaddocument/9783642332869-c1.pdf
+func L2Projection(grid []float64, f func(float64) float64) System {
+ n := len(grid)
+
+ // Assemble the symmetric banded mass matrix by iterating over all elements.
+ A := mat.NewSymBandDense(n, 1, nil)
+ for i := 0; i < n-1; i++ {
+ // h is the length of the i-th element.
+ h := grid[i+1] - grid[i]
+ // Add contribution from the i-th element.
+ A.SetSymBand(i, i, A.At(i, i)+h/3)
+ A.SetSymBand(i, i+1, h/6)
+ A.SetSymBand(i+1, i+1, A.At(i+1, i+1)+h/3)
+ }
+
+ // Assemble the load vector by iterating over all elements.
+ b := mat.NewVecDense(n, nil)
+ for i := 0; i < n-1; i++ {
+ h := grid[i+1] - grid[i]
+ b.SetVec(i, b.AtVec(i)+f(grid[i])*h/2)
+ b.SetVec(i+1, b.AtVec(i+1)+f(grid[i+1])*h/2)
+ }
+
+ return System{A, b}
+}
+
+func ExampleIterative() {
+ const (
+ n = 10
+ x0 = 0.0
+ x1 = 1.0
+ )
+ // Make a uniform grid.
+ grid := make([]float64, n+1)
+ floats.Span(grid, x0, x1)
+ sys := L2Projection(grid, func(x float64) float64 {
+ return x * math.Sin(x)
+ })
+
+ result, err := linsolve.Iterative(sys.A, sys.B, &linsolve.CG{}, nil)
+ if err != nil {
+ fmt.Println("Error:", err)
+ return
+ }
+
+ fmt.Printf("# iterations: %v\n", result.Stats.Iterations)
+ fmt.Printf("Final solution: %.6f\n", mat.Formatted(result.X.T()))
+
+ // Output:
+ // # iterations: 11
+ // Final solution: [-0.003339 0.006677 0.036530 0.085606 0.152981 0.237072 0.337006 0.447616 0.578244 0.682719 0.920847]
+}
diff --git a/linsolve/iterative_test.go b/linsolve/iterative_test.go
new file mode 100644
index 0000000..4456a43
--- /dev/null
+++ b/linsolve/iterative_test.go
@@ -0,0 +1,304 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "math"
+ "testing"
+
+ "golang.org/x/exp/rand"
+
+ "gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/mat"
+)
+
+func TestDefaultMethodDefaultSettings(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ newPDENonsymmetric(16, 16, rnd),
+ newPDEYang47(16, 16, rnd),
+ newPDEYang48(16, 16, rnd),
+ newPDEYang49(16, 16, rnd),
+ newPDEYang410(16, 16, rnd),
+ newPDEYang412(16, 16, rnd),
+ newPDEYang413(16, 16, rnd),
+ newPDEYang414(16, 16, rnd),
+ newPDEYang415(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ testMethodWithSettings(t, nil, nil, tc)
+ }
+}
+
+func TestCG(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ for _, tc := range testCases {
+ s := newTestSettings(rnd, tc)
+ testMethodWithSettings(t, &CG{}, s, tc)
+ }
+}
+
+func TestCGDefaultSettings(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ for _, tc := range testCases {
+ testMethodWithSettings(t, &CG{}, nil, tc)
+ }
+}
+
+func TestBiCG(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ newPDEYang47(16, 16, rnd),
+ newPDEYang48(16, 16, rnd),
+ newPDEYang410(16, 16, rnd),
+ newPDEYang413(16, 16, rnd),
+ newPDEYang414(16, 16, rnd),
+ newPDEYang415(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ s := newTestSettings(rnd, tc)
+ s.Tolerance = 1e-10
+ testMethodWithSettings(t, &BiCG{}, s, tc)
+ }
+}
+
+func TestBiCGDefaultSettings(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ newPDENonsymmetric(16, 16, rnd),
+ newPDEYang47(16, 16, rnd),
+ newPDEYang48(16, 16, rnd),
+ newPDEYang410(16, 16, rnd),
+ newPDEYang413(16, 16, rnd),
+ newPDEYang414(16, 16, rnd),
+ newPDEYang415(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ testMethodWithSettings(t, &BiCG{}, nil, tc)
+ }
+}
+
+func TestBiCGStab(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ s := newTestSettings(rnd, tc)
+ testMethodWithSettings(t, &BiCGStab{}, s, tc)
+ }
+}
+
+func TestBiCGStabDefaultSettings(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ testMethodWithSettings(t, &BiCGStab{}, nil, tc)
+ }
+}
+
+func TestGMRES(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ newPDENonsymmetric(16, 16, rnd),
+ newPDEYang47(16, 16, rnd),
+ newPDEYang48(16, 16, rnd),
+ newPDEYang49(16, 16, rnd),
+ newPDEYang410(16, 16, rnd),
+ newPDEYang412(16, 16, rnd),
+ newPDEYang413(16, 16, rnd),
+ newPDEYang414(16, 16, rnd),
+ newPDEYang415(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ s := newTestSettings(rnd, tc)
+ testMethodWithSettings(t, &GMRES{}, s, tc)
+ }
+}
+
+func TestGMRESDefaultSettings(t *testing.T) {
+ rnd := rand.New(rand.NewSource(1))
+
+ testCases := spdTestCases(rnd)
+ testCases = append(testCases,
+ nonsym3x3(),
+ nonsymTridiag(100),
+ newGreenbaum54(1, 1, rnd),
+ newGreenbaum54(1, 2, rnd),
+ newGreenbaum54(2, 4, rnd),
+ newGreenbaum54(10, 0, rnd),
+ newGreenbaum54(10, 20, rnd),
+ newGreenbaum54(50, 3, rnd),
+ newGreenbaum73(16, 16, rnd),
+ newPDENonsymmetric(16, 16, rnd),
+ newPDEYang47(16, 16, rnd),
+ newPDEYang48(16, 16, rnd),
+ newPDEYang49(16, 16, rnd),
+ newPDEYang410(16, 16, rnd),
+ newPDEYang412(16, 16, rnd),
+ newPDEYang413(16, 16, rnd),
+ newPDEYang414(16, 16, rnd),
+ newPDEYang415(16, 16, rnd),
+ )
+ for _, tc := range testCases {
+ testMethodWithSettings(t, &GMRES{}, nil, tc)
+ }
+}
+
+func newTestSettings(rnd *rand.Rand, tc testCase) *Settings {
+ n := len(tc.b)
+
+ // Initial guess is a random vector.
+ initX := mat.NewVecDense(n, nil)
+ for i := 0; i < initX.Len(); i++ {
+ initX.SetVec(i, rnd.NormFloat64())
+ }
+
+ // Allocate a destination vector and fill it with NaN.
+ dst := mat.NewVecDense(n, nil)
+ for i := 0; i < dst.Len(); i++ {
+ dst.SetVec(i, math.NaN())
+ }
+
+ // Preallocate a work context and fill it with NaN.
+ work := NewContext(n)
+ for i := 0; i < work.X.Len(); i++ {
+ work.X.SetVec(i, math.NaN())
+ work.Src.SetVec(i, math.NaN())
+ work.Dst.SetVec(i, math.NaN())
+ }
+ work.ResidualNorm = math.NaN()
+
+ return &Settings{
+ InitX: initX,
+ Dst: dst,
+ Tolerance: tc.tol,
+ MaxIterations: 5 * n,
+ PreconSolve: tc.PreconSolve,
+ Work: work,
+ }
+}
+
+func testMethodWithSettings(t *testing.T, m Method, s *Settings, tc testCase) {
+ wantTol := 1e-9
+ if s == nil {
+ // The default value of Settings.Tolerance is not as low as the tolerance in
+ // individual test cases, therefore we must use a higher tolerance for
+ // the expected accuracy of the computed solution.
+ wantTol = 1e-7
+ }
+
+ n := len(tc.b)
+
+ b := make([]float64, n)
+ copy(b, tc.b)
+ bVec := mat.NewVecDense(n, b)
+
+ result, err := Iterative(&tc, bVec, m, s)
+ if err != nil {
+ t.Errorf("%v: unexpected error %v", tc.name, err)
+ return
+ }
+
+ if !floats.Equal(tc.b, b) {
+ t.Errorf("%v: unexpected modification of b", tc.name)
+ }
+
+ wantVec := mat.NewVecDense(n, tc.want)
+ var diff mat.VecDense
+ diff.SubVec(result.X, wantVec)
+ dist := mat.Norm(&diff, 2) / mat.Norm(wantVec, 2)
+ if dist > wantTol {
+ t.Errorf("%v: unexpected solution, |want-got|/|want|=%v", tc.name, dist)
+ }
+
+ if s == nil {
+ return
+ }
+
+ if s.MaxIterations > 0 && result.Stats.Iterations > s.MaxIterations {
+ t.Errorf("%v: Result.Stats.Iterations greater than Settings.MaxIterations", tc.name)
+ }
+
+ if s.Dst != nil {
+ if !mat.Equal(s.Dst, result.X) {
+ t.Errorf("%v: Settings.Dst and Result.X not equal\n%v\n%v", tc.name, s.Dst, result.X)
+ }
+ result.X.SetVec(0, 123456.7)
+ if s.Dst.AtVec(0) != result.X.AtVec(0) {
+ t.Errorf("%v: Settings.Dst and Result.X are not the same vector", tc.name)
+ }
+ }
+}
diff --git a/linsolve/linsolve.go b/linsolve/linsolve.go
new file mode 100644
index 0000000..ca95043
--- /dev/null
+++ b/linsolve/linsolve.go
@@ -0,0 +1,155 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "fmt"
+
+ "gonum.org/v1/gonum/mat"
+)
+
+// BreakdownError signifies that a breakdown occured and the method cannot continue.
+type BreakdownError struct {
+ Value float64
+ Tolerance float64
+}
+
+func (e *BreakdownError) Error() string {
+ return fmt.Sprintf("linsolve: breakdown, value=%v tolerance=%v", e.Value, e.Tolerance)
+}
+
+// Method is an iterative method that produces a sequence of vectors that
+// converge to the solution of the system of linear equations
+// A * x = b,
+// where A is non-singular n×n matrix, and x and b are vectors of dimension n.
+//
+// Method uses a reverse-communication interface between the iterative algorithm
+// and the caller. Method acts as a client that commands the caller to perform
+// needed operations via an Operation returned from the Iterate method. This
+// provides independence of Method on representation of the matrix A, and
+// enables automation of common operations like checking for convergence and
+// maintaining statistics.
+type Method interface {
+ // Init initializes the method for solving an n×n
+ // linear system with an initial estimate x and
+ // the corresponding residual vector.
+ //
+ // Method will not retain x or residual.
+ Init(x, residual *mat.VecDense)
+
+ // Iterate performs a step in converging to the
+ // solution of a linear system.
+ //
+ // Iterate retrieves data from Context, updates it,
+ // and returns the next operation. The caller must
+ // perform the Operation using data in Context, and
+ // depending on the state call Iterate again.
+ Iterate(*Context) (Operation, error)
+}
+
+// Context mediates the communication between the Method and
+// the caller. The caller must not modify Context apart from
+// the commanded Operations.
+type Context struct {
+ // X will be set by Method to the current approximate
+ // solution when it commands ComputeResidual and MajorIteration.
+ X *mat.VecDense
+
+ // ResidualNorm is (an estimate of) a norm of
+ // the residual. Method will set it to the current
+ // value when it commands CheckResidualNorm.
+ ResidualNorm float64
+
+ // Converged indicates to Method whether ResidualNorm
+ // satisfies a stopping criterion as a result of
+ // CheckResidualNorm operation.
+ Converged bool
+
+ // Src and Dst are the source and destination vectors
+ // for various Operations. Src will be set by Method
+ // and the caller must store the result in Dst. See
+ // the Operation documentation for more information.
+ Src, Dst *mat.VecDense
+}
+
+// NewContext returns a new Context for work on problems of dimension n.
+// NewContext will panic if n is not positive.
+func NewContext(n int) *Context {
+ if n <= 0 {
+ panic("linsolve: context size is not positive")
+ }
+ return &Context{
+ X: mat.NewVecDense(n, nil),
+ Src: mat.NewVecDense(n, nil),
+ Dst: mat.NewVecDense(n, nil),
+ }
+}
+
+// Reset reinitializes the Context for work on problems of dimension n.
+// Reset will panic if n is not positive.
+func (ctx *Context) Reset(n int) {
+ if n <= 0 {
+ panic("linsolve: dimension not positive")
+ }
+ ctx.X.Reset()
+ ctx.X.ReuseAsVec(n)
+ ctx.Src.Reset()
+ ctx.Src.ReuseAsVec(n)
+ ctx.Dst.Reset()
+ ctx.Dst.ReuseAsVec(n)
+}
+
+// Operation specifies the type of operation.
+type Operation uint
+
+// Operations commanded by Method.Iterate.
+const (
+ NoOperation Operation = 0
+
+ // Compute A*x where x is stored in Context.Src. The
+ // result must be placed in Context.Dst.
+ MulVec Operation = 1 << (iota - 1)
+
+ // Perform a preconditioner solve
+ // M z = r
+ // where r is stored in Context.Src. The solution z
+ // must be placed in Context.Dst.
+ PreconSolve
+
+ // Trans indicates that MulVec or PreconSolve
+ // operation must be performed wih the transpose,
+ // that is, compute Aᵀ*x or solve Mᵀ*z = r. Method
+ // will command Trans only in bitwise OR combination
+ // with MulVec and PreconSolve.
+ Trans
+
+ // Compute b-A*x where x is stored in Context.X,
+ // and store the result in Context.Dst.
+ ComputeResidual
+
+ // Check convergence using (an estimate of) a
+ // residual norm in Context.ResidualNorm. Context.X
+ // does not need to be valid. The caller must set
+ // Context.Converged to indicate whether convergence
+ // has been determined, and then it must call
+ // Method.Iterate again.
+ CheckResidualNorm
+
+ // MajorIteration indicates that Method has finished
+ // what it considers to be one iteration. Method
+ // will make sure that Context.X is updated.
+ // If Context.Converged is true, the caller must
+ // terminate the iterative process, otherwise it
+ // should call Method.Iterate again.
+ MajorIteration
+)
+
+const (
+ // Machine epsilon.
+ eps = 1.0 / (1 << 53)
+
+ // Breakdown tolerance for BiCG and BiCGStab methods.
+ breakdownTol = eps * eps
+)
diff --git a/linsolve/linsolve_test.go b/linsolve/linsolve_test.go
new file mode 100644
index 0000000..f028bb8
--- /dev/null
+++ b/linsolve/linsolve_test.go
@@ -0,0 +1,762 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve
+
+import (
+ "fmt"
+ "math"
+
+ "golang.org/x/exp/rand"
+
+ "gonum.org/v1/gonum/lapack/testlapack"
+ "gonum.org/v1/gonum/linsolve/internal/triplet"
+ "gonum.org/v1/gonum/mat"
+)
+
+const defaultTol = 1e-13
+
+type testCase struct {
+ name string
+
+ mulVecTo func(*mat.VecDense, bool, mat.Vector) // Matrix-vector multiplication
+
+ b []float64 // Right-hand side vector
+ diag []float64 // Diagonal for the Jacobi preconditioner
+ tol float64 // Tolerance for the convergence criterion
+
+ want []float64 // Expected solution
+}
+
+func (tc *testCase) MulVecTo(dst *mat.VecDense, trans bool, x mat.Vector) {
+ tc.mulVecTo(dst, trans, x)
+}
+
+func (tc *testCase) PreconSolve(dst *mat.VecDense, trans bool, rhs mat.Vector) error {
+ if tc.diag == nil {
+ dst.CopyVec(rhs)
+ } else {
+ n := len(tc.diag)
+ diag := mat.NewVecDense(n, tc.diag)
+ dst.DivElemVec(rhs, diag)
+ }
+ return nil
+}
+
+func spdTestCases(rnd *rand.Rand) []testCase {
+ return []testCase{
+ newRandomSPD(1, rnd),
+ newRandomSPD(2, rnd),
+ newRandomSPD(3, rnd),
+ newRandomSPD(4, rnd),
+ newRandomSPD(5, rnd),
+ newRandomSPD(10, rnd),
+ newRandomSPD(20, rnd),
+ newRandomSPD(50, rnd),
+ newRandomDiagonal(2, rnd),
+ newRandomDiagonal(3, rnd),
+ newRandomDiagonal(4, rnd),
+ newRandomDiagonal(5, rnd),
+ newRandomDiagonal(10, rnd),
+ newRandomDiagonal(20, rnd),
+ newRandomDiagonal(50, rnd),
+ newGreenbaum41(24, 0.001, 1, 0.4, rnd),
+ newGreenbaum41(24, 0.001, 1, 0.6, rnd),
+ newGreenbaum41(24, 0.001, 1, 0.8, rnd),
+ newGreenbaum41(24, 0.001, 1, 1, rnd),
+ newPoisson1D(32, random(rnd)),
+ newPoisson2D(32, 32, one),
+ }
+}
+
+// newRandomSPD returns a test case with a random symmetric positive-definite
+// matrix of order n, and a random right-hand side.
+func newRandomSPD(n int, rnd *rand.Rand) testCase {
+ // Generate a random matrix C.
+ c := make([]float64, n*n)
+ for i := range c {
+ c[i] = rnd.NormFloat64()
+ }
+ C := mat.NewDense(n, n, c)
+ // Use C to generate a SPD matrix A.
+ var A mat.SymDense
+ A.SymOuterK(1, C)
+ for i := 0; i < n; i++ {
+ A.SetSym(i, i, A.At(i, i)+float64(n))
+ }
+ // Generate the right-hand side.
+ b := make([]float64, n)
+ for i := range b {
+ b[i] = 1 / math.Sqrt(float64(n))
+ }
+ // Compute the solution using the Cholesky factorization.
+ var chol mat.Cholesky
+ ok := chol.Factorize(&A)
+ if !ok {
+ panic("bad test matrix")
+ }
+ want := make([]float64, n)
+ chol.SolveVecTo(mat.NewVecDense(n, want), mat.NewVecDense(n, b))
+ // Matrix-vector multiplication.
+ mulVecTo := func(dst *mat.VecDense, _ bool, x mat.Vector) {
+ if dst.Len() != n || x.Len() != n {
+ panic("mismatched vector length")
+ }
+ dst.MulVec(&A, x)
+ }
+ // Store the diagonal for preconditioning.
+ diag := make([]float64, n)
+ for i := range diag {
+ diag[i] = A.At(i, i)
+ }
+ return testCase{
+ name: fmt.Sprintf("Random SPD n=%v", n),
+ mulVecTo: mulVecTo,
+ b: b,
+ tol: defaultTol,
+ diag: diag,
+ want: want,
+ }
+}
+
+// newRandomDiagonal returns a test case with a diagonal matrix with random positive elements,
+// a random right-hand side and a known solution.
+func newRandomDiagonal(n int, rnd *rand.Rand) testCase {
+ // Generate a diagonal matrix with random positive elements.
+ a := make([]float64, n)
+ diag := make([]float64, n)
+ for i := range a {
+ a[i] = 1 + 10*rnd.Float64()
+ diag[i] = a[i]
+ }
+ A := mat.NewDiagDense(n, a)
+ // Generate the right-hand side.
+ b := make([]float64, n)
+ for i := range b {
+ b[i] = 1 / math.Sqrt(float64(n))
+ }
+ // Compute the reference solution.
+ want := make([]float64, n)
+ for i := range want {
+ want[i] = b[i] / a[i]
+ }
+ // Matrix-vector multiplication.
+ mulVecTo := func(dst *mat.VecDense, _ bool, x mat.Vector) {
+ if dst.Len() != n || x.Len() != n {
+ panic("mismatched vector length")
+ }
+ dst.MulVec(A, x)
+ }
+ return testCase{
+ name: fmt.Sprintf("Random diagonal n=%v", n),
+ mulVecTo: mulVecTo,
+ b: b,
+ tol: defaultTol,
+ diag: diag,
+ want: want,
+ }
+}
+
+// newGreenbaum41 returns a test case with a symmetric positive definite matrix
+// A defined as
+// A = U * D * Uᵀ,
+// where U is a random orthogonal matrix and D is a diagonal matrix with entries
+// given by
+// d_i = d_1 + (i-1)/(n-1)*(d_n-d_1)*rho^{n-i}, i = 2,...,n-1.
+//
+// This matrix is described in Section 4.1 of
+// Greenbaum, A. (1997). Iterative Methods for Solving Linear Systems. SIAM.
+func newGreenbaum41(n int, d1, dn, rho float64, rnd *rand.Rand) testCase {
+ if n < 2 || dn < d1 {
+ panic("bad test")
+ }
+ // Generate the diagonal.
+ d := make([]float64, n)
+ d[0] = d1
+ d[n-1] = dn
+ for i := 1; i < n-1; i++ {
+ d[i] = d1 + float64(i)/float64(n-1)*(dn-d1)*math.Pow(rho, float64(n-i-1))
+ }
+ // Generate the SPD matrix A.
+ a := make([]float64, n*n)
+ testlapack.Dlagsy(n, 0, d, a, n, rnd, make([]float64, 2*n))
+ A := mat.NewSymDense(n, a)
+ // Generate the right-hand side.
+ b := make([]float64, n)
+ for i := range b {
+ b[i] = rnd.NormFloat64()
+ }
+ // Compute the solution using the Cholesky factorization.
+ var chol mat.Cholesky
+ ok := chol.Factorize(A)
+ if !ok {
+ panic("bad test matrix")
+ }
+ want := make([]float64, n)
+ chol.SolveVecTo(mat.NewVecDense(n, want), mat.NewVecDense(n, b))
+ // Matrix-vector multiplication.
+ mulVecTo := func(dst *mat.VecDense, _ bool, x mat.Vector) {
+ if dst.Len() != n || x.Len() != n {
+ panic("mismatched vector length")
+ }
+ dst.MulVec(A, x)
+ }
+ // Store the diagonal for preconditioning.
+ diag := make([]float64, n)
+ for i := range diag {
+ diag[i] = A.At(i, i)
+ }
+ return testCase{
+ name: fmt.Sprintf("Greenbaum 4.1 n=%v,d_1=%v,d_n=%v,rho=%v", n, d1, dn, rho),
+ mulVecTo: mulVecTo,
+ b: b,
+ tol: defaultTol,
+ diag: diag,
+ want: want,
+ }
+}
+
+func nonsym3x3() testCase {
+ return testCase{
+ name: "nonsym 3x3",
+ mulVecTo: func(dst *mat.VecDense, trans bool, x mat.Vector) {
+ A := mat.NewDense(3, 3, []float64{
+ 5, -1, 3,
+ -1, 2, -2,
+ 3, -2, 3,
+ })
+ if trans {
+ dst.MulVec(A.T(), x)
+ } else {
+ dst.MulVec(A, x)
+ }
+ },
+ b: []float64{7, -1, 4},
+ diag: []float64{5, 2, -3},
+ tol: defaultTol,
+ want: []float64{1, 1, 1},
+ }
+}
+
+func nonsymTridiag(n int) testCase {
+ A := triplet.NewMatrix(n, n)
+ for i := 0; i < n; i++ {
+ if i > 0 {
+ A.Append(i, i-1, -2)
+ }
+ A.Append(i, i, 4)
+ if i < n-1 {
+ A.Append(i, i+1, -1)
+ }
+ }
+ b := make([]float64, n)
+ for i := range b {
+ switch i {
+ case 0:
+ b[i] = 3
+ default:
+ b[i] = 1
+ case n - 1:
+ b[i] = 2
+ }
+ }
+ want := make([]float64, n)
+ for i := range want {
+ want[i] = 1
+ }
+ return testCase{
+ name: fmt.Sprintf("Nonsym tridiag n=%v", n),
+ mulVecTo: A.MulVecTo,
+ b: b,
+ tol: defaultTol,
+ want: want,
+ }
+}
+
+// newPoisson1D returns a test case that arises from a finite-difference discretization
+// of the Poisson equation
+// - ∂_x ∂_x u = f
+// on the interval [0,1].
+func newPoisson1D(nx int, f func(float64, float64) float64) testCase {
+ tc := newPDE(nx, 1, negOne, nil, zero, nil, zero, f)
+ tc.name = fmt.Sprintf("Poisson 1D nx=%v", nx)
+ return tc
+}
+
+// newPoisson2D returns a test case that arises from a finite-difference discretization
+// of the Poisson equation
+// - Δu = f
+// on the unit square [0,1]×[0,1].
+func newPoisson2D(nx, ny int, f func(float64, float64) float64) testCase {
+ tc := newPDE(nx, ny, negOne, negOne, zero, zero, zero, f)
+ tc.name = fmt.Sprintf("Poisson 2D nx=%v,ny=%v", nx, ny)
+ tc.tol = 1e-12
+ return tc
+}
+
+// newGreenbaum54 returns a test case with a general unsymmetric matrix
+// A defined as
+// A = V*D*V^{-1},
+// where V is a random matrix and D is a block-diagonal matrix with n1 complex
+// and n2 real eigenvalues.
+//
+// This matrix is described in Section 5.4 of
+// Greenbaum, A. (1997). Iterative Methods for Solving Linear Systems. SIAM.
+func newGreenbaum54(n1, n2 int, rnd *rand.Rand) testCase {
+ n := 2*n1 + n2
+ // Generate the block-diagonal matrix D.
+ d := make([]float64, n*3)
+ // Generate n1 blocks with complex eigenvalues.
+ for i := 0; i < 2*n1; i += 2 {
+ // The 2x2 block will have eigenvalues a±b*i.
+ // Real part of the eigenvalue is in [1,2).
+ a := rnd.Float64() + 1
+ // Imaginary part is in [-1,1).
+ b := 2*rnd.Float64() - 1
+ d[i*3+1] = a
+ d[i*3+2] = b
+ d[(i+1)*3] = -b
+ d[(i+1)*3+1] = a
+ }
+ // Generate n2 real eigenvalues.
+ for i := 2 * n1; i < n; i++ {
+ r := 9*rnd.Float64() + 1
+ if rnd.Intn(2) == 0 {
+ r *= -1
+ }
+ d[i*3+1] = r
+ }
+ D := mat.NewBandDense(n, n, 1, 1, d)
+ // Generate a random matrix V.
+ v := make([]float64, n*n)
+ for i := range v {
+ v[i] = rnd.NormFloat64()
+ }
+ V := mat.NewDense(n, n, v)
+ var luV mat.LU
+ luV.Factorize(V)
+ // Generate the right-hand side.
+ b := make([]float64, n)
+ for i := range b {
+ b[i] = rnd.NormFloat64()
+ }
+ // Compute V*D and (V*D)^{-1} for computing the reference solution and for the matrix-vector operation.
+ var VD mat.Dense
+ VD.Mul(V, D)
+ var luVD mat.LU
+ luVD.Factorize(&VD)
+ // Compute the solution of V*D*V^{-1}*x = b.
+ // First, compute the solution of V*D*y = b.
+ want := make([]float64, n)
+ wantVec := mat.NewVecDense(n, want)
+ err := luVD.SolveVecTo(wantVec, false, mat.NewVecDense(n, b))
+ if err != nil {
+ panic("luVD.SolveVecTo(notrans) failed")
+ }
+ // Second, compute the solution of V^{-1}*x = y, which amounts to just
+ // computing x = V*y.
+ wantVec.MulVec(V, wantVec)
+ // Matrix-vector multiplication.
+ mulVecTo := func(dst *mat.VecDense, trans bool, x mat.Vector) {
+ if dst.Len() != n || x.Len() != n {
+ panic("mismatched vector length")
+ }
+ if trans {
+ // Multiply (V*D*V^{-1})ᵀ * x which can be
+ // rewritten as V^{-1}ᵀ * (V*D)ᵀ * x.
+ dst.MulVec(VD.T(), x)
+ err := luV.SolveVecTo(dst, true, dst)
+ if err != nil {
+ panic("luV.SolveVecTo(trans) failed")
+ }
+ } else {
+ // Multiply V*D*V^{-1} * x.
+ err := luV.SolveVecTo(dst, false, x)
+ if err != nil {
+ panic("luV.SolveVecTo(notrans) failed")
+ }
+ dst.MulVec(&VD, dst)
+ }
+ }
+ return testCase{
+ name: fmt.Sprintf("Greenbaum 5.4 n=%v,n1=%v,n2=%v", n, n1, n2),
+ mulVecTo: mulVecTo,
+ b: b,
+ tol: defaultTol,
+ want: want,
+ }
+}
+
+// newGreenbaum73 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// - Δu + 40*(x*∂_x u + y*∂_y u) - 100*u = f
+//
+// This test problem is described in Section 7.3 of
+// Greenbaum, A. (1997). Iterative Methods for Solving Linear Systems. SIAM.
+func newGreenbaum73(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny,
+ negOne, negOne, // - ∂_x ∂_x u - ∂_y ∂_y u
+ func(x, _ float64) float64 { return 40 * x }, // 40 * x * ∂_x u
+ func(_, y float64) float64 { return 40 * y }, // 40 * y * ∂_y u
+ constant(-100), // -100 * u
+ random(rnd))
+ tc.name = fmt.Sprintf("Greenbaum 7.3 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDENonsymmetric returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu + henk*∂_x u + (∂_x henk/2)*u = f
+// where henk(x,y) := 20*exp(3.5*(x^2 + y^2))
+//
+// This test case is adapted from
+// http://www.netlib.org/templates/dftemplates.tgz
+func newPDENonsymmetric(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny, one, one, henk, zero, dhenkdx, random(rnd))
+ tc.name = fmt.Sprintf("PDE Nonsymmetric nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+func henk(x, y float64) float64 {
+ return 20 * math.Exp(3.5*(x*x+y*y))
+}
+
+func dhenkdx(x, y float64) float64 {
+ return 70 * x * math.Exp(3.5*(x*x+y*y))
+}
+
+// newPDEYang47 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu + 1000*∂_x u = f
+// The large coefficient of ∂_xu causes a loss of diagonal dominance in the system matrix.
+//
+// This test case corresponds to Eq. 4.7 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang47(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny, one, one, constant(1000), zero, zero, random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.7 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang48 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu + 1000*∂_y u = f
+// The large coefficient of ∂_yu causes a loss of diagonal dominance in the system matrix.
+//
+// This test case corresponds to Eq. 4.8 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang48(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny, one, one, zero, constant(1000), zero, random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.8 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang49 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu + 1000*exp(x*y)*(∂_x u - ∂_y u) = f
+//
+// This test case corresponds to Eq. 4.9 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang49(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(16, 16,
+ one, one,
+ func(x, y float64) float64 { return 1000 * math.Exp(x*y) },
+ func(x, y float64) float64 { return -1000 * math.Exp(x*y) },
+ zero, random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.9 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang410 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// ∂_x (exp(x*y) * ∂_x u) + ∂_y (exp(-x*y) * ∂_y u) - 500*(x + y)*∂_x u - (250*(1 + y) + 1/(1 + x + y))*u = f
+//
+// This test case corresponds to Eq. 4.10 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang410(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny,
+ func(x, y float64) float64 { return math.Exp(x * y) }, // ∂_x (exp(x*y) * ∂_x u)
+ func(x, y float64) float64 { return math.Exp(-x * y) }, // ∂_y (exp(-x*y) * ∂_y u)
+ func(x, y float64) float64 { return -500 * (x + y) }, // -500*(x + y)*∂_x u
+ zero,
+ func(x, y float64) float64 { return -250*(1+y) - 1/(1+x+y) }, // - (250*(1 + y) + 1/(1+x+y)) * u
+ random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.10 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang412 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu - 100000*x^2*(∂_x u + ∂_y u) = f
+//
+// This test case corresponds to Eq. 4.12 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang412(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny,
+ one, one,
+ func(x, _ float64) float64 { return -100000 * x * x },
+ func(x, _ float64) float64 { return -100000 * x * x },
+ zero, random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.12 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang413 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu - 1000*(1 + x^2)*∂_x u + 100*∂_y u = f
+//
+// This test case corresponds to Eq. 4.13 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang413(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny,
+ one, one,
+ func(x, _ float64) float64 { return -1000 * (1 + x*x) },
+ constant(100),
+ zero, random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.13 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang414 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu - 1000*x^2*∂_x u + 1000*u = f
+//
+// This test case corresponds to Eq. 4.14 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang414(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny,
+ one, one,
+ func(x, _ float64) float64 { return -1000 * x * x },
+ zero,
+ constant(1000), random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.14 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDEYang415 returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// Δu - 1000*(1 - 2*x)*∂_x u - 1000*(1 - 2*y)*∂_y u = f
+//
+// This test case corresponds to Eq. 4.15 in
+// Ulrike Meier Yang (1994), Preconditioned Conjugate Gradient-Like Methods for Nonsymmetric Linear Systems.
+func newPDEYang415(nx, ny int, rnd *rand.Rand) testCase {
+ tc := newPDE(nx, ny,
+ one, one,
+ func(x, _ float64) float64 { return -1000 * (1 - 2*x) },
+ func(_, y float64) float64 { return -1000 * (1 - 2*y) },
+ zero, random(rnd))
+ tc.name = fmt.Sprintf("PDE Yang, Eq. 4.15 nx=%v,ny=%v", nx, ny)
+ return tc
+}
+
+// newPDE returns a test case that arises from a finite-difference discretization
+// of the partial differential equation
+// ∂_x (a ∂_x u) + ∂_y (b ∂_y u) + c ∂_x u + d ∂_y u + e u = f
+// on the unit square [0,1]×[0,1] with zero Dirichlet boundary conditions.
+//
+// nx and ny must be positive. If ny is 1, a 1D variant of the equation on [0,1]×{0} will be used,
+// and b and d will not be referenced.
+func newPDE(nx, ny int, a, b, c, d, e, f func(float64, float64) float64) testCase {
+ if nx <= 0 || ny <= 0 {
+ panic("invalid mesh size")
+ }
+
+ var (
+ A *triplet.Matrix
+ rhs []float64
+ diag []float64
+ )
+ if ny == 1 {
+ A, rhs, diag = newPDESystem1D(nx, a, c, e, f)
+ } else {
+ A, rhs, diag = newPDESystem2D(nx, ny, a, b, c, d, e, f)
+ }
+ // Use a dense solver to obtain a reference solution.
+ Ad := A.DenseCopy()
+ var lu mat.LU
+ lu.Factorize(Ad)
+ n := len(rhs)
+ want := make([]float64, n)
+ err := lu.SolveVecTo(mat.NewVecDense(n, want), false, mat.NewVecDense(n, rhs))
+ if err != nil {
+ panic("lu.SolveVec failed")
+ }
+ return testCase{
+ mulVecTo: A.MulVecTo,
+ b: rhs,
+ tol: defaultTol,
+ diag: diag,
+ want: want,
+ }
+}
+
+// newPDESystem1D assembles and returns the matrix A, the right-hand side vector,
+// and the diagonal of A for a 1-dimensional PDE problem.
+func newPDESystem1D(nx int, a, b, c, f func(float64, float64) float64) (A *triplet.Matrix, rhs, diag []float64) {
+ h := 1 / float64(nx+1)
+ A = triplet.NewMatrix(nx, nx)
+ rhs = make([]float64, nx)
+ diag = make([]float64, nx)
+ var i int
+ for ix := 0; ix < nx; ix++ {
+ s := newStencil1D(ix, h, a, b, c)
+ // Add stencil elements to the system matrix, skipping boundary nodes.
+ if ix > 0 {
+ A.Append(i, i-1, s.left)
+ }
+ A.Append(i, i, s.center)
+ diag[i] = s.center
+ if ix < nx-1 {
+ A.Append(i, i+1, s.right)
+ }
+ // Assemble the right-hand side.
+ x := float64(ix+1) * h
+ rhs[i] = f(x, 0) * h * h
+ i++
+ }
+ return A, rhs, diag
+}
+
+type stencil1D struct {
+ left, right float64
+ center float64
+}
+
+// newStencil1D returns a finite difference stencil that approximates the differential operator
+// ∂_x (a ∂_x u) + b ∂_x u + c u
+// at point [(i+1)*h] using the formula
+// a(i+1/2,0)*(u(i+1) - u(i)) + a(i-1/2,0)*(u(i-1) - u(i))
+// + (h/2)*b(i,0)*(u(i+1) - u(i-1))
+// + h^2*c(i,0)*u(i)
+func newStencil1D(i int, h float64, a, b, c func(float64, float64) float64) (s stencil1D) {
+ x := float64(i+1) * h
+
+ coeff := a(x+0.5*h, 0)
+ s.center -= coeff
+ s.right = coeff
+ coeff = a(x-0.5*h, 0)
+ s.center -= coeff
+ s.left = coeff
+
+ coeff = b(x, 0)
+ s.right += 0.5 * h * coeff
+ s.left -= 0.5 * h * coeff
+
+ s.center += h * h * c(x, 0)
+
+ return s
+}
+
+// newPDESystem2D assembles and returns the matrix A, the right-hand side vector,
+// and the diagonal of A for a 2-dimensional PDE problem.
+func newPDESystem2D(nx, ny int, a, b, c, d, e, f func(float64, float64) float64) (A *triplet.Matrix, rhs, diag []float64) {
+ // Finite difference stencil:
+ // * (ix,iy+1)
+ // |
+ // | (ix,iy)
+ // * ---- * ---- *
+ // (ix-1,iy) | (ix+1,iy)
+ // |
+ // * (ix,iy-1)
+ // Node (ix,iy) is mapped to the index ix+iy*k.
+
+ h := 1 / float64(nx+1)
+ n := nx * ny
+ A = triplet.NewMatrix(n, n)
+ rhs = make([]float64, n)
+ diag = make([]float64, n)
+ var i int
+ for iy := 0; iy < ny; iy++ {
+ y := float64(iy+1) * h
+ for ix := 0; ix < nx; ix++ {
+ s := newStencil2D(ix, iy, h, a, b, c, d, e)
+ // Add the coefficients from the stencil to the system matrix, skipping boundary nodes.
+ if iy > 0 {
+ A.Append(i, i-nx, s.down)
+ }
+ if ix > 0 {
+ A.Append(i, i-1, s.left)
+ }
+ A.Append(i, i, s.center)
+ diag[i] = s.center
+ if ix < nx-1 {
+ A.Append(i, i+1, s.right)
+ }
+ if iy < ny-1 {
+ A.Append(i, i+nx, s.up)
+ }
+ // Assemble the right-hand side.
+ x := float64(ix+1) * h
+ rhs[i] = f(x, y) * h * h
+ i++
+ }
+ }
+ return A, rhs, diag
+}
+
+type stencil2D struct {
+ left, right float64
+ up, down float64
+ center float64
+}
+
+// newStencil2D returns a finite difference stencil that approximates the differential operator
+// ∂_x (a ∂_x u) + ∂_y (b ∂_y u) + c ∂_x u + d ∂_y u + e u
+// at point [(i+1)*h,(j+1)*h] using the formula
+// a(i+1/2,j)*(u(i+1,j) - u(i,j)) + a(i-1/2,j)*(u(i-1,j) - u(i,j))
+// + b(i,j+1/2)*(u(i,j+1) - u(i,j)) + b(i,j-1/2)*(u(i,j-1) - u(i,j))
+// + (h/2)*c(i,j)*(u(i+1,j) - u(i-1,j)) + (h/2)*d(i,j)*(u(i,j+1) - u(i,j-1))
+// + h^2*e(i,j)*u(i,j)
+func newStencil2D(i, j int, h float64, a, b, c, d, e func(float64, float64) float64) (s stencil2D) {
+ x := float64(i+1) * h
+ y := float64(j+1) * h
+
+ coeff := a(x+0.5*h, y)
+ s.center -= coeff
+ s.right = coeff
+ coeff = a(x-0.5*h, y)
+ s.center -= coeff
+ s.left = coeff
+ coeff = b(x+0.5*h, y)
+ s.center -= coeff
+ s.up = coeff
+ coeff = b(x-0.5*h, y)
+ s.center -= coeff
+ s.down = coeff
+
+ coeff = c(x, y)
+ s.right += 0.5 * h * coeff
+ s.left -= 0.5 * h * coeff
+ coeff = d(x, y)
+ s.up += 0.5 * h * coeff
+ s.down -= 0.5 * h * coeff
+
+ s.center += h * h * c(x, y)
+
+ return s
+}
+
+func zero(_, _ float64) float64 {
+ return 0
+}
+
+func one(_, _ float64) float64 {
+ return 1
+}
+
+func negOne(_, _ float64) float64 {
+ return -1
+}
+
+func constant(c float64) func(_, _ float64) float64 {
+ return func(_, _ float64) float64 {
+ return c
+ }
+}
+
+func random(rnd *rand.Rand) func(_, _ float64) float64 {
+ return func(_, _ float64) float64 {
+ return rnd.NormFloat64()
+ }
+}
diff --git a/linsolve/pde_example_test.go b/linsolve/pde_example_test.go
new file mode 100644
index 0000000..292c1ac
--- /dev/null
+++ b/linsolve/pde_example_test.go
@@ -0,0 +1,256 @@
+// Copyright ©2019 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package linsolve_test
+
+import (
+ "fmt"
+ "log"
+
+ "golang.org/x/exp/rand"
+
+ "gonum.org/v1/gonum/linsolve"
+ "gonum.org/v1/gonum/mat"
+ "gonum.org/v1/plot"
+ "gonum.org/v1/plot/plotter"
+ "gonum.org/v1/plot/plotutil"
+ "gonum.org/v1/plot/vg"
+)
+
+// AllenCahnFD implements a semi-implicit finite difference scheme for the
+// solution of the one-dimensional Allen-Cahn equation
+// u_t = u_xx - 1/ξ²·f'(u) in (0,L)×(0,T)
+// u_x = 0 on (0,T)
+// u(0) = u0 on (0,L)
+// where f is a double-well-shaped function with two minima at ±1
+// f(s) = 1/4·(s²-1)²
+//
+// The equation arises in materials science in the description of phase
+// transitions, e.g. solidification in crystal growth, but also in other areas
+// like image processing due to its connection to mean-curvature motion.
+// Starting the evolution from an initial distribution u0, the solution u
+// develops a thin steep layer, an interface between regions of the domain
+// (0,L) where u is constant and close to one of the minima of f.
+//
+// AllenCahnFD approximates derivatives by finite differences and the solution
+// is advanced in time by a semi-implicit Euler scheme where the nonlinear term
+// is taken from the previous time step. Therefore, at each time step a linear
+// system must be solved.
+type AllenCahnFD struct {
+ // Xi is the ξ parameter that determines the interface width.
+ Xi float64
+
+ // InitCond is the initial condition u0.
+ InitCond func(x float64) float64
+
+ h float64 // Spatial step size
+ tau float64 // Time step size
+
+ a *mat.SymBandDense
+ b *mat.VecDense
+ u *mat.VecDense
+
+ ls linsolve.Method
+ lssettings linsolve.Settings
+}
+
+// FPrime returns the value of the derivative of the double-well potential f at s.
+// f'(s) = s·(s²-1)
+func FPrime(s float64) float64 {
+ return s * (s*s - 1)
+}
+
+// Setup initializes the receiver for solving the Allen-Cahn equation on a
+// uniform grid with n+1 nodes on the spatial interval (0,L) and with the time
+// step size tau.
+func (ac *AllenCahnFD) Setup(n int, L float64, tau float64) {
+ ac.h = L / float64(n)
+ ac.tau = tau
+
+ // We solve this PDE numerically by replacing the derivatives with finite
+ // differences. For the spatial derivative, we use a central difference
+ // scheme, and for the time derivative derivative we use semi-implicit Euler
+ // integration where the non-linear term with f' is taken from the previous
+ // time step.
+ //
+ // After replacing the derivatives we get
+ // 1/tau*(u^{k+1}_i - u^k_i) = 1/h²*(u^{k+1}_{i-1} - 2*u^{k+1}_i + u^{k+1}_{i+1}) - 1/ξ²*f'(u^k_i)
+ // where tau is the time step size, h is the spatial step size, and u^k_i
+ // denotes the numerical solution that approximates u at time level k and
+ // grid node i, that is,
+ // u^k_i ≅ u(k*tau,i*h)
+ // Multiplying the equation by tau and collecting the terms from the same
+ // time level on each side gives
+ // -tau/h²*u^{k+1}_{i-1} + (1+2*tau/h²)*u^{k+1}_i - tau/h²*u^{k+1}_{i+1}) = u^k_i - tau/ξ²*f'(u^k_i)
+ // If we denote C:=tau/h² we can simplify this to
+ // -C*u^{k+1}_{i-1} + (1+2*C)*u^{k+1}_i - C*u^{k+1}_{i+1} = u^k_i - tau/ξ²*f'(u^k_i) (1)
+ // which must hold for all i=0,...,n.
+ //
+ // When i=0 or i=n, the expression (1) refers to values at nodes -1 and n+1
+ // which lie outside of the domain. We can eliminate them by using the fact
+ // that the first derivative is zero at the boundary. We approximate it by
+ // central difference:
+ // -1/h*(u^{k+1}_{-1} - u^{k+1}_1) = 0
+ // 1/h*(u^{k+1}_{n+1} - u^{k+1}_{n-1}) = 0
+ // which after simplifying gives
+ // u^{k+1}_{-1} = u^{k+1}_1
+ // u^{k+1}_{n+1} = u^{k+1}_{n-1}
+ // By substituting these two expressions into (1) at i=0 and i=n values at
+ // outside nodes do not appear in the expressions. If we further divide them
+ // by 2 (to obtain a symmetric matrix), we finally get
+ // (1/2+C)*u^{k+1}_0 - C*u^{k+1}_1 = 1/2*u^k_0 - 1/2*tau/ξ²*f'(u^k_0) (2)
+ // -C*u^{k+1}_{n-1} + (1/2+C)*u^{k+1}_n = 1/2*u^k_n - 1/2*tau/ξ²*f'(u^k_n) (3)
+ // Note that simply means that we treat values at the boundary nodes the
+ // same as the inner nodes.
+ //
+ // For a given k the equations (1),(2),(3) form a linear system for the
+ // unknown vector [u^{k+1}_i], i=0,...,n that must be solved at each step in
+ // order to advance the solution in time. The matrix of this system is
+ // tridiagonal and symmetric positive-definite:
+ // ⎛1/2+C -C 0 . . . . 0⎞
+ // ⎜ -C 1+2C -C .⎟
+ // ⎜ 0 -C 1+2C -C .⎟
+ // ⎜ . -C . . .⎟
+ // ⎜ . . . . .⎟
+ // ⎜ . . . -C 0⎟
+ // ⎜ . -C 1+2C -C⎟
+ // ⎝ 0 . . . . 0 -C 1/2+C⎠
+ // The right-hand side is:
+ // ⎛1/2*u^k_0 - 1/2*tau/ξ²*f'(u^k_0) ⎞
+ // ⎜ u^k_1 - tau/ξ²*f'(u^k_1) ⎟
+ // ⎜ u^k_2 - tau/ξ²*f'(u^k_2) ⎟
+ // ⎜ . ⎟
+ // ⎜ . ⎟
+ // ⎜ . ⎟
+ // ⎜ u^k_{n-1} - tau/ξ²*f'(u^k_{n-1})⎟
+ // ⎝1/2*u^k_n - 1/2*tau/ξ²*f'(u^k_n) ⎠
+
+ // Assemble the system matrix A based on the discretization scheme above.
+ // Since the matrix is symmetric, we only need to set elements in the upper
+ // triangle.
+ A := mat.NewSymBandDense(n+1, 1, nil)
+ c := ac.tau / ac.h / ac.h
+ // Boundary condition at the left node
+ A.SetSymBand(0, 0, 0.5+c)
+ A.SetSymBand(0, 1, -c)
+ // Interior nodes
+ for i := 1; i < n; i++ {
+ A.SetSymBand(i, i, 1+2*c)
+ A.SetSymBand(i, i+1, -c)
+ }
+ // Boundary condition at the right node
+ A.SetSymBand(n, n, 0.5+c)
+ ac.a = A
+
+ // Allocate the right-hand side b.
+ ac.b = mat.NewVecDense(n+1, nil)
+
+ // Allocate and set up the initial condition.
+ ac.u = mat.NewVecDense(n+1, nil)
+ for i := 0; i < ac.u.Len(); i++ {
+ ac.u.SetVec(i, ac.InitCond(float64(i)*ac.h))
+ }
+
+ // Allocate the linear solver and the settings.
+ ac.ls = &linsolve.CG{}
+ ac.lssettings = linsolve.Settings{
+ // Solution from the previous time step will be a good initial estimate.
+ InitX: ac.u,
+ // Store the solution into the existing vector.
+ Dst: ac.u,
+ // Provide context to reduce memory allocation and GC pressure.
+ Work: linsolve.NewContext(n + 1),
+ }
+}
+
+// Step advances the solution one step in time.
+func (ac *AllenCahnFD) Step() error {
+ // Assemble the right-hand side vector b.
+ tauXi2 := ac.tau / ac.Xi / ac.Xi
+ n := ac.u.Len()
+ for i := 0; i < ac.u.Len(); i++ {
+ ui := ac.u.AtVec(i)
+ bi := ui - tauXi2*FPrime(ui)
+ if i == 0 || i == n-1 {
+ bi *= 0.5
+ }
+ ac.b.SetVec(i, bi)
+ }
+ // Solve the system.
+ _, err := linsolve.Iterative(ac, ac.b, ac.ls, &ac.lssettings)
+ return err
+}
+
+// MulVecTo implements the MulVecToer interface.
+func (ac *AllenCahnFD) MulVecTo(dst *mat.VecDense, _ bool, x mat.Vector) {
+ dst.MulVec(ac.a, x)
+}
+
+func (ac *AllenCahnFD) Solution() *mat.VecDense {
+ return ac.u
+}
+
+func output(u mat.Vector, L float64, step int) error {
+ p, err := plot.New()
+ if err != nil {
+ return err
+ }
+
+ p.Title.Text = fmt.Sprintf("Step %d", step)
+ p.X.Label.Text = "x"
+ p.X.Min = 0
+ p.X.Max = L
+ p.Y.Min = -1.1
+ p.Y.Max = 1.1
+
+ n := u.Len()
+ h := L / float64(n-1)
+ pts := make(plotter.XYs, n)
+ for i := 0; i < u.Len(); i++ {
+ pts[i].X = float64(i) * h
+ pts[i].Y = u.AtVec(i)
+ }
+ err = plotutil.AddLines(p, "u", pts)
+ if err != nil {
+ return err
+ }
+ err = p.Save(20*vg.Centimeter, 10*vg.Centimeter, fmt.Sprintf("u%04d.png", step))
+ if err != nil {
+ return err
+ }
+ return nil
+}
+
+func ExampleIterative_evolutionPDE() {
+ const (
+ L = 10.0
+ nx = 1000
+ nt = 200
+ tau = 0.1 * L / nx
+ xi = 6.0 * L / nx
+ )
+ rnd := rand.New(rand.NewSource(1))
+ ac := AllenCahnFD{
+ Xi: xi,
+ InitCond: func(x float64) float64 {
+ // Initial condition is a perturbation of the (unstable) zero state
+ // (the peak in the double-well function f).
+ return 0.01 * rnd.NormFloat64()
+ },
+ }
+ ac.Setup(nx, L, tau)
+ for i := 1; i <= nt; i++ {
+ err := ac.Step()
+ if err != nil {
+ log.Fatal(err)
+ }
+ // Generate plots of u as PNG images that can be converted to video
+ // using for example
+ // ffmpeg -i u%04d.png output.webm
+ err = output(ac.Solution(), L, i)
+ if err != nil {
+ log.Fatal(err)
+ }
+ }
+}
