Need help sorting matix

[Dantecross12]

Hi I'm trying to write a program to rearrange a two different matrices such that the numbers below each number in the diagonal are smaller than the diagonal number {This is going to be a sub-routine for a linear equations solver). I was able to get the code to work for one array, but the other one is not sorting properly (It needs to line up with the first array). Any help or advice would be very much appreciated. Thanks

#include <iostream>

void main()
{
	double u[4][4] = { {1,6,3,11},{4,5,6,2},{7,-8,9,6},{12,2,-5,-12} };
	double b[4] = { 1,2,3,4 };

	for (int k = 0; k < 4; k++)
	{
		for (int i = k + 1; i < 4; i++)
		{
			if (u[k][k] < u[i][k])
			{
				for (int col = k; col < 4; col++)
				{
					std::swap(u[k][col], u[i][col]);
				}
				std::swap(b[k], b[i]);
				
			}
		}
	}
	for (int rowu = 0; rowu < 4; ++rowu)
	{
		for (int colu = 0; colu < 4; ++colu)
		{
			std::cout << u[rowu][colu] << " ";
		}
		std::cout << '\n';
	}

	std::cout << '\n';

	for (int i = 0; i < 4; i++)
	{
		std::cout << b[i] << " \n";
	}
}
 

[lastchance]

If you put your code in code tags we would be able to run it online and see what your problem is.

Your problem statement is unclear. If you simply want to line up columns nicely then you will have to format it: either by using setw() or by outputting a tab character ('\t') before each element output.

Your partial pivoting looks as if it would work (difficult to tell without the indentation that code tags would retain). However, you are doing an unnecessary number of swaps: find the largest subdiagonal element in the kth column and swap just that row. Usually, when this sort of swapping is done it is on the basis of the ABSOLUTE values of elements (i.e. disregarding signs) - not as you have at present. Incidentally, if this were a stand-alone subroutine, then it wouldn't help you with Gaussian elimination: in that algorithm this sort of swapping would be done once for each set of reduction operations, not as a stand-alone piece of sorting.

You are also making excessive use of the "magic number" 4: make it a const int that you define only in one place. Otherwise your code is only useful for 4x4 matrices.

BTW, it's int main(), not void main(), for standard C++.

[Dantecross12]

I've updated the post so that the code shows up properly. Sorry if I didn't explain this problem properly. I'm trying to set up a program that solves a set of equations using gaussian elimination with partial pivoting. I've gotten everything to work except for the partial pivoting. The program is able to reorganize array u properly, but messes up array b.

Thank you very much for your help. I've been struggling with this problem for days and I've slowly inched my way to the answer, but this last problem has left me stumped.

[lastchance]

No, it doesn't reorganise anything properly.

You don't understand how partial pivoting works within Gaussian elimination. You have to do both TOGETHER, not this row swapping separately. At the moment, lines 14-17 are swapping the tails of rows only - leaving the first part of the row in the wrong equation. This would be OK if you had already got zeros in this part by Gaussian elimination ... but you haven't.

So, do both together:

	for (int k = 0; k < KMAX; k++)               // k is the row you are working with; [k][k] is the diagonal element
	{

           // Partial pivoting - find the row i BELOW OR EQUAL k with the largest ABSOLUTE element in column k and swap rows i and k

           // Forward reduction - subtract multiples of the new row k from the rows below to get 0 in the rest of column k.
	}


	for (int k = KMAX-1; k >= 0; k++)
	{
             // Now do the back substitution
        }

If at any stage you get a zero diagonal element then the matrix is non-invertible.


Please make it int main(), then we can run your code as-is in CPP-shell.

[Dantecross12]

I'm sorry, I should have also included this. I'm trying to solve each step using simple matrices and then I incorporate them into the main program once I understand how it's done. I understand that the partial pivot step will need to be incorporated into the main loop for row reduction.

How do I get the program to locate the row below [k][k] with the largest element?

Thank you for all your help!

[lastchance]

