sacar - Código inverso de matriz 3x3 simple(C++)
matriz inversa (13)
¿Por qué no intentas codificarlo tú mismo? Tómalo como un desafío. :)
Para una matriz de 3 × 3
texto alternativo http://mathworld.wolfram.com/images/equations/MatrixInverse/NumberedEquation3.gif
la matriz inversa es
texto alternativo http://mathworld.wolfram.com/images/equations/MatrixInverse/NumberedEquation4.gif
Supongo que sabes cuál es el determinante de una matriz | A | es.
Imágenes (c) Wolfram | Alpha y mathworld.wolfram (06-11-09, 22.06)
¿Cuál es la forma más fácil de calcular una matriz inversa de 3x3?
Solo busco un pequeño fragmento de código que haga el truco para matrices no singulares, posiblemente usando la regla de Cramer. No necesita ser altamente optimizado. Preferiría simplicidad sobre velocidad. Prefiero no vincular en bibliotecas adicionales.
Acabo de crear una clase QMatrix. Utiliza el contenedor> vector integrado. QMatrix.h Utiliza el método de Jordan-Gauss para calcular el inverso de una matriz cuadrada.
Puedes usarlo de la siguiente manera:
#include "QMatrix.h"
#include <iostream>
int main(){
QMatrix<double> A(3,3,true);
QMatrix<double> Result = A.inverse()*A; //should give the idendity matrix
std::cout<<A.inverse()<<std::endl;
std::cout<<Result<<std::endl; // for checking
return 0;
}
La función inversa se implementa de la siguiente manera:
Dada una clase con los siguientes campos:
template<class T> class QMatrix{
public:
int rows, cols;
std::vector<std::vector<T> > A;
la función inversa ():
template<class T>
QMatrix<T> QMatrix<T>:: inverse(){
Identity<T> Id(rows); //the Identity Matrix as a subclass of QMatrix.
QMatrix<T> Result = *this; // making a copy and transforming it to the Identity matrix
T epsilon = 0.000001;
for(int i=0;i<rows;++i){
//check if Result(i,i)==0, if true, switch the row with another
for(int j=i;j<rows;++j){
if(std::abs(Result(j,j))<epsilon) { //uses Overloading()(int int) to extract element from Result Matrix
Result.replace_rows(i,j+1); //switches rows i with j+1
}
else break;
}
// main part, making a triangular matrix
Id(i)=Id(i)*(1.0/Result(i,i));
Result(i)=Result(i)*(1.0/Result(i,i)); // Using overloading ()(int) to get a row form the matrix
for(int j=i+1;j<rows;++j){
T temp = Result(j,i);
Result(j) = Result(j) - Result(i)*temp;
Id(j) = Id(j) - Id(i)*temp; //doing the same operations to the identity matrix
Result(j,i)=0; //not necessary, but looks nicer than 10^-15
}
}
// solving a triangular matrix
for(int i=rows-1;i>0;--i){
for(int j=i-1;j>=0;--j){
T temp = Result(j,i);
Id(j) = Id(j) - temp*Id(i);
Result(j)=Result(j)-temp*Result(i);
}
}
return Id;
}
Aquí hay una versión de la respuesta de Batty, pero esta calcula la inversa correcta . La versión de Batty calcula la transposición de la inversa.
// computes the inverse of a matrix m
double det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) -
m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +
m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
double invdet = 1 / det;
Matrix33d minv; // inverse of matrix m
minv(0, 0) = (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) * invdet;
minv(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * invdet;
minv(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * invdet;
minv(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * invdet;
minv(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * invdet;
minv(1, 2) = (m(1, 0) * m(0, 2) - m(0, 0) * m(1, 2)) * invdet;
minv(2, 0) = (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1)) * invdet;
minv(2, 1) = (m(2, 0) * m(0, 1) - m(0, 0) * m(2, 1)) * invdet;
minv(2, 2) = (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1)) * invdet;
Con todo el respeto debido a nuestro cartel desconocido (yahoo), miro el código así y me muero un poco por dentro. La sopa de alfabeto es tan terriblemente difícil de depurar. Un error tipográfico en cualquier lugar puede arruinar tu día. Tristemente, este ejemplo particular careció de variables con guiones bajos. Es mucho más divertido cuando tenemos a_b-c_d * e_f-g_h. Especialmente cuando se usa una fuente donde _ y - tienen la misma longitud de píxel.
Tomando Suvesh Pratapa en su sugerencia, observo:
Given 3x3 matrix:
y0x0 y0x1 y0x2
y1x0 y1x1 y1x2
y2x0 y2x1 y2x2
Declared as double matrix [/*Y=*/3] [/*X=*/3];
(A) Cuando tomamos a un menor de una matriz de 3x3, tenemos 4 valores de interés. El índice X / Y más bajo siempre es 0 o 1. El índice X / Y más alto siempre es 1 o 2. ¡Siempre! Por lo tanto:
double determinantOfMinor( int theRowHeightY,
int theColumnWidthX,
const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
int x1 = theColumnWidthX == 0 ? 1 : 0; /* always either 0 or 1 */
int x2 = theColumnWidthX == 2 ? 1 : 2; /* always either 1 or 2 */
int y1 = theRowHeightY == 0 ? 1 : 0; /* always either 0 or 1 */
int y2 = theRowHeightY == 2 ? 1 : 2; /* always either 1 or 2 */
return ( theMatrix [y1] [x1] * theMatrix [y2] [x2] )
- ( theMatrix [y1] [x2] * theMatrix [y2] [x1] );
}
(B) El determinante es ahora: (¡fíjate en el signo menos!)
double determinant( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
return ( theMatrix [0] [0] * determinantOfMinor( 0, 0, theMatrix ) )
- ( theMatrix [0] [1] * determinantOfMinor( 0, 1, theMatrix ) )
+ ( theMatrix [0] [2] * determinantOfMinor( 0, 2, theMatrix ) );
}
(C) Y lo inverso es ahora:
bool inverse( const double theMatrix [/*Y=*/3] [/*X=*/3],
double theOutput [/*Y=*/3] [/*X=*/3] )
{
double det = determinant( theMatrix );
/* Arbitrary for now. This should be something nicer... */
if ( ABS(det) < 1e-2 )
{
memset( theOutput, 0, sizeof theOutput );
return false;
}
double oneOverDeterminant = 1.0 / det;
for ( int y = 0; y < 3; y ++ )
for ( int x = 0; x < 3; x ++ )
{
/* Rule is inverse = 1/det * minor of the TRANSPOSE matrix. *
* Note (y,x) becomes (x,y) INTENTIONALLY here! */
theOutput [y] [x]
= determinantOfMinor( x, y, theMatrix ) * oneOverDeterminant;
/* (y0,x1) (y1,x0) (y1,x2) and (y2,x1) all need to be negated. */
if( 1 == ((x + y) % 2) )
theOutput [y] [x] = - theOutput [y] [x];
}
return true;
}
Y redondearlo con un código de prueba de menor calidad:
void printMatrix( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
for ( int y = 0; y < 3; y ++ )
{
cout << "[ ";
for ( int x = 0; x < 3; x ++ )
cout << theMatrix [y] [x] << " ";
cout << "]" << endl;
}
cout << endl;
}
void matrixMultiply( const double theMatrixA [/*Y=*/3] [/*X=*/3],
const double theMatrixB [/*Y=*/3] [/*X=*/3],
double theOutput [/*Y=*/3] [/*X=*/3] )
{
for ( int y = 0; y < 3; y ++ )
for ( int x = 0; x < 3; x ++ )
{
theOutput [y] [x] = 0;
for ( int i = 0; i < 3; i ++ )
theOutput [y] [x] += theMatrixA [y] [i] * theMatrixB [i] [x];
}
}
int
main(int argc, char **argv)
{
if ( argc > 1 )
SRANDOM( atoi( argv[1] ) );
double m[3][3] = { { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
{ RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
{ RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) } };
double o[3][3], mm[3][3];
if ( argc <= 2 )
cout << fixed << setprecision(3);
printMatrix(m);
cout << endl << endl;
SHOW( determinant(m) );
cout << endl << endl;
BOUT( inverse(m, o) );
printMatrix(m);
printMatrix(o);
cout << endl << endl;
matrixMultiply (m, o, mm );
printMatrix(m);
printMatrix(o);
printMatrix(mm);
cout << endl << endl;
}
Idea tardía:
También es posible que desee detectar determinantes muy grandes ya que los errores de redondeo afectarán su precisión.
Este fragmento de código calcula la inversa transpuesta de la matriz A:
double determinant = +A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2)) -A(0,1)*(A(1,0)*A(2,2)-A(1,2)*A(2,0)) +A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0)); double invdet = 1/determinant; result(0,0) = (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet; result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet; result(2,0) = (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet; result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet; result(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet; result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet; result(0,2) = (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet; result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet; result(2,2) = (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
Aunque la pregunta estipuló matrices no singulares, es posible que desee comprobar si el determinante es igual a cero (o muy cercano a cero) y marcarlo de alguna manera para estar seguro.
No trates de hacer esto tú mismo si eres serio acerca de cómo solucionar los casos de borde. Entonces, aunque muchos métodos ingenuos / simples son teóricamente exactos, pueden tener un comportamiento numérico desagradable para matrices casi singulares. En particular, puede obtener errores de cancelación / redondeo que hacen que obtenga arbitrariamente malos resultados.
Una forma "correcta" es la eliminación gaussiana con pivote de fila y columna para que siempre se divida por el valor numérico restante más grande. (Esto también es estable para matrices NxN). Tenga en cuenta que la fila pivotando sola no detecta todos los casos malos.
Sin embargo, IMO implementar este derecho y rápido no vale la pena: use una biblioteca bien probada y hay un montón de encabezados únicos.
Seguí adelante y lo escribí en Python ya que creo que es mucho más legible que en c ++ para un problema como este. El orden de las funciones está en orden de operaciones para resolver esto a mano a través de este video . Simplemente importe esto y llame a "print_invert" en su matriz.
def print_invert (matrix):
i_matrix = invert_matrix (matrix)
for line in i_matrix:
print (line)
return
def invert_matrix (matrix):
determinant = str (determinant_of_3x3 (matrix))
cofactor = make_cofactor (matrix)
trans_matrix = transpose_cofactor (cofactor)
trans_matrix[:] = [[str (element) +''/''+ determinant for element in row] for row in trans_matrix]
return trans_matrix
def determinant_of_3x3 (matrix):
multiplication = 1
neg_multiplication = 1
total = 0
for start_column in range (3):
for row in range (3):
multiplication *= matrix[row][(start_column+row)%3]
neg_multiplication *= matrix[row][(start_column-row)%3]
total += multiplication - neg_multiplication
multiplication = neg_multiplication = 1
if total == 0:
total = 1
return total
def make_cofactor (matrix):
cofactor = [[0,0,0],[0,0,0],[0,0,0]]
matrix_2x2 = [[0,0],[0,0]]
# For each element in matrix...
for row in range (3):
for column in range (3):
# ...make the 2x2 matrix in this inner loop
matrix_2x2 = make_2x2_from_spot_in_3x3 (row, column, matrix)
cofactor[row][column] = determinant_of_2x2 (matrix_2x2)
return flip_signs (cofactor)
def make_2x2_from_spot_in_3x3 (row, column, matrix):
c_count = 0
r_count = 0
matrix_2x2 = [[0,0],[0,0]]
# ...make the 2x2 matrix in this inner loop
for inner_row in range (3):
for inner_column in range (3):
if row is not inner_row and inner_column is not column:
matrix_2x2[r_count % 2][c_count % 2] = matrix[inner_row][inner_column]
c_count += 1
if row is not inner_row:
r_count += 1
return matrix_2x2
def determinant_of_2x2 (matrix):
total = matrix[0][0] * matrix [1][1]
return total - (matrix [1][0] * matrix [0][1])
def flip_signs (cofactor):
sign_pos = True
# For each element in matrix...
for row in range (3):
for column in range (3):
if sign_pos:
sign_pos = False
else:
cofactor[row][column] *= -1
sign_pos = True
return cofactor
def transpose_cofactor (cofactor):
new_cofactor = [[0,0,0],[0,0,0],[0,0,0]]
for row in range (3):
for column in range (3):
new_cofactor[column][row] = cofactor[row][column]
return new_cofactor
También recomendaría Ilmbase, que es parte de OpenEXR . Es un buen conjunto de rutinas modulares de 2,3,4 vectores y matrices.
Un archivo de encabezado bastante agradable (creo) que contiene macros para la mayoría de las operaciones de matriz 2x2, 3x3 y 4x4 ha estado disponible con la mayoría de los toolkits de OpenGL. No como estándar, pero lo he visto en varios lugares.
Puede verificarlo aquí. Al final del mismo, encontrarás tanto el inverso de 2x2, 3x3 y 4x4.
# include <conio.h>
# include<iostream.h>
const int size = 9;
int main()
{
char ch;
do
{
clrscr();
int i, j, x, y, z, det, a[size], b[size];
cout << " **** MATRIX OF 3x3 ORDER ****"
<< endl
<< endl
<< endl;
for (i = 0; i <= size; i++)
a[i]=0;
for (i = 0; i < size; i++)
{
cout << "Enter "
<< i + 1
<< " element of matrix=";
cin >> a[i];
cout << endl
<<endl;
}
clrscr();
cout << "your entered matrix is "
<< endl
<<endl;
for (i = 0; i < size; i += 3)
cout << a[i]
<< " "
<< a[i+1]
<< " "
<< a[i+2]
<< endl
<<endl;
cout << "Transpose of given matrix is"
<< endl
<< endl;
for (i = 0; i < 3; i++)
cout << a[i]
<< " "
<< a[i+3]
<< " "
<< a[i+6]
<< endl
<< endl;
cout << "Determinent of given matrix is = ";
x = a[0] * (a[4] * a[8] -a [5] * a[7]);
y = a[1] * (a[3] * a[8] -a [5] * a[6]);
z = a[2] * (a[3] * a[7] -a [4] * a[6]);
det = x - y + z;
cout << det
<< endl
<< endl
<< endl
<< endl;
if (det == 0)
{
cout << "As Determinent=0 so it is singular matrix and its inverse cannot exist"
<< endl
<< endl;
goto quit;
}
b[0] = a[4] * a[8] - a[5] * a[7];
b[1] = a[5] * a[6] - a[3] * a[8];
b[2] = a[3] * a[7] - a[4] * a[6];
b[3] = a[2] * a[7] - a[1] * a[8];
b[4] = a[0] * a[8] - a[2] * a[6];
b[5] = a[1] * a[6] - a[0] * a[7];
b[6] = a[1] * a[5] - a[2] * a[4];
b[7] = a[2] * a[3] - a[0] * a[5];
b[8] = a[0] * a[4] - a[1] * a[3];
cout << "Adjoint of given matrix is"
<< endl
<< endl;
for (i = 0; i < 3; i++)
{
cout << b[i]
<< " "
<< b[i+3]
<< " "
<< b[i+6]
<< endl
<<endl;
}
cout << endl
<<endl;
cout << "Inverse of given matrix is "
<< endl
<< endl
<< endl;
for (i = 0; i < 3; i++)
{
cout << b[i]
<< "/"
<< det
<< " "
<< b[i+3]
<< "/"
<< det
<< " "
<< b[i+6]
<< "/"
<< det
<< endl
<<endl;
}
quit:
cout << endl
<< endl;
cout << "Do You want to continue this again press (y/yes,n/no)";
cin >> ch;
cout << endl
<< endl;
} /* end do */
while (ch == ''y'');
getch ();
return 0;
}
#include <iostream>
using namespace std;
int main()
{
double A11, A12, A13;
double A21, A22, A23;
double A31, A32, A33;
double B11, B12, B13;
double B21, B22, B23;
double B31, B32, B33;
cout << "Enter all number from left to right, from top to bottom, and press enter after every number: ";
cin >> A11;
cin >> A12;
cin >> A13;
cin >> A21;
cin >> A22;
cin >> A23;
cin >> A31;
cin >> A32;
cin >> A33;
B11 = 1 / ((A22 * A33) - (A23 * A32));
B12 = 1 / ((A13 * A32) - (A12 * A33));
B13 = 1 / ((A12 * A23) - (A13 * A22));
B21 = 1 / ((A23 * A31) - (A21 * A33));
B22 = 1 / ((A11 * A33) - (A13 * A31));
B23 = 1 / ((A13 * A21) - (A11 * A23));
B31 = 1 / ((A21 * A32) - (A22 * A31));
B32 = 1 / ((A12 * A31) - (A11 * A32));
B33 = 1 / ((A11 * A22) - (A12 * A21));
cout << B11 << "/t" << B12 << "/t" << B13 << endl;
cout << B21 << "/t" << B22 << "/t" << B23 << endl;
cout << B31 << "/t" << B32 << "/t" << B33 << endl;
return 0;
}
//Function for inverse of the input square matrix ''J'' of dimension ''dim'':
vector<vector<double > > inverseVec33(vector<vector<double > > J, int dim)
{
//Matrix of Minors
vector<vector<double > > invJ(dim,vector<double > (dim));
for(int i=0; i<dim; i++)
{
for(int j=0; j<dim; j++)
{
invJ[i][j] = (J[(i+1)%dim][(j+1)%dim]*J[(i+2)%dim][(j+2)%dim] -
J[(i+2)%dim][(j+1)%dim]*J[(i+1)%dim][(j+2)%dim]);
}
}
//determinant of the matrix:
double detJ = 0.0;
for(int j=0; j<dim; j++)
{ detJ += J[0][j]*invJ[0][j];}
//Inverse of the given matrix.
vector<vector<double > > invJT(dim,vector<double > (dim));
for(int i=0; i<dim; i++)
{
for(int j=0; j<dim; j++)
{
invJT[i][j] = invJ[j][i]/detJ;
}
}
return invJT;
}
void main()
{
//given matrix:
vector<vector<double > > Jac(3,vector<double > (3));
Jac[0][0] = 1; Jac[0][1] = 2; Jac[0][2] = 6;
Jac[1][0] = -3; Jac[1][1] = 4; Jac[1][2] = 3;
Jac[2][0] = 5; Jac[2][1] = 1; Jac[2][2] = -4;`
//Inverse of the matrix Jac:
vector<vector<double > > JacI(3,vector<double > (3));
//call function and store inverse of J as JacI:
JacI = inverseVec33(Jac,3);
}
//Title: Matrix Header File
//Writer: Say OL
//This is a beginner code not an expert one
//No responsibilty for any errors
//Use for your own risk
using namespace std;
int row,col,Row,Col;
double Coefficient;
//Input Matrix
void Input(double Matrix[9][9],int Row,int Col)
{
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
{
cout<<"e["<<row<<"]["<<col<<"]=";
cin>>Matrix[row][col];
}
}
//Output Matrix
void Output(double Matrix[9][9],int Row,int Col)
{
for(row=1;row<=Row;row++)
{
for(col=1;col<=Col;col++)
cout<<Matrix[row][col]<<"/t";
cout<<endl;
}
}
//Copy Pointer to Matrix
void CopyPointer(double (*Pointer)[9],double Matrix[9][9],int Row,int Col)
{
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
Matrix[row][col]=Pointer[row][col];
}
//Copy Matrix to Matrix
void CopyMatrix(double MatrixInput[9][9],double MatrixTarget[9][9],int Row,int Col)
{
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
MatrixTarget[row][col]=MatrixInput[row][col];
}
//Transpose of Matrix
double MatrixTran[9][9];
double (*(Transpose)(double MatrixInput[9][9],int Row,int Col))[9]
{
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
MatrixTran[col][row]=MatrixInput[row][col];
return MatrixTran;
}
//Matrix Addition
double MatrixAdd[9][9];
double (*(Addition)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
MatrixAdd[row][col]=MatrixA[row][col]+MatrixB[row][col];
return MatrixAdd;
}
//Matrix Subtraction
double MatrixSub[9][9];
double (*(Subtraction)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
MatrixSub[row][col]=MatrixA[row][col]-MatrixB[row][col];
return MatrixSub;
}
//Matrix Multiplication
int mRow,nCol,pCol,kcol;
double MatrixMult[9][9];
double (*(Multiplication)(double MatrixA[9][9],double MatrixB[9][9],int mRow,int nCol,int pCol))[9]
{
for(row=1;row<=mRow;row++)
for(col=1;col<=pCol;col++)
{
MatrixMult[row][col]=0.0;
for(kcol=1;kcol<=nCol;kcol++)
MatrixMult[row][col]+=MatrixA[row][kcol]*MatrixB[kcol][col];
}
return MatrixMult;
}
//Interchange Two Rows
double RowTemp[9][9];
double MatrixInter[9][9];
double (*(InterchangeRow)(double MatrixInput[9][9],int Row,int Col,int iRow,int jRow))[9]
{
CopyMatrix(MatrixInput,MatrixInter,Row,Col);
for(col=1;col<=Col;col++)
{
RowTemp[iRow][col]=MatrixInter[iRow][col];
MatrixInter[iRow][col]=MatrixInter[jRow][col];
MatrixInter[jRow][col]=RowTemp[iRow][col];
}
return MatrixInter;
}
//Pivote Downward
double MatrixDown[9][9];
double (*(PivoteDown)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
CopyMatrix(MatrixInput,MatrixDown,Row,Col);
Coefficient=MatrixDown[tRow][tCol];
if(Coefficient!=1.0)
for(col=1;col<=Col;col++)
MatrixDown[tRow][col]/=Coefficient;
if(tRow<Row)
for(row=tRow+1;row<=Row;row++)
{
Coefficient=MatrixDown[row][tCol];
for(col=1;col<=Col;col++)
MatrixDown[row][col]-=Coefficient*MatrixDown[tRow][col];
}
return MatrixDown;
}
//Pivote Upward
double MatrixUp[9][9];
double (*(PivoteUp)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
CopyMatrix(MatrixInput,MatrixUp,Row,Col);
Coefficient=MatrixUp[tRow][tCol];
if(Coefficient!=1.0)
for(col=1;col<=Col;col++)
MatrixUp[tRow][col]/=Coefficient;
if(tRow>1)
for(row=tRow-1;row>=1;row--)
{
Coefficient=MatrixUp[row][tCol];
for(col=1;col<=Col;col++)
MatrixUp[row][col]-=Coefficient*MatrixUp[tRow][col];
}
return MatrixUp;
}
//Pivote in Determinant
double MatrixPiv[9][9];
double (*(Pivote)(double MatrixInput[9][9],int Dim,int pTarget))[9]
{
CopyMatrix(MatrixInput,MatrixPiv,Dim,Dim);
for(row=pTarget+1;row<=Dim;row++)
{
Coefficient=MatrixPiv[row][pTarget]/MatrixPiv[pTarget][pTarget];
for(col=1;col<=Dim;col++)
{
MatrixPiv[row][col]-=Coefficient*MatrixPiv[pTarget][col];
}
}
return MatrixPiv;
}
//Determinant of Square Matrix
int dCounter,dRow;
double Det;
double MatrixDet[9][9];
double Determinant(double MatrixInput[9][9],int Dim)
{
CopyMatrix(MatrixInput,MatrixDet,Dim,Dim);
Det=1.0;
if(Dim>1)
{
for(dRow=1;dRow<Dim;dRow++)
{
dCounter=dRow;
while((MatrixDet[dRow][dRow]==0.0)&(dCounter<=Dim))
{
dCounter++;
Det*=-1.0;
CopyPointer(InterchangeRow(MatrixDet,Dim,Dim,dRow,dCounter),MatrixDet,Dim,Dim);
}
if(MatrixDet[dRow][dRow]==0)
{
Det=0.0;
break;
}
else
{
Det*=MatrixDet[dRow][dRow];
CopyPointer(Pivote(MatrixDet,Dim,dRow),MatrixDet,Dim,Dim);
}
}
Det*=MatrixDet[Dim][Dim];
}
else Det=MatrixDet[1][1];
return Det;
}
//Matrix Identity
double MatrixIdent[9][9];
double (*(Identity)(int Dim))[9]
{
for(row=1;row<=Dim;row++)
for(col=1;col<=Dim;col++)
if(row==col)
MatrixIdent[row][col]=1.0;
else
MatrixIdent[row][col]=0.0;
return MatrixIdent;
}
//Join Matrix to be Augmented Matrix
double MatrixJoin[9][9];
double (*(JoinMatrix)(double MatrixA[9][9],double MatrixB[9][9],int Row,int ColA,int ColB))[9]
{
Col=ColA+ColB;
for(row=1;row<=Row;row++)
for(col=1;col<=Col;col++)
if(col<=ColA)
MatrixJoin[row][col]=MatrixA[row][col];
else
MatrixJoin[row][col]=MatrixB[row][col-ColA];
return MatrixJoin;
}
//Inverse of Matrix
double (*Pointer)[9];
double IdentMatrix[9][9];
int Counter;
double MatrixAug[9][9];
double MatrixInv[9][9];
double (*(Inverse)(double MatrixInput[9][9],int Dim))[9]
{
Row=Dim;
Col=Dim+Dim;
Pointer=Identity(Dim);
CopyPointer(Pointer,IdentMatrix,Dim,Dim);
Pointer=JoinMatrix(MatrixInput,IdentMatrix,Dim,Dim,Dim);
CopyPointer(Pointer,MatrixAug,Row,Col);
for(Counter=1;Counter<=Dim;Counter++)
{
Pointer=PivoteDown(MatrixAug,Row,Col,Counter,Counter);
CopyPointer(Pointer,MatrixAug,Row,Col);
}
for(Counter=Dim;Counter>1;Counter--)
{
Pointer=PivoteUp(MatrixAug,Row,Col,Counter,Counter);
CopyPointer(Pointer,MatrixAug,Row,Col);
}
for(row=1;row<=Dim;row++)
for(col=1;col<=Dim;col++)
MatrixInv[row][col]=MatrixAug[row][col+Dim];
return MatrixInv;
}
//Gauss-Jordan Elemination
double MatrixGJ[9][9];
double VectorGJ[9][9];
double (*(GaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim))[9]
{
Row=Dim;
Col=Dim+1;
Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,1);
CopyPointer(Pointer,MatrixGJ,Row,Col);
for(Counter=1;Counter<=Dim;Counter++)
{
Pointer=PivoteDown(MatrixGJ,Row,Col,Counter,Counter);
CopyPointer(Pointer,MatrixGJ,Row,Col);
}
for(Counter=Dim;Counter>1;Counter--)
{
Pointer=PivoteUp(MatrixGJ,Row,Col,Counter,Counter);
CopyPointer(Pointer,MatrixGJ,Row,Col);
}
for(row=1;row<=Dim;row++)
for(col=1;col<=1;col++)
VectorGJ[row][col]=MatrixGJ[row][col+Dim];
return VectorGJ;
}
//Generalized Gauss-Jordan Elemination
double MatrixGGJ[9][9];
double VectorGGJ[9][9];
double (*(GeneralizedGaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim,int vCol))[9]
{
Row=Dim;
Col=Dim+vCol;
Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,vCol);
CopyPointer(Pointer,MatrixGGJ,Row,Col);
for(Counter=1;Counter<=Dim;Counter++)
{
Pointer=PivoteDown(MatrixGGJ,Row,Col,Counter,Counter);
CopyPointer(Pointer,MatrixGGJ,Row,Col);
}
for(Counter=Dim;Counter>1;Counter--)
{
Pointer=PivoteUp(MatrixGGJ,Row,Col,Counter,Counter);
CopyPointer(Pointer,MatrixGGJ,Row,Col);
}
for(row=1;row<=Row;row++)
for(col=1;col<=vCol;col++)
VectorGGJ[row][col]=MatrixGGJ[row][col+Dim];
return VectorGGJ;
}
//Matrix Sparse, Three Diagonal Non-Zero Elements
double MatrixSpa[9][9];
double (*(Sparse)(int Dimension,double FirstElement,double SecondElement,double ThirdElement))[9]
{
MatrixSpa[1][1]=SecondElement;
MatrixSpa[1][2]=ThirdElement;
MatrixSpa[Dimension][Dimension-1]=FirstElement;
MatrixSpa[Dimension][Dimension]=SecondElement;
for(int Counter=2;Counter<Dimension;Counter++)
{
MatrixSpa[Counter][Counter-1]=FirstElement;
MatrixSpa[Counter][Counter]=SecondElement;
MatrixSpa[Counter][Counter+1]=ThirdElement;
}
return MatrixSpa;
}
Copie y guarde el código anterior como Matrix.h y luego pruebe el siguiente código:
#include<iostream>
#include<conio.h>
#include"Matrix.h"
int Dim;
double Matrix[9][9];
int main()
{
cout<<"Enter your matrix dimension: ";
cin>>Dim;
Input(Matrix,Dim,Dim);
cout<<"Your matrix:"<<endl;
Output(Matrix,Dim,Dim);
cout<<"The inverse:"<<endl;
Output(Inverse(Matrix,Dim),Dim,Dim);
getch();
}