Extract rows/columns of matrices in C++ vs MATLAV

JamieAl (236)
In MATLAB you can extract the first and last rows of a matrix by
1
2
 
  A(1,:);
  A(N,:);

Is there a similar way in C++ where we can do this?
For example,
1
2
3
4
5
6
7
 
   for (int i = 0; i < ny+1; i++){
		for (int j = 0; j < ny+1; j++){
            A[j + ny*i] = cos(acos(ygl[j]) *(i)); 
            cout << A[j + ny*i] << endl ;
		
		}
	}

How can I extract the first and last rows here and say replace them with all ones instead?
deleted account xyzzy (5768)
If you are using a 2D C++ container, such as std::vector, first and last row replacement can be done using operator[] indexing or iterators. A quick and dirty slap-dash bit of code to show how it might be done:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
 
#include <iostream>
#include <vector>

template <typename T>
std::ostream& operator<<(std::ostream&, const std::vector<std::vector<T>>&);

template <typename T>
std::ostream& operator<<(std::ostream&, const std::vector<T>&);

int main()
{
   
   // create a 2 dimensional int vector with known dimensions
   std::vector<std::vector<int>> vec { { 1, 2, 3, 4 },
                                       { 5, 6, 7, 8 },
                                       { 9, 10, 11, 12 }};

   std::cout << "Let's verify the sizes of the 2D vector....\n";
   std::cout << "The vector has " << vec.size() << " rows.\n";
   std::cout << "The vector has " << vec[0].size() << " columns.\n\n";

   // display the 2D vector
   std::cout << vec << '\n';

   // row element replacement using operator[] indexing
   for (size_t i { }; i < vec[0].size(); ++i)
   {
      size_t fst_row { };
      size_t lst_row { vec.size() - 1 };

      vec[fst_row][i] = 0;
      vec[lst_row][i] = 0;
   }

   std::cout << vec << '\n';

   // I'm sure there is a less messy way to do this.......

   // row element replacement using iterators
   for (auto itr { vec[0].begin() }; itr != vec[0].end(); ++itr)
   {
      *itr = 20;
   }

   for (auto itr { vec[vec.size() - 1].begin() }; itr != vec[vec.size() - 1].end(); ++itr)
   {
      *itr = 20;
   }

   std::cout << vec;
}

template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<std::vector<T>>& v)
{
   for (auto const& x : v)
   {
      os << x << '\n';
   }

   return os;
}


template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
{
   for (auto const& x : v)
   {
      os << x << ' ';
   }

   return os;
}
Let's verify the sizes of the 2D vector....
The vector has 3 rows.
The vector has 4 columns.

1 2 3 4
5 6 7 8
9 10 11 12

0 0 0 0
5 6 7 8
0 0 0 0

20 20 20 20
5 6 7 8
20 20 20 20
jonnin (11497)
depending on what you are doing, though, 2d vectors can cause memory access aggravations and other woes. For very small problems, it does not matter too much.

valarray's slice will get you any row, column, or stuff like the diagonals even if its a 1-d being used for 2d as you are already doing.
a derpy example:

1
2
3
4
5
6
7
8
9
10
11
12
 

int main()
{
 valarray<int> a{1,2,3,4,5,6,7,8,9,10,11,12};
 valarray<int> col = a[slice(0,4,3)];
 for(auto i:col)
 {cout << i << " ";}
 valarray<int> lastrow = a[slice(8,4,1)];
 cout << endl;
 for(auto i:lastrow)
 {cout << i << " ";}
}

this assumes a 3 rows x 4 cols entity. It takes a little getting used to getting the correct indices for the slices, but once you do, the amount of work you can do with near 0 effort is astonishing.
Last edited on
JamieAl (236)
I guess I forgot to mention, I am not using std:vector for this array, so I am initializing A as,
1
2
3
 
        double *A; 
	A= (double*) fftw_malloc((nx*ny)*sizeof(double));
	memset(A, 42, (nx*ny)* sizeof(double));


Which is different, I see what @George P has done and it makes sense but it feels a bit complex tbh.
JamieAl (236)
@jonnin I am using really large arrays that's why I am choosing to allocate memory, but can I initialize valarray to be arbitrary size of nx*ny?
jonnin (11497)
yes, valarray should be fine for a large matrix: it does the dynamic memory (new or malloc or whatever it uses) for you, and cleans the memory when it goes out of scope, all for you (vector does this too).

you can make a 2d valarray**, same as vector, but again, there are reasons to keep it 1d.
an arbitrary one?

