From: <rvaradhan_at_jhmi.edu>

Date: Thu, 18 Jun 2009 16:05:22 +0200 (CEST)

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat 20 Jun 2009 - 11:52:32 GMT

Date: Thu, 18 Jun 2009 16:05:22 +0200 (CEST)

Yes=2C Peter=2E I did look at it=2C but not carefully enought to catch =
that=2E

Thanks=2C

Ravi=2E

*=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=
*

=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=

=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F

Ravi Varadhan=2C Ph=2ED=2E

Assistant Professor=2C

Division of Geriatric Medicine and Gerontology
School of Medicine

Johns Hopkins University

Ph=2E (410) 502-2619

email=3A rvaradhan=40jhmi=2Eedu

- Original Message ----- From=3A Peter Dalgaard =3CP=2EDalgaard=40biostat=2Eku=2Edk=3E Date=3A Thursday=2C June 18=2C 2009 9=3A15 am Subject=3A Re=3A =5BRd=5D Inverting a square=2E=2E=2E (PR=2313762) To=3A rvaradhan=40jhmi=2Eedu Cc=3A R-bugs=40r-project=2Eorg

=3E Refiling this=2E The actual fix was slightly more complicated=2E Wi=

ll soon

=3E be committed to R-Patched (aka 2=2E9=2E1 beta)=2E

*=3E =
*

=3E -p

=3E =

*=3E rvaradhan=40jhmi=2Eedu wrote=3A
**=3E =3E Full=5FName=3A Ravi Varadhan
**=3E =3E Version=3A 2=2E8=2E1
**=3E =3E OS=3A Windows
*

=3E =3E Submission from=3A (NULL) (162=2E129=2E251=2E19)

=3E =3E =

=3E =3E =

=3E =3E Inverting a matrix with solve()=2C but using LAPACK=3DTRUE=2C g=

ives erroneous

=3E =3E results=3A

=3E =

=3E Thanks=2C but there seems to be a much easier fix=2E

*=3E =
*

=3E Inside coef=2Eqr=2C we have

=3E =

=3E coef=5Bqr=24pivot=2C =5D =3C-

=3E =2ECall(=22qr=5Fcoef=5Freal=22=2C qr=2C y=2C PACKAGE =3D =22base=22=

)=5Bseq=5Flen(p)=5D

=3E =

=3E which should be =5Bseq=5Flen(p)=2C=5D

=3E =

=3E (otherwise=2C in the matrix case=2C the RHS will recycle only the 1=

st p

=3E elements=2C i=2Ee=2E=2C the 1st column)=2E

=3E =

=3E =3E =

=3E =3E Here is an example=3A

*=3E =3E =
*

=3E =3E hilbert =3C- function(n) =7B i =3C- 1=3An=3B 1 / outer(i -=

1=2C i=2C =22+=22) =7D

*=3E =3E h5 =3C- hilbert(5)
*

=3E =3E hinv1 =3C- solve(qr(h5))

=3E =3E hinv2 =3C- solve(qr(h5=2C LAPACK=3DTRUE)) =

=3E =3E all=2Eequal(hinv1=2C hinv2) =23 They are not equal

=3E =3E =

=3E =3E Here is a function that I wrote to correct this problem=3A

*=3E =3E =
*

=3E =3E solve=2Elapack =3C- function(A=2C LAPACK=3DTRUE=2C tol=3D1=2Ee=

-07) =7B

=3E =3E =23 A function to invert a matrix using =22LAPACK=22 or =22LIN=

**PACK=22
**

*=3E =3E if (nrow(A) !=3D ncol(A)) stop(=22Matrix muxt be square=
**=22)
**=3E =3E qrA =3C- qr(A=2C LAPACK=3DLAPACK=2C tol=3Dtol)
*

=3E =3E if (LAPACK) =7B

=3E =3E apply(diag(1=2C ncol(A))=2C 2=2C function(x) solve(qrA=2C x)) =

*=3E =3E =7D else solve(qrA)
*

=3E =3E =7D

=3E =3E =

*=3E =3E hinv3 =3C- solve=2Elapack(h5)
*

=3E =3E all=2Eequal(hinv1=2C hinv3) =23 Now=2C they are equal

=3E =3E =

*=3E =3E =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=
**=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=
**=5F
*

=3E =3E R-devel=40r-project=2Eorg mailing list

=3E =3E =

=3E =

=3E =

=3E -- =

=3E O=5F=5F ---- Peter Dalgaard =D8ster Farimagsgade 5=2C=

Entr=2EB

=3E c/ /=27=5F --- Dept=2E of Biostatistics PO Box 2099=2C 1014 C=

ph=2E K

=3E (*) =5C(*) -- University of Copenhagen Denmark Ph=3A (+45)=

35327918

=3E =7E=7E=7E=7E=7E=7E=7E=7E=7E=7E - (p=2Edalgaard=40biostat=2Eku=2Edk)=

** FAX=3A (+45) 35327907
**

=3E =

*=3E
*

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat 20 Jun 2009 - 11:52:32 GMT

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 20 Jun 2009 - 13:31:00 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.
*