home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
vis-ftp.cs.umass.edu
/
vis-ftp.cs.umass.edu.tar
/
vis-ftp.cs.umass.edu
/
pub
/
Software
/
ASCENDER
/
NRroutines.c
< prev
next >
Wrap
C/C++ Source or Header
|
1996-04-25
|
11KB
|
583 lines
#include "../polygons.h"
#define TINY 1.0e-20;
void nrerror(error_text)
char error_text[];
{
void exit();
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(int nl,int nh)
{
float *v;
v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
if (!v) nrerror("allocation failure in vector()");
return v-nl;
}
int *ivector(int nl,int nh)
{
int *v;
v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
if (!v) nrerror("allocation failure in ivector()");
return v-nl;
}
double *dvector(int nl,int nh)
{
double *v;
v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
if (!v) nrerror("allocation failure in dvector()");
return v-nl;
}
float **matrix(int nrl,int nrh,int ncl,int nch)
{
int i;
float **m;
m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*));
if (!m) nrerror("allocation failure 1 in matrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
if (!m[i]) nrerror("allocation failure 2 in matrix()");
m[i] -= ncl;
}
return m;
}
double **dmatrix(int nrl,int nrh,int ncl,int nch)
{
int i;
double **m;
m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
m[i] -= ncl;
}
return m;
}
int **imatrix(int nrl,int nrh,int ncl,int nch)
{
int i,**m;
m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
if (!m) nrerror("allocation failure 1 in imatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
if (!m[i]) nrerror("allocation failure 2 in imatrix()");
m[i] -= ncl;
}
return m;
}
float **submatrix(float **a,int oldrl,int oldrh,int oldcl,
int oldch,int newrl,int newcl)
{
int i,j;
float **m;
m=(float **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(float*));
if (!m) nrerror("allocation failure in submatrix()");
m -= newrl;
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;
return m;
}
void free_vector(float *v,int nl,int nh)
{
free((char*) (v+nl));
}
void free_ivector(int *v,int nl,int nh)
{
free((char*) (v+nl));
}
void free_dvector(double *v,int nl,int nh)
{
free((char*) (v+nl));
}
void free_matrix(float **m,int nrl,int nrh,int ncl,int nch)
{
int i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}
void free_dmatrix(double **m,int nrl,int nrh,int ncl,int nch)
{
int i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}
void free_imatrix(int **m,int nrl,int nrh,int ncl,int nch)
{
int i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}
void free_submatrix(float **b,int nrl,int nrh,int ncl,int nch)
{
free((char*) (b+nrl));
}
float **convert_matrix(float *a,int nrl,int nrh,int ncl,int nch)
{
int i,j,nrow,ncol;
float **m;
nrow=nrh-nrl+1;
ncol=nch-ncl+1;
m = (float **) malloc((unsigned) (nrow)*sizeof(float*));
if (!m) nrerror("allocation failure in convert_matrix()");
m -= nrl;
for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
return m;
}
void free_convert_matrix(float **b,int nrl,int nrh,int ncl,int nch)
{
free((char*) (b+nrl));
}
void ludcmp(float **a,int n,int *indx,float *d)
{
int i,imax,j,k;
float big,dum,sum,temp;
float *vv,*vector();
void nrerror(),free_vector();
vv=vector(1,n);
*d=1.0;
for (i=1;i<=n;i++) {
big=0.0;
for (j=1;j<=n;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
vv[i]=1.0/big;
}
for (j=1;j<=n;j++) {
for (i=1;i<j;i++) {
sum=a[i][j];
for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<=n;i++) {
sum=a[i][j];
for (k=1;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big) {
big=dum;
imax=i;
}
}
if (j != imax) {
for (k=1;k<=n;k++) {
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d = -(*d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != n) {
dum=1.0/(a[j][j]);
for (i=j+1;i<=n;i++) a[i][j] *= dum;
}
}
free_vector(vv,1,n);
}
void bksub(int ne,int nb,int jf,int k1,int k2,float ***c)
{
int nbf,im,kp,k,j,i;
float xx;
nbf=ne-nb;
im=1;
for (k=k2;k>=k1;k--) {
if (k == k1) im=nbf+1;
kp=k+1;
for (j=1;j<=nbf;j++) {
xx=c[j][jf][kp];
for (i=im;i<=ne;i++)
c[i][jf][k] -= c[i][j][k]*xx;
}
}
for (k=k1;k<=k2;k++) {
kp=k+1;
for (i=1;i<=nb;i++) c[i][1][k]=c[i+nbf][jf][k];
for (i=1;i<=nbf;i++) c[i+nb][1][k]=c[i][jf][kp];
}
}
void lubksb(float **a,int n,int *indx,float b[])
{
int i,ii=0,ip,j;
float sum;
for (i=1;i<=n;i++) {
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
if (ii)
for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
else if (sum) ii=i;
b[i]=sum;
}
for (i=n;i>=1;i--) {
sum=b[i];
for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];
}
}
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
a[k][l]=h+s*(g-h*tau);
void jacobi(a,n,d,v,nrot)
float **a,d[],**v;
int n,*nrot;
{
int j,iq,ip,i;
float tresh,theta,tau,t,sm,s,h,g,c,*b,*z,*vector();
void nrerror(),free_vector();
b=vector(1,n);
z=vector(1,n);
for (ip=1;ip<=n;ip++) {
for (iq=1;iq<=n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}
for (ip=1;ip<=n;ip++) {
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;
for (i=1;i<=50;i++) {
sm=0.0;
for (ip=1;ip<=n-1;ip++) {
for (iq=ip+1;iq<=n;iq++)
sm += fabs(a[ip][iq]);
}
if (sm == 0.0) {
free_vector(z,1,n);
free_vector(b,1,n);
return;
}
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=1;ip<=n-1;ip++) {
for (iq=ip+1;iq<=n;iq++) {
g=100.0*fabs(a[ip][iq]);
if (i > 4 && fabs(d[ip])+g == fabs(d[ip])
&& fabs(d[iq])+g == fabs(d[iq]))
a[ip][iq]=0.0;
else if (fabs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if (fabs(h)+g == fabs(h))
t=(a[ip][iq])/h;
else {
theta=0.5*h/(a[ip][iq]);
t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=1;j<=ip-1;j++) {
ROTATE(a,j,ip,j,iq)
}
for (j=ip+1;j<=iq-1;j++) {
ROTATE(a,ip,j,j,iq)
}
for (j=iq+1;j<=n;j++) {
ROTATE(a,ip,j,iq,j)
}
for (j=1;j<=n;j++) {
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=1;ip<=n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
nrerror("Too many iterations in routine JACOBI");
}
#undef ROTATE
static float at,bt,ct;
#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
(ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
static float maxarg1,maxarg2;
#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
void svdcmp(a,m,n,w,v)
float **a,*w,**v;
int m,n;
{
int flag,i,its,j,jj,k,l,nm;
float c,f,h,s,x,y,z;
float anorm=0.0,g=0.0,scale=0.0;
float *rv1,*vector();
void nrerror(),free_vector();
if (m < n) nrerror("SVDCMP: You must augment A with extra zero rows");
rv1=vector(1,n);
for (i=1;i<=n;i++) {
l=i+1;
rv1[i]=scale*g;
g=s=scale=0.0;
if (i <= m) {
for (k=i;k<=m;k++) scale += fabs(a[k][i]);
if (scale) {
for (k=i;k<=m;k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
f=a[i][i];
g = -SIGN(sqrt(s),f);
h=f*g-s;
a[i][i]=f-g;
if (i != n) {
for (j=l;j<=n;j++) {
for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
f=s/h;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
}
for (k=i;k<=m;k++) a[k][i] *= scale;
}
}
w[i]=scale*g;
g=s=scale=0.0;
if (i <= m && i != n) {
for (k=l;k<=n;k++) scale += fabs(a[i][k]);
if (scale) {
for (k=l;k<=n;k++) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
f=a[i][l];
g = -SIGN(sqrt(s),f);
h=f*g-s;
a[i][l]=f-g;
for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
if (i != m) {
for (j=l;j<=m;j++) {
for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
}
}
for (k=l;k<=n;k++) a[i][k] *= scale;
}
}
anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
}
for (i=n;i>=1;i--) {
if (i < n) {
if (g) {
for (j=l;j<=n;j++)
v[j][i]=(a[i][j]/a[i][l])/g;
for (j=l;j<=n;j++) {
for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
}
}
for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0;
g=rv1[i];
l=i;
}
for (i=n;i>=1;i--) {
l=i+1;
g=w[i];
if (i < n)
for (j=l;j<=n;j++) a[i][j]=0.0;
if (g) {
g=1.0/g;
if (i != n) {
for (j=l;j<=n;j++) {
for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
}
for (j=i;j<=m;j++) a[j][i] *= g;
} else {
for (j=i;j<=m;j++) a[j][i]=0.0;
}
++a[i][i];
}
for (k=n;k>=1;k--) {
for (its=1;its<=30;its++) {
flag=1;
for (l=k;l>=1;l--) {
nm=l-1;
if (fabs(rv1[l])+anorm == anorm) {
flag=0;
break;
}
if (fabs(w[nm])+anorm == anorm) break;
}
if (flag) {
c=0.0;
s=1.0;
for (i=l;i<=k;i++) {
f=s*rv1[i];
if (fabs(f)+anorm != anorm) {
g=w[i];
h=PYTHAG(f,g);
w[i]=h;
h=1.0/h;
c=g*h;
s=(-f*h);
for (j=1;j<=m;j++) {
y=a[j][nm];
z=a[j][i];
a[j][nm]=y*c+z*s;
a[j][i]=z*c-y*s;
}
}
}
}
z=w[k];
if (l == k) {
if (z < 0.0) {
w[k] = -z;
for (j=1;j<=n;j++) v[j][k]=(-v[j][k]);
}
break;
}
if (its == 30) {
w[1] = 0.0;
return;
}
x=w[l];
nm=k-1;
y=w[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=PYTHAG(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
c=s=1.0;
for (j=l;j<=nm;j++) {
i=j+1;
g=rv1[i];
y=w[i];
h=s*g;
g=c*g;
z=PYTHAG(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g=g*c-x*s;
h=y*s;
y=y*c;
for (jj=1;jj<=n;jj++) {
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*c+z*s;
v[jj][i]=z*c-x*s;
}
z=PYTHAG(f,h);
w[j]=z;
if (z) {
z=1.0/z;
c=f*z;
s=h*z;
}
f=(c*g)+(s*y);
x=(c*y)-(s*g);
for (jj=1;jj<=m;jj++) {
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*c+z*s;
a[jj][i]=z*c-y*s;
}
}
rv1[l]=0.0;
rv1[k]=f;
w[k]=x;
}
}
free_vector(rv1,1,n);
}
#undef SIGN
#undef MAX
#undef PYTHAG