Discussion :

[Brian1780]

Bonsoir,

J'ai débuté en Fortran 90 il y a 2 ou 3 jours et mon objectif est de comparer des temps de calcul entre des codes écrits en Matlab, C et F90. Pour un premier test j'ai réalisé un petit code d'éléments finis en élasticité linéaire qui consiste en gros à construire une matrice K, un vecteur F et de résoudre le système linéaire K*U = F.

Mon code F90 étant beaucoup plus lent que ses copains en C et Matlab, j'aimerais avoir un avis sur ce petit code pour savoir si je réalise des mauvaises opérations, mauvaises initialisations des matrices etc. 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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
program linear_elasticity
 
!---------------------VARIABLES------------------------
implicit none
 
integer :: e, NE, N, ND, i, j, NLIN, NRHS, INFO
integer, parameter :: r=3, s=6
integer, dimension(3) :: nodes
integer, dimension(6) :: connect
integer, dimension(:, :), allocatable :: t
integer, dimension(:), allocatable :: fixeddofs
integer, dimension(:), allocatable :: IPIV
double precision, dimension(:, :), allocatable :: p, K
double precision, dimension(:), allocatable :: F
double precision, dimension(3, 6) :: Be, D
double precision, dimension(6, 3) :: BeT
double precision, dimension(6, 6) :: Ke
double precision, dimension(3, 3) :: L
double precision, dimension(3) :: x, y
double precision, dimension(:, :), allocatable :: IDENTITY
double precision :: A, young, nu, lambda, mu, t1, t2, ALPHA, BETA
 
!--------------------INSTRUCTIONS----------------------
call read_file(10, 'mesh.txt', NE, N, ND, p, t, fixeddofs)
 
call cpu_time (t1)
 
ALPHA = 1.0
BETA = 0.0
NLIN = 2*N
NRHS = 1
 
allocate (K(1:2*N,1:2*N))
allocate (F(1:2*N))
allocate (IDENTITY(1:ND, 1:ND))
allocate (IPIV(1:2*N))
 
IDENTITY = 0.0
DO i = 1,ND
   IDENTITY(i, i) = 1
ENDDO
 
young = 200e-2 !just for testing
nu = 0.3 
lambda = (young*nu)/((nu+1)*(1-2*nu))
mu = young/(2*(nu+1))
lambda = 2*lambda*mu/(lambda+2*mu)
L(1, :) = (/ lambda+2*mu, lambda, 0.0d0 /)
L(2, :) = (/ lambda, lambda+2*mu, 0.0d0 /)
L(3, :) = (/ 0.0d0, 0.0d0, mu /)
 
K = 0.0d0
 
DO e = 1,NE
 
   nodes = t(e, :)
   x = p(nodes, 1)
   y = p(nodes, 2)
 
   A = 0.5*((x(2)-x(1))*(y(3)-y(1)) - (y(2)-y(1))*(x(3)-x(1)))
 
   Be(1, :) = (/ (0.5/A)*(y(2)-y(3)), 0.0d0, (0.5/A)*(y(3)-y(1)), 0.0d0, (0.5/A)*(y(1)-y(2)), 0.0d0 /)
   Be(2, :) = (/ 0.0d0, (0.5/A)*(x(3)-x(2)), 0.0d0, (0.5/A)*(x(1)-x(3)), 0.0d0, (0.5/A)*(x(2)-x(1)) /)
   Be(3, :) = (/ (0.5/A)*(x(3)-x(2)), (0.5/A)*(y(2)-y(3)), (0.5/A)*(x(1)-x(3)), (0.5/A)*(y(3)-y(1)), &
        &        (0.5/A)*(x(2)-x(1)), (0.5/A)*(y(1)-y(2)) /)
 
   D = 0.0d0
   Ke = 0.0d0
 
   DO j = 1,r
      DO i = 1,s
         BeT(i, j) = Be(j, i)
      ENDDO
   ENDDO
 
   CALL DGEMM('N', 'N', r, s, r, ALPHA, L, r, Be, r, BETA, D, r)   
   CALL DGEMM('N', 'N', s, s, r, ALPHA, BeT, s, D, r, BETA, Ke, s) 
 
   connect = (/ 2*nodes(1)-1, 2*nodes(1), 2*nodes(2)-1, 2*nodes(2), 2*nodes(3)-1, 2*nodes(3) /)
 
   K(connect, connect) = K(connect, connect) + A*Ke
 
ENDDO
 
K(fixeddofs, :) = 0.0
K(fixeddofs, fixeddofs) = IDENTITY
F(fixeddofs(1:ND/2)) = p((1/2)*fixeddofs(1:ND/2)+1, 1)
 
