matrix multiplication

hello,

im currently writing a class which stores matrices and performs arithmic operations with them.
the matrices are stored as follows:

Code:

---

float *pfaMatrix = new float[nColumnCount + nRowCount];

---

to access the matrix easily, the following is used:

Code:

---

index = nRow * nColumnCount + nCol

---

using those methods, i can use the array as if it was a two-dimensional array but without the restriction that i have to do this

Code:

---

float *pfaMatrix = new float[nRowCount][CONST];

---

my problem is now, that i do not really know, how to write the function of the overloaded * operator to multiply two matrices. maybe i just need the right approach. in the beginning of the funct. i already made a new matrix of the correct size:

Code:

---

float *pfaResMatrix = new float[cMatrix1.m_nRowCount + cMatrix2.m_nColumnCount];

---

i hope someone could help me! thanks!

yours sincerely,
replax

to start with, your matrix isnt the right size

Quote:

---

Code:

---

float *pfaMatrix = new float[nColumnCount + nRowCount];

---

---

as for your question, do you know how to multiply matrices? can you convert that into a series of 'mechanical actions'?

thanks for your quick response!

why do you think that the size is wrong? because it works well for me (doing matrix additions etc.).

i know how to multiply matrices, but heres the problem better illustrated:

Code:

---

a b  *  y  =  n
c d      z      m

---

the program would have to do:
[code]n=a*y+b*z and m=c*y+d*z

i figured that it has to be somewhat like this:

for(int nRo = 0; nRow

the program would have to do:

Code:

---

n=a*y+b*z and m=c*y+d*z

i figured that it has to be somewhat like this:

for(int nRo = 0; nRow < resmatrix.m_nRowCount; nRow++)
{
for(int nCol = 0; nCol < resmatrix.m_nColumnCount; nCol++)
{
now i have acces to the matrix element holding the resolution(n)
resmatrix.pfaMatrix[nRow * resmatrix.m_nColumnCount + nCol] = /*What should the equation be here*/
}
}

---

sorry for the screwed up double post.. no edit button.

Quote:

---

why do you think that the size is wrong? because it works well for me (doing matrix additions etc.).

---

if your matrix is 3x3 how many elements are there? what is 3 + 3?

Two matrices, A(i,j) and B(k,l) can only be multiplied if the inner dimensions are the same, i.e. j==k. The resulting matrix will be of size i x l.

therefore

M(i,l) = A(i,j) * B (k,l)

the 'i,j th' element of M is the sum over r of (A_ir * B_rj).

so you have too loop over row and column for each element of M (as you have already got), and then to calculate the value, you need *another* loop.

thanks again.
You were of course right with the size, I always tested with 2x2 matrices.

as for

Quote:

---

Originally Posted by Amleto

M(i,l) = A(i,j) * B (k,l)
the 'i,j th' element of M is the sum over r of (A_ir * B_rj).

so you have too loop over row and column for each element of M (as you have already got), and then to calculate the value, you need *another* loop.

---

could you maybe explain the text in bold a little bit further? thanks

yours sincerely,

replax

Code:

---

A=
1 2 3
0 0 0
0 0 0

B =
0 3 0
0 2 0
0 1 0

---

it means something like multiply the elements of the ith row of A with the elements of the jth column of B to get the i-jth element of M (where M = AB).

e.g.

Code:

---

M_12 = A_11 * B_12 + A_12 * B_22 + A_13 * B_32
M_12 = (1    *  3)  + ( 2  *  2)  +  ( 3  *  1)

---

thanks, I got that part.
now I am a little further, still I am getting not total correct results

Code:

---

Matrix operator*(const Matrix &cObj1, const Matrix &cObj2)
{
        float *pfaResMat = new float[cObj1.m_nRowCount * cObj2.m_nColCount];

        int i, k, j;

        for(i = 0; i < cObj1.m_nRowCount; i++)
        {
                for(j = 0; j < cObj2.m_nColCount; j++)
                {
                        for(k = 0; k < cObj1.m_nColCount; k++)
                        {
                                pfaResMat[i * cObj2.m_nColCount + j] += cObj1.m_pfMatrix[i * cObj1.m_nColCount + k] * cObj2.m_pfMatrix[k * cObj2.m_nColCount + j];
                        }
                }
        }
        return Matrix(cObj1.m_nRowCount, cObj2.m_nColCount, pfaResMat);
}

---

with that code, when cObj1 and Obj2 are 2x2 matrices with 1 2 \n3 4 , the result is:
2 4
6 8 instead of

7 10
15 22

do you see the flaw in the code?

yours sincerely,

replax

Quote:

---

Originally Posted by replax

thanks, I got that part.
now I am a little further, still I am getting not total correct results

---

1) Are you using a debugger to solve your problems? If not, I recommend you do so.

