support integer arithmetic

[dvc94ch]

seems like it would be interesting to check if an algorithm was implemented correctly. my conjugate gradient implementation doesn't seem to converge very well and was wondering if it's just due to floating point error or something else.

[sarah-ek]

how would field operations work with integers?

[dvc94ch]

not sure what you mean. my understanding of a field is +/* and the inverses need to be defined? I guess thechnically finite fields are also fields even if they may not be what people expect? alternatively fixedpoint or something like [0] might work. I guess it doesn't need to be fast, just accurate.

• [0] https://docs.rs/malachite-float/latest/malachite_float/struct.Float.html

[sarah-ek]

data types are required to be Copy and Pod since that's used for many optimizations.

so heap allocated multiprecision floats are unsupported. you can implement stack allocated big floats, but that is out of scope for faer and should be done in a separate library.

regarding conjugate gradient, there is currently a working implementation in src/linop/conjugate_gradient.rs that you can use to compare with yours, or wait until the next version is released

[dvc94ch]

Oh, cool, will try that...

[dvc94ch]

tried cg and bicgstab. I know I'm probably doing something wrong, but I'm not very impressed with the results of krylov methods so far. (bicgstab doesn't converge even for the simplest matrices, and cg seems to take longer than the theoretical matrix dim). I guess I'll stick to direct methods for now as they're probably much faster and more reliable for small problem sizes even if the matrix is sparse. (the main reason I was interested in krylov methods is because it seemed like they could theoretically achieve a better accuracy than direct methods since they could be as accurate as computing the residual allows)

[sarah-ek]

can you share the code you tried?

[dvc94ch]

#[derive(Debug)]
struct Sparse(SparseColMat<usize, c64>);
impl LinOp<c64> for Sparse {
    fn nrows(&self) -> usize {
        self.0.nrows()
    }

    fn ncols(&self) -> usize {
        self.0.ncols()
    }

    fn apply_req(
        &self,
        _rhs_ncols: usize,
        _parallelism: Parallelism,
    ) -> Result<StackReq, SizeOverflow> {
        Ok(StackReq::empty())
    }

    fn apply(
        &self,
        out: MatMut<'_, c64>,
        rhs: MatRef<'_, c64>,
        parallelism: Parallelism,
        _stack: PodStack<'_>,
    ) {
        matmul::sparse_dense_matmul(
            out,
            self.0.as_ref(),
            rhs,
            None,
            c64::new(1., 0.),
            parallelism,
        );
    }

    fn conj_apply(
        &self,
        out: MatMut<'_, c64>,
        rhs: MatRef<'_, c64>,
        parallelism: Parallelism,
        _stack: PodStack<'_>,
    ) {
        matmul::sparse_dense_matmul(
            out,
            self.0.as_ref().adjoint().to_col_major().unwrap().as_ref(),
            rhs,
            None,
            c64::new(1., 0.),
            parallelism,
        );
    }
}

#[allow(unused)]
fn solve_bicgstab(a: SparseColMat<usize, c64>, y: Col<c64>) -> Result<(Col<c64>, SolutionInfo)> {
    println!("{a:?} x = {y:?}");
    let dim = y.nrows();
    let precond = IdentityPrecond { dim };
    let a = Sparse(a);
    let mut params = BicgParams::default();
    params.max_iters = dim;
    let mut pod_buffer = GlobalPodBuffer::new(
        bicgstab_req::<c64>(&precond, &precond, &a, 1, Parallelism::None).unwrap(),
    );
    let podstack = PodStack::new(&mut pod_buffer);
    let mut x = Col::zeros(dim);
    let info = bicgstab(
        x.as_2d_mut(),
        precond,
        precond,
        a,
        y.as_2d_ref(),
        params,
        Parallelism::None,
        podstack,
    )
    .map_err(|e| anyhow::anyhow!("{:?}", e))?;
    let info = SolutionInfo {
        abs_residual: info.abs_residual,
        rel_residual: info.rel_residual,
        iter_count: info.iter_count,
    };
    Ok((x, info))
}

[
  (0, 0, 1.0 + 0.0 * I), (1, 0, -1.0 + 0.0 * I), (2, 0, -1.0 - 0.0 * I),
  (0, 1, 0.0 + 0.0 * I), (1, 1, 1.0 + 0.0 * I), (2, 1, 1.0 + 0.0 * I), 
  (0, 2, 0.0 + 0.0 * I), (1, 2, 1.0 + 0.0 * I)
] x = [
  [0.0 + 0.0 * I],
  [0.0 + 0.0 * I],
  [1.0 + 0.0 * I],
]
thread 'main' panicked at src/main.rs:10:33:
called `Result::unwrap()` on an `Err` value: NoConvergence { abs_residual: 1.0, rel_residual: 1.0 }
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

[sarah-ek]

thanks! i'll look into it soon and report back

[sarah-ek]

@dvc94ch the code seems to be correct. the problem comes from the original bicgstab algorithm that seems to perform a division by zero in the first iteration due to the residual and its image by A being orthogonal. the same thing happens in the eigen and scipy implementations as well, which i guess makes bicgstab not that good of a solver in the general case

[dvc94ch]

Ah thought it was a general solver. What would work better? bicgstab(n) which if I understand correctly uses gmres every few iterations to stabilize? Or restarted gmres? I thought the main advantage over gmres was that it doesn't need more memory every iteration.

Was building a linear elasticity solver to modulate the infill in slicing of 3d parts, but got side tracked building a toy circuit simulator to practice some numerical methods. The above matrix is a modified nodal analysis matrix for a 1V source and 1ohm resistor. I think for linear elasticity cg works well, but for a circuit simulator need something more general.

[sarah-ek]

im still relatively new to iterative solvers so i don't know what would be better. but once i implement gmres i'll check if it also has any weird edge cases.

in any case if your system can be solved by a direct solver, that's what i would recommend, as they're both faster and usually more accurate