In your code:

	//Performing row reduction
	for (int k = 0; k < x - 1; k++)
	{

	  // ====> FIND THE MAXIMUM abs(a[][]) BELOW [k][k] AND, IF NECESSARY,  DO A SINGLE ROW SWAP HERE <=====


	   for (int i = k + 1; i < x; i++)
	   {

It's rather different from your formulation but you can see how I implemented it below. I've removed the checking routines for brevity. The matrix corresponds to the first of your posts, not the most recent one.

FWIW I also wrote something to do it with parallel computing recently. A lot more fun when different rows are on different processors (potentially on different computers):
http://www.cplusplus.com/forum/unices/253337/#msg1114925

#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;

using vec = vector<double>;
using matrix = vector<vec>;   

// Function prototypes
bool GaussElimination( matrix A, vec B, vec &X );
void readData( matrix &A, vec &B );

//======================================================================

int main()
{
   matrix A;
   vec B, X;

   // Get equations
   readData( A, B );

   // Solve
   if ( GaussElimination( A, B, X ) )
   {
      cout << "Solution: " ;
      for ( auto e : X ) cout << e << ' ';
      cout << '\n';
   }
   else 
   {
      cout << "Unable to solve\n";
   }
}


//======================================================================


bool GaussElimination( matrix A, vec B, vec &X )
//-------------------------------------------------------------------
// Solve linear system AX = B by Gaussian elimination
// Note: copies made of A and B by passing arguments by value
//-------------------------------------------------------------------
{
   const double SMALL = 1.0E-30;
   int N = A.size();

   // Row operations for i = 0, ,,,, N - 1 (N-1 only used for singularity check)
   for ( int i = 0; i < N; i++ )
   {
      // Pivot: find row r below with largest element in column i
      int r = i;
      double maxA = abs( A[i][i] );
      for ( int k = i + 1; k < N; k++ )
      {
         double val = abs( A[k][i] );
         if ( val > maxA )
         {
            r = k;
            maxA = val;
         }
      }
      if ( r != i )                    // Swap row (up to i-1 is zero anyway)
      {
         for ( int j = i; j < N; j++ ) swap( A[i][j], A[r][j] );
         swap( B[i], B[r] );
      }

      // Forward reduction: row operations to clear column below this diagonal
      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;            // Singular matrix

      for ( int r = i + 1; r < N; r++ )                    // Lower rows
      {
         double multiple = A[r][i] / pivot;                // Multiple of row i to clear element in ith column
         for ( int j = i; j < N; j++ ) A[r][j] -= multiple * A[i][j];
         B[r] -= multiple * B[i];
      }
   }

   // Back-substitute
   X = B;
   for ( int i = N - 1; i >= 0; i-- )
   {
      for ( int j = i + 1; j < N; j++ )  X[i] -= A[i][j] * X[j];
      X[i] /= A[i][i];
   }

   return true;
}


//====================================================================== 


void readData( matrix &A, vec &B )
{
// ifstream in( "in.txt" );
   stringstream in( "4"
                    " 1  6  3  11     1\n"
                    " 4  5  6   2     2\n"
                    " 7 -8  9   6     3\n"
                    "12  2 -5 -12     4\n" );

   int N;
   in >> N;
   A = matrix( N, vec( N ) );
   B = vec( N );
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) in >> A[i][j];
      in >> B[i];
   }
}


//====================================================================== 

[Dantecross12]

Thank you very much. I'll go over the code after work. As you can tell, I'm not a programmer. In fact, I just started learning c++ a few days ago. I didn't think it would be this hard, but I'm sure it will get easier. I'll update here when I get the last piece of the code to work.

[Dantecross12]

Hey, so I understand most of how you went about doing the pivoting. It never occurred to me to use a variable to track the largest row. I implemented it into my code, and it runs, but it's giving me the wrong answer.

For some reason it's also modifying the first row when it does row reduction. Also, in your code, line 66, why did you set maxA = val ? That's the only part I don't understand. Again, thank you very very much for all your help! Also, here is my latest attempt to get it to work.

[lastchance]

You aren't quite looping correctly. The structure should be:

for each row k:
{
   FIND LARGEST ABS ELEMENT AT OR BELOW [k][k]
   {
      ...
   }

   SWAP ROWS r AND k (if necessary)
   {
      ...
   }

   DO THE REDUCTION FOR EACH ROW i BELOW ROW k
   { row i loop
      { columns j of A loop
      } end of columns j of A loop
      reduction of b
   } end row i loop

}

I have marked up the central portion of your code.

        //Performing Row reduction 

        for (int k = 0; k < m - 1; k++)
        {

                // ===== FIND LARGEST ABS ELEMENT AT OR BELOW [k][k]
                int r = k; //Assume the first element k is the largest element in the column. Number will be updated
                double maxa = abs(a[k][k]);

                for (int i = k + 1; i < m; i++)
                {
                        double maxb = abs(a[i][k]);
                        if (maxb > maxa)
                        {
                                r = i;
                                maxa = maxb;
                        }
                }  // <==== END THIS LOOP HERE *****


                // ====== SWAP ROWS (if necessary)
                if (r != k)
                {
                        for (int j = k; j < n; j++)
                        {
                                std::swap(a[k][j], a[r][j]);         // <===== SHOULD BE r NOT i *****
                        }
                        std::swap(b[k], b[r]);                       // <===== THERE IS ONLY ONE COLUMN IN b (MOVED ***** )
                }

                // ====== NOW DO THE REDUCTION BELOW ROW k
                for (int i = k + 1; i < m; i++)              // <==== NEED TO START A NEW LOOP HERE *****
                {                                            // <==== *****
                        p = a[i][k] / a[k][k];
                        for (int j = 0; j < n; j++)
                        {
                                a[i][j] -= p * a[k][j];
                        }

                        b[i] -= p * b[k];
                }
        }

This then runs and produces the answer

-0.0504009 
-0.408935 
-0.261168 
-0.520046 

which is correct ... for this matrix.



Your method is then fine for non-singular matrices. To perfect it, you need to be able to stop gracefully if at any point your diagonal element a[k][k] is zero (or floating-point small). Otherwise you will get a divide-by-zero error. For example, your code would produce nan and/or inf for the simple (singular) matrix

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 0

If you look back at my code you will note that I have a simple check

      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;            // Singular matrix 

at which point I chose to return a boolean false from the solving routine.

I also had to make sure that this was checked for the last row (which doesn't actually have any rows below it to reduce). (I had to edit my serial code here - sorry!)

If you wanted to do this you would have to:
- put a check on a[k][k] after the partial pivoting and before you did the below-row reduction
- amend your outer loop to < m rather than < m-1. (The last element wouldn't play any role in row reduction, it would just be there for a final diagonal check.)

	//Performing Row reduction 
	for (int k = 0; k < m; k++)   // m, not m-1 if you want to check a[k][k] != 0