Implementación de Algoritmos de Álgebra Lineal en C++

Implementación de Algoritmos de Álgebra Lineal

1. Cálculo de Determinantes (Recursivo)

if(A.dim1()==0 || A.dim1()!=A.dim2()) return 0.;
if(A.dim1()==1) return A[0][0];
real determinante=0.;
for(int k=0;k<A.dim1();k++){
  Array2D< real > B(A.dim1()-1,A.dim1()-1);
  for(int i=0;i<B.dim1();i++){
    for(int j=0;j<B.dim1();j++){
      if(j<k) B[i][j]=A[i+1][j];
      else  B[i][j]=A[i+1][j+1];
    }
  }
  if(k%2==0) determinante+=A[0][k]*mn_determinante_recursivo(B);
  else determinante-=A[0][k]*mn_determinante_recursivo(B);
}
return determinante;

2. Resolución de Sistemas Lineales (Eliminación Gaussiana)

Array2D< real > A=A_original.copy();
Array1D< real > b=b_original.copy();
if(A.dim1()!=A.dim2() || A.dim1()!=b.dim() || b.dim()==0) return Array1D<real>();
for(int k=0;k<b.dim()-1;k++){
  int kmax=max_pos(A,k);
  if(A[kmax][k]==0.) return Array1D<real>();
  if(kmax!=k){
    for(int j=k;j<b.dim();j++){
      mn_pivotar(A[k][j],A[kmax][j]);
    }
    mn_pivotar(b[kmax],b[k]);
  }
  for(int j=k+1;j<b.dim();j++){
    real mul=-A[j][k]/A[k][k];
    A[j][k]=0.;
    for(int n=k+1;n<b.dim();n++) A[j][n]+=mul*A[k][n];
    b[j]+=mul*b[k];
  }
}
return mn_remonte(A,b);

3. Eliminación Gaussiana con Pivoteo

Array2D< real > A=A_original.copy();
Array1D< real > b=b_original.copy();
if(A.dim1()!=A.dim2() || A.dim1()!=b.dim() || b.dim()==0) return Array1D<real>();
Array1D< int >  piv(b.dim());
for(int k=0;k<b.dim();k++) piv[k]=k;
for(int k=0;k<b.dim()-1;k++){
  int kmax=max_pos(A,k,piv);
  if(A[piv[kmax]][k]==0.) return Array1D<real>();
  if(kmax!=k) mn_pivotar(piv[k],piv[kmax]);
  for(int j=k+1;j<b.dim();j++){
    real mul=-A[piv[j]][k]/A[piv[k]][k];
    A[piv[j]][k]=0.;
    for(int n=k+1;n<b.dim();n++) A[piv[j]][n]+=mul*A[piv[k]][n];
    b[piv[j]]+=mul*b[piv[k]];
  }
}
return mn_remonte(A,b,piv);

4. Inversión de Matrices

Array2D< real > A=A_original.copy();
if(A.dim1()!=A.dim2() || A.dim1()==0) return Array2D<real>();
Array2D< real > B(A.dim1(),A.dim2(),0.);
for(int k=0;k<A.dim1();k++) B[k][k]=1.;
for(int k=0;k<A.dim1();k++){
  int kmax=max_pos(A,k);
  if(A[kmax][k]==0.) return Array2D< real >();
  if(kmax!=k){
    for(int i=0;i<A.dim1();i++){
      mn_pivotar(A[k][i],A[kmax][i]);
      mn_pivotar(B[k][i],B[kmax][i]);
    }
  }
  for(int j=k+1;j<A.dim1();j++){
    real mul=-A[j][k]/A[k][k];
    A[j][k]=0.;
    for(int n=k+1;n<A.dim1();n++) A[j][n]+=mul*A[k][n];
    for(int n=0;n<A.dim1();n++) B[j][n]+=mul*B[k][n];
  }
}
return mn_remonte (A,B);

5. Determinante mediante Eliminación Gaussiana

