metis/
lib.rs

1//! This crate provides a thin but idiomatic API around libmetis.
2//!
3//! See [`Graph`] for a usage example.
4
5#![deny(missing_docs)]
6
7use crate::option::Opt;
8use metis_sys as m;
9use std::convert::TryFrom;
10use std::fmt;
11use std::mem;
12use std::os;
13use std::ptr;
14use std::result::Result as StdResult;
15use std::slice;
16
17pub mod option;
18
19#[cfg(target_pointer_width = "16")]
20compile_error!("METIS does not support 16-bit architectures");
21
22/// Integer type used by METIS, can either be an [`i32`] or an [`i64`].
23pub type Idx = m::idx_t;
24
25/// Floating-point type used by METIS, can either be an [`f32`] or an [`f64`].
26pub type Real = m::real_t;
27
28/// The length of the `options` array.
29///
30/// See [`Graph::set_options`] for an example.  It is also used in
31/// [`Mesh::set_options`].
32pub const NOPTIONS: usize = m::METIS_NOPTIONS as usize;
33
34/// Error type returned by METIS.
35#[derive(Debug, PartialEq, Eq)]
36pub enum Error {
37    /// Input is invalid.
38    ///
39    /// These bindings should check for most input errors, if not all.
40    Input,
41
42    /// METIS hit an out-of-memory error.
43    Memory,
44
45    /// METIS returned an error but its meaning is unknown.
46    Other,
47}
48
49impl std::error::Error for Error {}
50
51impl From<NewGraphError> for Error {
52    fn from(_: NewGraphError) -> Self {
53        Self::Input
54    }
55}
56
57impl From<NewMeshError> for Error {
58    fn from(_: NewMeshError) -> Self {
59        Self::Input
60    }
61}
62
63impl fmt::Display for Error {
64    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
65        match self {
66            Error::Input => write!(f, "invalid input"),
67            Error::Memory => write!(f, "out of memory"),
68            Error::Other => write!(f, "METIS returned an error"),
69        }
70    }
71}
72
73/// The result of a partitioning.
74pub type Result<T> = StdResult<T, Error>;
75
76trait ErrorCode {
77    /// Makes a [`Result`] from a return code (int) from METIS.
78    fn wrap(self) -> Result<()>;
79}
80
81impl ErrorCode for m::rstatus_et {
82    fn wrap(self) -> Result<()> {
83        match self {
84            m::rstatus_et_METIS_OK => Ok(()),
85            m::rstatus_et_METIS_ERROR_INPUT => Err(Error::Input),
86            m::rstatus_et_METIS_ERROR_MEMORY => Err(Error::Memory),
87            m::rstatus_et_METIS_ERROR => Err(Error::Other),
88            other => panic!("unexpected error code ({}) from METIS", other),
89        }
90    }
91}
92
93/// Error raised when the graph data fed to [`Graph::new`] cannot be safely
94/// passed to METIS.
95///
96/// Graph data must follow the format described in [`Graph::new`].
97#[derive(Debug)]
98pub struct InvalidGraphError {
99    msg: &'static str,
100}
101
102impl fmt::Display for InvalidGraphError {
103    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
104        self.msg.fmt(f)
105    }
106}
107
108/// Error type returned by [`Graph::new`].
109///
110/// Unlike [`Error`], this error originates from the Rust bindings.
111#[derive(Debug)]
112#[non_exhaustive]
113pub enum NewGraphError {
114    /// `ncon` must be greater than 1.
115    NoConstraints,
116
117    /// `nparts` must be greater than 1.
118    NoParts,
119
120    /// Graph is too large. One of the array's length doesn't fit into [`Idx`].
121    TooLarge,
122
123    /// The input arrays are malformed and cannot be safely passed to METIS.
124    ///
125    /// Note that these bindings do not check for all the invariants. Some might
126    /// be raised during [`Graph::part_recursive`] and [`Graph::part_kway`] as
127    /// [`Error::Input`].
128    InvalidGraph(InvalidGraphError),
129}
130
131impl fmt::Display for NewGraphError {
132    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
133        match self {
134            Self::NoConstraints => write!(f, "there must be at least one constraint"),
135            Self::NoParts => write!(f, "there must be at least one part"),
136            Self::TooLarge => write!(f, "graph is too large"),
137            Self::InvalidGraph(err) => write!(f, "invalid graph structure: {err}"),
138        }
139    }
140}
141
142impl std::error::Error for NewGraphError {}
143
144impl NewGraphError {
145    fn msg(msg: &'static str) -> Self {
146        Self::InvalidGraph(InvalidGraphError { msg })
147    }
148}
149
150/// Helper function to convert an immutable slice ref to a mutable pointer
151unsafe fn slice_to_mut_ptr<T>(slice: &[T]) -> *mut T {
152    slice.as_ptr() as *mut T
153}
154
155/// Builder structure to set up a graph partition computation.
156///
157/// This structure holds the required arguments for METIS to compute a
158/// partition.  It also offers methods to easily set any optional argument.
159///
160/// # Example
161///
162/// ```rust
163/// # fn main() -> Result<(), metis::Error> {
164/// # use metis::Graph;
165/// // Make a graph with two vertices and an edge between the two.
166/// let xadj = &mut [0, 1, 2];
167/// let adjncy = &mut [1, 0];
168///
169/// // Allocate the partition array which stores the partition of each vertex.
170/// let mut part = [0, 0];
171///
172/// // There are one constraint and two parts.  The partitioning algorithm used
173/// // is recursive bisection.  The k-way algorithm can also be used.
174/// Graph::new(1, 2, xadj, adjncy)?
175///     .part_recursive(&mut part)?;
176///
177/// // The two vertices are placed in different parts.
178/// assert_ne!(part[0], part[1]);
179/// # Ok(())
180/// # }
181/// ```
182#[derive(Debug, PartialEq)]
183pub struct Graph<'a> {
184    /// The number of balancing constrains.
185    ncon: Idx,
186
187    /// The number of parts to partition the graph.
188    nparts: Idx,
189
190    /// The adjency structure of the graph (part 1).
191    xadj: &'a [Idx],
192
193    /// The adjency structure of the graph (part 2).
194    ///
195    /// Required size: xadj.last()
196    adjncy: &'a [Idx],
197
198    /// The computational weights of the vertices.
199    ///
200    /// Required size: ncon * (xadj.len()-1)
201    vwgt: Option<&'a [Idx]>,
202
203    /// The communication weights of the vertices.
204    ///
205    /// Required size: xadj.len()-1
206    vsize: Option<&'a [Idx]>,
207
208    /// The weight of the edges.
209    ///
210    /// Required size: xadj.last()
211    adjwgt: Option<&'a [Idx]>,
212
213    /// The target partition weights of the vertices.
214    ///
215    /// If `None` then the graph is equally divided among the partitions.
216    ///
217    /// Required size: ncon * nparts
218    tpwgts: Option<&'a [Real]>,
219
220    /// Imbalance tolerances for each constraint.
221    ///
222    /// Required size: ncon
223    ubvec: Option<&'a [Real]>,
224
225    /// Fine-tuning parameters.
226    options: [Idx; NOPTIONS],
227}
228
229impl<'a> Graph<'a> {
230    /// Creates a new [`Graph`] object to be partitioned.
231    ///
232    /// - `ncon` is the number of constraints on each vertex (at least 1),
233    /// - `nparts` is the number of parts wanted in the graph partition.
234    ///
235    /// # Input format
236    ///
237    /// CSR (Compressed Sparse Row) is a data structure for representing sparse
238    /// matrices and is the primary data structure used by METIS. A CSR
239    /// formatted graph is represented with two slices: an adjacency list
240    /// (`adjcny`) and an index list (`xadj`). The nodes adjacent to node `n`
241    /// are `adjncy[xadj[n]..xadj[n + 1]]`. Additionally, metis requires that
242    /// graphs are undirected: if `(u, v)` is in the graph, then `(v, u)` must
243    /// also be in the graph.
244    ///
245    /// Consider translating this simple graph to CSR format:
246    /// ```rust
247    /// // 5 - 3 - 4 - 0
248    /// //     |   | /
249    /// //     2 - 1
250    /// let adjncy = [1, 4, 0, 2, 4, 1, 3, 2, 4, 5, 0, 1, 3, 3];
251    /// let xadj = [0, 2, 5, 7, 10, 13, 14];
252    ///
253    /// // iterate over adjacent nodes
254    /// let mut it = xadj
255    ///     .windows(2)
256    ///     .map(|x| &adjncy[x[0]..x[1]]);
257    ///
258    /// // node 0 is adjacent to nodes 1 and 4
259    /// assert_eq!(it.next().unwrap(), &[1, 4]);
260    ///
261    /// // node 1 is adjacent to nodes 0, 2, and 4
262    /// assert_eq!(it.next().unwrap(), &[0, 2, 4]);
263    ///
264    /// // node 2 is adjacent to nodes 1 and 3
265    /// assert_eq!(it.next().unwrap(), &[1, 3]);
266    ///
267    /// // node 3 is adjacent to nodes 2, 4, and 5
268    /// assert_eq!(it.next().unwrap(), &[2, 4, 5]);
269    ///
270    /// // node 4 is adjacent to nodes 0, 1, and 3
271    /// assert_eq!(it.next().unwrap(), &[0, 1, 3]);
272    ///
273    /// // node 5 is adjacent to node 3
274    /// assert_eq!(it.next().unwrap(), &[3]);
275    ///
276    /// assert!(it.next().is_none());
277    /// ```
278    ///
279    /// More info can be found at:
280    /// <https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format)>
281    ///
282    /// # Errors
283    ///
284    /// The following invariants must be held, otherwise this function returns
285    /// an error:
286    ///
287    /// - all the arrays have a length that can be held by an [`Idx`],
288    /// - `ncon` is strictly greater than zero,
289    /// - `nparts` is strictly greater than zero,
290    /// - `xadj` has at least one element (its length is the one more than the
291    ///   number of vertices),
292    /// - `xadj` is sorted,
293    /// - elements of `xadj` are positive,
294    /// - the last element of `xadj` is the length of `adjncy`,
295    /// - elements of `adjncy` are within zero and the number of vertices.
296    ///
297    /// # Mutability
298    ///
299    /// [`Graph::part_kway`] and [`Graph::part_recursive`] may mutate the
300    /// contents of `xadj` and `adjncy`, but should revert all changes before
301    /// returning.
302    pub fn new(
303        ncon: Idx,
304        nparts: Idx,
305        xadj: &'a [Idx],
306        adjncy: &'a [Idx],
307    ) -> StdResult<Graph<'a>, NewGraphError> {
308        if ncon <= 0 {
309            return Err(NewGraphError::NoConstraints);
310        }
311        if nparts <= 0 {
312            return Err(NewGraphError::NoParts);
313        }
314
315        let last_xadj = *xadj
316            .last()
317            .ok_or(NewGraphError::msg("index list is empty"))?;
318        let adjncy_len = Idx::try_from(adjncy.len()).map_err(|_| NewGraphError::TooLarge)?;
319        if last_xadj != adjncy_len {
320            return Err(NewGraphError::msg(
321                "length mismatch between index and adjacency lists",
322            ));
323        }
324
325        let nvtxs = match Idx::try_from(xadj.len()) {
326            Ok(xadj_len) => xadj_len - 1,
327            Err(_) => {
328                return Err(NewGraphError::TooLarge);
329            }
330        };
331
332        let mut prev = 0;
333        for x in xadj {
334            if prev > *x {
335                return Err(NewGraphError::msg("index list is not sorted"));
336            }
337            prev = *x;
338        }
339
340        for a in adjncy {
341            if *a < 0 || *a >= nvtxs {
342                return Err(NewGraphError::msg(
343                    "some values in the adjacency list are out of bounds",
344                ));
345            }
346        }
347
348        Ok(unsafe { Graph::new_unchecked(ncon, nparts, xadj, adjncy) })
349    }
350
351    /// Creates a new [`Graph`] object to be partitioned (unchecked version).
352    ///
353    /// - `ncon` is the number of constraints on each vertex (at least 1),
354    /// - `nparts` is the number of parts wanted in the graph partition.
355    ///
356    /// # Input format
357    ///
358    /// `xadj` and `adjncy` are the [CSR encoding][0] of the adjacency matrix
359    /// that represents the graph. `xadj` is the row index and `adjncy` is the
360    /// column index.
361    ///
362    /// [0]: https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format)
363    ///
364    /// # Safety
365    ///
366    /// This function still does some checks listed in "Panics" below. However,
367    /// the caller is reponsible for upholding all invariants listed in the
368    /// "Errors" section of [`Graph::new`]. Otherwise, the behavior of this
369    /// function is undefined.
370    ///
371    /// # Panics
372    ///
373    /// This function panics if:
374    /// - any of the arrays have a length that cannot be held by an [`Idx`], or
375    /// - `ncon` is not strictly greater than zero, or
376    /// - `nparts` is not strictly greater than zero, or
377    /// - `xadj` is empty, or
378    /// - the length of `adjncy` is different from the last element of `xadj`.
379    ///
380    /// # Mutability
381    ///
382    /// [`Graph::part_kway`] and [`Graph::part_recursive`] may mutate the
383    /// contents of `xadj` and `adjncy`, but should revert all changes before
384    /// returning.
385    pub unsafe fn new_unchecked(
386        ncon: Idx,
387        nparts: Idx,
388        xadj: &'a [Idx],
389        adjncy: &'a [Idx],
390    ) -> Graph<'a> {
391        assert!(0 < ncon, "ncon must be strictly greater than zero");
392        assert!(0 < nparts, "nparts must be strictly greater than zero");
393        let _ = Idx::try_from(xadj.len()).expect("xadj array larger than Idx::MAX");
394        assert_ne!(xadj.len(), 0);
395        let adjncy_len = Idx::try_from(adjncy.len()).expect("adjncy array larger than Idx::MAX");
396        assert_eq!(adjncy_len, *xadj.last().unwrap());
397
398        Graph {
399            ncon,
400            nparts,
401            xadj,
402            adjncy,
403            vwgt: None,
404            vsize: None,
405            adjwgt: None,
406            tpwgts: None,
407            ubvec: None,
408            options: [-1; NOPTIONS],
409        }
410    }
411
412    /// Sets the computational weights of the vertices.
413    ///
414    /// By default, all vertices have the same weight.
415    ///
416    /// The `ncon` weights of the `i`th vertex must be located in
417    /// `vwgt[i*ncon..(i+1)*ncon]`, and all elements of `vwgt` must be positive.
418    ///
419    /// # Panics
420    ///
421    /// This function panics if the length of `vwgt` is not `ncon` times the
422    /// number of vertices.
423    pub fn set_vwgt(mut self, vwgt: &'a [Idx]) -> Graph<'a> {
424        let vwgt_len = Idx::try_from(vwgt.len()).expect("vwgt array too large");
425        assert_eq!(vwgt_len, self.ncon * (self.xadj.len() as Idx - 1));
426        self.vwgt = Some(vwgt);
427        self
428    }
429
430    /// Sets the communication weights of the vertices.
431    ///
432    /// By default, all vertices have the same communication weight.
433    ///
434    /// Vertices can only have one communication weight. The length of `vsize`
435    /// does not depend on `ncon`.
436    ///
437    /// # Panics
438    ///
439    /// This function panics if the length of `vsize` is not the number of
440    /// vertices.
441    pub fn set_vsize(mut self, vsize: &'a [Idx]) -> Graph<'a> {
442        let vsize_len = Idx::try_from(vsize.len()).expect("vsize array too large");
443        assert_eq!(vsize_len, self.xadj.len() as Idx - 1);
444        self.vsize = Some(vsize);
445        self
446    }
447
448    /// Sets the weights of the edges.
449    ///
450    /// By default, all edges have the same weight.
451    ///
452    /// All elements of `adjwgt` must be positive.
453    ///
454    /// # Panics
455    ///
456    /// This function panics if the length of `adjwgt` is not equal to the
457    /// length of `adjncy`.
458    pub fn set_adjwgt(mut self, adjwgt: &'a [Idx]) -> Graph<'a> {
459        let adjwgt_len = Idx::try_from(adjwgt.len()).expect("adjwgt array too large");
460        assert_eq!(adjwgt_len, *self.xadj.last().unwrap());
461        self.adjwgt = Some(adjwgt);
462        self
463    }
464
465    /// Sets the target partition weights for each part and constraint.
466    ///
467    /// By default, the graph is divided equally.
468    ///
469    /// The target partition weight for the `i`th part and `j`th constraint is
470    /// specified at `tpwgts[i*ncon+j]`. For each constraint `j`, the sum of the
471    /// target partition weights must be 1.0. Meaning
472    /// `(0..nparts).map(|i| tpwgts[i*ncon+j]).sum() == 1.0`.
473    ///
474    /// # Panics
475    ///
476    /// This function panics if the length of `tpwgts` is not equal to `ncon`
477    /// times `nparts`.
478    pub fn set_tpwgts(mut self, tpwgts: &'a [Real]) -> Graph<'a> {
479        let tpwgts_len = Idx::try_from(tpwgts.len()).expect("tpwgts array too large");
480        assert_eq!(tpwgts_len, self.ncon * self.nparts);
481        self.tpwgts = Some(tpwgts);
482        self
483    }
484
485    /// Sets the load imbalance tolerance for each constraint.
486    ///
487    /// By default, it equals to 1.001 if `ncon` equals 1 and 1.01 otherwise.
488    ///
489    /// For the `i`th partition and `j`th constraint the allowed weight is the
490    /// `ubvec[j]*tpwgts[i*ncon+j]` fraction of the `j`th's constraint total
491    /// weight. The load imbalances must be greater than 1.0.
492    ///
493    /// # Panics
494    ///
495    /// This function panics if the length of `ubvec` is not equal to `ncon`.
496    pub fn set_ubvec(mut self, ubvec: &'a [Real]) -> Graph<'a> {
497        let ubvec_len = Idx::try_from(ubvec.len()).expect("ubvec array too large");
498        assert_eq!(ubvec_len, self.ncon);
499        self.ubvec = Some(ubvec);
500        self
501    }
502
503    /// Sets the fine-tuning parameters for this partitioning.
504    ///
505    /// When few options are to be set, [`Graph::set_option`] might be a
506    /// better fit.
507    ///
508    /// See the [option] module for the list of available parameters.  Note that
509    /// not all are applicable to a given partitioning method.  Refer to the
510    /// documentation of METIS ([link]) for more info on this.
511    ///
512    /// [link]: http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/manual.pdf
513    ///
514    /// # Example
515    ///
516    /// ```rust
517    /// # fn main() -> Result<(), metis::Error> {
518    /// # use metis::Graph;
519    /// use metis::option::Opt as _;
520    ///
521    /// let xadj = &[0, 1, 2];
522    /// let adjncy = &[1, 0];
523    /// let mut part = [0, 0];
524    ///
525    /// // -1 is the default value.
526    /// let mut options = [-1; metis::NOPTIONS];
527    ///
528    /// // four refinement iterations instead of the default 10.
529    /// options[metis::option::NIter::INDEX] = 4;
530    ///
531    /// Graph::new(1, 2, xadj, adjncy)?
532    ///     .set_options(&options)
533    ///     .part_recursive(&mut part)?;
534    ///
535    /// // The two vertices are placed in different parts.
536    /// assert_ne!(part[0], part[1]);
537    /// # Ok(())
538    /// # }
539    /// ```
540    pub fn set_options(mut self, options: &[Idx; NOPTIONS]) -> Graph<'a> {
541        self.options.copy_from_slice(options);
542        self
543    }
544
545    /// Sets a fine-tuning parameter for this partitioning.
546    ///
547    /// When options are to be set in batches, [`Graph::set_options`] might be a
548    /// better fit.
549    ///
550    /// See the [option] module for the list of available parameters.  Note that
551    /// not all are applicable to a given partitioning method.  Refer to the
552    /// documentation of METIS ([link]) for more info on this.
553    ///
554    /// [link]: http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/manual.pdf
555    ///
556    /// # Example
557    ///
558    /// ```rust
559    /// # fn main() -> Result<(), metis::Error> {
560    /// # use metis::Graph;
561    /// let xadj = &[0, 1, 2];
562    /// let adjncy = &[1, 0];
563    /// let mut part = [0, 0];
564    ///
565    /// Graph::new(1, 2, xadj, adjncy)?
566    ///     .set_option(metis::option::NIter(4))
567    ///     .part_recursive(&mut part)?;
568    ///
569    /// // The two vertices are placed in different parts.
570    /// assert_ne!(part[0], part[1]);
571    /// # Ok(())
572    /// # }
573    /// ```
574    pub fn set_option<O>(mut self, option: O) -> Graph<'a>
575    where
576        O: option::Opt,
577    {
578        self.options[O::INDEX] = option.value();
579        self
580    }
581
582    /// Partition the graph using multilevel recursive bisection.
583    ///
584    /// Returns the edge-cut, the total communication volume of the
585    /// partitioning solution.
586    ///
587    /// Equivalent of `METIS_PartGraphRecursive`.
588    ///
589    /// # Panics
590    ///
591    /// This function panics if the length of `part` is not the number of
592    /// vertices.
593    pub fn part_recursive(mut self, part: &mut [Idx]) -> Result<Idx> {
594        self.options[option::Numbering::INDEX] = option::Numbering::C.value();
595        let part_len = Idx::try_from(part.len()).expect("part array larger than Idx::MAX");
596        assert_eq!(
597            part_len,
598            self.xadj.len() as Idx - 1,
599            "part.len() must be equal to the number of vertices",
600        );
601
602        if self.nparts == 1 {
603            // METIS does not handle this case well.
604            part.fill(0);
605            return Ok(0);
606        }
607
608        let nvtxs = self.xadj.len() as Idx - 1;
609        let mut edgecut = mem::MaybeUninit::uninit();
610        let part = part.as_mut_ptr();
611        unsafe {
612            m::METIS_PartGraphRecursive(
613                &nvtxs as *const Idx as *mut Idx,
614                &self.ncon as *const Idx as *mut Idx,
615                slice_to_mut_ptr(self.xadj),
616                slice_to_mut_ptr(self.adjncy),
617                self.vwgt
618                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
619                self.vsize
620                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
621                self.adjwgt
622                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
623                &self.nparts as *const Idx as *mut Idx,
624                self.tpwgts
625                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
626                self.ubvec
627                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
628                slice_to_mut_ptr(&self.options),
629                edgecut.as_mut_ptr(),
630                part,
631            )
632            .wrap()?;
633            Ok(edgecut.assume_init())
634        }
635    }
636
637    /// Partition the graph using multilevel k-way partitioning.
638    ///
639    /// Returns the edge-cut, the total communication volume of the
640    /// partitioning solution.
641    ///
642    /// Equivalent of `METIS_PartGraphKway`.
643    ///
644    /// # Panics
645    ///
646    /// This function panics if the length of `part` is not the number of
647    /// vertices.
648    pub fn part_kway(self, part: &mut [Idx]) -> Result<Idx> {
649        let part_len = Idx::try_from(part.len()).expect("part array larger than Idx::MAX");
650        assert_eq!(
651            part_len,
652            self.xadj.len() as Idx - 1,
653            "part.len() must be equal to the number of vertices",
654        );
655
656        if self.nparts == 1 {
657            // METIS does not handle this case well.
658            part.fill(0);
659            return Ok(0);
660        }
661
662        let nvtxs = self.xadj.len() as Idx - 1;
663        let mut edgecut = mem::MaybeUninit::uninit();
664        let part = part.as_mut_ptr();
665        unsafe {
666            m::METIS_PartGraphKway(
667                &nvtxs as *const Idx as *mut Idx,
668                &self.ncon as *const Idx as *mut Idx,
669                slice_to_mut_ptr(self.xadj),
670                slice_to_mut_ptr(self.adjncy),
671                self.vwgt
672                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
673                self.vsize
674                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
675                self.adjwgt
676                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
677                &self.nparts as *const Idx as *mut Idx,
678                self.tpwgts
679                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
680                self.ubvec
681                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
682                slice_to_mut_ptr(&self.options),
683                edgecut.as_mut_ptr(),
684                part,
685            )
686            .wrap()?;
687            Ok(edgecut.assume_init())
688        }
689    }
690}
691
692/// Error raised when the mesh data fed to [`Mesh::new`] cannot be safely passed
693/// to METIS.
694///
695/// Mesh data must follow the format described in [`Mesh::new`].
696#[derive(Debug)]
697pub struct InvalidMeshError {
698    msg: &'static str,
699}
700
701impl fmt::Display for InvalidMeshError {
702    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
703        self.msg.fmt(f)
704    }
705}
706
707/// Error type returned by [`Mesh::new`].
708///
709/// Unlike [`Error`], this error originates from the Rust bindings.
710#[derive(Debug)]
711#[non_exhaustive]
712pub enum NewMeshError {
713    /// `nparts` must be greater than 1.
714    NoParts,
715
716    /// Mesh is too large. One of the array's length doesn't fit into [`Idx`].
717    TooLarge,
718
719    /// The input arrays are malformed and cannot be safely passed to METIS.
720    ///
721    /// Note that these bindings do not check for all the invariants. Some might
722    /// be raised during [`Mesh::part_dual`] and [`Mesh::part_nodal`] as
723    /// [`Error::Input`].
724    InvalidMesh(InvalidMeshError),
725}
726
727impl fmt::Display for NewMeshError {
728    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
729        match self {
730            Self::NoParts => write!(f, "there must be at least one part"),
731            Self::TooLarge => write!(f, "mesh is too large"),
732            Self::InvalidMesh(err) => write!(f, "invalid mesh structure: {err}"),
733        }
734    }
735}
736
737impl std::error::Error for NewMeshError {}
738
739impl NewMeshError {
740    fn msg(msg: &'static str) -> Self {
741        Self::InvalidMesh(InvalidMeshError { msg })
742    }
743}
744
745/// Returns the number of elements and the number of nodes in the mesh.
746fn check_mesh_structure(eptr: &[Idx], eind: &[Idx]) -> StdResult<(Idx, Idx), NewMeshError> {
747    let last_eptr = *eptr
748        .last()
749        .ok_or(NewMeshError::msg("element index is empty"))?;
750    let eind_len = Idx::try_from(eind.len()).map_err(|_| NewMeshError::TooLarge)?;
751    if last_eptr != eind_len {
752        return Err(NewMeshError::msg(
753            "length mismatch between element and node indices",
754        ));
755    }
756
757    let ne = Idx::try_from(eptr.len()).map_err(|_| NewMeshError::TooLarge)? - 1;
758
759    let mut prev = 0;
760    for x in eptr {
761        if prev > *x {
762            return Err(NewMeshError::msg("element index is not sorted"));
763        }
764        prev = *x;
765    }
766
767    let mut max_node = 0;
768    for a in eind {
769        if *a < 0 {
770            return Err(NewMeshError::msg(
771                "values in the node index are out of bounds",
772            ));
773        }
774        if *a > max_node {
775            max_node = *a;
776        }
777    }
778
779    Ok((ne, max_node + 1))
780}
781
782/// Builder structure to set up a mesh partition computation.
783///
784/// This structure holds the required arguments for METIS to compute a
785/// partition.  It also offers methods to easily set any optional argument.
786///
787/// # Example
788///
789/// Usage is fairly similar to [`Graph`].  Refer to its documentation for
790/// details.
791#[derive(Debug, PartialEq)]
792pub struct Mesh<'a> {
793    /// The number of nodes in the mesh.
794    nn: Idx,
795
796    /// The number of parts to partition the mesh.
797    nparts: Idx,
798
799    /// The number of nodes two elements must share for an edge to appear in the
800    /// dual graph.
801    ncommon: Idx,
802
803    eptr: &'a [Idx], // mesh representation
804    eind: &'a [Idx], // mesh repr
805
806    /// The computational weights of the elements.
807    ///
808    /// Required size: ne
809    vwgt: Option<&'a [Idx]>,
810
811    /// The communication weights of the elements.
812    ///
813    /// Required size: ne
814    vsize: Option<&'a [Idx]>,
815
816    /// The target partition weights of the elements.
817    ///
818    /// If `None` then the mesh is equally divided among the partitions.
819    ///
820    /// Required size: nparts
821    tpwgts: Option<&'a [Real]>,
822
823    /// Fine-tuning parameters.
824    options: [Idx; NOPTIONS],
825}
826
827impl<'a> Mesh<'a> {
828    /// Creates a new [`Mesh`] object to be partitioned.
829    ///
830    /// `nparts` is the number of parts wanted in the mesh partition.
831    ///
832    /// # Input format
833    ///
834    /// The length of `eptr` is `n + 1`, where `n` is the number of elements in
835    /// the mesh. The length of `eind` is the sum of the number of nodes in all
836    /// the elements of the mesh. The list of nodes belonging to the `i`th
837    /// element of the mesh are stored in consecutive locations of `eind`
838    /// starting at position `eptr[i]` up to (but not including) position
839    /// `eptr[i+1]`.
840    ///
841    /// # Errors
842    ///
843    /// The following invariants must be held, otherwise this function returns
844    /// an error:
845    ///
846    /// - `nparts` is strictly greater than zero,
847    /// - `eptr` has at least one element (its length is the one more than the
848    ///   number of mesh elements),
849    /// - `eptr` is sorted,
850    /// - elements of `eptr` are positive,
851    /// - the last element of `eptr` is the length of `eind`,
852    /// - all the arrays have a length that can be held by an [`Idx`].
853    ///
854    /// # Mutability
855    ///
856    /// [`Mesh::part_dual`] and [`Mesh::part_nodal`] may mutate the contents of
857    /// `eptr` and `eind`, but should revert all changes before returning.
858    pub fn new(nparts: Idx, eptr: &'a [Idx], eind: &'a [Idx]) -> StdResult<Mesh<'a>, NewMeshError> {
859        if nparts <= 0 {
860            return Err(NewMeshError::NoParts);
861        }
862        let (_ne, nn) = check_mesh_structure(eptr, eind)?;
863        Ok(unsafe { Mesh::new_unchecked(nn, nparts, eptr, eind) })
864    }
865
866    /// Creates a new [`Mesh`] object to be partitioned (unchecked version).
867    ///
868    /// - `nn` is the number of nodes in the mesh,
869    /// - `nparts` is the number of parts wanted in the mesh partition.
870    ///
871    /// # Input format
872    ///
873    /// See [`Mesh::new`].
874    ///
875    /// # Safety
876    ///
877    /// This function still does some checks listed in "Panics" below. However,
878    /// the caller is reponsible for upholding all invariants listed in the
879    /// "Errors" section of [`Mesh::new`]. Otherwise, the behavior of this
880    /// function is undefined.
881    ///
882    /// # Panics
883    ///
884    /// This function panics if:
885    /// - any of the arrays have a length that cannot be hold by an [`Idx`], or
886    /// - `nn` is not strictly greater than zero, or
887    /// - `nparts` is not strictly greater than zero, or
888    /// - `eptr` is empty, or
889    /// - the length of `eind` is different from the last element of `eptr`.
890    ///
891    /// # Mutability
892    ///
893    /// While nothing should be modified by the [`Mesh`] structure, METIS
894    /// doesn't specify any `const` modifier, so everything must be mutable on
895    /// Rust's side.
896    pub unsafe fn new_unchecked(
897        nn: Idx,
898        nparts: Idx,
899        eptr: &'a [Idx],
900        eind: &'a [Idx],
901    ) -> Mesh<'a> {
902        assert!(0 < nn, "nn must be strictly greater than zero");
903        assert!(0 < nparts, "nparts must be strictly greater than zero");
904        let _ = Idx::try_from(eptr.len()).expect("eptr array larger than Idx::MAX");
905        assert_ne!(eptr.len(), 0);
906        let eind_len = Idx::try_from(eind.len()).expect("eind array larger than Idx::MAX");
907        assert_eq!(eind_len, *eptr.last().unwrap());
908
909        Mesh {
910            nn,
911            nparts,
912            ncommon: 1,
913            eptr,
914            eind,
915            vwgt: None,
916            vsize: None,
917            tpwgts: None,
918            options: [-1; NOPTIONS],
919        }
920    }
921
922    /// Sets the computational weights of the elements.
923    ///
924    /// By default, all elements have the same weight.
925    ///
926    /// All elements of `vwgt` must be positive.
927    ///
928    /// # Panics
929    ///
930    /// This function panics if the length of `vwgt` is not the number of
931    /// elements.
932    pub fn set_vwgt(mut self, vwgt: &'a [Idx]) -> Mesh<'a> {
933        let vwgt_len = Idx::try_from(vwgt.len()).expect("vwgt array too large");
934        assert_eq!(vwgt_len, self.eptr.len() as Idx - 1);
935        self.vwgt = Some(vwgt);
936        self
937    }
938
939    /// Sets the communication weights of the elements.
940    ///
941    /// By default, all elements have the same communication weight.
942    ///
943    /// # Panics
944    ///
945    /// This function panics if the length of `vsize` is not the number of
946    /// elements.
947    pub fn set_vsize(mut self, vsize: &'a [Idx]) -> Mesh<'a> {
948        let vsize_len = Idx::try_from(vsize.len()).expect("vsize array too large");
949        assert_eq!(vsize_len, self.eptr.len() as Idx - 1);
950        self.vsize = Some(vsize);
951        self
952    }
953
954    /// Sets the target partition weights for each part.
955    ///
956    /// By default, the mesh is divided equally.
957    ///
958    /// The sum of the target partition weights must be 1.0.
959    ///
960    /// # Panics
961    ///
962    /// This function panics if the length of `tpwgts` is not equal to `nparts`.
963    pub fn set_tpwgts(mut self, tpwgts: &'a [Real]) -> Mesh<'a> {
964        let tpwgts_len = Idx::try_from(tpwgts.len()).expect("tpwgts array too large");
965        assert_eq!(tpwgts_len, self.nparts);
966        self.tpwgts = Some(tpwgts);
967        self
968    }
969
970    /// Sets the fine-tuning parameters for this partitioning.
971    ///
972    /// When few options are to be set, [`Mesh::set_option`] might be a
973    /// better fit.
974    ///
975    /// See the [option] module for the list of available parameters.  Note that
976    /// not all are applicable to a given partitioning method.  Refer to the
977    /// documentation of METIS ([link]) for more info on this.
978    ///
979    /// See [`Graph::set_options`] for a usage example.
980    ///
981    /// [link]: http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/manual.pdf
982    pub fn set_options(mut self, options: &[Idx; NOPTIONS]) -> Mesh<'a> {
983        self.options.copy_from_slice(options);
984        self
985    }
986
987    /// Sets a fine-tuning parameter for this partitioning.
988    ///
989    /// When options are to be set in batches, [`Mesh::set_options`] might be a
990    /// better fit.
991    ///
992    /// See the [option] module for the list of available parameters.  Note that
993    /// not all are applicable to a given partitioning method.  Refer to the
994    /// documentation of METIS ([link]) for more info on this.
995    ///
996    /// See [`Graph::set_option`] for a usage example.
997    ///
998    /// [link]: http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/manual.pdf
999    pub fn set_option<O>(mut self, option: O) -> Mesh<'a>
1000    where
1001        O: option::Opt,
1002    {
1003        self.options[O::INDEX] = option.value();
1004        self
1005    }
1006
1007    /// Partition the mesh using its dual graph.
1008    ///
1009    /// Returns the edge-cut, the total communication volume of the
1010    /// partitioning solution.
1011    ///
1012    /// Equivalent of `METIS_PartMeshDual`.
1013    ///
1014    /// # Panics
1015    ///
1016    /// This function panics if the length of `epart` is not the number of
1017    /// elements, or if `nparts`'s is not the number of nodes.
1018    pub fn part_dual(mut self, epart: &mut [Idx], npart: &mut [Idx]) -> Result<Idx> {
1019        self.options[option::Numbering::INDEX] = option::Numbering::C.value();
1020        let epart_len = Idx::try_from(epart.len()).expect("epart array larger than Idx::MAX");
1021        assert_eq!(
1022            epart_len,
1023            self.eptr.len() as Idx - 1,
1024            "epart.len() must be equal to the number of elements",
1025        );
1026        let npart_len = Idx::try_from(npart.len()).expect("npart array larger than Idx::MAX");
1027        assert_eq!(
1028            npart_len, self.nn,
1029            "npart.len() must be equal to the number of nodes",
1030        );
1031
1032        if self.nparts == 1 {
1033            // METIS does not handle this case well.
1034            epart.fill(0);
1035            npart.fill(0);
1036            return Ok(0);
1037        }
1038
1039        let ne = self.eptr.len() as Idx - 1;
1040        let mut edgecut = mem::MaybeUninit::uninit();
1041        unsafe {
1042            m::METIS_PartMeshDual(
1043                &ne as *const Idx as *mut Idx,
1044                &self.nn as *const Idx as *mut Idx,
1045                slice_to_mut_ptr(self.eptr),
1046                slice_to_mut_ptr(self.eind),
1047                self.vwgt
1048                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
1049                self.vsize
1050                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
1051                &self.ncommon as *const Idx as *mut Idx,
1052                &self.nparts as *const Idx as *mut Idx,
1053                self.tpwgts
1054                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
1055                slice_to_mut_ptr(&self.options),
1056                edgecut.as_mut_ptr(),
1057                epart.as_mut_ptr(),
1058                npart.as_mut_ptr(),
1059            )
1060            .wrap()?;
1061            Ok(edgecut.assume_init())
1062        }
1063    }
1064
1065    /// Partition the mesh using its nodal graph.
1066    ///
1067    /// Returns the edge-cut, the total communication volume of the
1068    /// partitioning solution.
1069    ///
1070    /// Previous settings of `ncommon` are not used by this function.
1071    ///
1072    /// Equivalent of `METIS_PartMeshNodal`.
1073    ///
1074    /// # Panics
1075    ///
1076    /// This function panics if the length of `epart` is not the number of
1077    /// elements, or if `nparts`'s is not the number of nodes.
1078    pub fn part_nodal(mut self, epart: &mut [Idx], npart: &mut [Idx]) -> Result<Idx> {
1079        self.options[option::Numbering::INDEX] = option::Numbering::C.value();
1080        let epart_len = Idx::try_from(epart.len()).expect("epart array larger than Idx::MAX");
1081        assert_eq!(
1082            epart_len,
1083            self.eptr.len() as Idx - 1,
1084            "epart.len() must be equal to the number of elements",
1085        );
1086        let npart_len = Idx::try_from(npart.len()).expect("npart array larger than Idx::MAX");
1087        assert_eq!(
1088            npart_len, self.nn,
1089            "npart.len() must be equal to the number of nodes",
1090        );
1091
1092        if self.nparts == 1 {
1093            // METIS does not handle this case well.
1094            epart.fill(0);
1095            npart.fill(0);
1096            return Ok(0);
1097        }
1098
1099        let ne = self.eptr.len() as Idx - 1;
1100        let mut edgecut = mem::MaybeUninit::uninit();
1101        unsafe {
1102            m::METIS_PartMeshNodal(
1103                &ne as *const Idx as *mut Idx,
1104                &self.nn as *const Idx as *mut Idx,
1105                slice_to_mut_ptr(self.eptr),
1106                slice_to_mut_ptr(self.eind),
1107                self.vwgt
1108                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
1109                self.vsize
1110                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
1111                &self.nparts as *const Idx as *mut Idx,
1112                self.tpwgts
1113                    .map_or_else(ptr::null_mut, |s| slice_to_mut_ptr(s)),
1114                slice_to_mut_ptr(&self.options),
1115                edgecut.as_mut_ptr(),
1116                epart.as_mut_ptr(),
1117                npart.as_mut_ptr(),
1118            )
1119            .wrap()?;
1120            Ok(edgecut.assume_init())
1121        }
1122    }
1123}
1124
1125/// The dual of a mesh.
1126///
1127/// Result of [`mesh_to_dual`].
1128#[derive(Debug, PartialEq, Eq)]
1129pub struct Dual {
1130    xadj: &'static mut [Idx],
1131    adjncy: &'static mut [Idx],
1132}
1133
1134impl Dual {
1135    /// The adjacency index array.
1136    pub fn xadj(&self) -> &[Idx] {
1137        self.xadj
1138    }
1139
1140    /// The adjacency array.
1141    pub fn adjncy(&self) -> &[Idx] {
1142        self.adjncy
1143    }
1144
1145    /// The adjacency index array, and the adjacency array as mutable slices.
1146    pub fn as_mut(&mut self) -> (&mut [Idx], &mut [Idx]) {
1147        (self.xadj, self.adjncy)
1148    }
1149}
1150
1151impl Drop for Dual {
1152    fn drop(&mut self) {
1153        unsafe {
1154            m::METIS_Free(self.xadj.as_mut_ptr() as *mut os::raw::c_void);
1155            m::METIS_Free(self.adjncy.as_mut_ptr() as *mut os::raw::c_void);
1156        }
1157    }
1158}
1159
1160/// Generate the dual graph of a mesh.
1161///
1162/// # Errors
1163///
1164/// This function returns an error if `eptr` and `eind` don't follow the mesh
1165/// format given in [`Mesh::new`].
1166pub fn mesh_to_dual(eptr: &[Idx], eind: &[Idx], ncommon: Idx) -> Result<Dual> {
1167    let (ne, nn) = check_mesh_structure(eptr, eind)?;
1168    let mut xadj = mem::MaybeUninit::uninit();
1169    let mut adjncy = mem::MaybeUninit::uninit();
1170    let numbering_flag = 0;
1171
1172    // SAFETY: METIS_MeshToDual allocates the xadj and adjncy arrays.
1173    // SAFETY: hopefully those arrays are of correct length.
1174    unsafe {
1175        m::METIS_MeshToDual(
1176            &ne as *const Idx as *mut Idx,
1177            &nn as *const Idx as *mut Idx,
1178            slice_to_mut_ptr(eptr),
1179            slice_to_mut_ptr(eind),
1180            &ncommon as *const Idx as *mut Idx,
1181            &numbering_flag as *const Idx as *mut Idx,
1182            xadj.as_mut_ptr(),
1183            adjncy.as_mut_ptr(),
1184        )
1185        .wrap()?;
1186        let xadj = xadj.assume_init();
1187        let xadj = slice::from_raw_parts_mut(xadj, eptr.len());
1188        let adjncy = adjncy.assume_init();
1189        let adjncy = slice::from_raw_parts_mut(adjncy, xadj[xadj.len() - 1] as usize);
1190        Ok(Dual { xadj, adjncy })
1191    }
1192}