Advanced Introduction to C++, Scientific Computing and Machine Learning




Claudius Gros, WS 2019/20

Institut für theoretische Physik
Goethe-University Frankfurt a.M.

Elimination

linear equations

$$ \left.\begin{array}{cccccc} 2x&-3 y& -z & + 2w& =&4\\ 4x&-4 y& -z & + 4w& =&4\\ 2x&-5 y& -3z& + 3w& =&9\\ \end{array}\right\} \qquad\quad A\,\vec x=\vec b \qquad\quad \vec x = \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) $$

Gauss elimination

$$ \left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 4&-4&-1& 4&4\\ 2&-5&-3& 3&9\\ \end{array} \right) $$ $$ \left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 0& 2& 1& 0&-4\\ 0&-2&-2& 1&5\\ \end{array} \right) $$ $$ \left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 0& 2& 1& 0&-4\\ 0& 0&-1& 1&1\\ \end{array} \right) $$ $$ \tilde{A}\,\vec x = \tilde{b} \qquad\quad \left( \begin{array}{rrrr} 2&-3&-1& 2\\ 0& 2& 1& 0\\ 0& 0&-1& 1\\ \end{array} \right) \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} 4 \\ -4 \\ 1 \end{array}\right) $$

back substitution

$$ \left( \begin{array}{rrrr} 2&-3&-1& 2\\ 0& 2& 1& 0\\ 0& 0&-1& 1\\ \end{array} \right) \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} 4 \\ -4 \\ 1 \end{array}\right) $$ $$ -z+w=1 $$ $$ \begin{array}{rcl} z&=&w-1\\ 2y+z=-4\qquad\quad y &=&\frac{-3-w}{2} \\ 2x-3y-z+2w=4\qquad\quad x &=&\frac{-3-5w}{4} \end{array} $$

with

$$ 2x = 4+3y+z-2w= 4-\frac{9}{2}-\frac{3}{2}w+w-1-2w =-\frac{3}{2} -\frac{5}{2}w $$ $$ \vec x= \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} -(3+5w)/4 \\ -(3-w)/2 \\ w-1 \\ w \end{array}\right) \qquad\qquad \mbox{for any}\ \ w $$

pivoting

$$ \left( \begin{array}{cccc} *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{array} \right) $$ $$ \left( \begin{array}{cccc} *&[*]&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&[*]\\ *&0&*&*\\ *&0&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&*\\ *&0&[*]&0\\ *&0&*&0\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&*\\ *&0&*&0\\ *&0&0&0\\ \end{array} \right) $$

swapping rows

order of equations irrelevant

Gauss elimination with pivoting

Copy Copy to clipboad
Downlaod Download
#include <iostream>
#include <stdio.h>     // for printf
#include <string>      // c_str()
#include <cmath>       // for abs
#include <fstream>     // file streams
using namespace std;
 
// --- Gauss elimination of Ax=b with pivoting.
// --- On input quadratic matrix A and vector b.
// --- On ouput solution x.
 
template<int N>
void gaussElimination(double (&A)[N][N], double (&b)[N], double (&x)[N]) 
{
 double temp;                       // temp variable
 const double smallNumber = 1e-10;  // a small number
//
//--- loop over all columns A[][col];  A[][col]
 for (int col=0; col<N; col++) 
   {             
//--- find pivot row; swap pivot with actualy row
   int pivotRow = col; 
   for (int row=col+1;row<N;row++) 
     if ( abs(A[row][col]) > abs(A[pivotRow][col]) ) 
       pivotRow = row;
 
//--- swap row A[col][] with A[pivotRow][]
   for (int i=0;i<N;i++)
     {
     temp = A[col][i];
     A[col][i] = A[pivotRow][i];
     A[pivotRow][i] = temp;
     }
 
//--- swap elements b[col] with b[pivotRow]
   temp = b[col];  
   b[col] = b[pivotRow]; 
   b[pivotRow] = temp;
 
//--- throw a warning if matrix is singular or nearly singular
   if (abs(A[col][col]) <= smallNumber) 
     cerr << "Matrix is singular or nearly singular" << endl;
 
//--- Gauss with pivot within A and b
   for (int row=col+1;row<N;row++) 
     {
     double alpha = A[row][col] / A[col][col];  // divide by pivot
     b[row] -= alpha*b[col];
     for (int j=col;j<N;j++) 
       A[row][j] -= alpha * A[col][j];
     }
      }  // end of loop over all columns
 
//--- back substitution
 for (int i=N-1;i>=0;i--) 
   {
   double sum = 0.0;
   for (int j=i+1;j<N;j++) 
      sum += A[i][j] * x[j];
   x[i] = (b[i] - sum) / A[i][i];   // pivots on diagonal after swapping
   }
}  // end of gaussElimination()
 
 
// ---
// --- print system
// ---
template<int N>
void printSystem(double A[N][N], double b[N], double x[], string message) 
{
 printf("\n---------------- %s\n",message.c_str());
 for (int col=0;col<N;col++)
   printf("%7d",col);
 printf(" |%7s   :: %7s\n","b","x");
//
 for (int row=0;row<N;row++)
   {
   for (int col=0;col<N;col++)
     printf("%7.3f",A[row][col]);
   printf(" |%7.3f   :: %7.3f\n",b[row],x[row]);
   }
} // end of printSystem()
 
// ---
// --- main
// ---
int main() 
{
/*
 const int N = 2;
 double A0[N][N], A[N][N], b0[N], b[N], x[N];
 A0[0][0] = 1.0;  A0[0][1] = 3.0;   b0[0] = 1;
 A0[1][0] = 2.0;  A0[1][1] = 4.0;   b0[1] = 1;
*/
 const int N = 3;
 double A0[N][N], A[N][N], b0[N], b[N], x[N];
 A0[0][0] = 2.0;  A0[0][1] =-3.0; A0[0][2] = -1.0;   b0[0] = 4;
 A0[1][0] = 4.0;  A0[1][1] =-4.0; A0[1][2] = -1.0;   b0[1] = 4;
 A0[2][0] = 2.0;  A0[2][1] =-5.0; A0[2][2] = -3.0;   b0[2] = 9;
//
//
 for (int row=0;row<N;row++)
   {
   b[row] = b0[row];
   for (int col=0;col<N;col++)
     A[row][col] =A0[row][col];
   }
//
 printSystem(A,b,x,"at the start");
 gaussElimination(A,b,x);
 printSystem(A,b,x,"after elimination (A,b)");
 printSystem(A0,b0,x,"after elimination (A0,B0)");
//
 printf("\n");
 for (int row=0;row<N;row++)  // testing the result
   {
   double sum = 0.0;
   for (int col=0;col<N;col++)
      sum += A0[row][col]*x[col];
   printf("# row %d %8.3f | %8.3f\n",row, sum, b0[row]);
   }
//
 return 1;
}  // end of main()

matrices

LR decomposition

$$ A = L\,R \qquad\qquad A\,\vec x =\vec b $$ $$ L\,R\,\vec x =\vec b \qquad\qquad L\,\vec y =\vec b \qquad\qquad \vec y = R\,\vec x $$ $$ L\,\vec y = \vec b \qquad\qquad L = \left( \begin{array}{cccc} *&0&0&0\\ *&*&0&0\\ *&*&*&0\\ *&*&*&*\\ \end{array} \right) $$ $$ R\,\vec x = \vec y \qquad\qquad R = \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{array} \right) $$

LR decomposition with Gauss elimination