1use std::io::{BufWriter, Write};
10use std::path::Path;
11
12use sprs::CsMat;
13
14use crate::{Error, Result};
15
16pub fn write_mtx(matrix: &CsMat<f64>, path: impl AsRef<Path>) -> Result<()> {
17 let path = path.as_ref();
18 if is_structurally_symmetric(matrix) {
19 write_symmetric_mtx(matrix, path)
20 } else {
21 sprs::io::write_matrix_market(path, matrix.view()).map_err(|e| Error::Mtx(e.to_string()))
22 }
23}
24
25fn write_symmetric_mtx(matrix: &CsMat<f64>, path: &Path) -> Result<()> {
26 let f = std::fs::File::create(path)?;
27 let mut w = BufWriter::new(f);
28 writeln!(w, "%%MatrixMarket matrix coordinate real symmetric")?;
29 writeln!(w, "% written by powerio")?;
30
31 let nnz = matrix
33 .iter()
34 .filter(|&(_, (i, j))| i >= j)
35 .filter(|&(&v, _)| v != 0.0)
36 .count();
37 writeln!(w, "{} {} {}", matrix.rows(), matrix.cols(), nnz)?;
38
39 for (&v, (i, j)) in matrix {
40 if i < j || v == 0.0 {
41 continue;
42 }
43 writeln!(w, "{} {} {:.16e}", i + 1, j + 1, v)?;
44 }
45 Ok(())
46}
47
48pub fn read_mtx(path: impl AsRef<Path>) -> Result<CsMat<f64>> {
50 let tri: sprs::TriMat<f64> =
51 sprs::io::read_matrix_market(path).map_err(|e| Error::Mtx(e.to_string()))?;
52 Ok(tri.to_csr())
53}
54
55pub fn read_vector_mtx(path: impl AsRef<Path>) -> Result<Vec<f64>> {
58 let text = std::fs::read_to_string(path)?;
59 let mut lines = text.lines().filter(|l| !l.starts_with('%'));
60 let header = lines
61 .next()
62 .ok_or_else(|| Error::Mtx("empty vector file".into()))?;
63 let len: usize = header
64 .split_whitespace()
65 .next()
66 .and_then(|t| t.parse().ok())
67 .ok_or_else(|| Error::Mtx(format!("bad vector dimensions line: {header:?}")))?;
68 let values = lines
69 .take(len)
70 .map(|l| {
71 l.trim()
72 .parse::<f64>()
73 .map_err(|_| Error::Mtx(format!("bad vector entry: {l:?}")))
74 })
75 .collect::<Result<Vec<_>>>()?;
76 if values.len() != len {
77 return Err(Error::Mtx(format!(
78 "expected {len} entries, got {}",
79 values.len()
80 )));
81 }
82 Ok(values)
83}
84
85pub fn write_vector_mtx(values: &[f64], path: impl AsRef<Path>) -> Result<()> {
87 let f = std::fs::File::create(path)?;
88 let mut w = BufWriter::new(f);
89 writeln!(w, "%%MatrixMarket matrix array real general")?;
90 writeln!(w, "% written by powerio")?;
91 writeln!(w, "{} 1", values.len())?;
92 for v in values {
93 writeln!(w, "{v:.16e}")?;
94 }
95 Ok(())
96}
97
98fn is_structurally_symmetric(a: &CsMat<f64>) -> bool {
99 if a.rows() != a.cols() {
100 return false;
101 }
102 for (i, row) in a.outer_iterator().enumerate() {
103 for (j, &v) in row.iter() {
104 let mirror = a.get(j, i).copied().unwrap_or(0.0);
105 if (v - mirror).abs() > 1e-12 {
106 return false;
107 }
108 }
109 }
110 true
111}