Re: high-precision eigenvalue solver

In article <439F7DFB.5BA3314E@xxxxxxxxxxxx>, Julian V. Noble
<jvn@xxxxxxxxxxxx> wrote:

> wakun@xxxxxxxxx wrote:
> >
> > Thank you for all your suggestions. I found a library in c++ supporting
> > quadruple precision and the problem is solved at last with the help of
> > the lib. However, it takes almost 10 times longer than that in double
> > precision. I am searching a faster package. A guy told me lapack is
> > great. However, I found, from http://www.netlib.org/lapack/, lapack
> > only support double precision.
>
> The reason it takes 10x longer is that everything is done in high-level
> for portability. If you know how to use assembler you can hand-code the
> specific quad-precision arithmetic functions you need (you might look
> at Knuth for clues) and link them to your high-level calling program.
>
> That's something of a pain, but it will reduce the running time about
> 3-fold. You can't do better than that, because a quad-precision multi-
> plication uses 3 double-precision ones.
>
> The double-precision (aka 64-bit) operations take the same time
> as the (32-bit) single-precision ones because they are done on a math
> chip in 80-bit wide registers. (Since the instructions are overlapped
> on superscalar processors like the Pentia, the extra memory-access
> time does not increase this factor 3 by more than a few percent, even
> on systems with 32-bit bus width. Systems with 64-bit buses notice no
> difference resulting from memory accesses.)

Well, 80-bit on Intel. (And maybe AMD, I don't remember.) Never count
on more than 64-bit unless you're CERTAIN the code will run on bigger
registers.

Rather than hand-coding the routines, I would use a package like Gnu
Multi-Precision Library, and recode the lapack routines to call the
higher-precision routines. Be warned, this is painful, since lapack is
natively FORTRAN, and the f2c translation is pretty ugly. (And the
traceback of which routines calls which means there's a lot more to be
done than at first appears.)

Also, parts of lapack have 64-bit reals hard-wired into them. (Someday
somebody's going to have to do something about that.)

But the swox folks (who maintain GMP) have put a lot of effort into
optimizing the hand-coded a/l parts for Pentium (and, to a lesser
extent, 64-bit AMD). Use a 64-bit processor if you possibly can, the
throughput for multiplication is 4 times faster. The swox folks
explain this on their website (there's also a very interesting paper
there comparing Intel and AMD instruction timings).

GMP uses integer operations to do its fp operations, incidentally.

I have an implementation of the Haeberly-Overton algorithm (to minimize
the maximal eigenvalue of a "pencil" of positive semi-definite
matrices). The algorithm is (theoretically, and usually practically)
quadratically convergent. I occasionally calculate things to 32768
bits of precision. Since it's quadratically convergent, this isn't as
bad as it sounds: it's only a few more iterations. I built an option
into the program to do staged precision (start at 64-bit, run until
convergence, reset to 128-bit, run, reset to 256-bit, etc.) but this
isn't always stable, and it sometimes takes MANY, MANY more iterates.
But it's a good example of what can be done. It's THOUSANDS of times
faster than what can be done in CAS's like Mathematica and Maple, so
the pain of coding was worth it.

--Ron Bruck
.

Relevant Pages

  • Re: Ambiguous interfaces
    ... > I have a question about ambiguous interfaces. ... > subroutine DoSomething_DefaultInteger ... if you combine two or more specific routines ... default real precision and at double precision. ...
    (comp.lang.fortran)
  • Re: qfloat extended precision in Cephes library
    ... for Mathematical Functions (Ellis Horwood Series in Mathematics and Its ... precision offering in particular. ... qfloat.h C++ wrapper makes it easier to write math code that looks like ... Some of Mr. Moshier's built-in routines, ...
    (sci.math.num-analysis)
  • Re: Maples linear algebra routine and higher precision arithmetic
    ... precision results with results generated by increasing levels of precision. ... with linear algebra routines --- maybe by using NAG routines ... Digits. ... meaning that hardware is used if Digits <= evalhf(i.e. ...
    (sci.math.symbolic)
  • Re: Maples linear algebra routine and higher precision arithmetic
    ... precision results with results generated by increasing levels of precision. ... Sometime ago I recall hearing that Maple was changing that way it dealt ... with linear algebra routines --- maybe by using NAG routines ... Digits. ...
    (sci.math.symbolic)
  • Re: high-precision eigenvalue solver
    ... The algorithm is ... Since it's quadratically convergent, this isn't as ... > into the program to do staged precision (start at 64-bit, ... > convergence, reset to 128-bit, run, reset to 256-bit, etc.) but this ...
    (sci.math.num-analysis)