Les paso unas rutinas que creé para el cálculo matricial en C++.
Código:
struct stMatrices
{
int rx;
int ry;
vector<float> Mat;
public:
//FUNCIONES
stMatrices(void); //Constructor sobrecargado
stMatrices(int x, int y); //Constructor
~stMatrices(void); //Destructor
float get(int x,int y); //Retorna el valor de la posición (x,y)
void set(int x, int y, float val); //Asigna el valor val a la posición (x,y)
void redim (int rx, int ry); //Redimensiona la matriz
int ubound(int col); //Retorna el máximo valor de x o y
};
//FUNCIONES CON VECTORES
void MMuestre(vector<float> Mat,int R1,int R2);
vector<float> MInversa(vector<float> Mat, int R1);
vector<float> MMult(vector<float> MatA, vector<float> MatB, int RAF, int RAC, int RBC);
vector<float> MTransp(vector<float> Mat, int RF, int RC);
vector<float> MkMult(vector<float> Mat, int RF, int RC,float factK);
vector<float> MSuma(vector<float> MatA, vector<float> MatB,int RF, int RC,float factK);
//FUNCIONES CON ESTRUCTURAS
void MMuestre(stMatrices Mat);
stMatrices MInversa(stMatrices Mat);
stMatrices MMult(stMatrices MatA, stMatrices MatB);
stMatrices MTransp(stMatrices Mat);
stMatrices MkMult(stMatrices Mat,float factK);
stMatrices MSuma(stMatrices MatA, stMatrices MatB,float factK);
Código:
//[O]
//Herramientas matriciales
#include <iostream>
#include <string>
#include <cstdio>
#include <vector>
#include "mat.h"
using namespace std;
//*********************************
// FUNCIONES DEL OBJETO
//*********************************
stMatrices::stMatrices(void) //Constructor sobrecargado
{
redim(0,0);
}
stMatrices::stMatrices(int x, int y) //Constructor
{
redim(x,y);
}
stMatrices::~stMatrices(void) //Destructor
{
}
float stMatrices::get(int x,int y)
{
if (x <= rx)
if (y <= ry)
return Mat[x*ry+y];
}
void stMatrices::set(int x, int y, float val)
{
if (x <= rx)
if (y <= ry)
Mat[x*ry+y] = val;
}
void stMatrices::redim (int x, int y)
{
rx = x;
ry = y;
Mat.resize(rx*ry);
}
int stMatrices::ubound(int col)
{
if (col == 1)
return rx;
else if (col == 2)
return ry;
else
return 0;
}
//************************************************
// FUNCIONES MATRICIALES CON VECTORES
//************************************************
void MMuestre (vector<float> Mat,int R1,int R2)
{
for(int i = 0;i < R1; i++)
{
for (int j = 0; j < R2; j++)
printf("%f ",Mat[i*R2+j]);
printf("\n");
}
}
vector<float> MInversa(vector<float> Mat, int R1)
{
//Inversión de la Matriz 'Mat', de rango R1, por el método de Gauss-Jordan
//int R1 = El rango de la Matriz
vector<float> MatM;
int R12 = R1*2;
MatM.resize(R1*R12);
//Creo la Matriz complementaria de la forma [M|I], donde M = Mat e I = Identidad
//1º la parte de MatM = [Mat|...]
for (int i = 0; i < R1; i++)
for (int j = 0; j < R1; j++)
MatM[i * R12 + j] = Mat[i * R1 + j];
//2º La parte de MatM = [...|I]
for (int i = 0; i < R1; i++)
for (int j = 0; j < R1; j++)
if (i == j)
MatM[i * R12 + R1 + j] = 1;
else
MatM[i * R12 + R1 + j] = 0;
//INVERSIÓN PROPIAMENTE DICHA
float Pibot=0; //Pibot
int iP = 0 ; //Índice del Pibot
float factK = 0; //Factor K
for (int i = 0; i < R1; i++)
{
iP = i; //Índice del Pibot
Pibot = MatM[iP * R12 + iP]; //Obtiene el Pibot
for (int k = iP; k < (R1 * 2 + 1); k++) //Reduce la FILA del pibot
MatM[i * R12 + k] = MatM[i * R12 + k] / Pibot;
//Fila n = Fila n -k.Fila n
//Obtiene el 0 delante del pibot
for (int h = 0; h < R1; h++)
if (h != iP) //La fila del pibot no se toca
{
factK = MatM[h * R12 + iP]; //Obtiene el factor k
for (int k = 0; k < R1 * 2 + 1; k++)
MatM[h*R12+k] = MatM[h*R12+k] - factK * MatM[iP*R12+k];
}
}
//FINAL PRESENTACIÓN DE LOS DATOS
vector<float> MatT;
MatT.resize(R1*R1);
for (int i = 0; i < R1; i++)
for (int j = 0; j < R1; j++)
MatT[i*R1+j] = MatM[i * R12 + j + R1];
return MatT;
}
vector<float> MMult(vector<float> MatA, vector<float> MatB, int RAF, int RAC,
int RBC){
//Multiplica MatA x MatB, siendo
//RAF = Cantidad de filas de MatA
//RAC = Cantidad de columnas de MatA
//RBC = Cantidad de columnas de MatB
//La cantidad de filas de MatB es la cantidad de columnas de MatA
//La respuesta tiene la cantidad de filas de MatA y de columnas de MatB
int RBF = RAC;
vector<float> MatR;
float val=0;
MatR.resize(RAF*RBC);
for (int i = 0; i < RAF; i++)
for (int j = 0; j < RBC; j++)
{
for (int k = 0; k < RAC; k++)
val += MatA[i * RAC + k] * MatB[k * RBC + j];
MatR[i * RBC + j] = val;
val = 0;
}
return MatR;
}
vector<float> MTransp(vector<float> Mat, int RF, int RC)
{
//Transpone la matriz MatA
vector<float> MatR;
MatR.resize(RF*RC);
for (int i = 0; i < RF; i++)
for (int j = 0; j < RC; j++)
MatR[j * RF + i] = Mat[i * RC + j];
return MatR;
}
vector<float> MkMult(vector<float> Mat, int RF, int RC,float factK)
{
//Multiplica una matriz por una constante
vector<float> MatR;
MatR.resize(RF*RC);
for (int i = 0; i < RF; i++)
for (int j = 0; j < RC; j++)
MatR[i*RC+j]= Mat[i*RC+j] * factK;
return MatR;
}
vector<float> MSuma(vector<float> MatA, vector<float> MatB,int RF, int RC,float
factK)
{
//Multiplica suma 2 matrices multiplicando a la segunda por un factor
//MatR = MatA + k.MatB
//Ej. para una resta de matrices, factK = -1
vector<float> MatR;
MatR.resize(RF*RC);
for (int i = 0; i < RF; i++)
for (int j = 0; j < RC; j++)
MatR[i*RC+j]= MatA[i*RC+j] + factK * MatB[i*RC+j];
return MatR;
}
//****************************************************
// FUNCIONES MATRICIALES CON ESTRUCTURAS
//****************************************************
void MMuestre (stMatrices Mat)
{
for(int i = 0;i < Mat.ubound(1) ; i++)
{
for (int j = 0; j < Mat.ubound(2); j++)
printf("%f ",Mat.get(i,j));
printf("\n");
}
}
stMatrices MInversa(stMatrices Mat)
{
//Inversión de la Matriz 'Mat', de rango R1, por el método de Gauss-Jordan
int R1 = Mat.ubound(1);
int R12 = R1*2;
stMatrices MatM(R1,R12);
//MatM.resize(R1,R12);
//Creo la Matriz complementaria de la forma [M|I], donde M = Mat e I = Identidad
//1º la parte de MatM = [Mat|...]
for (int i = 0; i < R1; i++)
for (int j = 0; j < R1; j++)
MatM.set(i,j, Mat.get(i,j));
//2º La parte de MatM = [...|I]
for (int i = 0; i < R1; i++)
for (int j = 0; j < R1; j++)
if (i == j)
MatM.set(i,j+R1,1);
else
MatM.set(i,j+R1,0);
//INVERSIÓN PROPIAMENTE DICHA
float Pibot=0; //Pibot
int iP = 0 ; //Índice del Pibot
float factK = 0; //Factor K
for (int i = 0; i < R1; i++)
{
iP = i; //Índice del Pibot
Pibot = MatM.get(iP,iP); //Obtiene el Pibot
for (int k = iP; k < (R1 * 2 + 1); k++) //Reduce la FILA del pibot
MatM.set(i,k, MatM.get(i,k) / Pibot);
//Fila n = Fila n -k.Fila n
//Obtiene el 0 delante del pibot
for (int h = 0; h < R1; h++)
if (h != iP) //La fila del pibot no se toca
{
factK = MatM.get(h,iP); //Obtiene el factor k
for (int k = 0; k < R1 * 2 + 1; k++)
MatM.set(h,k, MatM.get(h,k) - factK * MatM.get(iP,k));
}
}
//FINAL PRESENTACIÓN DE LOS DATOS
stMatrices MatT(R1,R1);
for (int i = 0; i < R1; i++)
for (int j = 0; j < R1; j++)
MatT.set(i,j, MatM.get(i,j+R1));
return MatT;
}
stMatrices MMult(stMatrices MatA, stMatrices MatB)
{
//Multiplica MatA x MatB, siendo
//RAF = Cantidad de filas de MatA
//RAC = Cantidad de columnas de MatA
//RBC = Cantidad de columnas de MatB
//La cantidad de filas de MatB es la cantidad de columnas de MatA
//La respuesta tiene la cantidad de filas de MatA y de columnas de MatB
int RAF = MatA.ubound(1);
int RAC = MatA.ubound(2);
int RBF = MatB.ubound(1);
int RBC = MatB.ubound(2);
// if (RBF == RAC)
// return;
stMatrices MatR(RAF,RBC);
float val = 0;
for (int i = 0; i < RAF; i++)
for (int j = 0; j < RBC; j++)
{
for (int k = 0; k < RAC; k++)
val += MatA.get(i,k) * MatB.get(k,j);
MatR.set(i,j,val);
val = 0;
}
return MatR;
}
stMatrices MTransp(stMatrices Mat)
{
//Transpone la matriz Mat
int RF = Mat.ubound(1);
int RC = Mat.ubound(2);
stMatrices MatR(RC,RF);
for (int i = 0; i < RF; i++)
for (int j = 0; j < RC; j++)
MatR.set(j,i, Mat.get(i,j));
return MatR;
}
stMatrices MkMult(stMatrices Mat,float factK)
{
//Multiplica una matriz por una constante
int RF = Mat.ubound(1);
int RC = Mat.ubound(2);
stMatrices MatR(RF,RC);
for (int i = 0; i < RF; i++)
for (int j = 0; j < RC; j++)
MatR.set(i,j, Mat.get(i,j)*factK);
return MatR;
}
stMatrices MSuma(stMatrices MatA, stMatrices MatB, float factK)
{
//Multiplica suma 2 matrices multiplicando a la segunda por un factor
//MatR = MatA + k.MatB
//Ej. para una resta de matrices, factK = -1
int RF = MatA.ubound(1);
int RC = MatA.ubound(2);
stMatrices MatR(RF,RC);
for (int i = 0; i < RF; i++)
for (int j = 0; j < RC; j++)
MatR.set(i,j, MatA.get(i,j) + factK * MatB.get(i,j));
return MatR;
}