signal: 11, SIGSEGV: invalid memory reference

[NickRoz1]

Hi,

I have this function:

pub fn calc_partial_transpose(
    snps: &mut Vec<f64>,
    partial_matrix: &mut Vec<f64>,
    ids_num: usize,
  ) -> () {
    extern crate cblas;
    extern crate openblas_src;

    let n = ids_num;
    let k = snps.len() / n;
    println!(
      "N = {:?}\nK = {:?}\nRES_MATR_SIZE = {}\nSNPS_SIZE = {}",
      n,
      k,
      partial_matrix.len(),
      snps.len()
    );
    // Only the upper triangular part of matrix will be calculated.
    unsafe {
      cblas::dsyrk(
        cblas::Layout::RowMajor,
        cblas::Part::Upper,
        cblas::Transpose::Ordinary,
        n as i32,
        k as i32,
        1.0,
        &snps[..],
        k as i32,
        0.0,
        &mut partial_matrix[..],
        n as i32,
      );
    }
  }

On run it prints:

N = 17
K = 1000
RES_MATR_SIZE = 289
SNPS_SIZE = 17000

And fails with (signal: 11, SIGSEGV: invalid memory reference).

It is called in parallel from different threads created by this function:

pub fn calc_transpose(&mut self, batch_size: usize) -> std::io::Result<Vec<f64>> {
  if batch_size < 1 {
    panic!("Batch size can't be less than 1.");
  }
  let ids_num = self.width_of_matrix;

  let common_matrix: Arc<Mutex<Vec<f64>>> =
    Arc::new(Mutex::new(vec![0.0; ids_num * ids_num]));

  let buf_size = ids_num * batch_size;
  // For each physical thread a buffer will be created.
  let buf_num = num_cpus::get();
  use std::sync::{Arc, Mutex};

  let mut read_bufs = Vec::<Arc<Mutex<Vec<f64>>>>::new();
  let mut part_res_bufs = Vec::<Arc<Mutex<Vec<f64>>>>::new();
  for _ in 0..buf_num {
    read_bufs.push(Arc::new(Mutex::new(vec![0.0; buf_size])));
    part_res_bufs.push(Arc::new(Mutex::new(vec![0.0; ids_num * ids_num])));
  }

  use std::sync::mpsc::channel;
  use std::thread;
  let (processor, buffer_filler) = channel::<usize>();
  // Fill the concurrent queue with the buffers numbers, so the buffer
  // filler can start parsing into them.
  for buf_idx in 0..buf_num {
    processor.send(buf_idx).unwrap();
  }
  let mut threads = Vec::<thread::JoinHandle<()>>::new();
  let file_reader = self.file_reader.get_mut();
  let mut line_iter = BufReader::new(file_reader).lines();
  let mut total_snps_read: usize = 0;
  loop {
    let freed_buffer_idx = buffer_filler.recv().unwrap();
    let read_buf = read_bufs[freed_buffer_idx].clone();
    let part_buf = part_res_bufs[freed_buffer_idx].clone();

    let read_line_amount = match Self::fill_buffer(
      &mut *read_buf.lock().unwrap(),
      &mut line_iter,
      self.markers.len(),
      self.hab_mapper.clone(),
    )? {
      0 => break,
      n => {
        total_snps_read += n;
        n
      }
    };

    // Reached EOF (there is not enough records to form a full batch.)
    if read_line_amount < batch_size {
      let buf = &mut *read_buf.lock().unwrap();
      buf.resize(read_line_amount * ids_num, 0.0);
    }
    // CBlas demands first dimension of snps matrix (ids_num) to be no less
    // than a row count, thus it is increased in cases when the read batch
    // happens to be smaller than it (ids_num is 20, but only 19 records
    // were left to read from file(last batch)).
    if read_line_amount < ids_num {
      let buf = &mut *read_buf.lock().unwrap();
      buf.resize(ids_num * ids_num, 0.0);
    }
    {
      let (read_buf_arc, part_buf_arc, res_matrix_arc, proc_sender, buf_idx) = (
        read_buf.clone(),
        part_buf.clone(),
        common_matrix.clone(),
        processor.clone(),
        freed_buffer_idx.clone(),
      );
      threads.push(std::thread::spawn(move || {
        let mut threads_read_buf = read_buf_arc.lock().unwrap();
        let mut threads_part_buf = part_buf_arc.lock().unwrap();
        calc_partial_transpose(&mut threads_read_buf, &mut threads_part_buf, ids_num);
        let mut res_matrix = res_matrix_arc.lock().unwrap();
        for (buf_elem, common_matrix_elem) in
          threads_part_buf.iter_mut().zip(res_matrix.iter_mut())
        {
          *common_matrix_elem += *buf_elem;
        }
        proc_sender.send(buf_idx).unwrap();
      }));
    }
  }

  threads.into_iter().for_each(|thread| {
    thread
      .join()
      .expect("The thread creating or execution failed !")
  });

  self.file_reader.seek(SeekFrom::Start(self.snp_pos_start))?;

  let mut res = Arc::try_unwrap(common_matrix)
    .expect("Arc uwrapping failed.")
    .into_inner()
    .expect("Mutex uwrapping failed.");

  // Mirror matrix, since only the upper part was calculated.
  for i in 0..ids_num {
    let row_length = ids_num;
    for j in 0..i {
      res[i * row_length + j] = res[j * row_length + i];
    }
  }
  Ok(res)
}

