Discussion :

[Nulysses]

Bonjour,

j'ai le code suivant pour inverser une matrice contenue dans un fichier .dat, mon problème est que la réponse qu'il me donne est 47. Je ne vois vraiment pas ou il y a erreure(s). Voici le code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
PROGRAM MULTISCALE_1
      IMPLICIT NONE
      INTEGER A,N,INDX,X
C     
C     n is the number of row minus 1 in the data file (the first line contains n)
C     
      OPEN (22,file='Cijkl.dat',form='formatted',status='old')
      READ (22,*) A
      CALL MIGS(A,N,X,INDX)
      PRINT *,X
      END
 
 
 
      SUBROUTINE MIGS(A,N,X,INDX)
C     
C     Subroutine to invert matrix A(N,N) with the inverse stored
C     in X(N,N) in the output.
 
C     
      DIMENSION A(N,N),X(N,N),INDX(N),B(N,N)
C     
 
      DO      20 I = 1, N
         DO    10 J = 1, N
            B(I,J) = 0.0
 10      CONTINUE
 20   CONTINUE
      DO      30 I = 1, N
         B(I,I) = 1.0
 30   CONTINUE
C     
      CALL ELGS(A,N,INDX)
C     
      DO     100 I = 1, N-1
         DO    90 J = I+1, N
            DO  80 K = 1, N
               B(INDX(J),K) = B(INDX(J),K)
     *              -A(INDX(J),I)*B(INDX(I),K)
 80         CONTINUE
 90      CONTINUE
 100  CONTINUE
C     
      DO     200 I = 1, N
         X(N,I) = B(INDX(N),I)/A(INDX(N),N)
         DO   190 J = N-1, 1, -1
            X(J,I) = B(INDX(J),I)
            DO 180 K = J+1, N
               X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
 180        CONTINUE
            X(J,I) =  X(J,I)/A(INDX(J),J)
 190     CONTINUE
 200  CONTINUE
C     
      RETURN
      END
C     
      SUBROUTINE ELGS(A,N,INDX)
C     
C     Subroutine to perform the partial-pivoting Gaussian elimination.
C     A(N,N) is the original matrix in the input and transformed
C     matrix plus the pivoting element ratios below the diagonal in
C     the output.  INDX(N) records the pivoting order.
 
C     
      DIMENSION A(N,N),INDX(N),C(N)
C     
C     Initialize the index
C     
      DO     50    I = 1, N
         INDX(I) = I
 50   CONTINUE
C     
C     Find the rescaling factors, one from each row
C     
      DO     100   I = 1, N
         C1= 0.0
         DO    90   J = 1, N
            C1 = AMAX1(C1,ABS(A(I,J)))
 90      CONTINUE
         C(I) = C1
 100  CONTINUE
C     
C     Search the pivoting (largest) element from each column
C     
      DO     200   J = 1, N-1
         PI1 = 0.0
         DO   150   I = J, N
            PI = ABS(A(INDX(I),J))/C(INDX(I))
            IF (PI.GT.PI1) THEN
               PI1 = PI
               K   = I
            ELSE
            ENDIF
 150     CONTINUE
C     
C     Interchange the rows via INDX(N) to record pivoting order
C     
         ITMP    = INDX(J)
         INDX(J) = INDX(K)
         INDX(K) = ITMP
         DO   170   I = J+1, N
            PJ  = A(INDX(I),J)/A(INDX(J),J)
C     
C     Record pivoting ratios below the diagonal
C     
            A(INDX(I),J) = PJ
C     
C     Modify other elements accordingly
C     
            DO 160   K = J+1, N
               A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
 160        CONTINUE
 170     CONTINUE
 200  CONTINUE
C     
      RETURN
      END

Merci d'avance pour votre aide...

Nulysses

[FR119492]

Salut!

inverser une matrice

A quoi cela sert-il? En d'autres termes, quand tu auras obtenu ta matrice inverse, qu'est-ce que tu en feras?
Jean-Marc Blanc

[Nulysses]

Bonjour,

la matrice inverse me servira à calculer une matrice de contraintes à partir d'une matrice de déformation et je n'ai que la matrice inverse de celle qui me servirait !

Merci,
Nulysses

[FR119492]

Salut!
Si je te comprends bien, tu as une matrice de déformation D, une certaine matrice A et tu veux calculer une matrice de contrainte C telle que A*C=D. Est-ce bien ça?
Jean-Marc Blanc

[flavien tetart]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

READ (22,*) A

J'ai l'impression que ton problème vient de la lecture de ta matrice. FORTRAN lit les fichiers ligne par ligne. Ton code ne doit lire que la première ligne de ta matrice.

[Ehouarn]

Bonjour,

Il y a clairement un problème de déclaration des variables dans le programme principal: A est déclaré comme un entier (la routine MIGS s'attends à recevoir un tableau de réels), idem pour X et N n'est pas initialisé (alors que c'est sensé donner la taille de la matrice A).