FACTOR Y SOLVE PARA RESOLUCION DE SISTEMAS DE ECUACIONES

triff
13 de Marzo del 2010
Hola a todos!! Estoy empezando a trabajar en c y tengo que
traducir de fortran a C los siguientes codigos. Para alguien
que controle un poco seguro que es algo muy facil y me ahorraria mucho tiempo. Y asi quedan escritos para el uso y disfrute de mas gente del foro, que veo que estan sedientos de metodos numericos. El codigo es el siguiente:

subroutine factor(a,n)
!donde a es una matriz nxn

real(8)::a(n,n),am
integer::n,i,j,k
!calcula la descomposición LU de la matriz A
!el algoritmo de Crout
!avanza por la diagonal
do j=1,n
!calcula l(j,j) a l(n,j)
do i=j,n
do k=1,j-1
a(i,j)=a(i,j)-a(i,k)*a(k,j)
end do
end do

!calcula u(j,j+1) a u (j,n)
am=1d0/a(j,j)
do i=j+1,n
do k=1,j-1
a(j,i)=a(j,i)-a(j,k)*a(k,i)
end do
a(j,i)=a(j,i)*am
end do
end do

return
end subroutine

!esta surutina permite descomponer cualquier matriz A !en dos matrices diagonales LU y almacenarlas en la !misma matriz A. Si partimos de un sistema A*x=b, !obtendriamos LUx=b. Resolvemos con la siguiente !subrutina haciendo Ly=b , y Ux=y.

subroutine solve(a,b,n)
real(8)::a(n,n),b(n)
integer::n,i,j
!resuelve el sistema Ax=b con la matriz A factorizada en !LU mediante el algoritmo de Crout

!resuelve Ly=b
do i=1,n
do j=1,i-1
b(i)=b(i)-a(i,j)*b(j)
end do
b(i)=b(i)/a(i,i)
end do

!resuelve Ux=y
do i=n,1,-1
do j=i+1,n
b(i)=b(i)-a(i,j)*b(i)
end do
end do

return
end subroutine