-2

Не уверен, что это правильное сообщество, чтобы спросить об этом, но здесь идет.

Один из моих коллег натолкнулся на этот код в ходе использования пакета Quantile Regression with R. Соответствующий код говорит, что он называет код Fortran. Мы нашли этот файл на Фортране, но не знаем, что это за другой файл. Это файл .f, но это не Фортран, потому что он использует!= если условие (Фортран использует NE), оно использует фигурные скобки, и, наконец, операторы if/do завершаются неправильно.

Тем не менее, он очень похож на fortran, настолько, что мои редакторы распознают большую часть кода и цветовых выражений соответственно. Что мы знаем:

-комментарии, указанные хэшем, #

-не равно, обозначается взрыва равных,!знак равно

подпрограммы, определенные как:

subroutine sub_name(arguments)
sub_code
return
end

-if/do операторы не имеют оконечных строк, но содержатся в фигурных скобках

-Если формат выписки:

if(dabs(s(i))<eps){                                             
        z(i)=dmax1(s(i), zero) + eps                        
        w(i)=dmax1(-s(i),zero) + eps
        }
    else {
        z(i)=dmax1(s(i), zero) 
        w(i)=dmax1(-s(i),zero)
        }

Формат заявления:

do i = 1,n{
        d(i) = one/(z(i)/x(i) + w(i)/s(i))
        ds(i)=z(i)-w(i)
        dz(i)=d(i)*ds(i)
        }

вот весь файл, если это поможет. rqfnb.f

subroutine rqfnb(n,p,a,y,rhs,d,u,beta,eps,wn,wp,nit,info)
integer n,p,info,nit(3)
double precision a(p,n),y(n),rhs(p),d(n),u(n),wn(n,9),wp(p,p+3)
double precision one,beta,eps
parameter( one = 1.0d0)
call lpfnb(n,p,a,y,rhs,d,u,beta,eps,wn(1,1),wn(1,2),
        wp(1,1),wn(1,3),wn(1,4),wn(1,5),wn(1,6),
        wp(1,2),wn(1,7),wn(1,8),wn(1,9),wp(1,3),wp(1,4),nit,info)
return
end


# This is a revised form of my primal-dual log barrier form of the
# interior point LP solver based on Lustig, Marsten and Shanno ORSA J Opt 1992.
# It is a projected Newton primal-dual logarithmic barrier method which uses
# the predictor-corrector approach of Mehrotra for the mu steps.
# For the sake of brevity we will call it a Frisch-Newton algorithm.
# Problem:
#   min c'x s.t. Ax=b, 0<=x<=u
# 
# Denote dx,dy,dw,ds,dz as the steps for the respective variables x,y,w,s,z
# 
subroutine lpfnb(n,p,a,c,b,d,u,beta,eps,x,s,y,z,w,
    dx,ds,dy,dz,dw,dr,rhs,ada,nit,info)

integer n,p,pp,i,info,nit(3),maxit
double precision a(p,n),c(n),b(p)
double precision zero,one,mone,big,ddot,dmax1,dmin1,dxdz,dsdw
double precision deltap,deltad,beta,eps,mu,gap,g
double precision x(n),u(n),s(n),y(p),z(n),w(n),d(n),rhs(p),ada(p,p)
double precision dx(n),ds(n),dy(p),dz(n),dw(n),dr(n)

parameter( zero  = 0.0d0)
parameter( one   = 1.0d0)
parameter( mone  = -1.0d0)
parameter( big  = 1.0d+20)
parameter( maxit  = 50)

# Initialization:  We follow the notation of LMS
# On input we require:
# 
#   c = n-vector of marginal costs (-y in the rq problem)
#   a = p by n matrix of linear constraints (x' in rq)
#   b = p-vector of rhs ((1-tau)x'e in rq)
#       u = upper bound vector ( e in rq)
#   beta = barrier parameter, LMS recommend .99995
#   eps = convergence tolerance, LMS recommend 10d-8
#   
#   the integer vector nit returns iteration counts
#   the integer info contains an error code from the Cholesky in stepy
#       info = 0 is fine
#       info < 0 invalid argument to dposv 
#       info > 0 singular matrix
nit(1)=0
nit(2)=0
nit(3)=n
pp=p*p
# Start at the OLS estimate for the dual vector y
call dgemv('N',p,n,one,a,p,c,1,zero,y,1)
do i=1,n
    d(i)=one
call stepy(n,p,a,d,y,ada,info)
if(info != 0) return
# put current residual vector in s (temporarily)
call dcopy(n,c,1,s,1)
call dgemv('T',p,n,mone,a,p,y,1,one,s,1)
# Initialize remaining variables
do i=1,n{
    if(dabs(s(i))<eps){                                             
        z(i)=dmax1(s(i), zero) + eps                        
        w(i)=dmax1(-s(i),zero) + eps
        }
    else {
        z(i)=dmax1(s(i), zero) 
        w(i)=dmax1(-s(i),zero)
        }
    s(i)=u(i)-x(i)
    }
