Eugeen
Member | Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
Уважаемые знатоки фортрана и С!
Есть подпрограмма, являющаяся критической в общем комплексе по времени выполнения:
Код: S U B R O U T I N E GLBQRT
C##################################################
C# ПOДПPOГPAMMA-ФУHKЦИЯ ДЛЯ PACЧETA #
C# ФУHKЦИИ Q(RO,TAU) И EE ЧACTHЫX #
C# ПPOИЗBOДHЫX ПO RO И TAU #
C# ( RO = R/1000 #
C# TAU = 1000/T , ГДE #
C# R-ПЛOTHOCTЬ , KГ/KУБ.M #
C# T-TEMПEPATУPA , ГPAД.K #
C# JDR-УMEHЬШEHHЫЙ HA 1 ПOPЯДOK #
C# ЧACTHOЙ ПPOИЗBOДHOЙ Q(RO,TAU) ПO RO #
C# JDT-УMEHЬШEHHЫЙ HA 1 ПOPЯДOK #
C# ЧACTHOЙ ПPOИЗBOДHOЙ Q(RO,TAU) ПO TAU ) #
C# Автор: Катковский Е.А. Москва, 1975 г. #
C##################################################
I M P L I C I T R E A L * 8 ( A - H , O - Z )
D I M E N S I O N AK( 36 ) , BK( 14 ) , DK( 7 )
C O M M O N /PROPTY/ VAR , TAU , RO , JDT , JDR
D A T A AK
* / +3.9636085D+0 , -5.1028070D+0
* , +2.7407416D+0 , -6.9048554D-1
* , +2.9568796D+1 , -2.9142470D+1
* , +1.5453061D+1 , -3.9661401D+0
* , +5.6323130D+1 , -4.7740374D+1
* , +2.6409282D+1 , -6.3972405D+0
* , +4.3125840D+0 , -9.2734289D+0
* , -7.2546108D-1 , -1.5641040D-1
* , -2.6181978D+1 , +6.5326396D+1
* , -2.6149751D+1 , +6.8335354D+0
* , +1.5597836D+2 , +1.3746153D+2
* , +1.2748742D+2 , -1.7731074D+2
* , -1.6254622D+1 , -3.3301902D+1
* , +7.7779182D+0 , -5.1985860D+0
* , +5.9728487D+0 , +1.5518535D+2
* , -2.4450042D+2 , +3.4218431D+2
* , -3.6093828D+2 , +2.7464632D+2
* , -1.3213917D+2 , +2.9492937D+1 /
D A T A BK
* / +7.1531353D+1 , +1.3041253D+1
* , +3.9917570D+2 , +7.9847970D+1
* , +6.4581880D+2 , +1.3687317D+2
* , +1.0401717D+1 , +6.7874983D+0
* , -7.3396848D+2 , -1.3746618D+2
* , -2.0988866D+2 , +3.3731180D+2
* , -4.1605860D+2 , -4.1030808D+2 /
C############################################
C# BЫЧИCЛEHИE ФУHKЦИИ Q(RO,TAU) #
C# ДЛЯ ЗAДAHHЫX ЗHAЧEHИЙ RO И TAU #
C############################################
TCR = TAU-1.D0/6.47286D-1
T25 = TAU- 2.5D0
RX1 = RO-1.D0
XE = DEXP ( -4.8D0*RO )
N = 1
M = 4
LABL = 0
DO 1 K = 1,5
X = 4.D0
IF ( JDR ) 9 , 12 , 14
1000 N = N+4
1 M = N+3
IF ( JDT ) 1001 , 1001 , 6
1001 M = M+4
LABL = 1
K = 6
X = 8.D0
IF ( JDR ) 9 , 12 , 14
2000 IF ( JDT ) 2001, 4 , 2001
2001 N = N+8
M = N+7
RX1 = RO-0.634D0
LABL = -1
K = 7
X = 8.D0
IF ( JDR ) 9 , 12 , 14
E N T R Y GLBQUT
C########################################
C# BXOД ДЛЯ ПOBTOPHOГO PACЧETA #
C# ФУHKЦИИ Q(RO,TAU) #
C# И EE ЧACTHЫX ПPOИЗBOДHЫX ПO TAU #
C# ПPИ HEИЗMEHHOM ЗHAЧEHИИ RO #
C########################################
T25 = TAU-2.5D0
TCR = TAU-1.D0/6.47286D-1
IF ( JDT ) 2 , 4 , 6
E N T R Y GLBQ1T
C########################################
C# BXOД ДЛЯ ПOBTOPHOГO PACЧETA #
C# ФУHKЦИИ Q(RO,TAU) #
C# ПPИ HEИЗMEHHЫX RO И TAU #
C########################################
2 VAR = ( ( ( ( ( DK( 1 )*T25
* +DK( 2 ) )*T25+DK( 3 ) )*T25
* +DK( 4 ) )*T25+DK( 5 ) )*T25
* +DK( 6 ) )*TCR+DK( 7 )
R E T U R N
E N T R Y GLBQ2T
C########################################
C# BXOД ДЛЯ ПOBTOPHOГO PACЧETA #
C# ПEPBOЙ ЧACTHOЙ ПPOИЗBOДHOЙ OT #
C# Q(RO,TAU) ПO TAU #
C# ПPИ HEИЗMEHHЫX RO И TAU #
C########################################
4 X = TCR+T25
VAR = ( ( ( ( ( ( 4.D0*TCR+X )*DK( 1 ) )*T25
* +( 3.D0*TCR+X )*DK( 2 ) )*T25
* +( TCR+TCR+X )*DK( 3 ) )*T25
* +( TCR+X )*DK( 4 ) )*T25
* +X*DK( 5 ) )+DK( 6 )
R E T U R N
E N T R Y GLBQ3T
C########################################
C# BXOД ДЛЯ ПOBTOPHOГO PACЧETA #
C# BTOPOЙ ЧACTHOЙ ПPOИЗBOДHOЙ OT #
C# Q(RO,TAU) ПO TAU #
C# ПPИ HEИЗMEHHЫX RO И TAU #
C########################################
6 X = T25+T25+TCR
VAR = ( ( ( ( 3.D0*TCR+X )*5.D0*DK( 1 ) )*T25
* +( TCR+TCR+X )*4.D0*DK( 2 ) )*T25
* +( TCR+X )*3.D0*DK( 3 ) )*T25
* +X*( DK( 4 )+DK( 4 ) )
* +DK( 5 )+DK( 5 )
R E T U R N
9 BX = AK( N )
JM = N+1
DO 10 J = JM,M
10 BX = BX*RX1+AK( J )
BJ = BK( 2*K-1 )*RO+BK( 2*K )
11 BJ = BJ*XE+BX
DK( K ) = BJ
IF ( LABL ) 2 , 1000 , 2000
12 JM = M-1
BX = 0.D0
DO 13 J = N,JM
X = X-1.D0
13 BX = BX*RX1+X*AK( J )
BJ = ( BK( 2*K-1 )*RO+BK( 2*K ) )*( -4.8D0 )+BK( 2*K-1 )
GO TO 11
14 JM = M-2
BX = 0.D0
DO 15 J = N,JM
X = X-1.D0
15 BX = BX*RX1+( X*X-X )*AK( J )
BJ = ( BK( 2*K-1 )*RO+BK( 2*K ) )*23.04D0
* -9.6D0*BK( 2*K-1 )
GO TO 11
E N D |
Мой сын сказал мне что я уже отстал от жизни и ".. Фортран - полный отстой", а он легко ускорит это в 2 раза на С. Он сделал:
Код: /*********************************************************************/
/* */
/* GLBQRT.C October,2th,1994 */
/* Version 1.10 */
/* Katkovsky Sergei */
/* */
/*********************************************************************/
/*********************************************************************/
/* */
/* Подпрограмма для расчета функции Q(RO,TAU) */
/* и ее частных производных по RO и TAU. */
/* RO=R/1000;TAU00/T; */
/* R-плотность , Кг/Куб.М;T-Температура , Град.К */
/* JDR-уменьшенный на 1 порядок частной */
/* производной Q(RO,TAU) по RO */
/* JDT-уменьшенный на 1 порядок частной */
/* производной Q(RO,TAU) по TAU */
/* */
/*********************************************************************/
#pragma optimize("acegltz",on)
#pragma loop_opt(on)
#pragma check_stack(off)
#pragma check_pointer(off)
#include"prop.h"
#include<math.h>
#include<stdio.h>
#pragma intrinsic(exp)
static void _fastcall f9(int,int,int);
static void _fastcall f12(int,int,int,double);
static void _fastcall f14(int,int,int,double);
static const double ak[36]={
3.9636085E+0 , -5.1028070E+0
, 2.7407416E+0 , -6.9048554E-1
, 2.9568796E+1 , -2.9142470E+1
, 1.5453061E+1 , -3.9661401E+0
, 5.6323130E+1 , -4.7740374E+1
, 2.6409282E+1 , -6.3972405E+0
, 4.3125840E+0 , -9.2734289E+0
,-7.2546108E-1 , -1.5641040E-1
,-2.6181978E+1 , 6.5326396E+1
,-2.6149751E+1 , 6.8335354E+0
, 1.5597836E+2 , 1.3746153E+2
, 1.2748742E+2 , -1.7731074E+2
,-1.6254622E+1 , -3.3301902E+1
, 7.7779182E+0 , -5.1985860E+0
, 5.9728487E+0 , 1.5518535E+2
,-2.4450042E+2 , 3.4218431E+2
,-3.6093828E+2 , 2.7464632E+2
,-1.3213917E+2 , 2.9492937E+1
};
static const double bk[14]={
7.1531353E+1 , 1.3041253E+1
, 3.9917570E+2 , 7.9847970E+1
, 6.4581880E+2 , 1.3687317E+2
, 1.0401717E+1 , 6.7874983E+0
,-7.3396848E+2 ,-1.3746618E+2
,-2.0988866E+2 , 3.3731180E+2
,-4.1605860E+2 ,-4.1030808E+2
};
static double dk[7],tcr,t25,rx1,xe;
/* Вычисление функции Q(RO,TAU) для заданных RO и TAU */
void fortran glbQrt()
{
double x;
int n=1,m=4,k;
tcr=propty.tau-1/0.647286;
t25=propty.tau-2.5;
rx1=propty.ro-1;
xe=exp(-4.8*propty.ro);
for(k=1;k<=5;k++){
x=4;
if(propty.jdr<0)f9(k,m,n);
if(!propty.jdr)f12(k,m,n,x);
if(propty.jdr>0)f14(k,m,n,x);
n=n+4;m=n+3;
}
if(propty.jdt>0){glbQ3t();return;}
m=m+4;k=6;x=8;
if(propty.jdr<0)f9(k,m,n);
if(!propty.jdr)f12(k,m,n,x);
if(propty.jdr>0)f14(k,m,n,x);
if(!propty.jdt){glbQ2t();return;}
n=n+8;m=n+7;
rx1=propty.ro-0.634;
k=7;x=8;
if(propty.jdr<0)f9(k,m,n);
if(!propty.jdr)f12(k,m,n,x);
if(propty.jdr>0)f14(k,m,n,x);
glbQ1t();
return;
}
/* Функция для повторного расчета функции Q(RO,TAU) и ее */
/* частных производных по TAU при неизменном значении RO */
void fortran glbQut()
{
t25=propty.tau-2.5;tcr=propty.tau-1/0.647286;
if(propty.jdt<0)glbQ1t();
if(!propty.jdt)glbQ2t();
if(propty.jdt>0)glbQ3t();
return;
}
/* Функция для повторного расчета Q(RO,TAU) при неизменных RO и TAU */
void fortran glbQ1t()
{
propty.var=(((((dk[0]*t25+dk[1])*t25
+dk[2])*t25+dk[3])*t25
+dk[4])*t25+dk[5])*tcr+dk[6];
return;
}
/* Функция для повторного расчета первой частной производной */
/* от Q(RO,TAU) по TAU при неизменных RO и TAU */
void fortran glbQ2t()
{
double x=tcrn;
propty.var=((((((4*tcr+x)*dk[0])*t25
+(3*tcr+x)*dk[1])*t25+(2*tcr+x)*dk[2])*t25
+(tcr+x)*dk[3])*t25+x*dk[4])+dk[5];
return;
}
/* Функция для повторного расчета второй частной производной */
/* от Q(RO,TAU) по TAU при неизменных RO и TAU */
void fortran glbQ3t()
{
double x=2*(t25К);
propty.var=((((2*tcr+x)*5*dk[0])*t25
+(tcr+x)*4*dk[1])*t25
+x*3*dk[2])*t25
+(x-tcr)*(2*dk[3])
+2*dk[4];
return;
}
static void _fastcall f9(int k,int m,int n)
{
double bx=ak[n-1],bj;
int j,jm=n+1;
for(j=jm;j<=m;j++)bx=bx*rx1+ak[j-1];
bj=bk[2*k-2]*propty.ro+bk[2*k-1];
bj=bj*xe+bx;
dk[k-1]=bj;
return;
}
static void _fastcall f12(int k,int m,int n,double x)
{
double bx=0,bj;
int j,jm=m-1;
for(j=n;j<=jm;j++){x--;bx=bx*rx1+x*ak[j-1];}
bj=-4.8*(bk[2*k-2]*propty.ro+bk[2*k-1])+bk[2*k-2];
bj=bj*xe+bx;
dk[k-1]=bj;
return;
}
static void _fastcall f14(int k,int m,int n,double x)
{
double bx=0,bj;
int j,jm=m-2;
for(j=n;j<=jm;j++){x--;bx=bx*rx1+(x*x-x)*ak[j-1];}
bj#.04*(bk[2*k-2]*propty.ro+bk[2*k-1])-9.6*bk[2*k-2];
bj=bj*xe+bx;
dk[k-1]=bj;
return;
} |
Но видимого ускорения не получилось.
Может быть кто-то скажет почему?
Плохой код на С? Или вообще нет смысла тягаться с Фортраном?
----------
С нами БОГ и Крестная сила! |
|
