如何在 quantreg 包中查看 .Fortran 函数

Question

我试图在quantreg 包中查看函数“crq.fit.pen”的源代码。出于这个原因，我尝试做以下工作。 一开始我只是输入了函数名，结果如下，

crq.fit.pen <-
function (x, y, cen, weights = NULL, grid, ctype = "right") 
{ 
    p <- ncol(x) 
    n <- length(y) 
    if (missing(grid)) 
        grid <- seq(1/n, 1 - 1/n, by = min(0.01, 1/(2 * length(y)^0.7))) 
    if (!is.numeric(grid)) 
        stop("Invalid grid") 
    if (any(grid < 0) || any(grid > 1)) 
        stop("Invalid grid") 
    m <- length(grid) 
    xbar <- apply(x, 2, mean) 
    if (length(weights)) { 
        if (any(weights < 0)) 
            stop("negative weights not allowed") 
        contr <- attr(x, "contrasts") 
        x <- x * weights 
        y <- y * weights 
    } 
    if (ctype == "left") 
        y <- -y 
    s <- rep(0, n) 
    u <- rep(1, n) 
    d <- rep(1, n) 
    r <- rep(1, p) 
    B <- matrix(0, p, m) 
    cc <- as.logical(cen) 
    y1 <- y[cc] 
    n1 <- length(y1) 
    x1 <- x[cc, ] 
    z <- .Fortran("crqfnb", as.integer(n), as.integer(p), a1 = as.double(t(as.matrix(x1))), 
        c1 = as.double(-y1), n1 = as.integer(n1), as.double(x), 
        as.double(y), as.double(cen), B = as.double(B), g = as.double(grid), 
        m = as.integer(m), as.double(r), as.double(s), as.double(d), 
        as.double(u), wn = double(n1 * 9), wp = double((p + 3) * 
            p), info = integer(1), PACKAGE = "quantreg") 
    J <- z$m - 1 
    B <- matrix(-z$B, p, m) 
    B <- B[, 1:J, drop = FALSE] 
    qhat <- t(xbar) %*% B 
    B <- rbind(grid[1:J], B, qhat) 
    dimnames(B) <- list(c("tau", dimnames(x)[[2]], "Qhat"), NULL) 
    if (ctype == "left") { 
        B[1, ] <- 1 - B[1, ] 
        B[-1, ] <- -B[-1, ] 
        B <- B[, ncol(B):1] 
    } 
    B <- list(sol = B, ctype = ctype) 
    class(B) <- "crq" 
    B 
} 
<environment: namespace:quantreg>

正如你在上面看到的，这个函数的主要工作是由另一个函数完成的，这个函数被引用：

z <- .Fortran("crqfnb", as.integer(n), as.integer(p), a1 = as.double(t(as.matrix(x1))), 
        c1 = as.double(-y1), n1 = as.integer(n1), as.double(x), 
        as.double(y), as.double(cen), B = as.double(B), g = as.double(grid), 
        m = as.integer(m), as.double(r), as.double(s), as.double(d), 
        as.double(u), wn = double(n1 * 9), wp = double((p + 3) * 
            p), info = integer(1), PACKAGE = "quantreg") 

我现在的问题是如何查看crqfnb Fortran 基函数？

之后，我执行了以下任务并获得了结果，但我看不到函数 crqfnb 的完整代码。

> untar(download.packages(pkgs = "quantreg", 
+                   destdir = ".", 
+                   type = "source")[,2]) 
trying URL 'http://cran.rstudio.com/src/contrib/quantreg_5.05.tar.gz'
Content type 'application/x-gzip' length 1636075 bytes (1.6 Mb) 
opened URL 
================================================== 
downloaded 1.6 Mb 

sh: /usr/bin/gnutar: No such file or directory 
gzip: error writing to output: Broken pipe 
gzip: ./quantreg_5.05.tar.gz: uncompress failed 
Warning message: 
In untar(download.packages(pkgs = "quantreg", destdir = ".", type = "source")[,  : 
  ‘/usr/bin/gzip -dc './quantreg_5.05.tar.gz' | /usr/bin/gnutar -xf '-'’ returned error code 127 

请帮我看看如何查看函数crqfnb的完整代码？

Answer 1

我从 CRAN 下载代码，进入 src 文件夹，打开 crqfnb.f 瞧？

C Output from Public domain Ratfor, version 1.0
      subroutine crqfnb(n,p,a1,c1,n1,x,y,c,b,g,m,r,s,d,u,wn,wp,info)
      integer n,p,n1,m,info,nit(3)
      double precision a1(p,n1),c1(n),x(n,p),y(n),c(n),b(p,m),g(m)
      double precision wn(n,9),wp(p,p+3),r(p),s(n),d(n),u(n)
      double precision zero,half,one,beta,eps,dh
      parameter( zero = 0.0d0)
      parameter( half = 0.5d0)
      parameter( one = 1.0d0)
      parameter( beta = 0.99995d0)
      parameter( eps = 1.0d-8)
      do23000 k = 2,m 
      dh = -log(one - g(k)) + log(one - g(k-1))
      do23002 i = 1,n 
      u(i) = one
      wn(i,1) = half
      if(d(i) .ge. zero)then
      s(i) = s(i) + dh
      endif
      d(i) = c(i) - s(i)
23002 continue
23003 continue
      call dgemv('T',n,p,one,x,n,d,1,zero,r,1)
      call rqfnb(n1,p,a1,c1,r,d,u,beta,eps,wn,wp,nit,info)
      if(info .ne. 0)then
      goto 23001
      endif
      call dcopy(p,wp,1,b(1,k-1),1)
      call dcopy(n,y,1,d,1)
      call dgemv('N',n,p,one,x,n,b(1,k-1),1,one,d,1)
23000 continue
23001 continue
      m = k-1
      return
      end