triff
Registrado: 13 Mar 2010 Mensajes: 1
| Publicado: 13/03/2010 2:57 pm | | | Título: subrutinas solve y factor |
| 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. Estas subrutinas permiten resolver sistemas de ecuaciones lineales, o el tipico "resuelve Ax=b". 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
he probado a intentar pasar una de las subrutinas a c, pero con mis conocimientos tan escasos de c...
void factor(int a[n][n],int n) { /* le estamos poniendo int para provar */ int am,i,j,k;
/*calcula la descomposicion LU de la matriz A*/ /*el algoritmo de Crout*/ /* avanza por la diagonal*/ /*calcula l(i,j) a l(n,j)*/ for (j=1;j<=n;j++) { for (i=j;i<=n;i++) { for (k=1;k<=j-1;k++) { a(i,j)=a(i,j)-a(i,k)*a(k,j); } } } /*calcula u(j,j+1) a u(j,n) */ am=1/a(j,j); for (i=j+1;i<=n;i++) { for (k=1;k<=j-1;k++) { a(j,i)=a(j,i)-a(j,k)*a(k,i); } } a(j,i)=a(j,i)*am; } } |
|