matrix problem
How can I define this algorithm to calculate the matrix that has uneven nr of rows and columns?
//==============================================================================
//
// Linear System Solution by Gauss method
//
//
//
//==============================================================================
#include <iostream>
#include <stdio.h>
#include <windows.h>
#include <cmath>
//==============================================================================
void VectorPrint(int nDim, double* pfVect)
{
int i;
printf("------------------------------------------------------------------\n");
for(i=0; i<nDim; i++)
{
printf("%9.2lf ", pfVect[i]);
}
printf("\n");
}
//==============================================================================
void MatrixPrint(int nDim, double* pfMatr)
{
int i,j;
printf("------------------------------------------------------------------\n");
for(i=0; i<nDim; i++)
{
for(j=0; j<nDim; j++)
{
printf("%9.2lf ", pfMatr[nDim*i + j]);
}
printf("\n");
}
}
//==============================================================================
// return 1 if system not solving
// nDim - system dimension
// pfMatr - matrix with coefficients
// pfVect - vector with free members
// pfSolution - vector with system solution
// pfMatr becames trianglular after function call
// pfVect changes after function call
//
//
//
//==============================================================================
int LinearEquationsSolving(int nDim, double* pfMatr, double* pfVect, double* pfSolution)
{
double fMaxElem;
double fAcc;
int i , j, k, m;
for(k=0; k<(nDim-1); k++) // base row of matrix
{
// search of line with max element
fMaxElem= fabs(pfMatr[k*nDim + k]);
m = k;
for(i=k+1; i<nDim; i++)
{
if(fMaxElem < fabs(pfMatr[i*nDim + k]) )
{
fMaxElem = pfMatr[i*nDim + k];
m = i;
}
}
// permutation of base line (index k) and max element line(index m)
if(m != k)
{
for(i=k; i<nDim; i++)
{
fAcc = pfMatr[k*nDim + i];
pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
pfMatr[m*nDim + i] = fAcc;
}
fAcc = pfVect[k];
pfVect[k] = pfVect[m];
pfVect[m] = fAcc;
}
if( pfMatr[k*nDim + k] == 0.) return 1; // needs improvement !!!
// triangulation of matrix with coefficients
for(j=(k+1); j<nDim; j++) // current row of matrix
{
fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
for(i=k; i<nDim; i++)
{
pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
}
pfVect[j] = pfVect[j] + fAcc*pfVect[k]; // free member recalculation
}
}
for(k=(nDim-1); k>=0; k--)
{
pfSolution[k] = pfVect[k];
for(i=(k+1); i<nDim; i++)
{
pfSolution[k] -= (pfMatr[k*nDim + i]*pfSolution[i]);
}
pfSolution[k] = pfSolution[k] / pfMatr[k*nDim + k];
}
return 0;
}
//==============================================================================
// testing of function
//==============================================================================
#define MATRIX_DIMENSION 11
int main(int nArgs, char** pArgs)
{
int nDim = MATRIX_DIMENSION;
double fMatr[MATRIX_DIMENSION*MATRIX_DIMENSION] =
{
9.00, 12.00, 11.00, 15.00, 27.00, 46.00, 34.00, 50.00, 64.00, 2.71, 35.00,
8.00, 10.45, 9.11, 12.1, 21.76, 39.42, 30.19, 39.00, 36.83, 0.00, 0.00,
4.00, 2.00, 2.00, 4.00, 8.00, 1.50, 1.50, 10.00, 39.00, 0.00, 0.00,
1.50, 1.50, 2.50, 5.50, 4.50, 6.00, 0.00, 30.00, 15.50, 0.00, 0.00,
33.80, 30.80, 27.10, 14.50, 23.48, 22.30, 13.00, 23.50, 29.40, 88.00, 18.14,
2.00, 2.00, 5.00, 10.50, 8.50, 10.50, 23.00, 0.00, 0.00, 0.00, 4.00,
0.26, 0.34, 0.41, 0.63, 0.60, 2.89, 1.21, 2.42, 5.05, 0.00, 0.95,
0.20, 0.28, 0.33, 0.46, 0.36, 2.56, 1.05, 1.90, 4.32, 0.00, 0.00,
};
double fVec[MATRIX_DIMENSION] = {};
double fSolution[MATRIX_DIMENSION];
int res;
int i;
res = LinearEquationsSolving(nDim, fMatr, fVec, fSolution); // !!!
if(res)
{
printf("No solution!\n");
return 1;
}
else
{
printf("Solution:\n");
VectorPrint(nDim, fSolution);
}
return 0;
}
//==============================================================================
//
// Linear System Solution by Gauss method
//
//
//
//==============================================================================
#include <iostream>
#include <stdio.h>
#include <windows.h>
#include <cmath>
//==============================================================================
void VectorPrint(int nDim, double* pfVect)
{
int i;
printf("------------------------------------------------------------------\n");
for(i=0; i<nDim; i++)
{
printf("%9.2lf ", pfVect[i]);
}
printf("\n");
}
//==============================================================================
void MatrixPrint(int nDim, double* pfMatr)
{
int i,j;
printf("------------------------------------------------------------------\n");
for(i=0; i<nDim; i++)
{
for(j=0; j<nDim; j++)
{
printf("%9.2lf ", pfMatr[nDim*i + j]);
}
printf("\n");
}
}
//==============================================================================
// return 1 if system not solving
// nDim - system dimension
// pfMatr - matrix with coefficients
// pfVect - vector with free members
// pfSolution - vector with system solution
// pfMatr becames trianglular after function call
// pfVect changes after function call
//
//
//
//==============================================================================
int LinearEquationsSolving(int nDim, double* pfMatr, double* pfVect, double* pfSolution)
{
double fMaxElem;
double fAcc;
int i , j, k, m;
for(k=0; k<(nDim-1); k++) // base row of matrix
{
// search of line with max element
fMaxElem= fabs(pfMatr[k*nDim + k]);
m = k;
for(i=k+1; i<nDim; i++)
{
if(fMaxElem < fabs(pfMatr[i*nDim + k]) )
{
fMaxElem = pfMatr[i*nDim + k];
m = i;
}
}
// permutation of base line (index k) and max element line(index m)
if(m != k)
{
for(i=k; i<nDim; i++)
{
fAcc = pfMatr[k*nDim + i];
pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
pfMatr[m*nDim + i] = fAcc;
}
fAcc = pfVect[k];
pfVect[k] = pfVect[m];
pfVect[m] = fAcc;
}
if( pfMatr[k*nDim + k] == 0.) return 1; // needs improvement !!!
// triangulation of matrix with coefficients
for(j=(k+1); j<nDim; j++) // current row of matrix
{
fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
for(i=k; i<nDim; i++)
{
pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
}
pfVect[j] = pfVect[j] + fAcc*pfVect[k]; // free member recalculation
}
}
for(k=(nDim-1); k>=0; k--)
{
pfSolution[k] = pfVect[k];
for(i=(k+1); i<nDim; i++)
{
pfSolution[k] -= (pfMatr[k*nDim + i]*pfSolution[i]);
}
pfSolution[k] = pfSolution[k] / pfMatr[k*nDim + k];
}
return 0;
}
//==============================================================================
// testing of function
//==============================================================================
#define MATRIX_DIMENSION 11
int main(int nArgs, char** pArgs)
{
int nDim = MATRIX_DIMENSION;
double fMatr[MATRIX_DIMENSION*MATRIX_DIMENSION] =
{
9.00, 12.00, 11.00, 15.00, 27.00, 46.00, 34.00, 50.00, 64.00, 2.71, 35.00,
8.00, 10.45, 9.11, 12.1, 21.76, 39.42, 30.19, 39.00, 36.83, 0.00, 0.00,
4.00, 2.00, 2.00, 4.00, 8.00, 1.50, 1.50, 10.00, 39.00, 0.00, 0.00,
1.50, 1.50, 2.50, 5.50, 4.50, 6.00, 0.00, 30.00, 15.50, 0.00, 0.00,
33.80, 30.80, 27.10, 14.50, 23.48, 22.30, 13.00, 23.50, 29.40, 88.00, 18.14,
2.00, 2.00, 5.00, 10.50, 8.50, 10.50, 23.00, 0.00, 0.00, 0.00, 4.00,
0.26, 0.34, 0.41, 0.63, 0.60, 2.89, 1.21, 2.42, 5.05, 0.00, 0.95,
0.20, 0.28, 0.33, 0.46, 0.36, 2.56, 1.05, 1.90, 4.32, 0.00, 0.00,
};
double fVec[MATRIX_DIMENSION] = {};
double fSolution[MATRIX_DIMENSION];
int res;
int i;
res = LinearEquationsSolving(nDim, fMatr, fVec, fSolution); // !!!
if(res)
{
printf("No solution!\n");
return 1;
}
else
{
printf("Solution:\n");
VectorPrint(nDim, fSolution);
}
return 0;
}
Not 100% sure what the question is, but I'll have a go:
First some nitpicking:
* Very difficult to read the code if it isn't in code tags.
* Why include <iostream> if you are only going to use printf?
* The values in fMatr don't fill the memory allocated for the matrix (8 x 11 instead of 11 x 11).
* No values for fVec were given (needs 11 values).
1) If the number of rows/columns is an odd number ("uneven"), Gaussian elimination still works (where the inverse of the matrix exists).
or, what I think you actually meant...
2) If the number of rows is different from the number of columns (ie, the matrix is not square) then consider the general problem:
Solve A x = b, for the unknown x, where A is an [m x n] matrix, x is a n-vector, and b is a m-vector.
There are 3 possibilities:
i) m == n: the number of equations is equal to the number of unknowns, and a unique solution may exist (depends if the matrix is nonsingular). Gaussian elimination or LU decomposition is appropriate. Example: Solving a set of linear simultaneous equations.
ii) m < n: the number of equations is less than the number of unknowns, in which case the system is underdetermined. There exists no unique solution, and a space of solutions exist, all which satisfy the original equation. Gaussian elimination will not solve this problem. QR decomposition or the Singular Value Decomposition may be used to find the range/null space of the problem, but again, no unique solution necessarily exists. Example: Describing the remaining degrees of freedom for a mechanical system, where constraints have been imposed.
iii) m > n: the number of equations is more than the number of unknowns, in which case the system is overdetermined. In general no solution exists, but is it common to find a solution that minimises the least squared error ||A x - b||_2^2. Gaussian elimination will not help here also. QR decomposition is the usual way to "solve" this problem. Example: Finding the coefficients of a polynomial given a set of data.
Cases ii and iii usually represent vastly different physical problems to case i, so I suggest you think carefully about the question you are trying to ask.
If you don't know what I'm talking about, I would suggest reading:
Matrix Computations 3rd Edition, Golub and van Loan, 1996.
This book describes how to solve all the various problems I listed above, and gives algorithms in a MATLAB-like pseudocode, which is relatively easy to translate into C/C++.
First some nitpicking:
* Very difficult to read the code if it isn't in code tags.
* Why include <iostream> if you are only going to use printf?
* The values in fMatr don't fill the memory allocated for the matrix (8 x 11 instead of 11 x 11).
* No values for fVec were given (needs 11 values).
1) If the number of rows/columns is an odd number ("uneven"), Gaussian elimination still works (where the inverse of the matrix exists).
or, what I think you actually meant...
2) If the number of rows is different from the number of columns (ie, the matrix is not square) then consider the general problem:
Solve A x = b, for the unknown x, where A is an [m x n] matrix, x is a n-vector, and b is a m-vector.
There are 3 possibilities:
i) m == n: the number of equations is equal to the number of unknowns, and a unique solution may exist (depends if the matrix is nonsingular). Gaussian elimination or LU decomposition is appropriate. Example: Solving a set of linear simultaneous equations.
ii) m < n: the number of equations is less than the number of unknowns, in which case the system is underdetermined. There exists no unique solution, and a space of solutions exist, all which satisfy the original equation. Gaussian elimination will not solve this problem. QR decomposition or the Singular Value Decomposition may be used to find the range/null space of the problem, but again, no unique solution necessarily exists. Example: Describing the remaining degrees of freedom for a mechanical system, where constraints have been imposed.
iii) m > n: the number of equations is more than the number of unknowns, in which case the system is overdetermined. In general no solution exists, but is it common to find a solution that minimises the least squared error ||A x - b||_2^2. Gaussian elimination will not help here also. QR decomposition is the usual way to "solve" this problem. Example: Finding the coefficients of a polynomial given a set of data.
Cases ii and iii usually represent vastly different physical problems to case i, so I suggest you think carefully about the question you are trying to ask.
If you don't know what I'm talking about, I would suggest reading:
Matrix Computations 3rd Edition, Golub and van Loan, 1996.
This book describes how to solve all the various problems I listed above, and gives algorithms in a MATLAB-like pseudocode, which is relatively easy to translate into C/C++.
Topic archived. No new replies allowed.