How does .Fortran "dqrls" work?

   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"))

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.