[Rd] How does .Fortran "dqrls" work?

[Duncan Murdoch]

On 12-04-27 8:56 AM, David L Lorenz wrote:
>    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.

I'm not sure what you mean by "theoretically".  R is open source, the 
source code is available.  Uwe Ligges wrote an article in R News a few 
years ago telling you where to look for it.

Are you thinking of things that aren't really "in R"?  R can load 
packages that aren't open source (though most of them are), and it can 
make calls to the run-time or operating system libraries, and with some 
compilers/operating systems, those may not be.

Duncan Murdoch

> 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"))
>> do.call(.Fortran,fargs)$coef
> [1] 11.176571 -0.883272 -1.262772
>>
>
> TIP: It often helps to use something like
>
>       debug(function.calling.Fortran)
>
> and then step thru the function till the call you want to study is
> invoked. Then inspect the inputs one-by-one and tinker with them and
> recall the function or save them via
>
>           dput( list(...) , file="fargs" )
>
> so you can later invoke the function as above.
>
> HTH,
>
> Chuck
>
>
>>
>> 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.
>>
>