1 /* qrecvc.c CCMATH mathematics library source code.
2 *
3 * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
4 * This code may be redistributed under the terms of the GNU library
5 * public license (LGPL). ( See the lgpl.license file for details.)
6 * ------------------------------------------------------------------------
7 */
8 #include "complex.h"
qrecvc(double * ev,Cpx * evec,double * dp,int n)9 void qrecvc(double *ev,Cpx *evec,double *dp,int n)
10 { double cc,sc,d,x,y,h,tzr=1.e-15;
11 int i,j,k,m,nqr=50*n;
12 Cpx *p;
13 for(j=0,m=n-1;j<nqr;++j){
14 while(1){
15 if(m<1) break;
16 k=m-1;
17 if(fabs(dp[k])<=fabs(ev[m])*tzr) --m;
18 else{ x=(ev[k]-ev[m])/2.; h=sqrt(x*x+dp[k]*dp[k]);
19 if(m>1 && fabs(dp[m-2])>fabs(ev[k])*tzr) break;
20 if((cc=sqrt((1.+x/h)/2.))!=0.) sc=dp[k]/(2.*cc*h); else sc=1.;
21 x+=ev[m]; ev[m--]=x-h; ev[m--]=x+h;
22 for(i=0,p=evec+n*(m+1); i<n ;++i,++p){
23 h=p[0].re; p[0].re=cc*h+sc*p[n].re;
24 p[n].re=cc*p[n].re-sc*h;
25 h=p[0].im; p[0].im=cc*h+sc*p[n].im;
26 p[n].im=cc*p[n].im-sc*h;
27 }
28 }
29 }
30 if(x>0.) d=ev[m]+x-h; else d=ev[m]+x+h;
31 cc=1.; y=0.; ev[0]-=d;
32 for(k=0; k<m ;++k){
33 x=ev[k]*cc-y; y=dp[k]*cc; h=sqrt(x*x+dp[k]*dp[k]);
34 if(k>0) dp[k-1]=sc*h;
35 ev[k]=cc*h; cc=x/h; sc=dp[k]/h; ev[k+1]-=d; y*=sc;
36 ev[k]=cc*(ev[k]+y)+ev[k+1]*sc*sc+d;
37 for(i=0,p=evec+n*k; i<n ;++i,++p){
38 h=p[0].re; p[0].re=cc*h+sc*p[n].re;
39 p[n].re=cc*p[n].re-sc*h;
40 h=p[0].im; p[0].im=cc*h+sc*p[n].im;
41 p[n].im=cc*p[n].im-sc*h;
42 }
43 }
44 ev[k]=ev[k]*cc-y; dp[k-1]=ev[k]*sc; ev[k]=ev[k]*cc+d;
45 }
46 }
47