2) Why doesn't your Matrix constructor generate the memory?

Code:

---

Matrix operator*(const Matrix &cObj1, const Matrix &cObj2)
{
    float *pfaResMat = new float[cObj1.m_nRowCount * cObj2.m_nColCount];  // This line.

---

I would have expected something like this:

Code:

---

Matrix operator*(const Matrix &cObj1, const Matrix &cObj2)
{
  Matrix temp(cObj1.m_nRowCount, cObj2.m_nColCount);  // creates a matrix

---

And then, you manipulate temp. Having the client generate the memory for the matrix is disjoint and bad design.

3) Stemming from this problem is the reason why your code doesn't work. Where do you initialize this "new" memory to 0?

Code:

---

float *pfaResMat = new float[cObj1.m_nRowCount * cObj2.m_nColCount];

---

This does not set the data to 0. Then you do this:

Code:

---

pfaResMat[i * cObj2.m_nColCount + j] += cObj1.m_pfMatrix[i * cObj1.m_nColCount + k] * cObj2.m_pfMatrix[k * cObj2.m_nColCount + j];

---

You are adding something to junk, giving you junk.

If you designed the matrix where it did its own memory handling, then you probably would have noticed this problem sooner.

Regards,

Paul McKenzie

thank you for your answer!

Quote:

---

Originally Posted by Paul McKenzie

2) Why doesn't your Matrix constructor generate the memory?

---

the problem is, if I use my Matrix constructor to construct a matrix, it will call the destructor before returning the matrix value. thus crashing the programm.

the constructor is:

Code:

---

Matrix::Matrix(int nRowCount, int nColCount = 1) : m_nColCount(nColCount), m_nRowCount(nRowCount)
{
    m_pfMatrix = new float[m_nColCount * m_nRowCount];
        for(int nRow = 0; nRow < m_nRowCount; nRow++)
        {
                for(int nCol = 0; nCol < m_nColCount; nCol++)
                        m_pfMatrix[nRow * m_nColCount + nCol] = 0;
        }
}

Matrix::Matrix(int nRowCount, int nColCount, float *fMatrix) : m_nColCount(nColCount), m_nRowCount(nRowCount)
{
        m_pfMatrix = new float[m_nColCount * m_nRowCount];
    for(int nRow = 0; nRow < m_nRowCount; nRow++)
    {
        for(int nCol = 0; nCol < m_nColCount; nCol++)
        {
            m_pfMatrix[nRow * m_nColCount + nCol] = fMatrix[nRow * m_nColCount + nCol];
        }
    }
}

---

and the destructor is:

Code:

---

Matrix::~Matrix()
{
    if(m_pfMatrix)
    {
        delete[] m_pfMatrix;
        m_pfMatrix = 0; //NULL
    }
}

---

i solved the problem, i did an unbelievably stupid error. in the main.cpp file where i did execute the arithmic expressions i typed '+' instead of '*', which of course gives me my matrix addition instead of multiplication.

thanks to all of you for your kind help!

regards,
replax

Quote:

---

Originally Posted by replax

the problem is, if I use my Matrix constructor to construct a matrix, it will call the destructor before returning the matrix value. thus crashing the programm.

