FORUM MP | Fortran : programme decomposition LU

COMMUNAUTE MP

Identifiez vous ...

DOMAINE MP

Bavardages

Langages Généraux

Langages Web

Langages DotNet

Autres langages

Dev. Jeux Video

Sécurité

Sys. Exploitation

Graphismes

Logiciels

Réseaux

Bases de données

Méthodologies

Emplois High-tech

Aide juridique

Articles juridiques

FORUM

CHAT MP IRC

Votre pseudo ...

Srv: irc.moteurprog.com

Chan: #MoteurProg

PARTICIPER

Plus de 3500 emplois.

Visiteur MP

programme decomposition LU
	Forum : FORTRAN Sous Catégorie : Aucune Type du sujet : Sujet Normale FAQ : FAQ FORTRAN
		SUIVI PAR MAIL INACTIF
		SUJET NON RESOLU
		SUJET ACTIF
		N'APPARTIENT PAS A LA FAQ

FORUM FORTRAN
	Page précedente	Page suivante

gibobo
Nouveau membre

Inscrit : 24/03/2007
Messages : 1

#128122 Posté le 24/03/07 à 18:11

le programme que j'ai fais ne donne pas les bon resultat quelqu'un peu m'aider:

fichier d'entree
3
3. 2. 7.
1. 3. 2. 13.
1. 3. 11.

* Methode de factorisation L.U (sans ‚change de lignes)
C ****************************************
REAL A(20,21),X(20),Y(20)
OPEN (1,FILE='LU.DON',ACCESS='SEQUENTIAL',STATUS='OLD')
READ(1,*) N
N1=N+1
READ(1,*) ((A(I,J),J=1,N1),I=1,N)
CLOSE(1,STATUS='KEEP')
OPEN (1,FILE='LU.DAT',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
WRITE(1,100)N
100 FORMAT(T15,'SystŠme lin‚aire d''ordre ',I4)
WRITE(1,*)
WRITE(1,*) 'Matrice'
DO 150 I=1,N
150 WRITE(1,400)(A(I,J),J=1,N)
400 FORMAT(4E17.7)
WRITE(1,*)
WRITE(1,*)' Second membre '
WRITE(1,400)(A(I,N1),I=1,N)
DO 1 K=1,N-1
C
C Elimination de X(k)
C
3 PIVOT=1./A(K,K)
DO 1 I=K+1,N
S=A(I,K)
IF(S.EQ.0)GOTO 1
S=PIVOT*S
A(I,K)=S
DO 9 J=K+1,N
9 A(I,J)=A(I,J)-S*A(K,J)
1 CONTINUE
C
DO 550 I=1,N
550 WRITE(1,400)(A(I,J),J=1,N)

C Matrices L et U de la factorisation
WRITE(1,*)
WRITE(1,*) 'Matrices L et U'
DO 151 I=1,N
DO 152 J=1,N
Y(J)=0
152 X(J)=0
X(I)=1.
IF(I.EQ.1)GOTO 157
DO 153 J=1,I-1
153 X(J)=A(I,J)
157 DO 154 J=I,N
154 Y(J)=A(I,J)
151 WRITE(1,500)(X(J),Y(J),J=1,N)
500 FORMAT(1X,6E13.6)
WRITE(1,*)
C
C R‚solution du systŠme triangulaire L.Y=B
C
IF(A(N,N).EQ.0)GOTO 1000
c WRITE(1,*)
c WRITE(1,*)'Solution'
Y(1)=A(1,N1)
DO 155 I=2,N
S=A(I,N1)
DO 156 J=1,I-1
156 S=S-A(I,J)*Y(J)
155 Y(I)=S
C
C R‚solution du systŠme U.X=Y
X(N)=A(N,N1)/A(N,N)
DO 7 K=N-1,1,-1
S=A(K,N1)
DO 8 J=K+1,N
8 S=S-A(K,J)*X(J)
7 X(K)=S/A(K,K)
WRITE(1,200)
200 FORMAT(T30,19('*'))
DO 4 I=1,N
4 WRITE(1,300) X(I)
300 FORMAT(T30,'*',T32,E14.7,T48,'*')
WRITE(1,200)
GOTO 2000
1000 WRITE(1,*)' SystŠme singulier det A = 0'
2000 CLOSE(1,STATUS='KEEP')
END

HAUT DE PAGE

Publicité

Inscrit : X
Messages : X

#Aucun

HAUT DE PAGE

solos
Nouveau membre

Inscrit : 11/01/2008
Messages : 1

#147472 Posté le 11/01/08 à 17:25

Hello!!
Voila qq subroutine pour la resolution de Ax = b en f90

! Decomposition A = LU avec pivot
subroutine DecompLU(A, n, PIV)
integer, intent(in) :: n
real, dimension(N,N), intent(inout) :: A
real, dimension(N), intent(out) :: PIV
integer :: i, j, amax, rep, k
real, dimension(n) :: v
real :: temp

do i=1,n
PIV(i)=i
v(i)=0
end do

do j=1,n-1
amax=0
rep=j
do i=j,n
if (abs(A(i,j))>amax) then
amax=abs(A(i,j))
rep=i
endif
enddo
v(:)=A(j,:)
A(j,:)=A(rep,:)
A(rep,:)=v(:)
temp=PIV(j)
PIV(j)=PIV(rep)
PIV(rep)=temp
do i=j+1,n
A(i,j)=A(i,j)/A(j,j)
do k=j+1,n
A(i,k)=A(i,k)-A(i,j)*A(j,k)
enddo
enddo
enddo
endsubroutine DecompLU

! La resolution inferieure
subroutine ResolvInf(A,n,x,b)
integer, intent(in) :: n
real, dimension(n,n), intent(in) :: A
real, dimension(n), intent(in) :: b
real, dimension(n), intent(out) :: x
integer :: i,j
real :: somme

x(1)=b(1)
do i=2,n
somme=0
do j=1,i-1
somme=somme+A(i,j)*x(j)
end do
x(i)=(b(i)-somme)
enddo
end subroutine ResolvInf

!!La resolution superieure
subroutine ResolvSup(A,n,x,b)
integer, intent(in) :: n
real, dimension(n,n), intent(in) :: A
real, dimension(n), intent(in) :: b
real, dimension(n), intent(out) :: x
integer :: i,j
real :: somme

x(n)=b(n)/A(n,n)
do i=n-1,1,-1
somme=0
do j=n,i+1,-1
somme=somme+A(i,j)*x(j)
end do
x(i)=(b(i)-somme)/A(i,i)
enddo
end subroutine ResolvSup

Il ne te reste plus qu'à faire le main ;)

Bon Courage

Las

HAUT DE PAGE

PAGE : [1]

.: Site Web développé par Yann Pingot et l'équipe MPWG avec www.evolvia-web.com :.