CALL DGESV(NLIN, NRHS, K, NLIN, IPIV, F, NLIN, INFO);
 
call cpu_time (t2)
print*, t2-t1
 
if (INFO==0) then
write(*,*)'INFO=0:exit successful'
else
write(*,*)'bad exit, INFO=',INFO,'system not solved'
endif
 
!---------------------SUBROUTINES----------------------
contains
 
!-------------SUBROUTINE FOR READING A FILE------------
subroutine read_file (UnitNum, FileName, NE, N, ND, p, t, fixeddofs)
 
integer, intent (in) :: UnitNum
integer, intent (out) :: Ne, N, ND
integer, dimension(:, :), intent (out),  allocatable :: t
integer, dimension(:), intent (out), allocatable :: fixeddofs
integer :: i
character (len=*), intent (in) :: FileName
double precision, dimension(:, :), allocatable :: p
 
open(unit=UnitNum, file=FileName, status='old', action='read')
 
read(10,*) NE
read(10,*) N
read(10,*) ND
 
allocate (p(1:N,1:2))
allocate (t(1:NE,1:3))
allocate (fixeddofs(1:ND))
 
do i=1,N
read(10,fmt='(2F12.9)')p(i, 1),p(i, 2)
enddo
 
do i=1,NE
read(10,*)t(i, 1),t(i, 2),t(i, 3)
enddo
 
do i=1,ND
read(10,*)fixeddofs(i)
enddo
 
close(10)
 
end subroutine read_file
 
!---------------------END-SUBROUTINES------------------
 
end program linear_elasticity

Merci beaucoup et bonne soirée !

[__dardanos__]

Bonsoir.
La question est intéressante. Je n'ai pas vu d'anomalie ou de maladresse. On pourrait factoriser le calcul des tableaux (mais le gain devrait être faible).
Pour aller plus loin, il faudrait réaliser le profil d'exécution avec des outils comme Gprof ou Valgrind. On connaîtrait les parties du programme qui consomment le temps d’exécution.
Peux tu donner la comparaison avec les autres langages ? Le compilateur et les options utilisées ?

[Brian1780]

Bonsoir. Alors voici quelques éléments de comparaison.

