librheolef-dev 6.5-1build1 / usr / include / rheolef / tqli.h

///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
#include "rheolef/compiler.h"
#include <iterator>
#include <iostream>
#include "rheolef/pythag.h"


template <class T>
inline T SIGN(T a, T b) { return b >= T(0) ? std::fabs(a) : -std::fabs(a); }

template <class Iterator1, class Iterator2, class Size>
void tqli (Iterator1 d, Iterator2 e, Size n)
{
	typedef typename std::iterator_traits<Iterator1>::value_type T;
	Size m,l,iter,i;
	T s,r,p,g,f,dd,c,b;
	for (i=2;i<=n;i++) e[i-1]=e[i];
	e[n]=0.0;
	for (l=1;l<=n;l++) {
	    iter=0;
	    do {
	      for (m=l;m<=n-1;m++) {
		dd=std::fabs(d[m])+std::fabs(d[m+1]);
		if ((T)(std::fabs(e[m])+dd) == dd) break;
	      }
	      if (m != l) {
		iter++;
		if (iter == 30) { std::cerr << "Too many iterations in tqli\n"; exit(1); }
		g=(d[l+1]-d[l])/(2.0*e[l]);
		r=pythag(g,T(1));
		g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
		s=c=1.0;
		p=0.0;
		for (i=m-1;i>=l;i--) {
			f=s*e[i];
			b=c*e[i];
			e[i+1]=(r=pythag(f,g));
			if (r == T(0)) {
				d[i+1] -= p;
				e[m]=0.0;
				break;
			}
			s=f/r;
			c=g/r;
			g=d[i+1]-p;
			r=(d[i]-g)*s+2.0*c*b;
			d[i+1]=g+(p=s*r);
			g=c*r-b;
		}
		if (r == T(0) && i >= l) continue;
		d[l] -= p;
		e[l]=g;
		e[m]=0.0;
	      }
	    } while (m != l);
	}
}