if(M.dim1()==0 || M.dim1()!=M.dim2()) return 0.;
Array2D< real > A=M.copy();
real determinante;
int Npiv=0;
for(int k=0;k<A.dim1();k++){
  int kmax=max_pos(A,k);
  if(A[kmax][k]==0.) return 0.;
  if(kmax!=k){
    Npiv++;
    for(int i=0;i<A.dim1();i++) mn_pivotar(A[k][i],A[kmax][i]);
  }
  for(int j=k+1;j<A.dim1();j++){
    real mul=-A[j][k]/A[k][k];
    A[j][k]=0.;
    for(int n=k+1;n<A.dim1();n++) A[j][n]+=mul*A[k][n];
  }
}
determinante=1;
for(int k=0;k<A.dim1();k++) determinante*=A[k][k];
if(Npiv%2==1) return(-determinante);
return determinante;

6. Factorización de Cholesky

if(A.dim1()!=A.dim2() ) return( Array2D< real >());
int N=A.dim1();
Array2D< real > B(N,N);
for (int i=0;i<N;i++){
  real Sum = 0.;
  for (int k=0;k<i;k++) Sum = Sum + B[i][k]*B[i][k];
  if (A[i][i]< Sum) return( Array2D< real >());
  B[i][i]=sqrt(A[i][i]-Sum);
  for  (int j=i+1;j<N;j++){
    Sum = 0;
    if (B[i][i]==0) return( Array2D< real >());
    for (int k=0;k<i;k++) Sum = Sum + B[j][k]*B[i][k];
    B[j][i]=(A[j][i]-Sum)/B[i][i];
    B[i][j]= B[j][i];
  }
}
return( B );

7. Resolución mediante Cholesky

int N=A.dim1();
Array2D< real > CH = mn_cholesky_factorization (A);
if (CH.dim1()==0) return (Array1D< real >() );
Array1D< real > z=mn_descenso (CH,b);
if (z.dim()==0) return (Array1D< real >() );
Array1D< real > u = mn_remonte (CH,z);
if (u.dim()==0) return (Array1D< real >() );
return(u);

8. Factorización LU

int N=A.dim1();
L=Array2D< real >(N,N,0.);
U=Array2D< real >(N,N,0.);
if(N!=A.dim2() || N!=L.dim1() || N!=L.dim2() || N!=U.dim1() || N!=U.dim2()) return(-1);
for(int i=0;i<N;i++){
  L[i][i]=1;
  real Sum = 0.;
  for(int k=0;k<=(i-1);k++) Sum = Sum + L[i][k]*U[k][i];
  U[i][i]=A[i][i]-Sum;
  if(U[i][i]==0.) return(-1);
  for(int j=i+1;j<N;j++){
    Sum = 0;
    for (int k=0;k<=(i-1);k++) Sum = Sum + L[i][k]*U[k][j];
    U[i][j]=A[i][j]-Sum;
    Sum = 0;
    for (int k=0;k<=(i-1);k++) Sum = Sum + L[j][k]*U[k][i];
    L[j][i]=(A[j][i]-Sum)/U[i][i];
  }
}
return(0);

9. Resolución mediante LU

Array1D< real > c=mn_descenso(L,b);
return mn_remonte(U,c);

10. Descomposición de Crout

if(l.dim()!=a.dim() || l.dim()!=(m.dim()+1) || u.dim()!=m.dim() || u.dim()!=b.dim() || u.dim()!=c.dim()) return(-1);
l[0]=a[0];
if(l[0]==0) return(-2);
u[0]=b[0]/l[0];
for(int i=1;i<(a.dim()-1);i++){
  m[i-1]=c[i-1];
  l[i]=a[i]-m[i-1]*u[i-1];
  if(l[i]==0) return(-2);
  u[i]=b[i]/l[i];
}
m[m.dim()-1]=c[c.dim()-1];
l[a.dim()-1]=a[a.dim()-1]-m[a.dim()-2]*u[a.dim()-2];
return(0);

11. Resolución mediante Crout

Array1D< real > l(a.dim());
Array1D< real > m(c.dim());
Array1D< real > u(b.dim());
int error=crout_descomposicion(a,b,c,l,m,u);
if(error<0) return(Array1D<real>());
Array1D< real > z=crout_descenso (l,m,t);
if(z.dim()==0) return(Array1D<real>());
return crout_remonte (u,z);