1use std::env;
2use std::error::Error;
3use std::fs::{self, File};
4use std::io::{BufRead, BufReader};
5use std::path::{Path, PathBuf};
6use std::time::Instant;
7
8use hybit::{
9 analyze_csr32, read_matrix_market, Csr32Matrix, HybitSolver, ParallelCsr32Operator,
10 RigidBodyAggregation, SolverOptions, StructuralOptions, StructuralPcgVectorPolicy,
11 StructuralPreconditionerPolicy, StructuralSpmvPolicy,
12};
13
14#[derive(Debug)]
15struct Args {
16 matrix: PathBuf,
17 coordinates: PathBuf,
18 rhs: Option<PathBuf>,
19 relative_tolerance: f64,
20 max_iterations: usize,
21 target_coarse_dimension: usize,
22 aggregation: RigidBodyAggregation,
23 spmv_policy: StructuralSpmvPolicy,
24 preconditioner_policy: StructuralPreconditionerPolicy,
25 pcg_vector_policy: StructuralPcgVectorPolicy,
26}
27
28impl Args {
29 fn parse() -> Result<Self, Box<dyn Error>> {
30 let mut matrix = None;
31 let mut coordinates = None;
32 let mut rhs = None;
33 let mut relative_tolerance: f64 = 1.0e-8;
34 let mut max_iterations = 3000usize;
35 let mut target_coarse_dimension = 1536usize;
36 let mut aggregation = RigidBodyAggregation::Auto;
37 let mut spmv_policy = StructuralSpmvPolicy::Auto;
38 let mut preconditioner_policy = StructuralPreconditionerPolicy::Auto;
39 let mut pcg_vector_policy = StructuralPcgVectorPolicy::Auto;
40 let mut it = env::args().skip(1);
41 while let Some(arg) = it.next() {
42 match arg.as_str() {
43 "--matrix" => matrix = Some(PathBuf::from(next_value(&mut it, "--matrix")?)),
44 "--coords" => coordinates = Some(PathBuf::from(next_value(&mut it, "--coords")?)),
45 "--rhs" => rhs = Some(PathBuf::from(next_value(&mut it, "--rhs")?)),
46 "--tol" => relative_tolerance = next_value(&mut it, "--tol")?.parse()?,
47 "--max-iters" => max_iterations = next_value(&mut it, "--max-iters")?.parse()?,
48 "--target-coarse-dim" => {
49 target_coarse_dimension = next_value(&mut it, "--target-coarse-dim")?.parse()?
50 }
51 "--aggregation" => {
52 aggregation = match next_value(&mut it, "--aggregation")?
53 .to_ascii_lowercase()
54 .as_str()
55 {
56 "auto" => RigidBodyAggregation::Auto,
57 "contiguous" => RigidBodyAggregation::Contiguous,
58 "graph" => RigidBodyAggregation::Graph,
59 other => {
60 return Err(format!(
61 "unknown aggregation '{other}'; use auto, contiguous, or graph"
62 )
63 .into())
64 }
65 }
66 }
67 "--spmv" => {
68 spmv_policy = match next_value(&mut it, "--spmv")?.to_ascii_lowercase().as_str()
69 {
70 "auto" => StructuralSpmvPolicy::Auto,
71 "serial" => StructuralSpmvPolicy::Serial,
72 "parallel" => StructuralSpmvPolicy::Parallel,
73 other => {
74 return Err(format!(
75 "unknown SpMV policy '{other}'; use auto, serial, or parallel"
76 )
77 .into())
78 }
79 }
80 }
81 "--precond" => {
82 let value = next_value(&mut it, "--precond")?;
83 preconditioner_policy = match value.to_ascii_lowercase().as_str() {
84 "auto" => StructuralPreconditionerPolicy::Auto,
85 "serial" => StructuralPreconditionerPolicy::Serial,
86 "parallel" => StructuralPreconditionerPolicy::Parallel,
87 other => {
88 return Err(format!(
89 "unknown preconditioner policy '{other}'; use auto, serial, or parallel"
90 )
91 .into())
92 }
93 };
94 }
95 "--pcg-vectors" => {
96 pcg_vector_policy = match next_value(&mut it, "--pcg-vectors")?
97 .to_ascii_lowercase()
98 .as_str()
99 {
100 "auto" => StructuralPcgVectorPolicy::Auto,
101 "serial" => StructuralPcgVectorPolicy::Serial,
102 "parallel" => StructuralPcgVectorPolicy::Parallel,
103 other => {
104 return Err(format!(
105 "unknown PCG vector policy '{other}'; use auto, serial, or parallel"
106 )
107 .into())
108 }
109 }
110 }
111 "-h" | "--help" => {
112 print_usage();
113 std::process::exit(0);
114 }
115 other if !other.starts_with('-') && matrix.is_none() => {
116 matrix = Some(PathBuf::from(other))
117 }
118 other => return Err(format!("unknown argument '{other}'").into()),
119 }
120 }
121 let matrix = matrix.ok_or("missing matrix path; use --matrix FILE.mtx")?;
122 let coordinates = coordinates.unwrap_or_else(|| matrix.with_extension("coords"));
123 if !relative_tolerance.is_finite() || relative_tolerance <= 0.0 {
124 return Err("--tol must be finite and > 0".into());
125 }
126 if max_iterations == 0 {
127 return Err("--max-iters must be > 0".into());
128 }
129 if target_coarse_dimension < 6 {
130 return Err("--target-coarse-dim must be >= 6".into());
131 }
132 Ok(Self {
133 matrix,
134 coordinates,
135 rhs,
136 relative_tolerance,
137 max_iterations,
138 target_coarse_dimension,
139 aggregation,
140 spmv_policy,
141 preconditioner_policy,
142 pcg_vector_policy,
143 })
144 }
145}
146
147fn next_value<I: Iterator<Item = String>>(
148 it: &mut I,
149 flag: &str,
150) -> Result<String, Box<dyn Error>> {
151 it.next()
152 .ok_or_else(|| format!("missing value after {flag}").into())
153}
154
155fn print_usage() {
156 println!("HyBIT structural-auto FEM benchmark");
157 println!("Usage: fem_structural_auto --matrix K.mtx [--coords K.coords] [--rhs b.txt] [--tol 1e-8] [--max-iters 3000] [--target-coarse-dim 1536] [--aggregation auto|contiguous|graph] [--spmv auto|serial|parallel] [--precond auto|serial|parallel] [--pcg-vectors auto|serial|parallel]");
158}
159
160fn read_coordinates(path: &Path) -> Result<Vec<[f64; 3]>, Box<dyn Error>> {
161 let file = File::open(path)?;
162 let reader = BufReader::new(file);
163 let mut expected = None::<usize>;
164 let mut coordinates = Vec::new();
165
166 for (line_no, line) in reader.lines().enumerate() {
167 let line = line?;
168 let text = line.trim();
169 if text.is_empty() || text.starts_with('#') {
170 continue;
171 }
172 if expected.is_none() {
173 expected = Some(text.parse::<usize>().map_err(|e| {
174 format!(
175 "{}:{}: invalid coordinate count: {e}",
176 path.display(),
177 line_no + 1
178 )
179 })?);
180 coordinates.reserve(expected.unwrap());
181 continue;
182 }
183 let fields: Vec<&str> = text.split_whitespace().collect();
184 if fields.len() != 3 {
185 return Err(format!(
186 "{}:{}: expected three coordinates",
187 path.display(),
188 line_no + 1
189 )
190 .into());
191 }
192 let x: f64 = fields[0].parse()?;
193 let y: f64 = fields[1].parse()?;
194 let z: f64 = fields[2].parse()?;
195 if !x.is_finite() || !y.is_finite() || !z.is_finite() {
196 return Err(
197 format!("{}:{}: non-finite coordinate", path.display(), line_no + 1).into(),
198 );
199 }
200 coordinates.push([x, y, z]);
201 }
202
203 let expected =
204 expected.ok_or_else(|| format!("{}: missing coordinate count", path.display()))?;
205 if coordinates.len() != expected {
206 return Err(format!(
207 "{}: coordinate count mismatch: header says {}, read {}",
208 path.display(),
209 expected,
210 coordinates.len()
211 )
212 .into());
213 }
214 Ok(coordinates)
215}
216
217fn load_rhs(path: &Path, n: usize) -> Result<Vec<f64>, Box<dyn Error>> {
218 let text = fs::read_to_string(path)?;
219 let mut values = Vec::with_capacity(n);
220
221 for (line_index, raw_line) in text.lines().enumerate() {
222 let line = raw_line.trim_start_matches('\u{feff}').trim();
225 if line.is_empty() || line.starts_with('#') || line.starts_with('%') {
226 continue;
227 }
228 for (token_index, token) in line.split_whitespace().enumerate() {
229 let value = token.parse::<f64>().map_err(|e| {
230 format!(
231 "{}: invalid RHS float at line {}, token {}: {:?} ({e})",
232 path.display(),
233 line_index + 1,
234 token_index + 1,
235 token
236 )
237 })?;
238 if !value.is_finite() {
239 return Err(format!(
240 "{}: non-finite RHS value at line {}, token {}: {:?}",
241 path.display(),
242 line_index + 1,
243 token_index + 1,
244 token
245 )
246 .into());
247 }
248 values.push(value);
249 }
250 }
251
252 if values.len() != n {
253 return Err(format!(
254 "{}: RHS length mismatch: expected {n}, got {}",
255 path.display(),
256 values.len()
257 )
258 .into());
259 }
260 Ok(values)
261}
262
263fn norm2(x: &[f64]) -> f64 {
264 x.iter().map(|v| v * v).sum::<f64>().sqrt()
265}
266
267fn verified_relative_residual(
268 a: &Csr32Matrix,
269 b: &[f64],
270 x: &[f64],
271) -> Result<f64, Box<dyn Error>> {
272 let ax = a.spmv(x)?;
273 let sum = b
274 .iter()
275 .zip(ax.iter())
276 .map(|(&bi, &ai)| {
277 let r = bi - ai;
278 r * r
279 })
280 .sum::<f64>();
281 let denom = norm2(b);
282 Ok(if denom == 0.0 {
283 sum.sqrt()
284 } else {
285 sum.sqrt() / denom
286 })
287}
288
289fn relative_error_to_ones(x: &[f64]) -> f64 {
290 let diff = x
291 .iter()
292 .map(|&xi| {
293 let d = xi - 1.0;
294 d * d
295 })
296 .sum::<f64>();
297 diff.sqrt() / (x.len() as f64).sqrt().max(f64::MIN_POSITIVE)
298}
299
300fn mib(bytes: usize) -> f64 {
301 bytes as f64 / (1024.0 * 1024.0)
302}
303
304fn main() -> Result<(), Box<dyn Error>> {
305 let args = Args::parse()?;
306 println!(
307 "HyBIT {} structural-auto FEM benchmark",
308 env!("CARGO_PKG_VERSION")
309 );
310 println!("matrix : {}", args.matrix.display());
311 println!("coordinates : {}", args.coordinates.display());
312
313 let load_start = Instant::now();
314 let (matrix, mm) = read_matrix_market(&args.matrix)?;
315 let matrix_load_seconds = load_start.elapsed().as_secs_f64();
316 let coord_start = Instant::now();
317 let coordinates = read_coordinates(&args.coordinates)?;
318 let coord_load_seconds = coord_start.elapsed().as_secs_f64();
319 let profile = analyze_csr32(&matrix)?;
320
321 println!(
322 "Matrix Market : {:?}, {} input entries -> {} CSR nnz",
323 mm.symmetry, mm.input_entries, mm.csr_nnz
324 );
325 println!("dimensions : {} x {}", profile.nrows, profile.ncols);
326 println!("nnz : {}", profile.nnz);
327 println!(
328 "CSR storage : {:.3} MiB",
329 mib(matrix.storage_bytes())
330 );
331 println!("matrix load : {:.3} ms", matrix_load_seconds * 1.0e3);
332 println!("coordinate nodes : {}", coordinates.len());
333 println!("coordinate load : {:.3} ms", coord_load_seconds * 1.0e3);
334
335 if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
336 return Err(
337 "rigid-body two-level PCG requires a square matrix with a complete positive diagonal"
338 .into(),
339 );
340 }
341 if matrix.nrows() != coordinates.len() * 3 {
342 return Err(format!(
343 "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
344 matrix.nrows(),
345 coordinates.len()
346 )
347 .into());
348 }
349
350 let generated_rhs = args.rhs.is_none();
351 let b = if let Some(path) = args.rhs.as_deref() {
352 println!("RHS : {}", path.display());
353 load_rhs(path, matrix.nrows())?
354 } else {
355 println!("RHS : generated as b=A*1 (known exact solution)");
356 let ones = vec![1.0; matrix.ncols()];
357 matrix.spmv(&ones)?
358 };
359 let mut solver = HybitSolver::new();
360 solver.set_options(SolverOptions {
361 relative_tolerance: args.relative_tolerance,
362 absolute_tolerance: 0.0,
363 max_iterations: args.max_iterations,
364 })?;
365 solver.set_structural_options(StructuralOptions {
366 target_coarse_dimension: args.target_coarse_dimension,
367 aggregation: args.aggregation,
368 spmv_policy: args.spmv_policy,
369 preconditioner_policy: args.preconditioner_policy,
370 pcg_vector_policy: args.pcg_vector_policy,
371 })?;
372
373 let analysis = solver.analyze_csr32(&matrix)?;
374 let mut prepared = solver.prepare_structural_csr32(&matrix, &analysis, &coordinates)?;
375 println!("policy : StructuralAuto/RigidBodyTwoLevel");
376 println!("aggregation : {:?}", prepared.aggregation());
377 println!("SpMV policy : {:?}", prepared.spmv_policy());
378 println!(
379 "precond policy : {:?}",
380 prepared.structural_preconditioner_policy()
381 );
382 println!("PCG vector policy : {:?}", prepared.pcg_vector_policy());
383 if prepared.parallel_spmv_enabled()
384 || prepared.parallel_preconditioner_enabled()
385 || prepared.parallel_pcg_vectors_enabled()
386 {
387 println!(
388 "Rayon threads : {}",
389 ParallelCsr32Operator::new(&matrix).rayon_threads()
390 );
391 }
392 if prepared.parallel_preconditioner_enabled() {
393 println!(
394 "parallel index : {:.3} MiB",
395 mib(prepared.parallel_preconditioner_index_bytes())
396 );
397 }
398 println!("target coarse dim : {}", args.target_coarse_dimension);
399 println!(
400 "aggregate nodes : {} (auto-selected)",
401 prepared.aggregate_nodes()
402 );
403 println!("fine block size : 3");
404 println!("aggregate count : {}", prepared.aggregate_count());
405 println!(
406 "aggregate min/max : {} / {} nodes",
407 prepared.min_aggregate_nodes(),
408 prepared.max_aggregate_nodes()
409 );
410 println!("modes/aggregate : 6");
411 println!("coarse dimension : {}", prepared.coarse_dimension());
412 println!(
413 "base factor : {:.3} MiB",
414 mib(prepared.base_factor_bytes())
415 );
416 println!(
417 "coarse factor : {:.3} MiB",
418 mib(prepared.coarse_factor_bytes())
419 );
420 println!(
421 "geometry storage : {:.3} MiB",
422 mib(prepared.geometry_bytes())
423 );
424 println!(
425 "total prec storage : {:.3} MiB",
426 mib(prepared.preconditioner_bytes())
427 );
428 println!(
429 "analysis : {:.3} ms",
430 prepared.analysis_seconds() * 1.0e3
431 );
432 println!(
433 "prepare : {:.3} ms",
434 prepared.prepare_seconds() * 1.0e3
435 );
436
437 let mut x = vec![0.0; matrix.ncols()];
438 let report = prepared.solve(&matrix, &b, &mut x)?;
439 let verified = verified_relative_residual(&matrix, &b, &x)?;
440
441 println!();
442 println!("Structural Auto API: 3x3 Block-Jacobi + six rigid-body coarse modes");
443 println!("status : {:?}", report.status);
444 println!("preconditioner : {:?}", report.preconditioner);
445 println!("iterations : {}", report.iterations);
446 println!("reported residual : {:.6e}", report.relative_residual);
447 println!("verified residual : {:.6e}", verified);
448 if generated_rhs {
449 println!("relative x error : {:.6e}", relative_error_to_ones(&x));
450 }
451 println!(
452 "solve time : {:.3} ms",
453 report.solve_seconds * 1.0e3
454 );
455 println!(
456 "total setup+solve : {:.3} ms",
457 (report.setup_seconds + report.solve_seconds) * 1.0e3
458 );
459
460 if !verified.is_finite() {
461 return Err("non-finite independently verified residual".into());
462 }
463 Ok(())
464}