Skip to content

eigen

19 messages · Paul Gilbert, Douglas Bates, Ross Ihaka +5 more

Original
Paul Gilbert
# 
Eigen seems to be giving different results in 0.64 :

eigen(matrix(
 c(0.00000000,  0.0000000,  0.00000000, 1.000000000,
   0.01262442, -0.1141633,  0.04989433, 0.006112607,
   0.04781087,  0.4177869, -0.15569656, 0.074242237,
   0.20167331,  0.2531109, -0.04973573, 0.034046467) ,4,4))$values

R:
 [1]  0.35519256 -0.29294606 -0.26637369 -0.03168619

S:
[1]  0.4689288205 -0.4246063863 -0.2810596564  0.0009238293

I am fairly sure that S and R previously gave the same results.

Paul Gilbert

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Peter Dalgaard
# 
Paul Gilbert <pgilbert@bank-banque-canada.ca> writes:
Indeedy. At least the pre 63.3 binary I had lying about did give the S
results. The eigenvectors are also bogus as far as I can tell.

The weird thing is that eigen() goes via .Fortran to eigen.c
(which was f2c'ed at some point) and neither eigen.c or eigen.R has
been changed in 1999 according to the cvs logs.
Ross Ihaka
#
On 10 Apr 1999, Peter Dalgaard BSA wrote:

            

      
      
      
    

        
Weirder still.  I get the problem results under Solaris with a clean
checkout of 0.65 but not under FreeBSD with a clean 0.65.
Solaris is using egcs and FreeBSD is using gcc 2.7.2.1.

I am beginning to suspect a compiler / library bug.

	Ross

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
I also get the problem results under 0.64.0 compiled with egcs on
Linux/Intel - both RedHat and Debian systems.  When I checked the
RedHat system under 0.63.4 (not sure of the compiler) I got the other
set of results.  That pattern was repeated on a SPARC Solaris 2.6
system.  0.64.0 compiled with egcs gave the incorrect results and
0.63.4 compiled with an unknown compiler gave the correct results.
Same again on an Intel/Solaris 2.6 system.

I think the compiler that I used for the 0.63.4 builds was egcs-1.0.3
but I'm not sure.  I'll go back and re-build 0.64.0 with egcs-1.0.3
instead of egcs-1.1.1 and see what happens.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

          

      
      
      
    

        
DMB> I think the compiler that I used for the 0.63.4 builds was
  DMB> egcs-1.0.3 but I'm not sure.  I'll go back and re-build 0.64.0
  DMB> with egcs-1.0.3 instead of egcs-1.1.1 and see what happens.

Changing to egcs-1.0.3 for a re-build of 0.64.0 didn't clear up the
problem.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Ross Ihaka
    

    

      
      #
    

  


  

      
I don't know what the problem was, but I think I have a fix.  I have
replaced src/appl/eigen.c by the original set of fortran eispack routines. 
This seems to fix the problem on my Solaris box. 

I suspect that this bug may have been there for a while.

I am still puzzled that it was not tripped on my FreeBSD box.  It must be
the advanced compiler technology; it silently finds and fixes the bugs in
your code as you compile it :-).

Anyway could someone check that this change does the job on Linux etc.
(I've committed in both branches).

	Ross

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On Sat, 10 Apr 1999, Ross Ihaka wrote:

        

            

      
      
      
    

        
The bug is that Fortran.h had <math.h>, now has "Mathlib.h", and soemthing
in there is trampling over error.c

            

      
      
      
    

        
No, it _is_ a 0.64 change, 0.63.3 is unaffected.

            

      
      
      
    

        
Look to header files in such circumstances. Mathlib.h includes

#include <errno.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>

and on Solaris the culprit is the last one of these.

Guess what: there a load of abs/fabs errors in eigen.c.  On some
systems abs is a macro, on others a function and on others it
depends on the optimization level.

            

      
      
      
    

        
(Do you want to revert?)

What worries me is what else Mathlib.h might be trampling over.
I think we need to check this over.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ross Ihaka
    

    

      
      #
    

  


  

      
On Sat, 10 Apr 1999, Prof Brian D Ripley wrote:

        

            

      
      
      
    

        
I think eigen.c dates from days before we tried to be pure C for
portability purposes.  Provided we can live with the namespace pollution
which comes with Fortran (all function names being global), I'd actually
prefer to live with the original code. We've reverted the Linpack and
optimizer code already.

(Yeah, Yeah so why did I do the nmath thing?  I thought it might be more
generally useful than just for R ...)

            

      
      
      
    

        
Definitely!

	Ross

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
Douglas Bates <bates@stat.wisc.edu> writes:

            

      
      
      
    

        
I also get the same result with different compilers (gcc 2.7.2.3 and
egcs-2.91.57) Note BTW that it *fails* the
consistency test from the example file

            

      
      
      
    

        
[,1]          [,2]          [,3]          [,4]
[1,]  2.775558e-17  1.387779e-17  2.775558e-17 -3.382711e-17
[2,]  4.117019e-01  2.872937e-01 -2.687366e-01  7.624831e-02
[3,]  4.615645e-02  3.220889e-02 -3.012842e-02  8.548301e-03
[4,] -2.688891e-01 -1.876361e-01  1.755162e-01 -4.979899e-02

- whereas the older version did

            

      
      
      
    

        
[,1]          [,2]          [,3]          [,4]
[1,] -8.326673e-17  0.000000e+00  0.000000e+00 -2.622414e-18
[2,] -5.551115e-17  1.110223e-16  0.000000e+00 -1.368263e-16
[3,]  3.469447e-18 -2.775558e-17 -5.551115e-17  3.957338e-17
[4,] -5.551115e-17 -1.110223e-16  5.028530e-16  1.502704e-16

[0.63.4 in progress (March 9, 1999) on my office machine]

The test on symmetric matrices does not fail, though. (Make check
would have caught this, save for the m<-m+t(m) line...)

Most likely, we're being bitten  either by a code 'cleanup' or some
sort of change in the linked libraries.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
Peter Dalgaard BSA <p.dalgaard@biostat.ku.dk> writes:

            

      
      
      
    

        
Oops. Didn't see Brian+Ross's exchange before sending. Never mind.
Martin may still want to fix the example, though?

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ross Ihaka
    

    

      
      #
    

  


  

      
The problem in eigen.c seems to be as follows.

1) f2c defines a macro abs() which does the standard abs() which
   works for all numeric types. This macro was inserted in Fortran.h.

