/*  compute eigenvalues and eigenvectors of a symmetric tridiagonal matrix */

/*  (c) 2003 by Gunter Winkler (guwi17@gmx.de) */

/*  This routine is free software; you can redistribute it and/or modify     */
/*  it under the terms of the GNU Lesser General Public License as published */
/*  by the Free Software Foundation; either version 2.1 of the License, or   */
/*  (at your option) any later version. */

/*  This program is distributed in the hope that it will be useful,    */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of     */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      */
/*  GNU General Public License for more details. */

/*  You should have received a copy of the GNU Lesser General Public    */
/*  License along with this program; if not, write to the Free Software */
/*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.           */

/*  important info: This routine was adapted from a numeric repository  */
/*  of the Chemnitz University (http://www.tu-chemnitz.de) Many people  */
/*  including Prof. Arndt Meyer and Dr. Matthias Pester developed this  */
/*  code. They were inspired by similar EISPACK routines.               */

/* ****************************************************************** */

/*     Documentation   */
/*     -------------   */

/* tridi_eigenvalues computes by using an implicite QZ-algorithm  */
/* the eigenvalues and eigenvectors of a tridiagonal symmetric matrix */

/*     Parameters:   

 *       nm      (I) :   leading dimension of ev[] 
 *       n       (I) :   dimension of diag[] and sub[] 
 *       diag[]  (I) :   diagonal elements of matrix 
 *               (O) :   eigenvalues 
 *       sub[]   (I) :   sub diagonal elements
 *       ev[]    (I) :   a few lines of a transformation matrix
 *                       column major storage assumed
 *                       set the first row of ev[] to (1,0,...,0) if you
 *                       need the first component of each eigenvector 
 *               (O) :   (first) nev elements of eigenvectors
 *       nev     (I) :   =k: compute only first nev elements of eigenvectors
 *       ierr    (O) :   error indicator 
 *                       number of the eigenvalue where the maximal
 *                       iteration count (30) was exceeded
 *                       only (ierr-1) eigenvalues are correct
 *       result      :   zero means ok, nonzero means error -> check ierr

 * ********************************************************************** */

#include "eigtridi.h"

/* declarations for incomplete C libraries */
#define max(x,y) (((x)>(y))?(x):(y))

int tridi_eigenvalues ( int nm, int n, double diag[], double sub[], 
                        double ev[], int nev, int* ierr )
{

  double  anorm, b, f, c, g, p, r, s, t, pa, psi, psj;
  int     i, ii, j, k, ki, l, m, mm1, mml, ki1, ka, ke;

  /* constants for double precision */
  static double EPS = 1e-16;
  static double EPSS = 1e-8;
  static double DINF = 1e20;
  static double DINF2 = 1e40;
  static double DMIN = 1e-20;

  *ierr = 0;

  anorm = 0.0;
  r = 0.0;

  /* move elements of sub one position and compute the largest 1-norm */
  for (i=1; i<n; ++i) {
    s=sub[i];
    sub[i-1]=s;
    s=fabs(s);
    p=fabs(diag[i-1])+r+s;
    if (p>anorm) anorm=p;
    r=s;
  }
  p = fabs(diag[n-1]) + r;
  if (p>anorm) anorm=p;
  sub[n-1]=0.0;

  for (l=0; l<n; ++l) {
    j = 0;  /* iteration counter */

    while (1) {
      /* ********** LOOK FOR SMALL SUBDIAGONAL ELEMENT ********** */
      for (m=l; m<n; ++m) {
        if (m==n-1) break;
        if (fabs(sub[m]) <= EPS*(fabs(diag[m])+fabs(diag[m+1]))) break;
        if (fabs(sub[m]) <= DMIN*anorm) break;
      }
      p = diag[l];
      mm1 = m-1;
      if (m==l) break; /* exit while loop */
      if (30==j) {
        *ierr=l;
        return -1; /* error: maximal number of iterations reached */
      }
      ++j;
      /* ********** FORM SHIFT ********* */
      g = (diag[l+1]-p) / (2.0*sub[l]);
      r = sqrt(g*g + 1.0);
      s = p - sub[l]/(g + ((g>=0)?(fabs(r)):(-fabs(r))) );
      if (m!=l+1) {
        t=s;
        r=max( fabs(s), anorm/n );
        for (i=1; i<=6; ++i) {
          pa = p;
          psi = diag[m] - t;
          psj = -1.0;
          for (k=mm1; k>=l; --k) {
            if (fabs(psi) >= EPS * fabs(sub[k])) {
              p = sub[k]/psi;
              psi = diag[k] - t - p*sub[k];
              psj = p*p*psj - 1.0;
            } else {
              psi = DINF;
              psj = DINF2;
            }
          }
          if (fabs(psj) <= DMIN) break;
          p = psi/psj;
          c = p;
          if (fabs(p) > 0.5*r) c=((p>0)?(fabs(r)):(-fabs(r)));
          t = t - c;
          if (fabs(p) <= EPSS*r) {
            s = t;
            break;
          }
        } /* for i */
      } /* if (m != l+1) */
      g = diag[m]-s;
      s=1.0;
      c=1.0;
      p=0.0;
      ki1 = m*nm;
      ka = ki1 - nm;
      ke = ka + nev - 1;
      /* compute c, s for rotation matrix */
      for (i=mm1; i>=l; --i) {
        f = s * sub[i];
        b = c * sub[i];
        if (fabs(f) <= fabs(g)) {
          s = f/g;
          r = sqrt(1.0 + s*s);
          sub[i+1] = g*r;
          c = 1.0/r;
          s = s*c;
        } else {
          c = g/f;
          r = sqrt(1.0+c*c);
          sub[i+1] = f*r;
          s = 1.0/r;
          c = c*s;
        }
        /* apply rotation to matrix */
        g = diag[i+1]-p;
        r = (diag[i]-g)*s + 2.0*c*b;
        p = s*r;
        diag[i+1] = g+p;
        g = c*r-b;
        /* apply rotation to eigenvectors */
        if (nev>0) {
          for (ki=ka; ki<=ke; ++ki,++ki1) {
            f = ev[ki1];
            b = ev[ki];
            ev[ki1] = s*b+c*f;
            ev[ki]  = c*b-s*f;
          }
          ki1=ka;
          ka=ka-nm;
          ke=ke-nm;
        }
      } /* for i */
      diag[l] = diag[l]-p;
      sub[l] = g;
      sub[m] = 0.0;
    }
  }
  return 0;
}