nr.c
Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 #include "eus.h"
00026 #include "nr.h"
00027 #include <math.h>
00028 
00029 static eusfloat_t sqrarg;
00030 #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
00031 
00032 static eusfloat_t maxarg1, maxarg2;
00033 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
00034 
00035 static int iminarg1, iminarg2;
00036 #define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))
00037 
00038 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
00039 
00040 #define NR_END 1
00041 #define FREE_ARG char*
00042 
00043 void nrerror(char error_text[])
00044 {
00045   fprintf(stderr,"Numerical Recipes run-time error...\n");
00046   fprintf(stderr,"%s\n",error_text);
00047   fprintf(stderr,"...now existing to system...\n");
00048 }
00049 
00050 eusfloat_t *nr_vector(int nl, int nh)
00051 {
00052   eusfloat_t *v;
00053   v =(eusfloat_t *)malloc((size_t)((nh-nl+1+NR_END)*sizeof(eusfloat_t)));
00054   if (!v) {nrerror("allocation failure in nr_vector()"); return (eusfloat_t*)NULL;}
00055   return v-nl+NR_END;
00056 }
00057 
00058 eusfloat_t **nr_matrix(int nrl, int nrh, int ncl, int nch)
00059 {
00060   int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
00061   eusfloat_t **m;
00062   
00063   m = (eusfloat_t **)malloc((size_t)((nrow+NR_END)*sizeof(eusfloat_t *)));
00064   if (!m) {nrerror("allocation failure 1 in nr_matrix()"); return (eusfloat_t**)NULL;}
00065   m += NR_END;
00066   m -= nrl;
00067 
00068   m[nrl]=(eusfloat_t *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(eusfloat_t)));
00069   if (!m[nrl]) {nrerror("allocation failure 2 in nr_matrix()"); return (eusfloat_t**)NULL;}
00070   m[nrl] += NR_END;
00071   m[nrl] -= ncl;
00072 
00073   for (i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
00074 
00075   return m;
00076 }
00077 
00078 void free_nr_vector(eusfloat_t *v, int nl, int nh)
00079 {
00080   free((FREE_ARG)(v+nl-NR_END));
00081 }
00082 
00083 void free_nr_matrix(eusfloat_t **m, int nrl, int nrh, int ncl, int nch)
00084 {
00085   free((FREE_ARG)(m[nrl]+ncl-NR_END));
00086   free((FREE_ARG)(m+nrl-NR_END));
00087 }
00088 
00089 void lubksb(eusfloat_t **a, int n, int *indx, eusfloat_t b[]) {
00090   int i,ii=0,ip,j;
00091   eusfloat_t sum;
00092 
00093   for (i=1;i<=n;i++) {
00094     ip=indx[i];
00095     sum=b[ip];
00096     b[ip]=b[i];
00097     if (ii)
00098       for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
00099     else if (sum) ii=i;
00100     b[i]=sum;
00101   }
00102   for (i=n;i>=1;i--) {
00103     sum=b[i];
00104     for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
00105     b[i]=sum/a[i][i];
00106   }
00107 }
00108 
00109 int ludcmp(eusfloat_t **a, int n, int *indx, eusfloat_t *d) {
00110   int i,imax,j,k;
00111   eusfloat_t big,dum,sum,temp;
00112   eusfloat_t *vv;
00113 
00114   vv=nr_vector(1,n);
00115   *d=1.0;
00116   for (i=1;i<=n;i++) {
00117     big=0.0;
00118     for (j=1;j<=n;j++)
00119       if ((temp=fabs(a[i][j])) > big) big=temp;
00120     if (big == 0.0) { free_nr_vector(vv,1,n); return -1; }
00121     vv[i]=1.0/big;
00122   }
00123   for (j=1;j<=n;j++) {
00124     for (i=1;i<j;i++) {
00125       sum=a[i][j];
00126       for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
00127       a[i][j]=sum;
00128     }
00129     big=0.0;
00130     for (i=j;i<=n;i++) {
00131       sum=a[i][j];
00132       for (k=1;k<j;k++)
00133         sum -= a[i][k]*a[k][j];
00134       a[i][j]=sum;
00135       if ( (dum=vv[i]*fabs(sum)) >= big) {
00136         big=dum;
00137         imax=i;
00138       }
00139     }
00140     if (j != imax) {
00141       for (k=1;k<=n;k++) {
00142         dum=a[imax][k];
00143         a[imax][k]=a[j][k];
00144         a[j][k]=dum;
00145       }
00146       *d = -(*d);
00147       vv[imax]=vv[j];
00148     }
00149     indx[j]=imax;
00150     if (a[j][j] == 0.0) a[j][j]=TINY;
00151     if (j != n) {
00152       dum=1.0/(a[j][j]);
00153       for (i=j+1;i<=n;i++) a[i][j] *= dum;
00154     }
00155   }
00156   free_nr_vector(vv,1,n);
00157   return 0;
00158 }
00159 
00160 int svdsolve(eusfloat_t **a, int m, int n, eusfloat_t *b, eusfloat_t *x)
00161 {
00162   int j;
00163   eusfloat_t **v, *w, wmax, wmin;
00164   v = nr_matrix(1,n,1,n);
00165   w = nr_vector(1,n);
00166   if ( svdcmp(a,m,n,w,v) < 0 ) {
00167     free_nr_vector(w,1,n);
00168     free_nr_matrix(v,1,n,1,n);
00169     return -1;
00170   }
00171   wmax = 0.0;
00172   for (j=1; j<=n; j++) if (w[j] > wmax) wmax = w[j];
00173   wmin = wmax*1.0e-6;
00174   for (j=1; j<=n; j++) if (w[j] < wmin) w[j] = 0.0;
00175   svbksb(a,w,v,m,n,b,x);
00176   free_nr_vector(w,1,n);
00177   free_nr_matrix(v,1,n,1,n);
00178   return 1;
00179 }
00180 
00181 void svbksb(eusfloat_t **u, eusfloat_t w[], eusfloat_t **v, int m, int n, eusfloat_t b[], eusfloat_t x[])
00182 {
00183   int jj,j,i;
00184   eusfloat_t s, *tmp;
00185 
00186   tmp = nr_vector(1,n);
00187   for (j=1;j<=n;j++){
00188     s=0.0;
00189     if (w[j]){
00190       for (i=1;i<=m;i++) s += u[i][j]*b[i];
00191       s/= w[j];
00192     }
00193     tmp[j] = s;
00194   }
00195   for (j=1;j<=n;j++){
00196     s=0.0;
00197     for (jj=1;jj<=n;jj++) s+=v[j][jj]*tmp[jj];
00198     x[j]=s;
00199   }
00200   free_nr_vector(tmp,1,n);
00201 }
00202 
00203 int svdcmp(eusfloat_t **a, int m, int n, eusfloat_t w[], eusfloat_t **v)
00204 {
00205   eusfloat_t pythag(eusfloat_t a, eusfloat_t b);
00206   int flag,i,its,j,jj,k,l,nm;
00207   eusfloat_t anorm,c,f,g,h,s,scale,x,y,z,*rv1;
00208   
00209   rv1=nr_vector(1,n);
00210   g=scale=anorm=0.0;
00211   for (i=1;i<=n;i++){
00212     l=i+1;
00213     rv1[i]=scale*g;
00214     g=s=scale=0.0;
00215     if (i<=m){
00216       for (k=i;k<=m;k++) scale += fabs(a[k][i]);
00217       if (scale) {
00218         for (k=i;k<=m;k++){
00219           a[k][i] /= scale;
00220           s+=a[k][i]*a[k][i];
00221         }
00222         f=a[i][i];
00223         g = -SIGN(sqrt(s), f);
00224         h=f*g-s;
00225         a[i][i]=f-g;
00226         for (j=l;j<=n;j++){
00227           for (s=0.0,k=i;k<=m;k++) s+=a[k][i]*a[k][j];
00228           f = s/h;
00229           for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
00230         }
00231         for (k=i;k<=m;k++) a[k][i] *= scale;
00232       }
00233     }
00234     w[i]=scale*g;
00235     g=s=scale=0.0;
00236     if (i<=m && i!=n){
00237       for (k=l;k<=n;k++) scale += fabs(a[i][k]);
00238       if (scale){
00239         for (k=l;k<=n;k++){
00240           a[i][k] /= scale;
00241           s += a[i][k]*a[i][k];
00242         }
00243         f = a[i][l];
00244         g = -SIGN(sqrt(s), f);
00245         h=f*g-s;
00246         a[i][l]=f-g;
00247         for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
00248         for (j=l;j<=m;j++){
00249           for (s=0.0,k=l;k<=n;k++) s+=a[j][k]*a[i][k];
00250           for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
00251         }
00252         for (k=l;k<=n;k++) a[i][k] *= scale;
00253       }
00254     }
00255     anorm=FMAX(anorm, (fabs(w[i])+fabs(rv1[i])));
00256   }
00257   for (i=n;i>=1;i--){
00258     if (i<n){
00259       if (g){
00260         for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
00261         for (j=l;j<=n;j++){
00262           for (s=0.0,k=l;k<=n;k++) s+=a[i][k]*v[k][j];
00263           for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
00264         }
00265       }
00266       for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
00267     }
00268     v[i][i]=1.0;
00269     g=rv1[i];
00270     l=i;
00271   }
00272   for (i=IMIN(m,n);i>=1;i--){
00273     l=i+1;
00274     g=w[i];
00275     for (j=l;j<=n;j++) a[i][j]=0.0;
00276     if (g){
00277       g=1.0/g;
00278       for (j=l;j<=n;j++){
00279         for (s=0.0,k=l;k<=m;k++) s+= a[k][i]*a[k][j];
00280         f = (s/a[i][i])*g;
00281         for (k=i;k<=m;k++) a[k][j]+=f*a[k][i];
00282       }
00283       for (j=i;j<=m;j++) a[j][i] *= g;
00284     }else for (j=i;j<=m;j++) a[j][i]=0.0;
00285     ++a[i][i];
00286   }
00287   for (k=n;k>=1;k--){
00288     for (its=1;its<=30;its++){
00289       flag =1;
00290       for (l=k;l>=1;l--){
00291         nm=l-1;
00292         if ((eusfloat_t)(fabs(rv1[l])+anorm) == anorm){
00293           flag=0;
00294           break;
00295         }
00296         if ((eusfloat_t)(fabs(w[nm])+anorm)==anorm) break;
00297       }
00298       if (flag){
00299         c=0.0;
00300         s=1.0;
00301         for (i=l;i<=k;i++){
00302           f=s*rv1[i];
00303           rv1[i]=c*rv1[i];
00304           if ((eusfloat_t)(fabs(f)+anorm)==anorm) break;
00305           g=w[i];
00306           h=pythag(f,g);
00307           w[i]=h;
00308           h=1.0/h;
00309           c=g*h;
00310           s = -f*h;
00311           for (j=1;j<=m;j++){
00312             y=a[j][nm];
00313             z=a[j][i];
00314             a[j][nm]=y*c+z*s;
00315             a[j][i]=z*c-y*s;
00316           }
00317         }
00318       }
00319       z=w[k];
00320       if (l==k){
00321         if (z<0.0){
00322           w[k] = -z;
00323           for (j=1;j<=n;j++) v[j][k] = -v[j][k];
00324         }
00325         break;
00326       }
00327       if (its == 30) {nrerror("no convergence in 30 svdcmp iterations"); return -1;}
00328       x=w[l];
00329       nm=k-1;
00330       y=w[nm];
00331       g=rv1[nm];
00332       h=rv1[k];
00333       f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
00334       g=pythag(f, 1.0);
00335       f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
00336       c=s=1.0;
00337       for (j=l;j<=nm;j++){
00338         i=j+1;
00339         g=rv1[i];
00340         y=w[i];
00341         h=s*g;
00342         g=c*g;
00343         z=pythag(f,h);
00344         rv1[j]=z;
00345         c=f/z;
00346         s=h/z;
00347         f=x*c+g*s;
00348         g=g*c-x*s;
00349         h=y*s;
00350         y*=c;
00351         for (jj=1;jj<=n;jj++){
00352           x=v[jj][j];
00353           z=v[jj][i];
00354           v[jj][j]=x*c+z*s;
00355           v[jj][i]=z*c-x*s;
00356         }
00357         z=pythag(f,h);
00358         w[j]=z;
00359         if (z) {
00360           z=1.0/z;
00361           c=f*z;
00362           s=h*z;
00363         }
00364         f=c*g+s*y;
00365         x=c*y-s*g;
00366         for (jj=1;jj<=m;jj++){
00367           y=a[jj][j];
00368           z=a[jj][i];
00369           a[jj][j]=y*c+z*s;
00370           a[jj][i]=z*c-y*s;
00371         }
00372       }
00373       rv1[l]=0.0;
00374       rv1[k]=f;
00375       w[k]=x;
00376     }
00377   }
00378   free_nr_vector(rv1,1,n);
00379   return 1;
00380 }
00381 
00382 void tred2(eusfloat_t **a, int n, eusfloat_t d[], eusfloat_t e[])
00383 {
00384   int l,k,j,i;
00385   eusfloat_t scale,hh,h,g,f;
00386   for (i=n;i>=2;i--) {
00387     l=i-1;
00388     h=scale=0.0;
00389     if (l > 1) {
00390       for (k=1;k<=l;k++)
00391         scale += fabs(a[i][k]);
00392       if (scale == 0.0) // Skip transformation.
00393         e[i]=a[i][l];
00394       else {
00395         for (k=1;k<=l;k++) {
00396           a[i][k] /= scale; // Use scaled a's for transformation.
00397           h += a[i][k]*a[i][k]; // Form mu in h.
00398         }
00399         f=a[i][l];
00400         g=(f >= 0.0 ? -sqrt(h) : sqrt(h));
00401         e[i]=scale*g;
00402         h -= f*g; // Now h is equation (11.2.4).
00403         a[i][l]=f-g; // Store u in the ith row of a.
00404         f=0.0;
00405         for (j=1;j<=l;j++) {
00406           /* Next statement can be omitted if eigenvectors not wanted */
00407           a[j][i]=a[i][j]/h; // Store u/H in ith column of a.
00408           g=0.0; // Form an element of A u in g.
00409           for (k=1;k<=j;k++)
00410             g += a[j][k]*a[i][k];
00411           for (k=j+1;k<=l;k++)
00412             g += a[k][j]*a[i][k];
00413           e[j]=g/h; // Form element of p in temporarily unused element of e.
00414           f += e[j]*a[i][j];
00415         }
00416         hh=f/(h+h); // Form K, equation (11.2.11).
00417         for (j=1;j<=l;j++) { // Form q and store in e overwriting p.
00418           f=a[i][j];
00419           e[j]=g=e[j]-hh*f;
00420           for (k=1;k<=j;k++) // Reduce a, equation (11.2.13).
00421             a[j][k] -= (f*e[k]+g*a[i][k]);
00422         }
00423       }
00424     } else
00425       e[i]=a[i][l];
00426     d[i]=h;
00427   }
00428   /* Next statement can be omitted if eigenvectors not wanted */
00429   d[1]=0.0;
00430   e[1]=0.0;
00431   /* Contents of this loop can be omitted if eigenvectors not
00432      wanted except for statement d[i]=a[i][i]; */
00433   for (i=1;i<=n;i++) { // Begin accumulation of transformation matrices.
00434     l=i-1;
00435     if (d[i]) { // This block skipped when i=1.
00436       for (j=1;j<=l;j++) {
00437         g=0.0;
00438         for (k=1;k<=l;k++) // Use u and u/H stored in a to form PQ.
00439           g += a[i][k]*a[k][j];
00440         for (k=1;k<=l;k++)
00441           a[k][j] -= g*a[k][i];
00442       }
00443     }
00444     d[i]=a[i][i]; // This statement remains.
00445     a[i][i]=1.0;  // Reset row and column of a to identity  matrix for next iteration
00446     for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;
00447   }
00448 }
00449 
00450 int tqli(eusfloat_t d[], eusfloat_t e[], int n, eusfloat_t **z)
00451 {
00452   eusfloat_t pythag(eusfloat_t a, eusfloat_t b);
00453   int m,l,iter,i,k;
00454   eusfloat_t s,r,p,g,f,dd,c,b;
00455 
00456   for (i=2;i<=n;i++) e[i-1]=e[i]; // Convenient to renumber the elements of e. 
00457   e[n]=0.0;
00458   for (l=1;l<=n;l++) {
00459     iter=0;
00460     do {
00461       for (m=l;m<=n-1;m++) { // Look for a single small subdiagonal element to split the matrix.
00462         dd=fabs(d[m])+fabs(d[m+1]);
00463         if ((eusfloat_t)(fabs(e[m])+dd) == dd) break;
00464       }
00465       if (m != l) {
00466         if (iter++ == 30) {nrerror("Too many iterations in tqli"); return -1;}
00467         g=(d[l+1]-d[l])/(2.0*e[l]); // Form shift.
00468         r=pythag(g,1.0);
00469         g=d[m]-d[l]+e[l]/(g+SIGN(r,g)); // This is dm . ks.
00470         s=c=1.0;
00471         p=0.0;
00472         for (i=m-1;i>=l;i--) { // A plane rotation as in the original QL, followed by Givens rotations to restore tridiagonal form.
00473           f=s*e[i];
00474           b=c*e[i];
00475           e[i+1]=(r=pythag(f,g));
00476           if (r == 0.0) { // Recover from underflow.
00477             d[i+1] -= p;
00478             e[m]=0.0;
00479             break;
00480           }
00481           s=f/r;
00482           c=g/r;
00483           g=d[i+1]-p;
00484           r=(d[i]-g)*s+2.0*c*b;
00485           d[i+1]=g+(p=s*r);
00486           g=c*r-b;
00487           /* Next loop can be omitted if eigenvectors not wanted*/
00488           for (k=1;k<=n;k++) { // Form eigenvectors.
00489             f=z[k][i+1];
00490             z[k][i+1]=s*z[k][i]+c*f;
00491             z[k][i]=c*z[k][i]-s*f;
00492           }
00493         }
00494         if (r == 0.0 && i >= l) continue;
00495         d[l] -= p;
00496         e[l]=g;
00497         e[m]=0.0;
00498       }
00499     } while (m != l);
00500   }
00501   return 1;
00502 }
00503 
00504 eusfloat_t pythag(eusfloat_t a, eusfloat_t b)
00505 {
00506   eusfloat_t absa, absb;
00507   absa=fabs(a);
00508   absb=fabs(b);
00509   if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));
00510   else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
00511 }
00512 


jskeus
Author(s): JSK Alumnis
autogenerated on Fri Aug 28 2015 11:15:08