How does .Fortran "dqrls" work?
Locked
 
|
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 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. |
|
yangleicq <[hidden email]> 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. > -- 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
Re: 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: <[hidden email]> To: <[hidden email]> Date: 04/27/2012 06:28 AM Subject: Re: [Rd] How does .Fortran "dqrls" work? Sent by: [hidden email] yangleicq <[hidden email]> 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 > 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
Re: How does .Fortran "dqrls" work?
|
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: > <[hidden email]> > To: > <[hidden email]> > Date: > 04/27/2012 06:28 AM > Subject: > Re: [Rd] How does .Fortran "dqrls" work? > Sent by: > [hidden email] > > > > yangleicq<[hidden email]> 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. >> > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
Re: How does .Fortran "dqrls" work?
|
In reply to this post by David L Lorenz
On 27-04-2012, at 14:56, 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. > Dave Huh? Just get the source code of R. Do a bit of digging and you'll find it. @yang: you should have set n <- 10L p <- 2L and see reply by Chuck. In your case you should have done lm(y ~ qr + 0) Berend > > > From: > <[hidden email]> > To: > <[hidden email]> > Date: > 04/27/2012 06:28 AM > Subject: > Re: [Rd] How does .Fortran "dqrls" work? > Sent by: > [hidden email] > > > > yangleicq <[hidden email]> 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. >> > > -- > 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 > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > > > > [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
«
Return to R devel
|
1 view|%1 views
Loading...
| Powered by Nabble | Edit this page |