matrix_base.hpp

Go to the documentation of this file.

/*
Copyright 2020 Peter Rakyta, Ph.D.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

*/

#ifndef matrix_BASE_H
#define matrix_BASE_H

#include "QGDTypes.h"
#include <atomic>
#include <cstring>
#include <iostream>
#include <tbb/scalable_allocator.h>



template<typename scalar>
class matrix_base {

public:
  int rows;
  int cols;
  int stride;
  scalar* data;

protected:

  bool conjugated;
  bool transposed;
  bool owner;
  std::atomic<int64_t>* references;



public:

matrix_base() {

  // The number of rows
  rows = 0;
  // The number of columns
  cols = 0;
  // The column stride of the matrix
  stride = 0;
  // pointer to the stored data
  data = NULL;
  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  conjugated = false;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  transposed = false;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  owner = false;
  references = new std::atomic<int64_t>(1);
}


matrix_base( scalar* data_in, int rows_in, int cols_in) {

  // The number of rows
  rows = rows_in;
  // The number of columns
  cols = cols_in;
  // The column stride of the matrix
  stride = cols_in;
  // pointer to the stored data
  data = data_in;
  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  conjugated = false;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  transposed = false;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  owner = false;
  references = new std::atomic<int64_t>(1);
}



matrix_base( scalar* data_in, int rows_in, int cols_in, int stride_in) {

  // The number of rows
  rows = rows_in;
  // The number of columns
  cols = cols_in;
  // The column stride of the matrix
  stride = stride_in;
  // pointer to the stored data
  data = data_in;
  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  conjugated = false;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  transposed = false;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  owner = false;
  references = new std::atomic<int64_t>(1);
}



matrix_base( int rows_in, int cols_in) {

  // The number of rows
  rows = rows_in;
  // The number of columns
  cols = cols_in;
  // The column stride of the matrix
  stride = cols_in;
  // pointer to the stored data
  data = (scalar*)scalable_aligned_malloc( rows*cols*sizeof(scalar), CACHELINE);
  assert(data);
  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  conjugated = false;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  transposed = false;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  owner = true;
  references = new std::atomic<int64_t>(1);


}


matrix_base( int rows_in, int cols_in, int stride_in) {

  // The number of rows
  rows = rows_in;
  // The number of columns
  cols = cols_in;
  // The column stride of the matrix
  stride = stride_in;
  // pointer to the stored data
  data = (scalar*)scalable_aligned_malloc( rows*stride*sizeof(scalar), CACHELINE);
#ifdef DEBUG
  if (rows > 0 && cols>0) assert(data);
#endif
  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  conjugated = false;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  transposed = false;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  owner = true;
  references = new std::atomic<int64_t>(1);


}


matrix_base(const matrix_base<scalar> &in) {

    data = in.data;
    rows = in.rows;
    cols = in.cols;
    stride = in.stride;
    transposed = in.transposed;
    conjugated = in.conjugated;
    owner = in.owner;

    references = in.references;

    if (references != NULL) {
      references->fetch_add(1, std::memory_order_relaxed);
    }



}



~matrix_base() {
  release_data();
}

bool is_conjugated() const {
  return conjugated;
}

void conjugate() {

  conjugated = !conjugated;

}


bool is_transposed() const {

  return transposed;

}


void transpose()  {

  transposed = !transposed;

}




scalar* get_data() const {

  return data;

}


void replace_data( scalar* data_in, bool owner_in) {

    release_data();
    data = data_in;
    owner = owner_in;

    references = new std::atomic<int64_t>(1);

}


void release_data() {

    if (references == NULL) return;

    std::atomic<int64_t>* references_loc = references;
    scalar* data_loc = data;
    const bool owner_loc = owner;

    data = NULL;
    references = NULL;

    if (references_loc->fetch_sub(1, std::memory_order_acq_rel) == 1) {
      if (owner_loc) {
        scalable_aligned_free(data_loc);
      }
      delete references_loc;
    }

}



void set_owner( bool owner_in)  {

    owner=owner_in;

}

void operator= (const matrix_base& mtx ) {

  // releasing the containing data
  release_data();

  // The number of rows
  rows = mtx.rows;
  // The number of columns
  cols = mtx.cols;
  // The column stride of the matrix
  stride = mtx.stride;
  // pointer to the stored data
  data = mtx.data;
  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  conjugated = mtx.conjugated;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  transposed = mtx.transposed;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  owner = mtx.owner;

  references = mtx.references;

    if (references != NULL) {
      references->fetch_add(1, std::memory_order_relaxed);
  }

}


scalar& operator[](int idx) const {

#ifdef DEBUG
    if ( idx >= rows*stride || idx < 0) {
        std::cout << "Accessing element out of bonds. Exiting" << std::endl;
        exit(-1);
    }
#endif

    return data[idx];
}

bool operator<(const matrix_base<scalar>& other) const {
    if (rows != other.rows) return rows < other.rows;
    if (cols != other.cols) return cols < other.cols;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            if (data[i * stride + j] != other.data[i * stride + j])
                return data[i * stride + j] < other.data[i * stride + j];
        }
    }
    return false; // they are equal
}

matrix_base<scalar> copy() const {

  matrix_base<scalar> ret = matrix_base<scalar>(rows, cols, stride);

  // logical variable indicating whether the matrix needs to be conjugated in CBLAS operations
  ret.conjugated = conjugated;
  // logical variable indicating whether the matrix needs to be transposed in CBLAS operations
  ret.transposed = transposed;
  // logical value indicating whether the class instance is the owner of the stored data or not. (If true, the data array is released in the destructor)
  ret.owner = true;

  memcpy( ret.data, data, rows*stride*sizeof(scalar));

  return ret;

}

void copy_to(matrix_base<scalar>& target) const {

  if (target.rows != rows || target.cols != cols || target.stride != stride) {
    target = matrix_base<scalar>(rows, cols, stride);
  }

  if (data != NULL && target.data != data) {
    memcpy(target.data, data, rows*stride*sizeof(scalar));
  }

  target.conjugated = conjugated;
  target.transposed = transposed;

}

void ensure_aligned() {
    if (((uintptr_t)(void*)data & (CACHELINE-1)) == 0) return; //CACHELINE must be power of 2, 16 sufficient, 64 is okay though
    scalar* newdata = (scalar*)scalable_aligned_malloc( rows*stride*sizeof(scalar), CACHELINE);
    memcpy( newdata, data, rows*stride*sizeof(scalar));
    replace_data(newdata, true);
}



int size() const {

  return rows*cols;

}


void print_matrix() const {
    std::cout << std::endl << "The stored matrix:" << std::endl;
    for ( int row_idx=0; row_idx < rows; row_idx++ ) {
        for ( int col_idx=0; col_idx < cols; col_idx++ ) {
            int element_idx = row_idx*stride + col_idx;
              std::cout << " " << data[element_idx];
        }
        std::cout << std::endl;
    }
    std::cout << std::endl << std::endl << std::endl;

}







}; //matrix_base





#endif