2) recently it was recognised that there is an abs() in stdlib and
   the macro in Fortran.h was commented out.

3) unfortunately the abs() in stdlib is an integer valued function
   of an integer argument.

4) end of world.

I think the lesson is that partial translations of fortran into C and the
use of macros to get the rest of the way is VERY dangerous. 

Sometimes it is unavoidable, because we want to do some bit-fiddling or
platform dependent manipulation.  The lesson I think is to keep the use of
externally defined macros in the translations to a minimum. 

I am going to go through appl/ to look carefully at all the translations
for problems of this kind.

	Ross

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Ross Ihaka
    

    

      
      #
    

  


  

      
To follow up the "eigen" issue.  What should our position on translating
Fortran into C be?  Should we leave as much code as possible in Fortran or
should we move it to C as much as possible? 

I can see arguments both ways.

	Ross

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

          

      
      
      
    

        
Ross> To follow up the "eigen" issue.  What should our position on
  Ross> translating Fortran into C be?  Should we leave as much code
  Ross> as possible in Fortran or should we move it to C as much as
  Ross> possible?

I can think of three different sets of circumstances:
 a) creating new code: This one is obvious, just say "no" to Fortran.
 b) enhancing existing code: For something small do the automatic
      translation to C with f2c then modify the C sources and regard
      those as the definitive version.  For a big Fortran routine, try
      to write a C wrapper around it and make modifications in the
      wrapper. 
 c) using established code without modification (e.g. Eispack, Linpack):
      Regard the Fortran sources as the definitive version and don't
      translate.

Now that g77 is available and reasonably stable, a lot of the
motivation for using f2c is gone.  There are still the messy linkage
issues when calling Fortran from C but most of those are now resolved
by the configure script.  Passing character strings, etc. will
continue to be problematic but most of the legacy Fortran code we
would use will be strictly numerical routines with, we hope,
comparatively simple calling sequences.  (Lapack is an exception in
that it uses strings to declare options but we aren't using Lapack.)

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On 10 Apr 1999, Douglas Bates wrote:

        

            

      
      
      
    

        
[A long reply: I am much-bitten here!]

1) In my experience there are problems with error reporting via strings in
such Fortran code, as I found in loess (and found more easily in less
complicated packages). I wrote a two-step interface (in Fortran find the
string length, pass that and string back to C), and perhaps we should
set up warning and error subroutines to be called from Fortran to
standardize this that package writers/porters could use.


2) My motivation would be to use Fortran wherever it already exists. In the
world of scientific computing Fortran has been king and the effort in
compiler development has gone into efficient Fortran compilers with godd
optimizers. So on a Sun or DEC Alpha (and probably SGI) the original
Fortran used to run appreciably faster than the C translation (using the
native compilers), and although the gap has closed I think it is still
quite noticeable.  In my experience R can't afford that price. This may not
be an issue on Linux/Windows, but I expect is on Sun/DEC/SGI/HP.

BTW, I think that means that Doug's `just say "no" to Fortran' is not
obviously true if performance matters. This is likely to especially true if
single-precision calculations will do the job. (And, yes, I know about
Intel hardware but I don't see the need to slow a Sparc or Alpha FP unit
down to Intel speeds, especially older ones. I also know how to make C do
efficient float calculations using float libraries on Sparc, but that sort
of thing is a nightmare to configure portably, and also to link into S or 
R. It happens for free in Fortran.)