valarray<double> mat;
mat.resize(nx*ny); //ready to go! or resize(nx*ny,42.0) if you have a default value for all cells.

vector, valarray, and std::array are all objects /classes that replace a basic C array (or pointer equivalent).

you should look at valarray. it can do block work, like adding up a row or column for you, or take the sin() of all the elements, and many things of that nature that require your own loops for other containers. The code is greatly simplified, and they were designed for performance.

**of sorts. you may have to make a vector of valarrays or something goofy if you really want to fool with that approach.
Last edited on
deleted account xyzzy (5768)
I don't routinely deal with large gobs of data, manipulating the data as a whole, so std::valarray is not something that pops up as a "first use" option.

https://en.cppreference.com/w/cpp/numeric/valarray/slice shows a Matrix class using std::valarray and how it can be manipulated.
mbozzi (3944)
It's probably wrong to use plain valarray<double> or vector<double> here.

JamieAL is using the library FFTW3, whose interface requires memory allocated with fftw_malloc. Whereas std::vector<double> allocates memory with std::allocator. And std::valarray<double> does not allow the user to control its memory allocation at all.

It's critical to call the right function, because fftw_malloc might do more than simply forward on to malloc. For example it could return an over-aligned buffer, or modify some library-internal data structures.

Here is an example of how to use std::vector with a "special" allocation function like fftw_malloc:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
 
#include <memory_resource>
#include <stdexcept>
#include <cstdint>

#include <stdlib.h>
#include <unistd.h>

// Just for example:
void* fftw_malloc(std::size_t n) 
{ 
  static long const page_size = sysconf(_SC_PAGESIZE);
  if (page_size <= -1) 
    return nullptr;
    
  void* result;
  if (::posix_memalign(&result, page_size, n) != 0)
    return nullptr;
    
  return result;
}

void fftw_free(void* p) { return ::free(p); } 

struct bad_fftw_malloc: std::bad_alloc 
{ char const* what() { return "bad fftw_malloc"; } };  

struct fftw_resource: std::pmr::memory_resource
{
  void* do_allocate(std::size_t bytes, std::size_t const a) override
  {
    void* const ptr_orig = fftw_malloc(bytes);
    if (ptr_orig == nullptr) 
      throw bad_fftw_malloc{};
    
    // Ensure the alignment requirement is already satisfied by ptr
    void* ptr_copy = ptr_orig;
    if (ptr_orig != std::align(a, bytes, ptr_copy, bytes)) 
      throw bad_fftw_malloc{}; 
    
    return ptr_orig;
  }
  
  void do_deallocate(void* p, std::size_t, std::size_t) override 
  { 
    if (p) fftw_free(p); 
  } 
  
  bool do_is_equal(std::pmr::memory_resource const&r ) const noexcept override
  { 
    return this == &r; 
  }
};

std::pmr::memory_resource* fftw_malloc_resource() 
{ static fftw_resource res; return &res; } 

#include <iostream>
#include <vector>

int main()
{  
  // Now `my_vector` is suitable for use with fftw3
  std::pmr::vector<double> my_vector(fftw_malloc_resource());  
  
  int constexpr n = 10;
  for (int i = 0; i < n; ++i) my_vector.push_back(i + 0.5);
  for (int i = 0; i < n; ++i) std::cout << my_vector[i] << '\n'; 
}
Last edited on
JamieAl (236)
Okay this sounds way too messy for just a simple one line of code in MATLAB. I am pretty confused and shocked how the alternative in C++ is waaay harder. Mm I tried some of your suggestions and I have something primitive-ish:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 

    valarray<double> VGL1;
    VGL1.resize((ny+1)*(ny+1),0.0); 

    valarray<double> dVGL1;
    dVGL1.resize((ny+1)*(ny+1),0.0);

    valarray<double> firstrow;
    firstrow.resize((ny+1)*(ny+1),0.0);
    //

     for (int i = 0; i < ny+1; i++){
		for (int j = 0; j < ny+1; j++){
            VGL1[j + ny*i] = cos(acos(ygl[j]) *(i));
            dVGL1[j + ny*i] = dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j + ny*i];
            dVGL1[0] = pow(-1,(i))*(pow(i,2)); //first row. This is correct
            dVGL1[ny+1] = 1. * pow(i,2); //last row. This is correct
            //firstrow[j + ny*i] = dVGL1[slice()]; // this won't work
            
        }
    }

I guess I shouldn't use valarray here for the reasons @mbozzi mentioned, but as a simple test now I will leave it.

