/* ------------------------------------------------------------- -------------------- S V D -------------------- This is a tandem translation of a cdc 6600 fortran program to ibm 360 fortran iv and then to 'C'. The original routine used short precision arithmetic; thus, I have used single precision arithmetic. original cdc 6600 programmer= R. C. Singleton 360 version by= J. G. Lewis 'C' version by k.t. kilty latest revision October 31, 1988 ------------------------------------------------------------------- additional routine needed: rotate ------------------------------------------------------------------- This subroutine computes the singular value decomposition of a real m*n matrix a, i.e. it computes matrices u,s,and v such that a = u * s * vt where u is an m*m matrix and ut*u = i, (ut=transpose of u), v is an n*n matrix and vt*v = i, (vt=transpose of v), and s is an n*n diagonal matrix. description of parameters a= real array. a contains the matrix to be decomposed. the original data are lost. If withv=TRUE, then the matrix u is computed and stored in the array a. a is organized row major, and should work with 'C' two-D vectors. m,n = integer variables. The number of rows and columns in the matrix stored in a. If m0, then p columns n,..., n+p-1 of a are assumed to contain the columns of an m*p matrix b. This matrix is multiplied by ut, and upon exit, a contains in these same columns the n*p matrix ut*b. (p>=0) withu, withv = integers. If withu=TRUE, then the matrix u is computed and stored in the array a. If withv=TRUE, then the matrix v is computed and stored in the array v. s = real array. s[0], ..., s[n-1] contain the diagonal elements of the matrix s ordered so that s[i]>=s[i+1], i=0, ..., n-1. The s' are the singular values of the matrix and satisfy the shifted eigenvalue equations a.v=su a.u=sv or the s-squared are eigen values of the matrix (at).a. v = real array. v contains the matrix v. If withu and withv are not both = TRUE, then the actual parameter corresponding to a and v may be the same. This original routine is a version of a fortran subroutine by Businger and Golub, algorithm 358= singular value decomposition of a complex matrix, comm. acm, v. 12, no. 10, pp. 564-565 (Oct. 1969). with revisions by RC Singleton, May 1972. ------------------------------------------------------------------- */ #define BEGIN { #define END } #define TRUE 1 #define FALSE 0 #include #include #include void rotate(float *x,int offset1,int offset2,float cosine,float sine,int count,int cols) BEGIN int j,bias1,bias2; float xx; for(j=0; j=tolerance) BEGIN g=(float)sqrt((double)z); f=a[k*np+k]; if(f>=0.0) g=-g; s[k]=g; h=g*(f-g); a[k*np+k]=f-g; if (k<(np-1)) BEGIN for(j=l; j=tolerance */ /* elimination of a[k,j], j=k+2, ..., n */ temp=(float)(fabs((double)s[k])+fabs((double)t[k])); if (epsilon=tolerance) BEGIN g=(float)sqrt((double)z); f=a[k*np+l]; if(f>=0.0) g=-g; h=g*(f-g); a[k*np+l]=f-g; for(j=l; j=tolerance */ END /* if k-1) BEGIN if(t[l]!=0.0) BEGIN h=a[k*np+l]*t[l]; for(j=l; j-1) */ END /* if withv==TRUE */ if(withu==TRUE) BEGIN k=n-1; g=s[k]; if(g!=0.0) g=1.0/g; for(j=k; j-1) BEGIN for(j=l; j-1) */ END /* if withu==TRUE */ /* qr diagonalization */ k=n-1; /* test for split */ do BEGIN /* end while */ l=k; while((float)fabs((double)t[l])>epsilon) BEGIN l--; if ((float)fabs((double)s[l])<=epsilon) break; END /* cancellation */ if((float)fabs((double)s[l])<=epsilon) BEGIN sine=1.0; cosine=0.0; l1=l++; for(i=l; i<=k; i++) BEGIN f=sine*t[i]; t[i]*=cosine; if((float)fabs((double)f)>epsilon) BEGIN h=s[i]; w=(float)sqrt((double)(f*f+h*h)); s[i]=w; cosine=h/w; sine=-f/w; if(withu==TRUE) rotate(a,l1,i,cosine,sine,m,np); if(np>n) BEGIN for(j=n; jepsilon */ else break; END /* loop over i */ END /* if abs(s[l])-1); /* sort singular values */ for(k=0; k=g) BEGIN g=s[i]; j=i; END if(j!=k) BEGIN g=s[j]; s[j]=s[k]; s[k]=g; if(withv==TRUE) BEGIN for(i=0; in) BEGIN for(i=n; in */ END /* if j!=k */ END /* loop over k */ free(t); return(0); END /* svd */ main() /* To test svd routine */ BEGIN FILE *fp; float static x[10]=BEGIN 0,.5,1,1.5,2,2.5,3,3.5,4,4.5 END; float static t[10]=BEGIN 95,64,50,36,30,0,0,0,0,0 END; float static a[10*4]; float sum; float s[4],v[26]; int i,j,k,tindex,index,ROWS=5,COLS=4; fp = fopen("svdtest3.prn","w"); for(j=0; j