problem with my row echelon algorithm

[Sprechen sie von C++]

Seems to be a problem with my code, can anyone see a problem?


template <typename base> matrix<base> reduced_row_echelon_form(matrix<base>
&m) {
matrix<base> a;
base pivot, scale;
int pivotr, r;
int rows = (int)m.data.size();
int cols = (int)m.data[0].size();
a.resize(rows, cols);
a = m;
pivot = 0;
r = 0;
for (int j = 0; j < cols; j++) {
if (r >= rows) return a;
pivot = (base)0.0;
for (int i = r; i < rows; i++)
if (abs(a[j]) > pivot) {
pivot = abs(a[j]);
pivotr = i;
}
if (pivot <= (base)NZERO) {
for (int i = r; i < rows; i++) a[j] = (base)0;
continue; // try the next pivot candidate
}
std::swap(a[pivotr], a[r]);
// Eliminate elements below the pivot
scale = a[r][j];
a[r][j] = (base)1;
a[r] = a[r] / scale;
for (int i = r + 1; r < rows; i++) {
scale = a[j];
a[j] = (base)0;
a = a - a[r] * scale;
}
r++;
}
return a;
}

 

[red floyd]

Seems to be a problem with my code, can anyone see a problem?
[redacted]

Care to be a bit more specific? What are you seeing that you
believe to be a problem?

 

[Salt_Peter]

Seems to be a problem with my code, can anyone see a problem?

template <typename base> matrix<base> reduced_row_echelon_form(matrix<base>
&m) {
        matrix<base> a;
        base pivot, scale;
        int pivotr, r;
        int rows = (int)m.data.size();
        int cols = (int)m.data[0].size();
        a.resize(rows, cols);
        a = m;
        pivot = 0;
        r = 0;
        for (int j = 0; j < cols; j++) {
                if (r >= rows) return a;
                pivot = (base)0.0;
                for (int i = r; i < rows; i++)
                        if (abs(a[j]) > pivot) {
                                pivot = abs(a[j]);
                                pivotr = i;
            }
                if (pivot <= (base)NZERO) {
                        for (int i = r; i < rows; i++) a[j] = (base)0;
                        continue;  // try the next pivot candidate
                }
                std::swap(a[pivotr], a[r]);
                // Eliminate elements below the pivot
                scale = a[r][j];
                a[r][j] = (base)1;
                a[r] = a[r] / scale;
                for (int i = r + 1; r < rows; i++) {
                        scale = a[j];
                        a[j] = (base)0;
                        a = a - a[r] * scale;
         }
      r++;
        }
        return a;

}

Yes, i see pathetic code. How pathetic? Its so bad that even you can't
dissect what that function should be doing. If you can't provide a
compilable snippet that at least gives us a small clue as to what this
function is suppose to do, then go figure it out yourself.

From the point of view of constantness, casting, dependant types and
the 100 questions that would need to be asked about those particulars
we know nothing about in your code, an answer could easily require
1000 lines. C++ programmers are terrible at guessing and assuming.

Old-time (casts) are a no-no
A function parameter that isn't modified should be const
dependent types need keyword typename

the rest is noise

 

[Stuart Redmann]

Seems to be a problem with my code, can anyone see a problem?

Apart from what red_floyd and Salt_Peter have said, there seem to be
some algorithmic problems in your code (see below).

template <typename base> matrix<base> reduced_row_echelon_form(matrix<base>
&m) {
        matrix<base> a;
        base pivot, scale;
        int pivotr, r;
        int rows = (int)m.data.size();
        int cols = (int)m.data[0].size();
        a.resize(rows, cols);
        a = m;
        pivot = 0;
        r = 0;

Is r the same as the loop index j? If yes, drop this superfluous
variable or rename the loop variable to r. Everything else is just
confusing.

        for (int j = 0; j < cols; j++) {
                if (r >= rows) return a;
                pivot = (base)0.0;
                for (int i = r; i < rows; i++)
                        if (abs(a[j]) > pivot) {
                                pivot = abs(a[j]);
                                pivotr = i;
            }
                if (pivot <= (base)NZERO) {
                        for (int i = r; i < rows; i++) a[j] = (base)0;
                        continue;  // try the next pivot candidate
                }
                std::swap(a[pivotr], a[r]);
                // Eliminate elements below the pivot
                scale = a[r][j];

You want to multiply the pivot row by some factor that the pivot
element becomes one. However, you first set only the pivot element to
one ("a[r][j] = (base)1;") and then perform the multiplication ("a[r]
= a[r] / scale;"). So if your pivot element was x before these two
lines, it should then be 1/x (instead of one).

                a[r][j] = (base)1;
                a[r] = a[r] / scale;
                for (int i = r + 1; r < rows; i++) {
                        scale = a[j];
                        a[j] = (base)0;
                        a = a - a[r] * scale;
         }
      r++;
        }
        return a;



}>

Apart from that: Please DO document your code more extensively. As a
basic rule of programming you should first state what you want to do
with your own words (the comment), and only then add the code that you
think does what you just said. For example a comment like "// Find the
pivot element in the current column (the element that lies furthest
away from zero)." will make your code readable for everyone (even
those who cannot immediately remember which element should be the
pivot element). I once had to tutor students who had to write small
programmes in order to learn some programming language. When I
reviewed their homework, most of them could not remember what they had
written, and that after only one week!

