From pgilbert at bank-banque-canada.ca  Wed Apr 11 16:51:22 2007
From: pgilbert at bank-banque-canada.ca (Paul Gilbert)
Date: Wed, 11 Apr 2007 10:51:22 -0400
Subject: [Rd] eigen in beta
In-Reply-To: <Pine.LNX.4.64.0704110734210.1007@gannet.stats.ox.ac.uk>
References: <461BCE69.2020409@bank-banque-canada.ca>
	<Pine.LNX.4.64.0704101854560.32417@gannet.stats.ox.ac.uk>
	<461BD7E9.1050200@bank-banque-canada.ca>
	<461BEDCB.7090805@biostat.ku.dk>
	<Pine.LNX.4.64.0704110734210.1007@gannet.stats.ox.ac.uk>
Message-ID: <461CF5EA.4050001@bank-banque-canada.ca>

Hmmm.  It is a bit disconcerting that make check passes and I can get a 
fairly seriously wrong answer. Perhaps this test could be added to make 
check. Whether or not the one answer is correct may be questionable, but 
there is no question that

  prod(eigen(z, symmetric = FALSE, only.values = TRUE)$values ) *
  prod(eigen(solve(z), symmetric = FALSE, only.values = TRUE)$values )
   [1] 1.01677-0i

is wrong. (The product of the determinants should equal the determinant 
of the product, and the determinant of I is 1.)

