scirs2_sparse/linalg/
minres.rs

1use crate::error::{SparseError, SparseResult};
2use crate::linalg::interface::LinearOperator;
3use scirs2_core::numeric::{Float, NumAssign, SparseElement};
4use std::iter::Sum;
5
6/// Result of MINRES solver
7#[derive(Debug, Clone)]
8pub struct MINRESResult<F> {
9    pub x: Vec<F>,
10    pub iterations: usize,
11    pub residual_norm: F,
12    pub converged: bool,
13    pub message: String,
14}
15
16/// Options for MINRES solver
17pub struct MINRESOptions<F> {
18    pub max_iter: usize,
19    pub rtol: F,
20    pub atol: F,
21    pub x0: Option<Vec<F>>,
22    pub preconditioner: Option<Box<dyn LinearOperator<F>>>,
23}
24
25impl<F: Float> Default for MINRESOptions<F> {
26    fn default() -> Self {
27        Self {
28            max_iter: 1000,
29            rtol: F::from(1e-8).unwrap(),
30            atol: F::from(1e-12).unwrap(),
31            x0: None,
32            preconditioner: None,
33        }
34    }
35}
36
37/// MINRES solver for symmetric indefinite systems
38///
39/// Solves Ax = b where A is symmetric (possibly indefinite)
40/// using the Minimal Residual method.
41///
42/// Based on:
43/// C. C. Paige and M. A. Saunders (1975),
44/// "Solution of sparse indefinite systems of linear equations",
45/// SIAM J. Numer. Anal. 12(4), pp. 617-629.
46///
47/// # Arguments
48/// * `a` - The linear operator representing the matrix A
49/// * `b` - The right-hand side vector
50/// * `options` - Solver options including maximum iterations and tolerances
51///
52/// # Returns
53/// A `MINRESResult` containing the solution vector, iteration count, residual norm,
54/// and convergence status.
55///
56/// # Example
57/// ```
58/// use scirs2_sparse::linalg::{minres, MINRESOptions, IdentityOperator};
59///
60/// let identity = IdentityOperator::<f64>::new(3);
61/// let b = vec![1.0, 2.0, 3.0];
62/// let options = MINRESOptions::default();
63/// let result = minres(&identity, &b, options).unwrap();
64/// assert!(result.converged);
65/// ```
66#[allow(dead_code)]
67pub fn minres<F>(
68    a: &dyn LinearOperator<F>,
69    b: &[F],
70    options: MINRESOptions<F>,
71) -> SparseResult<MINRESResult<F>>
72where
73    F: Float + NumAssign + Sum + SparseElement + 'static,
74{
75    let (m, n) = a.shape();
76    if m != n {
77        return Err(SparseError::DimensionMismatch {
78            expected: m,
79            found: n,
80        });
81    }
82
83    if b.len() != n {
84        return Err(SparseError::DimensionMismatch {
85            expected: n,
86            found: b.len(),
87        });
88    }
89
90    // Initial guess
91    let mut x = options
92        .x0
93        .clone()
94        .unwrap_or_else(|| vec![F::sparse_zero(); n]);
95    let bnorm = norm2(b);
96
97    if bnorm < options.atol {
98        return Ok(MINRESResult {
99            x,
100            iterations: 0,
101            residual_norm: F::sparse_zero(),
102            converged: true,
103            message: "System has zero right-hand side".to_string(),
104        });
105    }
106
107    // Compute initial residual
108    let r1 = compute_residual(a, b, &x)?;
109
110    // Apply preconditioner if available
111    let y = if let Some(ref prec) = options.preconditioner {
112        prec.matvec(&r1)?
113    } else {
114        r1.clone()
115    };
116
117    // beta1 = sqrt(r1' * M^{-1} * r1)
118    let beta1_sq = inner(&r1, &y);
119
120    if beta1_sq < F::sparse_zero() {
121        return Err(SparseError::ComputationError(
122            "Indefinite preconditioner".to_string(),
123        ));
124    }
125
126    if beta1_sq == F::sparse_zero() {
127        return Ok(MINRESResult {
128            x,
129            iterations: 0,
130            residual_norm: F::sparse_zero(),
131            converged: true,
132            message: "Initial residual is zero".to_string(),
133        });
134    }
135
136    let beta1 = beta1_sq.sqrt();
137
138    // Initialize quantities for the iteration
139    let mut oldb = F::sparse_zero();
140    let mut beta = beta1;
141    let mut dbar = F::sparse_zero();
142    let mut epsln = F::sparse_zero();
143    let mut phibar = beta1;
144    let mut cs = F::from(-1.0).unwrap();
145    let mut sn = F::sparse_zero();
146    let mut w = vec![F::sparse_zero(); n];
147    let mut w2 = vec![F::sparse_zero(); n];
148    let mut r2 = r1.clone();
149    let mut v = vec![F::sparse_zero(); n];
150    let mut r1 = vec![F::sparse_zero(); n];
151    let mut y_vec = y;
152
153    let mut gmax = F::sparse_zero();
154    let mut gmin = F::max_value();
155    let mut tnorm2 = F::sparse_zero();
156    let mut qrnorm = beta1;
157
158    let eps = F::epsilon();
159
160    for itn in 0..options.max_iter {
161        // Lanczos iteration
162        let s = F::sparse_one() / beta;
163        for i in 0..n {
164            v[i] = s * y_vec[i];
165        }
166
167        let mut y_new = a.matvec(&v)?;
168
169        if itn >= 1 {
170            for i in 0..n {
171                y_new[i] -= (beta / oldb) * r1[i];
172            }
173        }
174
175        let alfa = inner(&v, &y_new);
176
177        for i in 0..n {
178            y_new[i] -= (alfa / beta) * r2[i];
179        }
180
181        r1 = r2;
182        r2 = y_new;
183
184        y_vec = if let Some(ref prec) = options.preconditioner {
185            prec.matvec(&r2)?
186        } else {
187            r2.clone()
188        };
189
190        oldb = beta;
191        let beta_sq = inner(&r2, &y_vec);
192
193        if beta_sq < F::sparse_zero() {
194            return Err(SparseError::ComputationError(
195                "Non-symmetric matrix".to_string(),
196            ));
197        }
198
199        beta = beta_sq.sqrt();
200        tnorm2 = tnorm2 + alfa * alfa + oldb * oldb + beta * beta;
201
202        // Apply previous rotation
203        let oldeps = epsln;
204        let delta = cs * dbar + sn * alfa;
205        let gbar = sn * dbar - cs * alfa;
206        epsln = sn * beta;
207        dbar = -cs * beta;
208
209        // Compute next rotation
210        let gamma = (gbar * gbar + beta * beta).sqrt();
211        let gamma_clamped = if gamma > eps { gamma } else { eps };
212        cs = gbar / gamma_clamped;
213        sn = beta / gamma_clamped;
214        let phi = cs * phibar;
215        phibar = sn * phibar;
216
217        // Update solution
218        let denom = F::sparse_one() / gamma_clamped;
219        let w1 = w2;
220        w2 = w;
221        w = vec![F::sparse_zero(); n];
222        for i in 0..n {
223            w[i] = (v[i] - oldeps * w1[i] - delta * w2[i]) * denom;
224            x[i] += phi * w[i];
225        }
226
227        // Update residual norm
228        gmax = gmax.max(gamma);
229        gmin = gmin.min(gamma);
230        qrnorm = phibar;
231        let rnorm = qrnorm;
232
233        // Check convergence
234        let anorm = tnorm2.sqrt();
235        let ynorm = norm2(&x);
236
237        let test1 = if ynorm == F::sparse_zero() || anorm == F::sparse_zero() {
238            F::infinity()
239        } else {
240            rnorm / (anorm * ynorm)
241        };
242
243        if test1 <= options.rtol || rnorm <= options.atol {
244            return Ok(MINRESResult {
245                x,
246                iterations: itn + 1,
247                residual_norm: rnorm,
248                converged: true,
249                message: "Converged".to_string(),
250            });
251        }
252    }
253
254    Ok(MINRESResult {
255        x,
256        iterations: options.max_iter,
257        residual_norm: qrnorm,
258        converged: false,
259        message: "Maximum iterations reached".to_string(),
260    })
261}
262
263/// Compute residual r = b - Ax
264#[allow(dead_code)]
265fn compute_residual<F>(a: &dyn LinearOperator<F>, b: &[F], x: &[F]) -> SparseResult<Vec<F>>
266where
267    F: Float + NumAssign,
268{
269    let ax = a.matvec(x)?;
270    Ok(b.iter()
271        .zip(ax.iter())
272        .map(|(&bi, &axi)| bi - axi)
273        .collect())
274}
275
276/// Compute 2-norm of a vector
277#[allow(dead_code)]
278fn norm2<F: Float + Sum>(v: &[F]) -> F {
279    v.iter().map(|&x| x * x).sum::<F>().sqrt()
280}
281
282/// Compute inner product of two vectors
283#[allow(dead_code)]
284fn inner<F: Float + Sum>(a: &[F], b: &[F]) -> F {
285    a.iter().zip(b.iter()).map(|(&ai, &bi)| ai * bi).sum()
286}
287
288#[cfg(test)]
289mod tests {
290    use super::*;
291    use crate::csr::CsrMatrix;
292    use crate::linalg::interface::{AsLinearOperator, DiagonalOperator, IdentityOperator};
293
294    #[test]
295    fn test_minres_identity() {
296        // Test MINRES on identity matrix: I * x = b => x = b
297        let identity = IdentityOperator::<f64>::new(3);
298        let b = vec![1.0, 2.0, 3.0];
299        let options = MINRESOptions::default();
300        let result = minres(&identity, &b, options).unwrap();
301
302        println!("Identity test result: {:?}", result);
303        println!("Expected x: {:?}", b);
304        println!("Actual x: {:?}", result.x);
305
306        // Verify solution
307        let ax = identity.matvec(&result.x).unwrap();
308        println!("Ax = {:?}", ax);
309        println!("b  = {:?}", b);
310
311        assert!(result.converged);
312        for (i, (xi, bi)) in result.x.iter().zip(&b).enumerate() {
313            let diff = (xi - bi).abs();
314            println!("x[{}]: expected {}, got {}, diff {}", i, bi, xi, diff);
315            assert!(
316                diff < 1e-10,
317                "x[{}]: expected {}, got {}, diff {}",
318                i,
319                bi,
320                xi,
321                diff
322            );
323        }
324    }
325
326    #[test]
327    fn test_minres_diagonal() {
328        // Test MINRES on diagonal matrix
329        let diag = DiagonalOperator::new(vec![2.0, -3.0, 4.0]); // Note: can have negative values
330        let b = vec![2.0, -6.0, 12.0];
331        let options = MINRESOptions::default();
332        let result = minres(&diag, &b, options).unwrap();
333
334        println!("Result: {:?}", result);
335        println!("Solution x: {:?}", result.x);
336        println!("Iterations: {}", result.iterations);
337        println!("Converged: {}", result.converged);
338        println!("Residual norm: {}", result.residual_norm);
339
340        // Verify solution
341        let ax = diag.matvec(&result.x).unwrap();
342        println!("Ax = {:?}", ax);
343        println!("b  = {:?}", b);
344
345        assert!(result.converged);
346        let expected = vec![1.0, 2.0, 3.0];
347        for (xi, ei) in result.x.iter().zip(&expected) {
348            assert!((xi - ei).abs() < 1e-9, "Expected {}, got {}", ei, xi);
349        }
350    }
351
352    #[test]
353    fn test_minres_symmetric_indefinite() {
354        // Test MINRES on a symmetric indefinite matrix
355        // A = [ 4  -1  0]
356        //     [-1   2 -1]
357        //     [ 0  -1  4]
358        // This is symmetric but not positive definite
359        let rows = vec![0, 0, 1, 1, 1, 2, 2];
360        let cols = vec![0, 1, 0, 1, 2, 1, 2];
361        let data = vec![4.0, -1.0, -1.0, 2.0, -1.0, -1.0, 4.0];
362        let shape = (3, 3);
363
364        let matrix = CsrMatrix::new(data, rows, cols, shape).unwrap();
365        let op = matrix.as_linear_operator();
366
367        let b = vec![1.0, 0.0, 1.0];
368        let options = MINRESOptions::default();
369        let result = minres(op.as_ref(), &b, options).unwrap();
370
371        assert!(result.converged);
372        // Verify solution by checking Ax = b
373        let ax = op.matvec(&result.x).unwrap();
374        for (axi, bi) in ax.iter().zip(&b) {
375            assert!((axi - bi).abs() < 1e-9);
376        }
377    }
378
379    #[test]
380    fn test_minres_preconditioned() {
381        // Test MINRES with a preconditioner
382        let rows = vec![0, 0, 1, 1, 1, 2, 2];
383        let cols = vec![0, 1, 0, 1, 2, 1, 2];
384        let data = vec![4.0, -1.0, -1.0, 4.0, -1.0, -1.0, 4.0];
385        let shape = (3, 3);
386
387        let matrix = CsrMatrix::new(data, rows, cols, shape).unwrap();
388        let op = matrix.as_linear_operator();
389
390        let b = vec![1.0, 2.0, 3.0];
391
392        // Use Jacobi preconditioner (diagonal matrix)
393        let preconditioner = Box::new(DiagonalOperator::new(vec![1.0 / 4.0, 1.0 / 4.0, 1.0 / 4.0]));
394
395        let options = MINRESOptions {
396            preconditioner: Some(preconditioner),
397            ..Default::default()
398        };
399
400        let result = minres(op.as_ref(), &b, options).unwrap();
401
402        assert!(result.converged);
403        // Verify solution by checking Ax = b
404        let ax = op.matvec(&result.x).unwrap();
405        for (axi, bi) in ax.iter().zip(&b) {
406            assert!((axi - bi).abs() < 1e-9);
407        }
408    }
409}
scirs2_sparse/linalg/minres.rs

scirs2_sparse/linalg/
minres.rs