---

That does not really make sense. The destructor will not normally be called from the constructor.

On a hunch: did you implement the copy constructor and copy assignment operator?

no, I did not implent the copy constructor.
So when doing

Code:

---

return cResMatrix;

---

return does copy the matrix and then returns the copy because the original matrix goes out of scope, right? and when the matrix goes out of scope (after the curly braces) it calles the destructor, correct?

so I would need to write a copy constructor, so that return knows how to copy a matrix? buy why do I not need to do that when I use:

Code:

---

return Matrix(nRowCount, nColCount, pfaResMat);

---

because the returned type is a matrix, too.

Or did I get something essential completly wrong?

regards,

replax

Quote:

---

Originally Posted by replax

no, I did not implent the copy constructor.
So when doing

Code:

---

return cResMatrix;

---

return does copy the matrix and then returns the copy because the original matrix goes out of scope, right? and when the matrix goes out of scope (after the curly braces) it calles the destructor, correct?

so I would need to write a copy constructor, so that return knows how to copy a matrix? buy why do I not need to do that when I use:

Code:

---

return Matrix(nRowCount, nColCount, pfaResMat);

---

because the returned type is a matrix, too.

Or did I get something essential completly wrong?

regards,

replax

---

The principle of the matter is robust code, that means that you should minimise the number of places that you allocate and deallocate memory in your class - when you get on to the subject of exception saftey you will appreciate this point more. A proven pattern for a safe/robust way of allocation and deallocation of resources is the RAII (Resource Acquisition Is Initialization) idiom. This would be a good thing for you to look up, as well as the Rule of Three - i.e. if you write one of the following for your class, then you should probably write all three:

- destructor
- copy constructor
- copy assignment operator

The point in hand, is that you could save yourself a lot of bother and bugs if you adhere to the advice that you are being given with respect to memory management.

Here is an example pattern that you could follow:

Code:

---

class matrix
{
public:
  matrix()
    :rows_(0)
    ,columns_(0)
    ,data_(0)
  {}

  matrix(size_t rows, size_t columns)
    :rows_(rows)
    ,columns_(columns)
    ,data_(new float[rows_*columns_])
  {
    for(size_t i = 0; i<rows_*columns_; i++)
    {
      data_[i] = 0;
    }
  }

  matrix(const matrix& other)
    :rows_(other.rows_)
    ,columns_(other.columns_)
    ,data_(new float[rows_*columns_])
  {
    size_t size = rows_*columns_;
    std::copy(other.data_, other.data_+size, data_);
  }

  void swap(matrix& other)
  {
    std::swap(rows_, other.rows_);
    std::swap(columns_, other.columns_);
    std::swap(data_, other.data_);
  }

  matrix& operator=(const matrix& rhs)
  {
    matrix tmp(rhs);
    swap(tmp);
    return *this;
  }

  ~matrix()
  {
    delete [] data_;
  }

  float& operator()(size_t row, size_t column)
  {
    return data_[row*columns_+column];
  }

  const float& operator()(size_t row, size_t column) const
  {
    return data_[row*columns_+column];
  }

private:
  size_t rows_;
  size_t columns_;
  float* data_;
};

---

Of course, you could save yourself even more if you use std::vector:

Code:

---

#include <vector>

class matrix
{
public:
  matrix()
    :rows_(0)
    ,columns_(0)
    ,data_()
  {}

  matrix(size_t rows, size_t columns)
    :rows_(rows)
    ,columns_(columns)
    ,data_(rows_*columns_)
  {}

  float& operator()(size_t row, size_t column)
  {
    return data_[row*columns_+column];
  }

  const float& operator()(size_t row, size_t column) const
  {
    return data_[row*columns_+column];
  }

private:
  size_t rows_;
  size_t columns_;
  std::vector<float> data_;
};

---

std::vector will handle the allocation and deallocation for you - you also do not need to write a destructor, a copy assignment operator or a copy constructor, since the default generated ones will be perfectly fine.