[NickRoz1]

The SIGSEGV happens once amount of elements in input data grows from 2822 to 2839 (+ 1 row (17 elements per row)).

I.e.

Fine:

N = 17
K = 166
RES_MATR_SIZE = 289
SNPS_SIZE = 2822

Crash:

N = 17
K = 167
RES_MATR_SIZE = 289
SNPS_SIZE = 2839

[NickRoz1]

Stack trace:

Thread 13 "check" received signal SIGSEGV, Segmentation fault.
(gdb) bt
#0  0x00007ffff5ebdee5 in dgemm_otcopy_PRESCOTT () from /usr/lib/x86_64-linux-gnu/libopenblas.so.0
#1  0x00007fffd2020100 in ?? ()
#2  0x0000000000000011 in ?? ()
#3  0x0000000000000011 in ?? ()
#4  0x00007ffff7fc9670 in ?? ()
#5  0x0000000000000053 in ?? ()
#6  0x0000000000000011 in ?? ()
#7  0x00007ffff5d07eb5 in dsyrk_LN () from /usr/lib/x86_64-linux-gnu/libopenblas.so.0
#8  0x00007ffff5d16ccc in dsyrk_thread_LN () from /usr/lib/x86_64-linux-gnu/libopenblas.so.0
#9  0x00007ffff5c248fb in cblas_dsyrk () from /usr/lib/x86_64-linux-gnu/libopenblas.so.0
#10 0x0000555555578862 in cblas::dsyrk (layout=cblas::Layout::RowMajor, uplo=cblas::Part::Upper, 
    trans=cblas::Transpose::Ordinary, n=17, k=167, alpha=1, a=..., lda=167, beta=0, c=..., ldc=17)

[IvanUkhov]

Thanks for the issue! My spontaneous reaction is that the function is misused. Some of the input parameter probably do not make sense together. Another idea is to try to run the code without multithreading to exclude that factor.

[NickRoz1]

Hi,

I tried to disable multithreading in BLAS:

$ ldconfig -p | grep openblas
	libopenblas.so.0 (libc6,x86-64) => /usr/lib/x86_64-linux-gnu/libopenblas.so.0
	libopenblas.so (libc6,x86-64) => /usr/lib/x86_64-linux-gnu/libopenblas.so

https://github.com/xianyi/OpenBLAS/blob/ce3651516f12079f3ca2418aa85b9ad571c3a391/driver/others/openblas_set_num_threads.c

      extern crate libc;
      use libc::{c_int, c_void};

      extern "C" {
        fn openblas_set_num_threads(num_threads: c_int) -> c_void;
      }

      unsafe {
        openblas_set_num_threads(1);
      }

Not sure that I'm doing it correctly, though.

Exporting variable (export OPENBLAS_NUM_THREADS=1) as suggested here doesn't seem to help also.

Still:
(signal: 11, SIGSEGV: invalid memory reference).

[NickRoz1]