So I am able to create each of the first and last row separately but not sure how to put it back in the array. I tried using slice in the loop but it wasn't really working (probably my syntax is off). I know for a standard 2D matrix in C++ I can do something like,
1
2
 
for (int i=0; i<3; i++)
    Matrix[0][i] = 4;

to replace my first row elements with 4 for example, but I am struggling to do the same with my code. I just feel like there is a straight forward solution that I am not seeing.
mbozzi (3944)
I am pretty confused and shocked how the alternative in C++ is waaay harder.

It takes a lot of effort to go from bits and bytes to an interface that acts like a matrix. The authors of MATLAB have already done the job. But there's nothing like that in the C++ standard library (valarray is the best it's got).

I guess I shouldn't use valarray here for the reasons @mbozzi mentioned, but as a simple test now I will leave it.

It's fine, but perhaps only on your computer on Friday. Or, it might be totally okay - who knows?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
 
#include <valarray>
#include <numeric>
#include <iostream>
#include <iomanip>

// Print `a` as a matrix of dimensions `width` by `height` to standard output
// Assuming row-major storage order.
void print_as_matrix(std::valarray<double> a, int width, int height)
{
  for (int y = 0; y < height; ++y)
  {
     std::cout << "[";
     for (int x = 0; x < width; ++x)
     {
       std::cout << std::fixed << a[y * width + x] << (x < width - 1? ", ": ""); 
     }
     std::cout << "]\n";
  }
}

std::slice column_slice(int i, int width, int height) { return{i, height, width}; }
std::slice row_slice(int i, int width, int) { return{width * i, width, 1}; }

int main()
{
  int constexpr width = 8;
  int constexpr height = 4;
  
  // allocate a width-by-height matrix in a single valarray
  // arbitrarily fill it with data
  std::valarray<double> a(width * height);      
  std::iota(std::begin(a), std::end(a), 1.0); 
  std::sin(a /= (width / std::atan(1) * 4.0));
  
  print_as_matrix(a, width, height);
  std::cout << '\n';
  
  // Change column 3 to all zeros
  a[column_slice(3, width, height)] = 0.0;
  // Change the last row to all ones
  a[row_slice(height - 1, width, height)] = 1.0;
  
  print_as_matrix(a, width, height);
  std::cout << '\n';
}


The a[column_slice(3, width, height)] has type std::slice_array
Assigning to it modifies the subset of a described by the slice.
https://en.cppreference.com/w/cpp/numeric/valarray/slice_array
Last edited on
jonnin (11497)
is this a good time to mention that there are many matrix capable libraries, one of the best being eigen? https://eigen.tuxfamily.org/index.php?title=Main_Page

or, to do you one better, it may cost like 1/2 a million bucks by now, but matlab has, had, something a tool to export c++ (that depended upon a matlab library file). It wasn't the fastest, but it works, and it saves translation across 2 very, very different languages.
Last edited on
deleted account xyzzy (5768)
Thanks for the link, jonnin. :)

I notice using the Eigen library is supported on most compilers, and doesn't require another 3rd party library to use.

I also notice it is available via vcpkg; to make downloading, compiling and integrating into the compiler's header/lib pathing much less of a headache.

There's even a header-only version, called Spectra:
https://spectralib.org/
Last edited on
JamieAl (236)
Well, spent the past 24 hrs setting up Eigen with my VS code (was not easy) and now I am trying to rewrite some of my my arrays including std::vector to work with eigen and will see if that works better. My only concern if Eigen works with fftw3 library in C++ since I have already written a lot of functions with fftw3. Thanks all.
Last edited on
deleted account xyzzy (5768)
JamieAl wrote:
Well, spent the past 24 hrs setting up Eigen with my VS code (was not easy)

There is a somewhat easy way to manage 3rd party library package integration with Visual Studio, vcpkg.
https://vcpkg.io/en/index.html

It is a command-line tool, but once you get it working installation of a library is a single command away. You can search for packages at either the website or in a command prompt with vcpkg.

