Fortran TD 8 Résolution d'un système linéaire Routine issue by hcw25539

VIEWS: 328 PAGES: 7

									                                 Fortran TD 8

                        Résolution d'un système linéaire

                            Routine issue de LAPACK

________________________________________________________________________________

Ecrire un programme principal qui appelle la routine DGESV


      Voici le code du main, notez la déclaration de n et NRHS en tant que
paramètres afin que la matrice utilisée par la routine DGESV ait bien la même
dimension. En effet si la dimension utilisée est inférieure à la dimension
déclarée, le compilateur rempli les lignes et les colonnes par des 0 et la
matrice appelée par DGESV est automatiquement singulière et le système ne
peut pas être résolu.
On effectue une sauvegarde de la matrice A et du second membre B dans deux
matrices A1 et B1, car DGESV les modifie lors de la résolution du système. A
est décomposée en LU et B est le vecteur solution en sortie.
      Pour s'assurer de la précision de la solution nous avons calculé la norme
de Ax-B qui tends vers 0 lorsque la méthode converge. La norme utilisée est la
norme 1, plus simple à calculer et à manipuler.
      L'écriture des résultats dépend de la variable info en sortie de la routine
DGESV, en effet il peut arriver que la matrice soit singulière et que dans ce
cas la décomposition LU ne marche pas. En conséquence il est important
d'informer l'utilisateur que son système n'a pas pu être résolu.

    program main

c*******************************************************************************
c declaration des variables
c*******************************************************************************

    parameter (n=3,NRHS=1)

    integer info
    double precision A(n,n),B(n,NRHS),A1(n,n)
    double precision B1(n,NRHS),IPIV(n), R(n), norme

c*******************************************************************************
c initialisation de la matrice A et du second membre
c*******************************************************************************

    do j=1,n
      read(*,*) (A(i,j),i=1,n)
    enddo
    do j=1,NRHS
      read(*,*) (B(i,j),i=1,n)
    enddo


c*******************************************************************************
c sauvegarde de la matrice A et de la matrice B dans A1 et B1
c*******************************************************************************

    do i=1,n
      do j=1,n
      A1(i,j)=A(i,j)
      enddo
    B1(i,NRHS)=B(i,NRHS)
    enddo


c*******************************************************************************
c resolution du systeme
c*******************************************************************************

    call dgesv(n,NRHS,A,n,IPIV,B,n,info)

c*******************************************************************************
c calcul du reste puis de la norme
c*******************************************************************************

    norme =0

    do i=1,n

     R(i)=0

       do j=1,n
           R(i)=R(i)+A1(i,j)*B(j,NRHS)
       enddo

     R(i)=R(i)-B1(i,NRHS)
     norme = norme +abs(R(i))

    enddo


c*******************************************************************************
c ecriture des resultats
c*******************************************************************************

     write(*,*) 'resolution du systeme de AX=B'
     write(*,*) 'matrice A'
     write(*,*)
     do i=1,n
   write(*,*) (A1(i,j),j=1,n)
   write(*,*)
 enddo

 write(*,*) 'matrice B'
 write(*,*)
  do i=1,n
   write(*,*) (B1(i,j),j=1,NRHS)
   write(*,*)
 enddo



 if(info.eq.0)then
    write(*,*) 'le systeme a ete resolu'
    write(*,*)
    write(*,*) 'vecteur X solution'
    write(*,*)
     do i=1,n
       write(*,*) (B(i,j),j=1,NRHS)
       write(*,*)
     enddo

   write(*,*) 'norme de AX-B (verification de la validite
cdu resultat)'
   write(*,*)
   write(*,*) norme
 endif


 if(info.gt.0)then
    write(*,*) 'le systeme n''a pas ete resolu'
    write(*,*) 'la matrice est singuliere'
 endif


 if(info.lt.0)then
    write(*,*) 'DGESV a ete appele avec un mauvais parametre'
    write(*,*) 'verifiez votre fichier don'
 endif


stop
end
Ecrire un Makefile qui utilise les bibliothèques liblapack.a et libatlas.a


#FC = fc
#FC = fort77
FC = f77
#FC = g77

OPT = -O
#OPT = -g

main.o:main.f
     $(FC) $(OPT) main.f -c

exe: main.o liblapack.a libatlas.a
      $(FC) $(OPT) main.o liblapack.a libatlas.a -o exe


Test de validité de la méthode

On a pris soin de rentrer les matrices à l'aide d'un fichier don et de sortir les
résultats dans un fichier out. Nous avons trois tests :

                                      Test 1

contenu du fichier don1 :

0
1
1
1
0
1
1
1
0
7
8
9

contenu du fichier out1 :

resolution du systeme de AX=B
matrice A

    0.   1.0000000000000    1.0000000000000

     1.0000000000000 0.     1.0000000000000

     1.0000000000000     1.0000000000000 0.
matrice B

  7.0000000000000

  8.0000000000000

  9.0000000000000

le systeme a ete resolu

vecteur X solution

  5.0000000000000

  4.0000000000000

  3.0000000000000

norme de AX-B (verification de la validite du resultat)

 0.
                                    Test 2

contenu du fichier don2 :

4
1
3
5
2
1
4
5
1
2
3
4
5

contenu du fichier out2 :

resolution du systeme de AX=B
matrice A

    4.0000000000000       5.0000000000000    4.0000000000000

    1.0000000000000       2.0000000000000    5.0000000000000

    3.0000000000000       1.0000000000000    1.0000000000000

matrice B

    2.0000000000000

    3.0000000000000

    4.0000000000000

le systeme a ete resolu

vecteur X solution

    1.5526315789474

    -1.5789473684211

    0.92105263157895

norme de AX-B (verification de la validite    du resultat)

    1.3322676295502D-15
                                   Test 3

contenu du fichier don3 :

1
2
3
4
5
6
7
8
9
1
1
1

contenu du fichier out3 :

resolution du systeme de AX=B
matrice A

    1.0000000000000     4.0000000000000     7.0000000000000

    2.0000000000000     5.0000000000000     8.0000000000000

    3.0000000000000     6.0000000000000     9.0000000000000

matrice B

    1.0000000000000

    1.0000000000000

    1.0000000000000

le systeme n'a pas ete resolu
la matrice est singuliere

On constate que le programme est opérationnel et que la routine DGESV
résoud bien des systèmes linéaires AX=B, avec une précision satisfaisante.

								
To top