In order to use the routines PLUFactorization and SolvingPhase, you will need to dynamically allocate the matrix A and the vectors b, x, and ipivot. The beginning of your main program should be:

#include

#include

#include

#include "GESolver.h"int main(int argc, char *argv[])

{

FILE *input = fopen(argv[1], "r");

FILE *output = fopen(argv[2], "w");

GESolver gauss;float detA;

int iflag;

int i;

int j;

int n;

fscanf(input,"%d",&n);

float **A = new float*[n];

for (i = 0; i

{

A[i] = new float[n];

}

float *b= new

float[n];

float *x= new

float[n];

int *ipivot = new int[n];

for (i = 0; i

{

for (j = 0; j

{

fscanf(input,"%f", &A[i][j]);

}

fscanf(input,"%f", &b[i]);

}

fprintf(output,"A, b: ");

for (i = 0; i

{

for (j = 0; j

{

fprintf(output,"%18.9f", A[i][j]);

}

fprintf(output,"%18.9f ", b[i]);

}fprintf(output," ");

where n is the order of the system. After you have solved a given system, you must free the dynamically allocated memory and close the files. This can be done using

for (i = 0; i

{

delete[] A[i];

}

delete[] A;

delete[] b;

delete[] x;

delete[] ipivot;

fclose(input);

fclose(output);

return (0);

}

at the end of your program. Run your program on the four systems given on the website.

If your program name is ms572hw1

you can run your program by typing

ms572hw1 inputFilename outputFilename

For example, using ms572hw1 system1.dat system1.out

will read the input file system1.dat and create the outputsystem1.out.

// GESolver.h

#ifndef _GESOLVER_

#define _GESOLVER_

class GESolver

{

public: GESolver(); // Constructor

~GESolver(); // Destructor

void PLUFactorization(float **w,

int n,

int *ipivot,

int &iflag);

void SolvingPhase(float **w,

int *ipivot,

float *b,

int n,

float *x);

private:

const float EPSILON = 0.000001f;

};

#endif

// GESolver.cpp

#include

#include

#include "GESolver.h"

GESolver::GESolver() // Constructor

{

};

GESolver::~GESolver() // Destructor

{

};

// This function performs the triangular factorization of a

// matrix using Gaussian elimination and scaled partial pivoting.

// It keeps track of the row interchanges in the array ipivot.

//

// NOTE: This is a modified version of the FACTOR subroutine

// on page 165 of Elementary Numerical Analysis by

// S. D. Conte and Carl de Boor written in C++. void GESolver::PLUFactorization(float **w,

// An NxN matrix which initially contains the matrix

// A. On output it contains the LU-factorization. int n,

// The order of the matrix A. int *ipivot,

// An integer array which will contain the permutation

// of the rows. int &iflag)

// A flag which will be -1 if the permutation ipivot

// is odd, 0 if Gaussian elimination fails, and 1 if

// ipivot an even permutation.

{

float *d = new float[n]; n--; // Adjust n for zero indexing.

iflag = 1;

for (int i = 0; i

ipivot[i] = i; // Initialize ipivot.

// Determine the absolute value of the maximum element

// of row i and store in d[i].

float rowmax = 0.0;

for (int j = 0; j

{

rowmax = std::max(rowmax, (float(fabs(double(w[i][j])))));

}

// If rowmax = 0.0, then row i is all zeros, failure.

if (rowmax

{

iflag = 0;

rowmax = 1.0;

}

d[i] = rowmax;

}

if (n > 0) // Elimination Steps.

{

for (int k = 0; k

{

// Find the equation which has the largest

// scaled coefficient of x[k].

float colmax = (float(fabs(double(w[k][k])))) / d[k]; int istar = k;

for (int i = k + 1; i

{

float awikod = (float(fabs(double(w[i][k])))) / d[i]; if (awikod > colmax)

{

colmax = awikod; istar = i;

}

}

if (colmax

{

// If colmax = 0.0, then matrix is singlar, failure.

iflag = 0;

}

else

{

if (istar > k)

{

// If the equation is not the K th, interchange rows

// and corresponding data.

iflag = -iflag;

int i = ipivot[istar];

ipivot[istar] = ipivot[k];

ipivot[k] = i;

float temp = d[istar];

d[istar] = d[k];

d[k] = temp;

for (int j = 0; j

{

temp = w[istar][j];

w[istar][j] = w[k][j];

w[k][j] = temp;

}

}

// Eliminate x[k] from the new equations k+1, ... , n.

for (int i = k+1; i

{

float ratio = w[i][k] / w[k][k];

for (int j = k+1; j

{

w[i][j] -= ratio * w[k][j];

}

w[i][k] = ratio;

}

}

}

if ((float(fabs(double(w[n][n]))))

{

// If the last element is zero, matrix is singular, failure. iflag = 0;

}

}

delete[] d;

}

// This function performs the solving phase of solving Ax = b

// using the triangular factorization of A.

//

// NOTE: This is a modified version of the SUBST subroutine

// on page 164 of Elementary Numerical Analysis by

// S. D. Conte and Carl de Boor written in C++.

void GESolver::SolvingPhase(float **w, // An NxN matrix containing the triangular factorization.

int *ipivot, // An integer array containing the permutation of the rows.

float *b, // The right-hand-side of Ax = b, an N-vector.

int n, // The order of the matrix A. float *x)

// The output, solution to Ax = b.

{

n--; // Adjust n for zero indexing.

if (n == 0) { x[0] = b[0] / w[0][0];

}

else

{

int ip = ipivot[0];

// Forward Substitution Step. x[0] = b[ip];

for (int i = 1; i

{

float sum = 0.0;

for (int j = 0; j

{

sum += w[i][j] * x[j];

}

ip = ipivot[i]; x[i] = b[ip] - sum;

}

// Back Substitution Step.

x[n] /= w[n][n];

for (int i = n - 1; i >= 0; i--)

{

float sum = 0.0;

for (int j = i + 1; j

{

sum += w[i][j] * x[j];

}

x[i] = (x[i] - sum) / w[i][i];

}

}

}

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

system1.dat

3

1.998 2.001 4.133 4.136

-4.125 3.333 -3.234 4.224

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

system2.dat

51.000000 0.5000000 0.3333333 0.2500000 0.2000000 1.000000

0.500000 0.3333333 0.2500000 0.2000000 0.1666667 0.000000

0.333333 0.2500000 0.2000000 0.1666667 0.1428571 0.000000

0.2500000 0.2000000 0.1666667 0.1428571 0.1250000 0.000000

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

system3.dat

36.0 2.0 2.0 -2.0

2.0 0.6666667 0.3333333 1.0

1.0 2.0 -1.0 0.0

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

system4.dat

4

1.0 2.0 3.0 4.0 1.0

1.0 2.0 3.0 4.0 1.0

3.0 3.0 3.0 4.0 1.0

4.0 4.0 4.0 4.0 1.0

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>.

Please write your program along with the outputs of each system run.

Write a program to use the C++ class GESolver that contains the routines PLUFactorization and SolvingPhase (see Canvas for code) to solve the linear system Ax = b to obtain the approximate solution . (If this is possible. A may be noninvertible. Remember to check before calling SolvingPhase) You should also compute det(A) using the triangular factorization produced by PLUFactorization. The first record of each system contains the order n of the system. The next n records each contain a row of A and the corresponding element of b. OUTPUT: For each system, write the following at the top of a new page. A.b: ain 011 021 a12 022 bi b2 a2n : : : anl an2 ann bn W: Win W11 W21 W12 W22 W2n : : Wnl Wn2 Wnn det(A): det A X: i 2 n NOTE: Use a %18.9f format for all output. Write a program to use the C++ class GESolver that contains the routines PLUFactorization and SolvingPhase (see Canvas for code) to solve the linear system Ax = b to obtain the approximate solution . (If this is possible. A may be noninvertible. Remember to check before calling SolvingPhase) You should also compute det(A) using the triangular factorization produced by PLUFactorization. The first record of each system contains the order n of the system. The next n records each contain a row of A and the corresponding element of b. OUTPUT: For each system, write the following at the top of a new page. A.b: ain 011 021 a12 022 bi b2 a2n : : : anl an2 ann bn W: Win W11 W21 W12 W22 W2n : : Wnl Wn2 Wnn det(A): det A X: i 2 n NOTE: Use a %18.9f format for all output