On this machine I am using GNU Fortran (GCC 3.2.3 20030502 (Red Hat 
Linux 3.2.3-54)) 3.2.3 20030502 (Red Hat Linux 3.2.3-54),  which is a 
bit old but not really old. (I am using gcc 4.1.1 on my home machine, 
but the failing machine is suppose to be fairly new and "supported." 
(There is an OS upgrade planned.)  Should I really be thinking of this 
as an old compiler?  Is the compiler the most likely problem or is it 
possible I have a bad BLAS configuration, or something else? Previous 
versions of R have compiled without problems on this machine. (I am 
never very sure where to find all the information to report for a 
problem like this. Is there a simple way to get all the relevant 
information?)

Paul Gilbert

Prof Brian Ripley wrote:
> All the systems I tried this on give the 'correct' answer, including
>
> x86_64 Linux, FC5 (gcc 4.1.1)
> i686 Linux, FC5
> ix86 Windows (both gcc 3.4.5 and gcc pre-4.3.0)
> Sparc Solaris, with gcc3, gcc4 and SunPro compilers.
>
> Mainly with R 2.5.0 beta, some with R-devel (where the code is 
> unchanged).
>
> We have seen problems specific to RHEL's Fortran compilers on x86_64 
> several times before. I would strongly recommend compiler updates.
>
>
> On Tue, 10 Apr 2007, Peter Dalgaard wrote:
>
>> Paul Gilbert wrote:
>>> Here is the example. Pehaps others could check on other platforms. 
>>> It is
>>> only the first eigenvalue that is different.  I am relatively sure the
>>> old values are correct, since I compare with an alternate calculation
>>> using the expansion of a polynomial determinant.
>>>
>>>
>>> z <- t(matrix(c(
>>>  0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0064083373167516857,
>>> -0.14786612501440565826,  0.368411802235074137,
>>>  0, 0, 0, 0, 0, 0, 0, 0, 0,  0.0568624483195125444,
>>> 0.08575928008564302762, -0.101993668348446601,
>>>  0, 0, 0, 0, 0, 0, 0, 0, 0,  0.0039684327579889069,
>>> -0.00002857482925046247,  0.202241897806646448,
>>>  1, 0, 0, 0, 0, 0, 0, 0, 0, -0.0222834092601282285,
>>> -0.09126708346036176145,  0.644249961695308682,
>>>  0, 1, 0, 0, 0, 0, 0, 0, 0, -0.0032676036920228878,
>>> 0.16985862929849462888,  0.057282326361118636,
>>>  0, 0, 1, 0, 0, 0, 0, 0, 0,  0.0148488735227452068,
>>> -0.06175528918915401677,  0.109566197834008949,
>>>  0, 0, 0, 1, 0, 0, 0, 0, 0, -0.0392756265125193960,
>>> 0.04921079262665441212,  0.078176878215115805,
>>>  0, 0, 0, 0, 1, 0, 0, 0, 0, -0.0013937451966661973,
>>> 0.02009823693764142133, -0.207228935136287512,
>>>  0, 0, 0, 0, 0, 1, 0, 0, 0,  0.0273358858605219357,
>>> 0.03830466468488327725,  0.224426004034737836,
>>>  0, 0, 0, 0, 0, 0, 1, 0, 0, -0.1456426235151105919,
>>> 0.28688029213315069388,  0.326933845656016908,
>>>  0, 0, 0, 0, 0, 0, 0, 1, 0,  0.0164670122082246559,
>>> -0.21966261349875662590,  0.036404179329694988,
>>>  0, 0, 0, 0, 0, 0, 0, 0, 1,  0.0146156940584119890,
>>> 0.07505490943478997090,  0.077660578370038813
>>>  ), 12, 12))
>>>
>>>
>>> R-2.5.0 gives
>>> >  eigen(z, symmetric = FALSE, only.values = TRUE)$values
>>>  [1]  0.8465266+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>>>  [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>>>  [7]  0.6177419+0.0000000i -0.5604582+0.1958709i -0.5604582-0.1958709i
>>> [10]  0.1458799+0.4909300i  0.1458799-0.4909300i  0.3378356+0.0000000i
>>>
>>> R-2.4.1 and many, many previous versions gave
>>> >  eigen(z, symmetric = FALSE, only.values = TRUE)$values
>>>  [1]  0.8794798+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>>>  [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>>>  [7] -0.5604582+0.1958709i -0.5604582-0.1958709i  0.5847887+0.0000000i
>>> [10]  0.1458799+0.4909300i  0.1458799-0.4909300i  0.3378356+0.0000000i
>>>
>>> Sys.info()
>>>                              sysname                              
>>> release
>>>                              "Linux"                    
>>> "2.4.21-40.ELsmp"
>>>                              version                             
>>> nodename
>>> "#1 SMP Thu Feb 2 22:13:55 EST 2006"                           
>>> "mfa04559"
>>>                              machine
>>>                             "x86_64"
>>>
>>> Paul Gilbert
>> Hmm, I don't get that
>>
>>> version$version.string
>> [1] "R version 2.5.0 beta (2007-04-10 r41105)"
>>> eigen(z, symmetric = FALSE, only.values = TRUE)$values
>> [1]  0.8794798+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>> [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>> [7] -0.5604582+0.1958709i -0.5604582-0.1958709i  0.5847887+0.0000000i
>> [10]  0.1458799+0.4909300i  0.1458799-0.4909300i  0.3378356+0.0000000i
>>> Sys.info()
>>                              sysname                               
>> release
>>                              "Linux"                   
>> "2.6.20-1.2933.fc6"
>>                              version                              
>> nodename
>> "#1 SMP Mon Mar 19 11:38:26 EDT 2007"              
>> "titmouse2.kubism.ku.dk"
>>                              machine                                 
>> login
>>                               "i686"                                  
>> "pd"
>>                                 user
>>                                 "pd"
>>
>>
>>
>> And
>>
>>> version$version.string
>> [1] "R version 2.5.0 beta (2007-04-09 r41098)"
>>>  eigen(z, symmetric = FALSE, only.values = TRUE)$values
>> [1]  0.8794798+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>> [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>> [7] -0.5604582+0.1958709i -0.5604582-0.1958709i  0.5847887+0.0000000i
>> [10]  0.1458799+0.4909300i  0.1458799-0.4909300i  0.3378356+0.0000000i
>>> Sys.info()
>>                             sysname                              release
>>                             "Linux"               "2.6.18.8-0.1-default"
>>                             version                             nodename
>> "#1 SMP Fri Mar 2 13:51:59 UTC 2007"                              
>> "viggo"
>>                             machine                                login
>>                            "x86_64"                                 "pd"
>>                                user
>>                                "pd"
>>
>>
>>
>> The latter should be the actual build used in the current beta tarball
>> (which is what you used, right?).
>>
>> I would suspect one of the following:
>>
>> - RHEL compilers
>> - over-optimizing compiler settings
>> - system blas/lapack libraries
>> - system glibc
>>
>> but Brian probably has more concrete information.
>>
>>> Prof Brian Ripley wrote:
>>>
>>>
>>>> We are only aware of better behaviour from LAPACK 3.1 (which is what I
>>>> suppose you are talking about, that is R compiled with its internal
>>>> LAPACK).
>>>>
>>>> But in at least one case that means finding a complex set of
>>>> eigenvalues where previously a real one was found.
>>>>
>>>> On Tue, 10 Apr 2007, Paul Gilbert wrote:
>>>>
>>>>
>>>>> I am having some trouble with a case where  eigen in R-beta  gives a
>>>>> different largest value than in previous versions of R. Other values
>>>>> seem to be the same. Before I spend too much time, is anyone aware 
>>>>> of a
>>>>> problem (symmetric = FALSE, only.values = TRUE).
>>>>>
>>>>> Paul Gilbert
>>>>> ==================================================================================== 
>>>>>
>>>>>
>>>>>
>>>>> La version fran?aise suit le texte anglais.
>>>>>
>>>>> ------------------------------------------------------------------------------------ 
>>>>>
>>>>>
>>>>>
>>>>> This email may contain privileged and/or confidential
>>>>> inform...{{dropped}}
>>>>>
>>>>> ______________________________________________
>>>>> R-devel at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>>
>>>>>
>>> ==================================================================================== 
>>>
>>>
>>> La version fran?aise suit le texte anglais.
>>>
>>> ------------------------------------------------------------------------------------ 
>>>
>>>
>>> This email may contain privileged and/or confidential 
>>> inform...{{dropped}}
>>>
>>> ______________________________________________
>>> R-devel at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
====================================================================================

La version fran?aise suit le texte anglais.

------------------------------------------------------------------------------------

This email may contain privileged and/or confidential inform...{{dropped}}



[Rd] eigen in beta
Paul Gilbert