Eigen is one of the available packages:
D:\Programming\vcpkg>vcpkg search eigen
ceres[eigensparse]                        Use of Eigen as a sparse linear algebra library in Ceres
eigen3                   3.4.0            C++ template library for linear algebra: matrices, vectors, numerical solv...
highfive[eigen3]                          Enable Eigen testing
libigl                   2.3.0#1          libigl is a simple C++ geometry processing library. We have a wide functio...
libigl[embree]                            Build with embree
libigl[glfw]                              Build with glfw
libigl[imgui]                             Build with imgui
libigl[opengl]                            Build with opengl
libigl[xml]                               Build with libxml
magnum-integration[eigen]                 EigenIntegration library
opencv[eigen]                             Eigen support for opencv
opencv2[eigen]                            Eigen support for opencv
opencv3[eigen]                            Eigen support for opencv
opencv4[eigen]                            Eigen support for opencv
pangolin[eigen]                           Build support for Eigen matrix types
spectra                  1.0.0            A header-only C++ library for large scale eigenvalue problems
uvatlas[eigen]                            Use Eigen & Spectra for eigen-value computations
The result may be outdated. Run `git pull` to get the latest results.

If your port is not listed, please open an issue at and/or consider making a pull request:
    https://github.com/Microsoft/vcpkg/issues

The packages I have installed at the moment (and somewhat use):
D:\Programming\vcpkg>vcpkg list
7zip:x64-windows                                   21.07#1          Library for archiving file with a high compressi...
7zip:x86-windows                                   21.07#1          Library for archiving file with a high compressi...
bigint:x64-windows                                 2010.04.30#7     C++ Big Integer Library
bigint:x86-windows                                 2010.04.30#7     C++ Big Integer Library
catch2:x64-windows                                 3.0.1#1          A modern, header-only test framework for unit te...
catch2:x86-windows                                 3.0.1#1          A modern, header-only test framework for unit te...
cppunit:x64-windows                                1.15.1#3         Unit testing framework module for the C++ progra...
cppunit:x86-windows                                1.15.1#3         Unit testing framework module for the C++ progra...
fmt:x64-windows                                    8.1.1#1          Formatting library for C++. It can be used as a ...
fmt:x86-windows                                    8.1.1#1          Formatting library for C++. It can be used as a ...
gsl-lite:x64-windows                               0.40.0           A single-file header-only implementation of ISO ...
gsl-lite:x86-windows                               0.40.0           A single-file header-only implementation of ISO ...
gtest:x64-windows                                  1.11.0#5         GoogleTest and GoogleMock testing frameworks
gtest:x86-windows                                  1.11.0#5         GoogleTest and GoogleMock testing frameworks
ms-gsl:x64-windows                                 4.0.0            Microsoft implementation of the Guidelines Suppo...
ms-gsl:x86-windows                                 4.0.0            Microsoft implementation of the Guidelines Suppo...
pdcurses:x64-windows                               3.9#3            Public Domain Curses - a curses library for envi...
pdcurses:x86-windows                               3.9#3            Public Domain Curses - a curses library for envi...
vcpkg-cmake-config:x64-windows                     2022-02-06#1
vcpkg-cmake-config:x86-windows                     2022-02-06#1
vcpkg-cmake:x64-windows                            2022-05-10#1
vcpkg-cmake:x86-windows                            2022-05-10#1

(using the output tags works so much better than the quote tags)
Last edited on
MikeyBoy (5631)
Okay this sounds way too messy for just a simple one line of code in MATLAB. I am pretty confused and shocked how the alternative in C++ is waaay harder.

C++ is designed to be a powerful, general-purpose programming language that allows developers to express an enormous variety of concepts while still maintaining backward-compatibility with older versions of the language, and C, as much as possible.

Matlab is designed specifically to enable mathematical experts to express, as easily as possible, complicated mathematical expressions and models.

What would be shocking and confusing would be if it were just as easy to express mathematical expressions in C++ as it were in Matlab, because that would mean one of those things wasn't doing its job properly.
Last edited on
deleted account xyzzy (5768)
There is a 3rd party library available for reading and writing MATLAB binary files.
https://github.com/tbeu/matio

Available via vcpkg.
deleted account xyzzy (5768)
Also available via vcpkg is the FFTW3 library:
https://vcpkg.io/en/packages.html
jonnin (11497)
I spent about 3 years making a tool to convert matlab to c++ but that was in the late 90s, the code is horribly out of date now. And it wasn't that much stuff, basic matrix and vector math, eigenvalues, inverse, solver, LU decomp, lyap, kroneker (sp?), a few dozen other common items.

but it started with lyap.
one matlab line, one function call .. probably 20 pages of C++ just for that one, and 3 different ways to do it before I found one that was fast enough.
JamieAl (236)
Thanks everyone for the incredible feedback/info! I am working with Eigen now and have managed to get few things working. I think, I will have to make a new post regarding new questions. Thanks again.
Topic archived. No new replies allowed.