From: David L Lorenz <lorenz_at_usgs.gov>

Date: Fri, 27 Apr 2012 07:56:12 -0500

*> do.call(.Fortran,fargs)$coef
*

[1] 11.176571 -0.883272 -1.262772

*>
*

Date: Fri, 27 Apr 2012 07:56:12 -0500

Of course, what you could do is Google dqrls and get the source and documentation. That is because it is in the publically available linpack. If it were not publically available that would not work. Theoretically, all FORTRAN or C code in R should be publically available. Dave

From:

<cberry_at_tajo.ucsd.edu>

To:

<r-devel_at_stat.math.ethz.ch>

Date:

04/27/2012 06:28 AM

Subject:

Re: [Rd] How does .Fortran "dqrls" work?
Sent by:

r-devel-bounces_at_r-project.org

yangleicq <yanglei_cq_at_126.com> writes:

> Hi, all.

*> I want to write some functions like glm() so i studied it.
**> In glm.fit(), it calls a fortran subroutine named "dqrfit" to compute
*

least

> squares solutions

*> to the system
**> x * b = y
**>
**> To learn how "dqrfit" works, I just follow how glm() calls "dqrfit" by
*

my

*> own example, my codes are given below:
**>
**>> qr <-
**>>
*

matrix(c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14,2.3,1.7,1.3,1.7,1.7,1.6,1,1.7,1.7,1.7),ncol=2)

>> qr

*> [,1] [,2]
**> [1,] 4.17 2.3
**> [2,] 5.58 1.7
**> [3,] 5.18 1.3
**> [4,] 6.11 1.7
**> [5,] 4.50 1.7
**> [6,] 4.61 1.6
**> [7,] 5.17 1.0
**> [8,] 4.53 1.7
**> [9,] 5.33 1.7
**> [10,] 5.14 1.7
**>> n=10
**>> p=2
**>> y <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
**>> ny=1L
**>> tol=1e-07
**>> coefficients=double(p)
**>> residuals=double(n)
**>> effects=double(n)
**>> rank=integer(1L)
**>> pivot=1:n
**>> qraux=double(n)
**>> work=double(2*n)
**>>
**>>
**>>
**>> fittt<-.Fortran("dqrls", qr =qr, n = n,
**> + p = p, y = y, ny = ny, tol = tol,
**> coefficients=coefficients,
**> + residuals = residuals, effects = effects,
**> + rank = rank, pivot = pivot, qraux = qraux,
**> + work = work, PACKAGE = "base")
**>>
**>> fittt$coefficients
**> [1] 0 0
*

You have the args for .Fortran wrong. Try:

*> fargs <- structure(list("dqrls", qr = structure(c(1, 1, 1, 1, 1, 1, 1,
*

+ 1, 1, 1, 4.17, 5.58, 5.18, 6.11, 4.5, 4.61, 5.17, 4.53, 5.33, + 5.14, 2.3, 1.7, 1.3, 1.7, 1.7, 1.6, 1, 1.7, 1.7, 1.7), .Dim = c(10L, + 3L)), n = 10L, p = 3L, y = c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, + 6.03, 4.89, 4.32, 4.69), ny = 1L, tol = 1e-11, coefficients = c(0, + 0, 0), residuals = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), effects = c(0, + 0, 0, 0, 0, 0, 0, 0, 0, 0), rank = 0L, pivot = 1:3, qraux = c(0, + 0, 0), work = c(0, 0, 0, 0, 0, 0), PACKAGE = "base"), .Names = c("", + "qr", "n", "p", "y", "ny", "tol", "coefficients", "residuals", + "effects", "rank", "pivot", "qraux", "work", "PACKAGE"))

[1] 11.176571 -0.883272 -1.262772

TIP: It often helps to use something like

debug(function.calling.Fortran)

dput( list(...) , file="fargs" )

so you can later invoke the function as above.

*>
**> but when i use lm() which also calls "dqrls" internally to fit this
*

model,

*> it gives reasonable result.
**>
*

>> lm(y~qr)

*>
**> Call:
**> lm(formula = y ~ qr)
**>
**> Coefficients:
**> (Intercept) qr1 qr2
**> 11.1766 -0.8833 -1.2628
**>
**>
**> when I change the coefficients to be c(1,1), the output from "dqrls",
**> fittt$coefficients also equals to c(1,1). That means the
*

.Fortran("dqrls",

*> qr=qr,n=n,p=p,...) did nothing to the coefficients! I don't know why, is
*

> there anything I did wrong or missed? How can I get the result from

"dqrls"

*> as what lm() or glm() gets from "dqrls"?
**>
*

> Thanks in advance. Best Regards.

*>
**>
**>
**> --
**> View this message in context:
*

http://r.789695.n4.nabble.com/How-does-Fortran-dqrls-work-tp4588973p4588973.html

*> Sent from the R devel mailing list archive at Nabble.com.
**>
*

-- Charles C. Berry Dept of Family/Preventive Medicine cberry at ucsd edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 ______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel [[alternative HTML version deleted]] ______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-develReceived on Fri 27 Apr 2012 - 12:59:05 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

*
Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.
Archive generated by hypermail 2.2.0, at Sat 28 Apr 2012 - 00:10:52 GMT.
*

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*