Fast Diag. Routines: RESPONSES

 Here is a summary of the response to my mesage about possible
 fast diagonalization methods, posted to the Computational Chemistry List
 and sci.math, sci.numerical-analysis, sci.physics, and sci.chem.
 Concerning 'fast' algorithms (ie faster than O(n^3) ) there don't seem
 to be any for non-sparse matrices. The routines most recommended are
 tred2 (Householder reduction of a symmetric matrix to a tridiagonal
 matrix) followed by tqli (finding eigenvalues, and eigenvectors if
 desired, of the tridiagonal matrix). These routines are given in
 Press et al, _Numerical Recipes_, available in FORTRAN, C, and Pascal
 versions, to my knowledge. (2nd ed.) One can purchase the programs
 on disk, and they may already be on whatever machine you currently use
 (but I would check to see if they are 1st ed.  and if there have been
 any problems with the 1st ed programs; i have no idea. Probably
 could just type in the necessary changes from the 2nd ed. book).
 I especially recommend that one read  Chapter 11 on Eigensystems, if they
 have not already done so.
 The eispack library has versions of these routines, and lapack also
 has versions of these routines. For access to these libraries, the
 full bibliography for Press, and other quips and references, see below.
 There are also one or two routines offered to me, which I have
 included. I have not tried any of these, so I would check with the
 person who sent the code (name and e-mail are provided).
 Keith Ball
 From: "M. J. Saltzman" <mjs()at()>
 In sci.math you write:
 >  I am wondering if anyone knows of any references to fast algorithms
 >for diagonalization (i.e. finding the eigenvalues and eigenvectors)
 >of a symmetric matrix.
 >  The problem I am using this for consists of finding stationary points
 >(specifically saddle points) on the 3*N-dimensional Cartesian potential
 >surface for a system of N particles interacting via pairwise interactions.
 >The matrix in question is the Hessian second-derivative matrix.
 >Any references or suggestions would be greatly appreciated!
 >Please send your responses directly by e-mail.
 Check out LAPACK at Netlib.  Send mail to netlib()at() with
 the message "send index".  --
 		Matthew Saltzman
 		Clemson University Math Sciences
 From: David Heisterberg <djh()at()>
 For the general case I don't know of anything that beats the old
 workhorse of Householder reduction followed by QL with implicit
 shifts, as implemented in Eispack (and now Lapack).  This is
 assuming you want all eigenvalues and vectors.
 If you are repeatedly solving for eigenvalues and vectors of matrices
 that change only a small amount at each step, then maybe an inverse
 iteration technique would work.
 David J. Heisterberg (djh()at()      Gee, it's so beautiful, I gotta
 The Ohio Supercomputer Center           give somebody a sock in the jaw.
 Columbus, Ohio                          -- Little Skippy (Percy Crosby)
 From: sling()at() (Song Ling)
 Hi Keith.
 I suggest that you look into existing program packages, one is
 called IMSL and another NAG, I think they are commercial.
 If you want to try it on your own, there are some fortran codes
 in a book called "Numerical Recipes" by Press et al, it is a very
 popular book, and I suggest that you find it at the bookstore of
 U of Chicago.
 		Song Ling
 From: hinsenk()at()ERE.UMontreal.CA (Hinsen Konrad)
 All linear algebra libraries provide routines for that, and unless
 your matrix has some exotic properties and/or you are looking for
 only a small subset of eigenvalues, there aren't any faster routines
 than the standard ones. My personal preference are the routines
 from the LAPACK library, which are free and in my experience more
 stable and faster than many others. You can get them from many
 sources, e.g. from AT&T's netlib (ftp to
 Konrad Hinsen                     | E-Mail: hinsenk()at()
 Departement de Chimie             | Tel.: +1-514-343-6111 ext. 3953
 Universite de Montreal            | Fax:  +1-514-343-7586
 C.P. 6128, succ. A                | Deutsch/Esperanto/English/Nederlands/
 Montreal (QC) H3C 3J7             | Francais (phase experimentale)
 From: mike()at()mycenae.CChem.Berkeley.EDU (Mike L. Greenfield)
 For calculating eigenvalues and eigenvectors, I recommend the
 combination of tred2 and tql2 from the EISPACK library.  You can
 obtain the Fortran source code by anonymous ftp to
 There are probably also LAPACK routines to do this; I don't know their
 names.  I think these two routines are also discussed in the
 eigenvalues section of Numerical Recipes.
 I assume you've looked at Jon Baker's paper (and references therein)
 on finding transition states?  He discusses a Cerjan-Miller type
 algorithm in detail:
   author = 	 "Jon Baker",
   title = 	 "An Algorithm for the Location of Transition States",
   journal = 	 "J. Comput. Chem.",
   year = 	 1986,
   volume = 	 7,
   pages = 	 "385--395"
 Good luck with your calculations.
 Mike Greenfield
 Dept. of Chemical Engineering
 UC Berkeley
 From: cooksj()at() (Stephen J. Cook)
 I'm not sure how fast they are, but have you looked at the
 netlib archives?  These can be found at either   (or)
 One of these names is obsolete, but I know that one of them
 works.  These machines hold linpack,lapack,eispack and other
 numerical algorithm suites that may be helpful to you.  They
 are all available in Fortran.
 Good luck.
 Steve Cook
 * Steve Cook                                cooksj()at()   *
 * Air Products and Chemicals, Inc.          Tel. (610) 481-2135     *
 * 7201 Hamilton Blvd.                       FAX  (610) 481-2446     *
 * Allentown, PA 18195                                               *
 * USA                                                               *
 *             Emacs - the choice of a GNU generation                *
 * Disclaimer: The opinions expressed here are those of the author.  *
 *             Any resemblance between my opinions and those of Air  *
 *             Products is purely coincidental...                    *
 From: bittner()at() (bittner)
 I use the LAPACK equlvalents to the EISPACK  tred2 and tql2 routines.
 All of which I know are on rainbow.uchicago.
 BTW, will their be a 2nd Keith Ball symposium for mathematical
 physics?  I forgot the subject for the first symposium...but it was
 great idea.
     Eric R. Bittner                  phone:  (512)-471-1092
     Dept. of Chemistry                 fax:  (512)-471-8696
     Univ. of Texas
   (formerly U of C)
    Wenn is das Nunstuck git und Slottermeyer? Ja!
  ...Beiherhund das Oder die Flipperwaldt gersput!
 From: Ryszard Czerminski X 217 <ryszard()at()>
 Hi Keith,
 I am not quite sure if this one will be fast enough for you,
 but it served me well over past ~10 years.
 Ryszard Czerminski
       IMPLICIT REAL*8 (A-H,O-Z)
 C     INPUT:
 C        Z  - Matrix to be diagonalized
 C        NZ - first dimension in calling program
 C        N  - order of Z
 C     OUTPUT:
 C        D - eigenvalues in decreasing order
 C        Z - genvectors
 C        ERROR=0 - normal solution
 C        ERROR=J - failure in J eigenvalue after 30 iterations
 C       MACHEP - relative machine precision
 C       SMALL  - smallest floating point number
       REAL*8 MACHEP
       real*8 D(N),E(N),Z(NZ,N)
 	real*8 b,c,f,g,h,p,r,s,small,tol
 c       DATA SMALL/5.4D-79/      IBM 370
 c       MACHEP=2.D0**(-56)       IBM 370
 C       DATA SMALL/2.225D-308/   titan, IBM PC real*8 ms-fortran 4.0
 C                                = 2**(-1022)
 C       MACHEP=2.D0**(-52)       titan, IBM pc real*8 ms-fortran 4.0
 C       DATA SMALL/0.495D-323/   sun : 2**(-1074)
 C       DATA SMALL/0.495D-323/   sun : 2**(-1074)
         DATA SMALL/0.495D-323/
 C       DABS(X)=DDABS(X)
       IF (N .EQ. 1) GO TO 320
       DO 300 II = 2, N
          I = N + 2 - II
          L = I - 2
          F = Z(I-1,I)
          G = 0.0E0
       IF (L .LT. 1) GO TO 140
       DO 120 K = 1, L
   120    G = G + Z(K,I) * Z(K,I)
   140    H = G + F * F
          IF (G .GT. TOL) GO TO 160
       E(I) = F
       H = 0.0E0
       GO TO 280
   160    L = L + 1
          G = -DSIGN(DSQRT(H),F)
          E(I) = G
          H = H - F * G
          Z(I-1,I) = F - G
          F = 0.0E0
          DO 240 J = 1, L
             Z(I,J) = Z(J,I)
             Z(J,I) = Z(J,I) / H
             G = 0.0E0
             DO 180 K = 1, J
   180       G = G + Z(K,J) * Z(I,K)
             JP1 = J + 1
             IF (L .LT. JP1) GO TO 220
             DO 200 K = JP1, L
   200      G = G + Z(J,K) * Z(K,I)
   220       E(J) = G / H
             F = F + G * Z(J,I)
   240      CONTINUE
             HH = F / (H + H)
             DO 260 J = 1, L
               F = Z(I,J)
               G = E(J) - HH * F
             E(J) = G
             DO 260 K = 1, J
                Z(K,J) = Z(K,J) - F * E(K) - G * Z(I,K)
   260      CONTINUE
   280      D(I) = H
   320 D(1) = 0.0E0
       E(1) = 0.0E0
       DO 500 I = 1, N
          L = I - 1
          IF (D(I) .EQ. 0.0E0) GO TO 380
          DO 360 J = 1, L
             G = 0.0E0
             DO 340 K = 1, L
   340       G = G + Z(I,K) * Z(K,J)
             DO 360 K = 1, L
                Z(K,J) = Z(K,J) - G * Z(K,I)
   360    CONTINUE
   380     D(I) = Z(I,I)
          Z(I,I) = 1.0E0
          IF (L .LT. 1) GO TO 500
          DO 400 J = 1, L
             Z(I,J) = 0.0E0
             Z(J,I) = 0.0E0
   400    CONTINUE
       ERROR = 0
       IF (N .EQ. 1) GO TO 1500
       DO 1100 I = 2, N
  1100 E(I-1) = E(I)
       F = 0.0E0
       B = 0.0E0
       E(N) = 0.0E0
       DO 1240 L = 1, N
          J = 0
          H = MACHEP * (DABS(D(L)) + DABS(E(L)))
          IF (B .LT. H) B = H
          DO 1110 M = L, N
             IF (DABS(E(M)) .LE. B) GO TO 1120
  1110      CONTINUE
  1120      IF (M .EQ. L) GO TO 1220
  1130      IF (J .EQ. 30) GO TO 1400
          J = J + 1
           P = (D(L+1) - D(L)) / (2.0E0 * E(L))
           R = DSQRT(P * P + 1.0E0)
           H = D(L) - E(L) / (P + DSIGN(R,P))
           DO 1140 I = L, N
  1140    D(I) = D(I) - H
           F = F + H
           P = D(M)
           C = 1.0E0
           S = 0.0E0
           MML = M - L
           DO 1200 II = 1, MML
             I = M - II
             G = C * E(I)
       H = C * P
       IF (DABS(P) .LT. DABS(E(I))) GO TO 1150
             C = E(I) / P
       R = DSQRT(C * C + 1.0E0)
             E(I+1) = S * P * R
             S = C / R
             C = 1.0E0 / R
             GO TO 1160
  1150        C = P / E(I)
             R = DSQRT(C * C + 1.0E0)
             E(I+1) = S * E(I) * R
             S = 1.0E0 / R
             C = C / R
  1160        P = C * D(I) - S * G
             D(I+1) = H + S * (C * G + S * D(I))
             DO 1180 K = 1, N
                H = Z(K,I+1)
                Z(K,I+1) = S * Z(K,I) + C * H
                Z(K,I) = C * Z(K,I) - S * H
  1180        CONTINUE
  1200    CONTINUE
          E(L) = S * P
          D(L) = C * P
          IF (DABS(E(L)) .GT. B) GO TO 1130
  1220    D(L) = D(L) + F
       NM1 = N - 1
       DO 1300 I = 1, NM1
          K = I
          P = D(I)
             IP1 = I + 1
          DO 1260 J = IP1, N
             IF (D(J) .LE. P) GO TO 1260
             K = J
             P = D(J)
  1260    CONTINUE
           IF (K .EQ. I) GO TO 1300
           D(K) = D(I)
           D(I) = P
           DO 1280 J = 1, N
             P = Z(J,I)
             Z(J,I) = Z(J,K)
             Z(J,K) = P
  1280    CONTINUE
       GO TO 1500
  1400 ERROR = L
  1500 RETURN
 	subroutine horder(a,e,n,ch)
 	character*1 ch
 	real*8        a(n,n),e(n),x
 c       ch = 'i' - increasing order of eigenvalues
 c       ch = 'd' - decreasing order of eigenvalues
 	do 30 i=1,n-1
 	do 20 j=i+1,n
 	if(e(i).gt.e(j)) then
 		if(ch.eq.'i') then
 			do 10 k=1,n
 			x      = a(k,i)
 			a(k,i) = a(k,j)
 			a(k,j) = x
 			x      = e(i)
 			e(i)   = e(j)
 			e(j)   = x
 10			continue
 20	continue
 30	continue
 From: dickson()at() (Ross M. Dickson)
 Have you looked into EISPACK or its successor, LAPACK?  They're
 public-domain FORTRAN libraries for eigenproblems.  Although I'm not
 intimately familiar with the algorithm available therein for real
 symmetric matrices, I know there is code for that special case.  The
 developers tried to choose the most efficient yet reliable algorithms
 possible for each matrix type covered.
 They can be obtained as a whole or in pieces through NETLIB:  For info,
 send a message consisting solely of the word 'help' to netlib()at()
 Ross Dickson,  dickson()at()
 Department of Chemistry, The University of Calgary, Alberta, Canada
 From: ucacsam <ucacsam()at()>
 The Householder method is among the best you can find.
 See texts by J H Wilkinson, L Fox and many others for descriptions.
 Basically the method calculates orthogonal matrices that are used to
 transform your original matrix to a (symmetric) tridiagonal form, which for
 an n x n matrix can be done in about 2/3 of the cost of one complete
 matrix multiplication (of two n x n matrices - that's very clever!!). From
 then on the process uses Laguerre's method for finding the eigenvalues and
 you can the vectors as well. All rapidly convergent and numericaly stable.
 A simpler method, which takes rather longer is known as Jacobi's method.
 Have fun.
 Paul Samet
 From: saj()at() (Stephen Jacobs)
 How big?  Up to about order 1000 brute force (and maybe some sparse
 storage techniques) seems adequate.  Have you looked at "Numerical
 Recipes"?  It gives references to the stuff it leaves out.  In the
 most general case, you're screwed anyway because the problem gets
 so badly ill-conditioned.  My personal favorite matrix book is
 Golub and Van Loan (title "Matrix Computations"); but it doesn't go
 into extreme practical details.
 From: rossi()at() (Ivan Rossi)
 Keith, I think the absolute best around are the LAPACK routines from
 J. Dongarra They use BLAS level 3 wherever possible, and if you use a
 workstation or a Cray BLAS come optimized for maximum performance
 (libblas on IBM, SGI; libsci on Cray C90) .  You can get the source
 from netlib send a message to netlib()at() or saying
 send index from lapack
 or search the WWW site at
 happy computing
 Ivan Rossi                                |
 Department of Chemistry                   |   EXPERT : (n)
 University of Minnesota                   |   someone who avoids all the
 207 Pleasant St. SE,Minneapolis, MN 55455 |   little errors, going straight
 Tel. +1-612-624-6099                      |   towards the catastrophe
 e-mail : rossi()at()            |
 				          |   Intel Inside : the warning label
 From: Bob Funchess <bobf()at()>
 Dear Keith,
    There's a fairly fast method for sparse matrices called the "Lanczos
 Algorithm"; I'm not at my reference papers right now so I can't give you
 the exact reference.  The original method had some numerical instabilities,
 but those are correctable... I know that the book "Biological Magnetic
 Resonance, vol. 8: spin labeling: theory and applications" by academic
 press (L. Berliner, ed.; pub. date 1988 or 1989) contains an EPR simulation
 program which uses this method.  The book comes with a DOS disk that has
 the source code to the program, and the chapter which discusses the program
 contains the full references to the original algorithm and the correction
           Bob Funchess
 Dr. Robert B. Funchess                    bobf()at()
 Senior Support Scientist                  Tel (617) 229-9800 x202
 Molecular Simulations Inc.                FAX (617) 229-9899
 16 New England Executive Park   
 Burlington, MA 01803
 From: Robert Fraczkiewicz <robert()at()>
 	I use Fortran routines DSYEV.F and DSPEV.F from LAPACK package
 available from NETLIB (a repository of mathematical software). If you
 have a WWW client try the following URL:
 Eventually you can ftp anonymously to; directory netlib/lapack.
 Some machines have specifically optimized LAPACK library; using it is much
 better than transferring source code. On UNIX machine you would type command
 "man lapack" to see if your sysadm installed one.
 I am also interested in high-speed diagonalization; if you get better answers,
 please, let me know and/or post the summary to CCL. Thanks in advance !
 Happy computing,
 +                               |\+_+/\+_+/|                                 +
 +   just Robert Fraczkiewicz    |=\ /||\ /=|    robert()at()
 +    Deparment of Chemistry     |--\-||-\--|      chem86()at()
 +    University of Houston      |=/_\||/_\=|         (713) 743-3236          +
 +                               |/+ +\/+ +\|                                 +
 +                                                                            +
 +   "It is also a good rule not to put too much confidence in experimental
 +    results until they have been confirmed by Theory"
 +                                                                            +
 +                                                     Sir Arthur Eddington   +
 +                                                                            +