Single threaded code produces same error (signal: 11, SIGSEGV: invalid memory reference):

    pub fn calc_transpose(&mut self, batch_size: usize) -> std::io::Result<Vec<f64>> {
      if batch_size < 1 {
        panic!("Batch size can't be less than 1.");
      }
      let ids_num = self.markers.len();

      let mut common_matrix: Vec<f64> = vec![0.0; ids_num * ids_num];

      // This amount of snps will be parsed and processed on each iteration.
      let buf_size = ids_num * batch_size;

      let mut read_buf: Vec<f64> = vec![0.0; buf_size];

      let file_reader = self.file_reader.get_mut();
      let mut line_iter = BufReader::new(file_reader).lines();
      let mut total_snps_read: usize = 0;
      loop {
        let read_line_amount = match Self::fill_buffer(
          &mut read_buf,
          &mut line_iter,
          self.markers.len(),
          self.hab_mapper.clone(),
        )? {
          0 => break,
          n => {
            total_snps_read += n;
            n
          }
        };
        // CBlas demands first dimension of snps matrix (ids_num) to be no less
        // than a row count, thus it is increased in cases when the read batch
        // happens to be smaller than it (ids_num is 20, but only 19 records
        // were left to read from file(last batch)).
        if read_line_amount < ids_num {
          let buf = &mut read_buf;
          buf.resize(ids_num * ids_num, 0.0);
        }
        calc_partial_transpose(&mut read_buf, &mut common_kinship_matrix, ids_num);
      }

      self.file_reader.seek(SeekFrom::Start(self.snp_pos_start))?;

      let mut res = common_matrix;

      // Mirror matrix
      for i in 0..ids_num {
        let row_length = ids_num;
        for j in 0..i + 1 {
          res[j * row_length + i] /= total_snps_read as f64;
          res[i * row_length + j] = res[j * row_length + i];
        }
      }
      Ok(res)
    }

  pub fn calc_partial_transpose(
    snps: &mut Vec<f64>,
    partial_matrix: &mut Vec<f64>,
    ids_num: usize,
  ) -> () {
    extern crate cblas;
    extern crate openblas_src;

    let n = ids_num;
    let k = snps.len() / n;

    // Only the upper triangular part of matrix will be calculated.
    unsafe {
      cblas::dsyrk(
        cblas::Layout::RowMajor,
        cblas::Part::Upper,
        cblas::Transpose::Ordinary,
        n as i32,
        k as i32,
        1.0,
        &snps[..],
        k as i32,
        1.0,
        &mut partial_matrix[..],
        n as i32,
      );
    }
  }

[NickRoz1]

By the way, as a temporary solution I manually reimplemented Fortran dsyrk subroutine and it works fine in multithreaded environment and gives an expected result (hence the input parameters must be right, I guess?).

  pub fn calc_partial_transpose(
    snps: &mut Vec<f64>,
    partial_matrix: &mut Vec<f64>,
    ids_num: usize,
  ) -> () {
    let n = ids_num;
    let k = snps.len() / n;
    // Algorithm from BLAS dsyrk:
    // http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gae0ba56279ae3fa27c75fefbc4cc73ddf.html#gae0ba56279ae3fa27c75fefbc4cc73ddf
    //
    // The BLAS Fortran stores array in a column-major format, but the this data
    // stored in a row-major format, so this algorithm corresponds to the branch
    // for non transposed, lower triangular part version, however in fact it
    // performs transposed, upper triangular part multiplication (G.T*G).
    //
    // This algorithm branch (Lower, Non transposed) chosen based on CBLAS
    // http://www.netlib.org/blas/blast-forum/cblas.tgz code for dsyrk
    // (cblas_dsyrk.c), which transforms options (Upper, Transposed) to these
    // arguments when called for row-major matrixes.
    //
    // When the matrix stored in row-major way read in column major way,
    // obtained data is a transpose of this matrix:
    // https://en.wikipedia.org/wiki/Row-_and_column-major_order#Transposition
    //
    // Since this is an exact copy of Fortran code, and Fortran utilizes double
    // index (i,j) to operate over single dimension array (which represents 2D
    // array), the code below performs index flattening for column-major
    // storages exactly how Fortran does. Normally, to flatten index in
    // row-major languages we will multiply row index i by row width and add
    // column index j, here, since this is a direct copy of Fortran code which
    // is a colum-major language, we flatten it as column index j *
    // column height + row index i.
    for j in 0..n {
      for l in 0..k {
        for i in j..n {
          partial_matrix[j * ids_num + i] += snps[l * ids_num + j] * snps[l * ids_num + i];
        }
      }
    }
  }

[IvanUkhov]

Thanks for investigating! You might have implemented the function according to your understanding of how the original one works, and hence it is correct. The original one might be working not as you expect.

If you are positive the inputs make sense, the next step would be to try to eliminate this crate (define the extern C function and do the linking manually). Then one can even try eliminate Rust all together (implement in C). It will then tell where the problem lies.