J'utilise un maillage (données stockées dans le mesh.txt) qui a les caractéristiques suivantes:
* NE = 1112 (nombre d'éléments du maillage)
* N = 612 (nombre de points dans le maillage)

Les temps de calcul sont en moyenne (estimés avec 100 réalisations):

* Matlab : 0.1 sec
* C : 0.05 sec
* F90 : 0.13 sec (!)

Maillage plus fin (NE = 3620, N = 1879)

* Matlab : 1.1 sec
* C : 0.95 sec
* F90 : 3.4 sec (!)

Pour compiler en F90 je fais :

gfortran -o LS LS.f90 -Wall -fcheck=a -lblas -llapack

Mon code C est un MEX-file en réalité, ce qui me permet de l'exécuter via Matlab et je le compile de la manière suivante :

mex linear_solver.c -lmwblas -lmwlapack

A titre d'information je vous mets les codes Matlab et C ci-dessous. Je ne connais pas Gprof et Valgrind, je vais y jeter un oeil. Sachant que les matrices peuvent avoir une taille important, je devrais peut-être tester une "vectorisation".

Merci, bonne soirée.

MATLAB :

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
E = 200E-2; nu = 0.3; 
lambda = (E*nu)/((nu+1)*(1-2*nu)); mu = E/(2*(nu+1));
Lambda = 2*lambda*mu/(lambda+2*mu); 
L = [Lambda+2*mu, Lambda, 0; Lambda, Lambda+2*mu, 0; 0, 0, mu];
Ne = size(t, 1); Nn = size(p, 1); ND = length(FIXEDDOFS);
 
K = zeros(2*Nn, 2*Nn);
F = zeros(2*Nn, 1);
Be = zeros(3, 6);
 
for e=1:Ne
 
    nodes = t(e, :);
    x = p(nodes, 1);
    y = p(nodes, 2);
 
    A = 0.5*( (x(2)-x(1))*(y(3)-y(1)) - (y(2)-y(1))*(x(3)-x(1)) );
    Be(1, 1:2:end) = (0.5/A)*[y(2)-y(3), y(3)-y(1), y(1)-y(2)];
    Be(2, 2:2:end) = (0.5/A)*[x(3)-x(2), x(1)-x(3), x(2)-x(1)];
    Be(3, 1:2:end) = Be(2, 2:2:end);
    Be(3, 2:2:end) = Be(1, 1:2:end);
 
    Ke = A*transpose(Be)*L*Be;
 
    connec = [2*nodes(1)-1, 2*nodes(1), 2*nodes(2)-1, 2*nodes(2), ...
        2*nodes(3)-1, 2*nodes(3)];
 
    K(connec, connec) = K(connec, connec) + Ke;
 
end
 
K(FIXEDDOFS, :) = zeros(length(FIXEDDOFS), 2*Nn);
K(FIXEDDOFS, FIXEDDOFS) = eye(length(FIXEDDOFS));
 
F(FIXEDDOFS, 1) = [p(0.5*(FIXEDDOFS(1:ND/2, 1)+1), 1); zeros(ND/2, 1)];
 
U = K\F;

C MEX-FILE :

Code c :

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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mex.h"
#include "lapack.h"
#include "blas.h"
 
#define t2D(A, i, j, N) (A[ (j)*(N) + (i) ])
 
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
 
    double *p, *t, *fixeddofs, *output, *Ke, *K, *F, *D;
    double x[3], y[3], Be[18], BeT[18], L[9];
    double A, Lambda, lambda, mu, E, nu;
    double one = 1.0, zero = 0.0;
    int nodes[3], connect[6];
    int e, I, J, II, JJ, INFO;
    int *IPIV;
 
    char *chn = "N";
 
    size_t N, NE, ND, nL, nB, NLIN, NRHS;
 
    if(nrhs!=3) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Three inputs required : p, t, fixeddofs");
    }
    if(nlhs!=2) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","Two outputs required : K and F");
    }
 
    p = mxGetPr(prhs[0]);
    t = mxGetPr(prhs[1]);
    fixeddofs = mxGetPr(prhs[2]);
 
    N = mxGetM(prhs[0]);
    NE = mxGetM(prhs[1]);
    ND = mxGetM(prhs[2]);
 
    nL = 3;
    nB = 6;
    NLIN = 2*N;
    NRHS = 1;
 
    D = mxCalloc( 3*6, sizeof(double) );
    Ke = mxCalloc( 6*6, sizeof(double) );
    K = mxCalloc( NLIN*NLIN, sizeof(double) );
    F = mxCalloc( NLIN*NRHS, sizeof(double) );
    /* IIV = malloc( NLIN * sizeof(int) ); */
 
    plhs[0] = mxCreateDoubleMatrix(NLIN, NLIN, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(NLIN, 1, mxREAL);
    K = mxGetPr(plhs[0]);
    F = mxGetPr(plhs[1]);
 
    /* computational routine */
 
    E = 200E-2; nu = 0.3; 
    Lambda = (E*nu)/((nu+1)*(1-2*nu)); mu = E/(2*(nu+1));
    lambda = 2*Lambda*mu/(Lambda+2*mu);
 
    L[0] = lambda+2*mu;
    L[1] = lambda;
    L[3] = lambda;
    L[4] = lambda+2*mu;
    L[8] = mu;
 
    for (e=0;e<NE;e++)
        {
        nodes[0] = t[e+NE*0]-1;
        nodes[1] = t[e+NE*1]-1;
        nodes[2] = t[e+NE*2]-1;
 
        x[0] = p[nodes[0]]; x[1] = p[nodes[1]]; x[2] = p[nodes[2]];
        y[0] = p[nodes[0]+N]; y[1] = p[nodes[1]+N]; y[2] = p[nodes[2]+N]; 
 
        A = 0.5*( (x[1]-x[0])*(y[2]-y[0]) - (y[1]-y[0])*(x[2]-x[0]) );
 
        t2D(Be, 0, 0, 3) = (0.5/A)*(y[1]-y[2]); 
        t2D(Be, 0, 2, 3) = (0.5/A)*(y[2]-y[0]);
        t2D(Be, 0, 4, 3) = (0.5/A)*(y[0]-y[1]);
 
        t2D(Be, 1, 1, 3) = (0.5/A)*(x[2]-x[1]);
        t2D(Be, 1, 3, 3) = (0.5/A)*(x[0]-x[2]);
        t2D(Be, 1, 5, 3) = (0.5/A)*(x[1]-x[0]);
 
        t2D(Be, 2, 0, 3) = (0.5/A)*(x[2]-x[1]);
        t2D(Be, 2, 1, 3) = (0.5/A)*(y[1]-y[2]);
        t2D(Be, 2, 2, 3) = (0.5/A)*(x[0]-x[2]);
        t2D(Be, 2, 3, 3) = (0.5/A)*(y[2]-y[0]);
        t2D(Be, 2, 4, 3) = (0.5/A)*(x[1]-x[0]);
        t2D(Be, 2, 5, 3) = (0.5/A)*(y[0]-y[1]);
 
        for (I=0;I<6;I++){ for (J=0;J<3;J++) { 
            t2D(BeT, I, J, 6) = A*t2D(Be, J, I, 3); } }
 
        dgemm(chn, chn, &nL, &nB, &nL, &one, L, &nL, Be, &nL, &zero, D, &nL);
        dgemm(chn, chn, &nB, &nB, &nL, &one, BeT, &nB, D, &nL, &zero, Ke, &nB);
 
        connect[0] = 2*nodes[0];
        connect[1] = 2*nodes[0]+1;   
        connect[2] = 2*nodes[1];
        connect[3] = 2*nodes[1]+1;
        connect[4] = 2*nodes[2];
        connect[5] = 2*nodes[2]+1;
 
        for(I=0;I<6;I++) { for(J=0;J<6;J++)
                {
                t2D(K, connect[I], connect[J], 2*N) = t2D(K, connect[I], connect[J], 2*N) + t2D(Ke, I, J, 6); 
                } } 
 
        }
 
for (I=0;I<ND;I++) { for (J=0;J<2*N;J++) { II = fixeddofs[I]-1;
        t2D(K, II, J, 2*N) = 0; } }
 
for (I=0;I<ND;I++) { II = fixeddofs[I]-1; 
    t2D(K, II, II, 2*N) = 1;
    }
 
for (I=0;I<ND/2;I++) { II = fixeddofs[I]-1; JJ = 0.5*(fixeddofs[I]-1+1); 
    F[II] = p[JJ];
    }
 
dgesv(&NLIN, &NRHS, K, &NLIN, IPIV, F, &NLIN, &INFO);
 
/* printf(" INFO = %d : ", INFO);
if (INFO==0)
    {
    printf("  exit successful  \n ");
    }
else
    {
    printf(" bad exit, system not solved \n ");
    }  */
}

