Crate iterative_solvers

IterativeSolvers

English | 简体中文

This library provides Rust implementations of iterative algorithms for solving linear system, drawing heavy inspiration from the Julia package IterativeSolvers.jl.

Iterative Algorithms

• Conjugate Gradient (CG)

Supported Linear Algebra Libraries

• nalgebra (default)
• faer

You can choose your preferred backend using Cargo features:

[dependencies]
iterative-solvers = { version = "0.1", default-features = false, features = ["faer"] }

Usage

Consider the following differential equation:

-\frac{d^2 u}{dx^2} = \pi^2 \sin(\pi x) \quad \text{for} \quad x \in (0, 1)

with the boundary conditions:

u(0) = 0, \quad u(1) = 0.

The solution to this problem is:

u(x) = \sin(\pi x).

Now we use the central difference method to discretize the differential equation:

-\frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} = \pi^2 \sin(\pi x_i) \quad \text{for} \quad i = 1, \ldots, N-1,

where $h = 1/N$ is the step size, and $u_i$ is the approximation of $u(x_i)$. The boundary conditions are:

u_0 = 0, \quad u_N = 0.

The coefficient matrix of this linear system is a symmetric tridiagonal matrix, whose diagonal elements are $2/h^2$ and the sub-diagonal elements are $-1/h^2$.

We can solve this linear system using the Conjugate Gradient (CG) method.

use iterative_solvers::{cg, CG, utils::sparse::symmetric_tridiagonal_csc};
use nalgebra::DVector;
use std::f64::consts::PI;

fn main() {
    let n = 1024;
    let h = 1.0 / 1024.0;
    let a = vec![2.0 / (h * h); n - 1];
    let b = vec![-1.0 / (h * h); n - 2];
    // Store the symmetric tridiagonal matrix in CSC format
    let mat = symmetric_tridiagonal_csc(&a, &b).unwrap();
    // Generate the right-hand side vector
    let rhs: Vec<_> = (1..n)
        .map(|i| PI * PI * (i as f64 * h * PI).sin())
        .collect();
    // Generate the exact solution
    let solution: Vec<_> = (1..n).map(|i| (i as f64 * h * PI).sin()).collect();
    let solution = DVector::from_vec(solution);
    let rhs = DVector::from_vec(rhs);
    // Solve the linear system using the CG method
    let result = cg(&mat, &rhs, 1e-10, 1e-8).unwrap();
    // Calculate the error
    let e = (solution - result.solution()).norm();
    println!("error: {}", e);
}

If you want to know the residual, the approximate solution and the conjugate direction at each iteration, the iterator will help you.

let abstol = 1e-10;
let reltol = 1e-8;
let mut solver = CG::new(&mat, &rhs, abstol, reltol).unwrap();
while let Some(residual) = solver.next() {
    println!("residual: {residual}");
    println!("solution: {:#?}", solver.solution());
    println!("conjugate direction: {:#?}", solver.conjugate_direction());
}
let e = (solution - solver.solution()).norm();
println!("error: {}", e);

License

Licensed under either of:

• Apache License, Version 2.0, (LICENSE-APACHE or https://www.apache.org/licenses/LICENSE-2.0)
• MIT license (LICENSE-MIT or https://opensource.org/licenses/MIT)

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

Re-exports

pub use error::*;

Modules

cg

Conjugate Gradient (CG) method.

error

Error and Result type for iterative solvers.

ops

utils

Utility functions for creating sparse and dense matrices.

Structs

CG

Conjugate Gradient (CG) method for solving linear systems Ax = b.

Traits

MatrixOp

Functions

cg

Solve a linear system Ax = b using the Conjugate Gradient method.

cg_with_initial_guess

Solve a linear system Ax = b using the Conjugate Gradient method.