Regards,
Stuart

 

[Jorgen Grahn]

....

Apart from that: Please DO document your code more extensively. As a
basic rule of programming you should first state what you want to do
with your own words (the comment), and only then add the code that you
think does what you just said. For example a comment like "// Find the
pivot element in the current column (the element that lies furthest
away from zero)." will make your code readable for everyone (even
those who cannot immediately remember which element should be the
pivot element).

Some people recommend that, but I don't believe it's a "basic rule of
programming" except maybe in assembly programming.

(I don't like it as a general rule, because IME people who comment
inline a lot describing *what* rather than *why* tend not to read
their own comments, and it pretty soon goes out of sync, and just
misleads the reader and takes up screen space. Plus, it takes focus
away from writing readable code, choosing readable names and so on.)

I noted though that the comment *at the top* was absent. That's
usually a really bad idea. But then I assumed that everyone who cares
already knows exactly what "reduced_row_echelon_form" means, just like
I know the concepts in *my* usual problem domains.

/Jorgen

 

[Sprechen sie von C++]

Here is the output from a test run:

Original matrix
1 2 1 3 3
2 4 0 4 4
1 2 3 5 5
2 4 0 4 7


reduced row echelon matrix program
0.5 2 0 2 2
-0.166667 0 0.333333 1 1
-0.333333 0 0 0 0.333333
-0.333333 0 -0.333333 0 0


manually solved reduced row echelon matrix
1 2 0 2 0
0 0 1 1 0
0 0 0 0 1
0 0 0 0 0


here is the code, with a few more comments

template <typename base> matrix<base> reduced_row_echelon_form(matrix<base>
&m) {
matrix<base> a;
base pivot, scale;
int pivotr, r;
int rows = (int)m.size();
int cols = (int)m[0].size();
a.resize(rows, cols);
a = m;
pivot = 0;
r = 0;
for (int j = 0; j < cols; j++) { // sweep by columns first
if (r >= rows) return a;
pivot = (base)0.0;
for (int i = r; i < rows; i++) // find the pivot
if (abs(a[j]) > pivot) {
pivot = abs(a[j]);
pivotr = i;
}
if (pivot <= (base)NZERO) {
for (int i = r; i < rows; i++) a[j] = (base)0;
continue; // try the next pivot candidate
}
std::swap(a[pivotr], a[r]);
// Eliminate elements below the pivot
scale = a[r][j];
a[r][j] = (base)1;
a[r] = a[r] / scale;
for (int i = r + 1; i < rows; i++) {
scale = a[j];
a[j] = (base)0;
a = a - a[r] * scale;
}
r++;
}
return a;
}

 

[Stuart Redmann]

Some people recommend that, but I don't believe it's a "basic rule of
programming" except maybe in assembly programming.

(I don't like it as a general rule, because IME people who comment
inline a lot describing *what* rather than *why* tend not to read
their own comments, and it pretty soon goes out of sync, and just
misleads the reader and takes up screen space. Plus, it takes focus
away from writing readable code, choosing readable names and so on.)

Point taken. I know from experience that even I fail to update the
comments, especially if the change had to be made quickly. Luckily,
I'll spot these inconsistency when I check in my code, since this only
happens when things cool off a bit: my code quality (that includes the
comments in my eyes) has definitely improved since I use CVS.

I noted though that the comment *at the top* was absent. That's
usually a really bad idea. But then I assumed that everyone who cares
already knows exactly what "reduced_row_echelon_form" means, just like
I know the concepts in *my* usual problem domains.

Agreed. Anyone reading your code should have at least a basic
knowledge of the problem domain at hand. Telling the user that the
pivot element is the element that lies furthest away from zero is
probably redundant, I just wanted to get over my point that the
comment should give as much information as necessary.

Sometimes I play with the idea to use the Knuth's WEB programming
language, but Knuth's book about Tex (which is effectively Tex's
source code) was not a pleasing sight to my eyes.

Regards,
Stuart

 

[Jorgen Grahn]

Point taken. I know from experience that even I fail to update the
comments, especially if the change had to be made quickly.

I thought I always updated my (sparse) comments, but recently found
one source file at work where the comments were way off, and had been
for quite a while. A humbling experience.

Luckily,
I'll spot these inconsistency when I check in my code, since this only
happens when things cool off a bit: my code quality (that includes the
comments in my eyes) has definitely improved since I use CVS.

Yes; that's one of the most powerful tools that version control offers
to you. I probably diff once every five minutes or so on average when
I'm coding, and check in/commit several times per hour.

/Jorgen

 

[Öö Tiib]

Point taken. I know from experience that even I fail to update the
comments, especially if the change had to be made quickly. Luckily,
I'll spot these inconsistency when I check in my code, since this only
happens when things cool off a bit: my code quality (that includes the
comments in my eyes) has definitely improved since I use CVS.

One trick is to use comments that answer to "what?" only in doxygen
headers. There they are in place. Elsewhere use only "why?" comments.

Sometimes you feel need for a "what?" comment within code. That proves
the code feels not readable enough for you.
There may be two reasons:
A) You are not familiar with the problem domain. That is best studied
from books and not comments.
B) Code is worth to be refactored to more readable. Adding comments
does not usually accomplish it.