[__dardanos__]

L'écart absolu n'est pas si important que ça. Il faudrait connaitre l'endroit où on passe le plus de temps.
Sinon gfortran a une option "O" pour optimiser la compilation.

[Brian1780]

Bonjour, je n'ai pas réussi à faire fonctionner gprof ou valgrind sur ma version de mac os x donc je me suis tourné vers le "Time Profiler" de Xcode. J'ai les résultats dans un fichier qui ne s'ouvre qu'avec Xcode alors j'ai pris un screenshot :



Je vais essayer de comprendre mais pour l'instant j'ai un peu du mal à interpreter ce qu'il me raconte

Bonne journée.

[Brian1780]

Bonjour,
En regardant manuellement avec des call_cpu time, j'ai trouvé que c'est la fin du code qui prend un temps fou !

Entre les lignes :

53 et 84 : 2.5480000000000017E-003
84 et 88 : 2.0440000000000041E-003

et le DGESV : 3.3610699999999998 !

Je vais vérifier si je ne fais pas d'erreur dans ma manière d'appeler DGESV.

Bonne journée.

[Ehouarn]

Bonjour,

Si une bonne partie du travail est fait par des routines de la librairie LAPACK, alors il est évident que la façon dont cette dernière a été compilé (degré d'optimisation?) est prépondérante.

[Brian1780]

Bonjour,

J'ai essayé avec un -O3 de la manière suivante :

gfortran -g -pg -o LS LS.f90 -Wall -O3 -fcheck=all -lblas -llapack ,

mais ça ne change rien du tout, c'est bizarre. J'ai la version suivante de gfortran : GNU Fortran (GCC) 5.2.0 (sur mac os x 10.11.2), peut-être qu'il existe une plus récente, je vais voir. Mon code C utilise les mêmes routines de blas et lapack, et pourtant il n'est pas si lent !

Bonne journée.

[Ehouarn]

Bonjour,

Mon code C utilise les mêmes routines de blas et lapack, et pourtant il n'est pas si lent !
Bonne journée.

Les routines oui, mais pas forcément les mêmes librairies (qui peuvent avoir été compilées avec des niveaux d'optimisation différents).

Tu disais auparavant

Mon code C est un MEX-file en réalité, ce qui me permet de l'exécuter via Matlab et je le compile de la manière suivante :
mex linear_solver.c -lmwblas -lmwlapack

Ce qui donne plutôt l'impression que ce ne sont pas les mêmes librairies qui sont utilisées quand tu compiles tes programmes C (via/dans Matlab) et Fortran.

[Brian1780]

Bonsoir,

D'accord, je ne savais pas qu'il y avait une différence ! Je vais essayer de passer mon code Fortran en MEX file par curiosité, ou peut-être mettre à jour mes librairies blas et lapack (installées par mac ports).

Merci, bonne soirée.