Note that the operator()(size_t rows, size_t columns) allows you to access the relevent matrix element.

Code:

---

int main()
{
  matrix m(3,3); //create a 3x3 matrix

  m(1,1) = 5.2f; //set the central element to 5.2
}

---

Quote:

---

Originally Posted by replax

thank you for your answer!

the problem is, if I use my Matrix constructor to construct a matrix, it will call the destructor before returning the matrix value. thus crashing the programm.

---

That is because you didn't write an appropriate copy constructor and assignment operator. How many C++ classes have you used, where your code needs to allocate the memory for the class? I bet none, except only yours.

For example std::string -- why is it that we don't need to allocate memory for the string data? What is the magic that std::string does that your Matrix class is not doing? The answer is that the std::string has implemented proper copy semantics, and your Matrix class hasn't.

Code:

---

#include <string>

std::string DoStuff()
{
  std::string s = "abc123";
  s += "456";
  s = s.substr(0,1);
  return s;
}

---

There is no crash, because on return, the std::string copy constructor makes a copy of s and returns it. Without a user-defined copy constructor, the std::string would suffer the same problem that your Matrix class suffers from.

To correct the problem, I could start naively allocating memory all over the place to make sure the crash doesn't happen, or I give std::string some intelligence and write functions (copy/assignment operator) that tells std::string what to do when making copies of itself.

When you create an object, that object is supposed to handle its own memory, and clean up its own memory. The user of the object shouldn't be doing this work -- it the user has to do this work, then that class is basically useless in a program. The user has to jump through hoops keeping the object valid, doing work the class should be doing. That defeats the whole purpose of using the class to begin with.

Quote:

---

i solved the problem, i did an unbelievably stupid error. in the main.cpp file where i did execute the arithmic expressions i typed '+' instead of '*', which of course gives me my matrix addition instead of multiplication.

---

That doesn't mention the problem I pointed out to you, and that is you are allocating memory that is not initialized to 0. Did you fix that problem?

Regards,

Paul McKenzie

thanks for the very constructive input!

as for now, I got it working as far as I can see it. Now my constructors are the same, I just fixed the allocation that was not initialized to 0 and added the copy constructor:

Code:

---

Matrix::Matrix(Matrix &cSource)
{
        m_nRowCount = cSource.m_nRowCount;
        m_nColCount = cSource.m_nColCount;
        m_pfMatrix = new float[m_nRowCount * m_nColCount];
        for(int nRow = 0; nRow < m_nRowCount; nRow++)
        {
                for(int nCol = 0; nCol < m_nColCount; nCol++)
                        m_pfMatrix[nRow * m_nColCount + nCol] = cSource(nRow, nCol);
        }
}

---

and also the overloaded assingment operator:

Code:

---

Matrix& Matrix::operator=(Matrix &cSource)
{
        if(this == &cSource)
                return *this;

    m_nRowCount = cSource.m_nRowCount;
    m_nColCount = cSource.m_nColCount;
    for(int nRow = 0; nRow < m_nRowCount; nRow++)
    {
        for(int nCol = 0; nCol < m_nColCount; nCol++)
            m_pfMatrix[nRow * m_nColCount + nCol] = cSource(nRow, nCol);
    }
    return *this;
}

---

so that my multiplication algorithm now looks like this:

Code:

---

Matrix operator*(Matrix &cObj1, Matrix &cObj2)
{
        Matrix pfaResMat(cObj1.m_nRowCount, cObj2.m_nColCount);
        int i, k, j;

        for(i = 0; i < pfaResMat.m_nRowCount; i++)
        {
                for(j = 0; j < pfaResMat.m_nColCount; j++)
                {
                        pfaResMat(i, j) = 0;
                        for(k = 0; k < cObj1.m_nColCount; k++)
                                pfaResMat(i, j) += cObj1(i, k) * cObj2(k, j);
                }
        }
        return pfaResMat;
}

---

and it works like a charm as far as I can tell.
all other overloaded operators( (), *(by int or float), +) also work perfectly well.

