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

From: David L Lorenz <lorenz_at_usgs.gov>
Date: Fri, 27 Apr 2012 07:56:12 -0500

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
>  0 0

> 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

 11.176571 -0.883272 -1.262772
>

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.

>
> 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-devel
```
Received 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.