3) The problem is finding the external calls from Fortran code. In the base
code this should be not problem now we use the Fortran front-end to link.  
I have problems with quite a few packages that use Fortran, and have
translated some code myself just to avoid those. (It was a big problem in
the first release of MASS software for S-PLUS, and we removed all the 
Fortran to save ourselves the hastle of `it won't work on my XYZ system'
reports. We never made the Fortran work on HPs, whose Fortran name space
clashes with the C one, and it seemed that S had already grabbed some
critical symbols as C ones.)

I think in R we can do better than at present, if only I knew
exactly how. One idea is to have (effectively) a `hints' file that
configure uses to add a base set of libraries to FLIBS, but on Solaris
at least, using the Fortran command (f77 or g77) to link a dynamic library
does the job.  For linking with ld the base FLIBS would need to be

-lF77 -lM77 -lV77 -lsunmath -lm

_if_ f77 is used.  How viable is `use FC to link' on other platforms?

And to come back to beating a drum, if the dlopen option were RTLD_NOW
not RTLD_LAZY, these things would be easier to check. People keep sending
packages to CRAN that are incomplete (mclust seems to be now) and with
RTLD_NOW this is _much_ more likely to be noticed. Can we re-consider?

Brian

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ross Ihaka
    

    

      
      #
    

  


  

      
On Sun, 11 Apr 1999, Prof Brian D Ripley wrote:

        

            

      
      
      
    

        
I've changed in the devel (non) branch.  Let's see if it bites anyone 
(HPUX?). 

	Ross

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Kurt Hornik
    

    

      
      #
    

  


  

          

      
      
      
    

        
Ross> To follow up the "eigen" issue.  What should our position on
Ross> translating Fortran into C be?  Should we leave as much code
Ross> as possible in Fortran or should we move it to C as much as
Ross> possible?

            

      
      
      
    

        
Are you suggesting using the Fortran compiler as SHLIBLD?

Please help me understanding this a bit more.  The above `-lF77 ...'
would be the default FLIBS to use when compiling/linking Fortran?

Short term fix should be to use platform-specific overrides in basic
configure.  Another idea is to have configure scripts supported for add
ons.

(If I fully understood the problem, I think I could fix it ...)

=k
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On Mon, 12 Apr 1999, Kurt Hornik wrote:

        

            

      
      
      
    

        
Yes, at least on platforms where this is known to be safe.

            

      
      
      
    

        
On Solaris.

            

      
      
      
    

        
Yes!

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Martin Maechler
    

    

      
      #
    

  


  

          

      
      
      
    

        

        
BDR> On Sat, 10 Apr 1999, Ross Ihaka wrote:

        
>> I don't know what the problem was, but I think I have a fix.  I
      >> have replaced src/appl/eigen.c by the original set of fortran
      >> eispack routines.  This seems to fix the problem on my Solaris
      >> box.

    BDR> The bug is that Fortran.h had <math.h>, now has "Mathlib.h", and
    BDR> something in there is trampling over error.c

    >> I suspect that this bug may have been there for a while.

    BDR> No, it _is_ a 0.64 change, 0.63.3 is unaffected.

    >> 
    >> I am still puzzled that it was not tripped on my FreeBSD box.  It must be
    >> the advanced compiler technology; it silently finds and fixes the bugs in
    >> your code as you compile it :-).

    BDR> Look to header files in such circumstances. Mathlib.h includes

    BDR> #include <errno.h>
    BDR> #include <limits.h>
    BDR> #include <float.h>
    BDR> #include <math.h>
    BDR> #include <stdlib.h>

    BDR> and on Solaris the culprit is the last one of these.
(grrreat)

    BDR> Guess what: there a load of abs/fabs errors in eigen.c.  On some
    BDR> systems abs is a macro, on others a function and on others it
    BDR> depends on the optimization level.

We must use a Mathlib.h -- src/nmath/  version of this

    >> Anyway could someone check that this change does the job on Linux etc.
    >> (I've committed in both branches).

    BDR> (Do you want to revert?)

Yes! -- The eigen.c has been massaged etc to be quite a bit more readable
     than its Fortran counterpart..


    BDR> What worries me is what else Mathlib.h might be trampling over.
    BDR> I think we need to check this over.

certainly
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Albrecht Gebhardt
    

    

      
      #
    

  


  

      
On Sun, 11 Apr 1999, Prof Brian D Ripley wrote:

        

            

      
      
      
    

        
I'm not sure if it does matter here, but the only way to build R on my
Alphas currently is using Fortran as LDCMD:

in config.site:
LDCMD="f77 -g -nofor_main -fpe3 "

SHLIBLD is "ld -shared -expect_unresolved '*'".



Albrecht

-------------------------------------------------------------------------------
Albrecht Gebhardt                   email   : albrecht.gebhardt@uni-klu.ac.at
Institut fuer Mathematik            Tel.    : (++43 463) 2700/837
Universitaet Klagenfurt             Fax     : (++43 463) 2700/834
Villacher Str. 161
A-9020 Klagenfurt, Austria
-------------------------------------------------------------------------------


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._