if you see some flaws, leaks or anything that should be inprooved, I would be very happy to know :)

regards,

replax

Your copy constructor and copy assignment operator should have the parameter as const Matrix&. After all, it should make sense to copy a matrix that is constant.

I notice that the copy assignment operator copies the row count and column count, i.e., the number of rows and columns of the source matrix could differ from the destination matrix. However, it makes no attempt to allocate space if the source matrix has more rows or columns than the destination matrix.

Thanks for your information.

Initially, I had them as const Matrix&, however, as I am using the overloaded () operator inside the constructor and assignment operator, I am not allowed to use const.
compilier error:

Code:

---

error: passing 'const Matrix' as 'this' argument of 'float& Matrix::operator()(int, int)' discards qualifiers

---

maybe it's just that my () operator is incorectly overloaded:

Code:

---

float& Matrix::operator()(const int nRow, const int nCol)
{
    return m_pfMatrix[nRow * m_nColCount + nCol];
}

---

as for the copy assignment operator, I'm gonna fix it :)
thanks for your input!

regards,

replax

It sounds like you need to const overload operator(), i.e.,

Code:

---

float& Matrix::operator()(const int nRow, const int nCol)
{
    return m_pfMatrix[nRow * m_nColCount + nCol];
}

const float& Matrix::operator()(const int nRow, const int nCol) const
{
    return m_pfMatrix[nRow * m_nColCount + nCol];
}

---

I thought of that, but then I would not be able to modify the matrix using

Code:

---

Matrix(nRow, nCol) = X;

---

is there another way of doing this?

Quote:

---

Originally Posted by replax

I thought of that, but then I would not be able to modify the matrix using

---

Notice that there are two versions of operator().

Ah ok I thought you ment to replace it. Anyway, I figured that it works like that:

Code:

---

float& Matrix::operator()(const int nRow, const int nCol) const
{
    return m_pfMatrix[nRow * m_nColCount + nCol];
}

---

thank you very much for your input.

replax

Yeah, though it may be better to have the const version return a const reference.

In which case would it get called?

When the Matrix object is const, or when you are accessing the Matrix object through a const reference or a pointer to a const Matrix.

As laserlight pointed out, the operator= needs to reallocate memory when the matrix size is changing. One easy way to do this is using the copy-and-swap trick that Paul demonstrated earlier.

It bears repeating that std::vector should be preferred to dynamic arrays whenever possible, because it's both safer and easier to work with.

One other point worth mentioning----you should not assume that returning a matrix from a function by value will invoke the copy constructor to create a copy, and the destructor of the matrix that's local to the function. It might; that's the behavior to expect, certainly; but don't assume it must do that, because one common form of compiler optimization tries to remove such unnecessary copies using return value optimization. Whatever it does end up doing will be equivalent to the above, but only if your copy constructor and operator= really do create logical copies. So don't try to get clever with those functions.

This is the almost final program, thanks for your help!

http://rapidshare.com/files/369676516/matrix.zip

@Lindley, thanks for the information on return values and copy constructors!
as for the vector, it was a kind of "requirment" to do it with arrays. But I agree that using a std::vector would be better!

Quote:

---

Originally Posted by replax

This is the almost final program, thanks for your help!

http://rapidshare.com/files/369676516/matrix.zip

@Lindley, thanks for the information on return values and copy constructors!
as for the vector, it was a kind of "requirment" to do it with arrays. But I agree that using a std::vector would be better!

---

Yes, usage of vector removes the need to write a copy ctor/assignment op (and destructor).

Also, as Lindley pointed out, your copies should be logically equivalent. In other words, if you replace the original with the copy, the copy must behave exactly the same as the original, regardless of the path of execution your code will take.

A source of bugs is writing copy/assignment functions that really are not making true copies of the object. This is called coding with side-effects. This type of coding is to be avoided when it comes to copy/assignment, as you don't really know when a copy will be made, due to the optimizations a compiler may use.

Regards,

Paul McKenzie