gap = ddot(n,z,1,x,1)+ddot(n,w,1,s,1)                           
while(gap > eps && nit(1)<maxit) {
    nit(1)=nit(1)+1
    do i = 1,n{
        d(i) = one/(z(i)/x(i) + w(i)/s(i))
        ds(i)=z(i)-w(i)
        dz(i)=d(i)*ds(i)
        }
    call dcopy(p,b,1,dy,1)#save rhs
    call dgemv('N',p,n,mone,a,p,x,1,one,dy,1) 
    call dgemv('N',p,n,one,a,p,dz,1,one,dy,1) 
    call dcopy(p,dy,1,rhs,1)#save rhs
    call stepy(n,p,a,d,dy,ada,info)
    if(info != 0) return
    call dgemv('T',p,n,one,a,p,dy,1,mone,ds,1) #ds -> A'dy - ds
    deltap=big
    deltad=big
    do i=1,n{
        dx(i)=d(i)*ds(i)
        ds(i)=-dx(i)
        dz(i)=-z(i)*(dx(i)/x(i) + one)
        dw(i)=-w(i)*(ds(i)/s(i) + one)
        if(dx(i)<0)deltap=dmin1(deltap,-x(i)/dx(i))                                 
        if(ds(i)<0)deltap=dmin1(deltap,-s(i)/ds(i))
        if(dz(i)<0)deltad=dmin1(deltad,-z(i)/dz(i))
        if(dw(i)<0)deltad=dmin1(deltad,-w(i)/dw(i))
        }
    deltap=dmin1(beta*deltap,one)
    deltad=dmin1(beta*deltad,one)
    if(min(deltap,deltad) < one){
        nit(2)=nit(2)+1
        # Update mu
        mu = ddot(n,x,1,z,1)+ddot(n,s,1,w,1)
        g = mu + deltap*ddot(n,dx,1,z,1)+
            deltad*ddot(n,dz,1,x,1) +
            deltap*deltad*ddot(n,dz,1,dx,1)+
            deltap*ddot(n,ds,1,w,1)+
            deltad*ddot(n,dw,1,s,1) +
            deltap*deltad*ddot(n,ds,1,dw,1)
        mu = mu * ((g/mu)**3) /dfloat(2*n)
        # Compute modified step
        do i=1,n{
            dr(i)=d(i)*(mu*(1/s(i)-1/x(i))+
                dx(i)*dz(i)/x(i)-ds(i)*dw(i)/s(i))
            }
        call dswap(p,rhs,1,dy,1)
        call dgemv('N',p,n,one,a,p,dr,1,one,dy,1)# new rhs
        call dpotrs('U',p,1,ada,p,dy,p,info)# backsolve for dy
        call dgemv('T',p,n,one,a,p,dy,1,zero,u,1)#ds=A'ddy
        deltap=big
        deltad=big
        do i=1,n{
            dxdz =  dx(i)*dz(i)
            dsdw =  ds(i)*dw(i)
            dx(i)=  d(i)*(u(i)-z(i)+w(i))-dr(i)
            ds(i)= -dx(i)
            dz(i)= -z(i)+(mu - z(i)*dx(i) - dxdz)/x(i)
            dw(i)= -w(i)+(mu - w(i)*ds(i) - dsdw)/s(i)
            if(dx(i)<0)deltap=dmin1(deltap,-x(i)/dx(i))
            if(ds(i)<0)deltap=dmin1(deltap,-s(i)/ds(i))
            if(dz(i)<0)deltad=dmin1(deltad,-z(i)/dz(i))
            if(dw(i)<0)deltad=dmin1(deltad,-w(i)/dw(i))
            }
        deltap=dmin1(beta*deltap,one)
        deltad=dmin1(beta*deltad,one)
        }
    call daxpy(n,deltap,dx,1,x,1)
    call daxpy(n,deltap,ds,1,s,1)
    call daxpy(p,deltad,dy,1,y,1)
    call daxpy(n,deltad,dz,1,z,1)
    call daxpy(n,deltad,dw,1,w,1)
    gap = ddot(n,z,1,x,1)+ddot(n,w,1,s,1)
    }
# return residuals in the vector x
call daxpy(n,mone,w,1,z,1)
call dswap(n,z,1,x,1)
return
end


subroutine stepy(n,p,a,d,b,ada,info)
integer n,p,pp,i,info
double precision a(p,n),b(p),d(n),ada(p,p),zero
parameter( zero = 0.0d0)
# Solve the linear system ada'x=b by Choleski -- d is diagonal
# Note that a isn't altered, and on output ada returns the upper
# triangle Choleski factor, which can be reused, eg with blas dtrtrs
pp=p*p
do j=1,p
    do k=1,p
        ada(j,k)=zero
do i=1,n
    call dsyr('U',p,d(i),a(1,i),1,ada,p)
call dposv('U',p,1,ada,p,b,p,info)
return
end

1 ответ1

0

В качестве примера кода вы указали Ratfor (Rational Fortran).

https://en.wikipedia.org/wiki/Ratfor

Оригинальный файл цитируется как Ratfor на этой веб-странице.

http://www.econ.uiuc.edu/~roger/research/rqn/rqn.html

Всё ещё ищете ответ? Посмотрите другие вопросы с метками .