diff --git a/doc/Programs/SVD/Cpp/svdcmp.c b/doc/Programs/SVD/Cpp/svdcmp.c new file mode 100755 index 000000000..d3611254c --- /dev/null +++ b/doc/Programs/SVD/Cpp/svdcmp.c @@ -0,0 +1,183 @@ +#include +#define NRANSI +#include "nrutil.h" + +void svdcmp(float **a, int m, int n, float w[], float **v) +{ + float pythag(float a, float b); + int flag,i,its,j,jj,k,l,nm; + float anorm,c,f,g,h,s,scale,x,y,z,*rv1; + + rv1=vector(1,n); + g=scale=anorm=0.0; + 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; + 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; + 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=FMAX(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=IMIN(m,n);i>=1;i--) { + l=i+1; + g=w[i]; + for (j=l;j<=n;j++) a[i][j]=0.0; + if (g) { + g=1.0/g; + 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 ((float)(fabs(rv1[l])+anorm) == anorm) { + flag=0; + break; + } + if ((float)(fabs(w[nm])+anorm) == anorm) break; + } + if (flag) { + c=0.0; + s=1.0; + for (i=l;i<=k;i++) { + f=s*rv1[i]; + rv1[i]=c*rv1[i]; + if ((float)(fabs(f)+anorm) == anorm) break; + 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) nrerror("no convergence in 30 svdcmp iterations"); + 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 *= 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 NRANSI diff --git a/doc/Programs/SVD/Cpp/svdfit.c b/doc/Programs/SVD/Cpp/svdfit.c new file mode 100755 index 000000000..47f3b74d8 --- /dev/null +++ b/doc/Programs/SVD/Cpp/svdfit.c @@ -0,0 +1,41 @@ +#define NRANSI +#include "nrutil.h" +#define TOL 1.0e-5 + +void svdfit(float x[], float y[], float sig[], int ndata, float a[], int ma, + float **u, float **v, float w[], float *chisq, + void (*funcs)(float, float [], int)) +{ + void svbksb(float **u, float w[], float **v, int m, int n, float b[], + float x[]); + void svdcmp(float **a, int m, int n, float w[], float **v); + int j,i; + float wmax,tmp,thresh,sum,*b,*afunc; + + b=vector(1,ndata); + afunc=vector(1,ma); + for (i=1;i<=ndata;i++) { + (*funcs)(x[i],afunc,ma); + tmp=1.0/sig[i]; + for (j=1;j<=ma;j++) u[i][j]=afunc[j]*tmp; + b[i]=y[i]*tmp; + } + svdcmp(u,ndata,ma,w,v); + wmax=0.0; + for (j=1;j<=ma;j++) + if (w[j] > wmax) wmax=w[j]; + thresh=TOL*wmax; + for (j=1;j<=ma;j++) + if (w[j] < thresh) w[j]=0.0; + svbksb(u,w,v,ndata,ma,b,a); + *chisq=0.0; + for (i=1;i<=ndata;i++) { + (*funcs)(x[i],afunc,ma); + for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j]; + *chisq += (tmp=(y[i]-sum)/sig[i],tmp*tmp); + } + free_vector(afunc,1,ma); + free_vector(b,1,ndata); +} +#undef TOL +#undef NRANSI diff --git a/doc/Programs/SVD/Cpp/svdvar.c b/doc/Programs/SVD/Cpp/svdvar.c new file mode 100755 index 000000000..23c711fa3 --- /dev/null +++ b/doc/Programs/SVD/Cpp/svdvar.c @@ -0,0 +1,22 @@ +#define NRANSI +#include "nrutil.h" + +void svdvar(float **v, int ma, float w[], float **cvm) +{ + int k,j,i; + float sum,*wti; + + wti=vector(1,ma); + for (i=1;i<=ma;i++) { + wti[i]=0.0; + if (w[i]) wti[i]=1.0/(w[i]*w[i]); + } + for (i=1;i<=ma;i++) { + for (j=1;j<=i;j++) { + for (sum=0.0,k=1;k<=ma;k++) sum += v[i][k]*v[j][k]*wti[k]; + cvm[j][i]=cvm[i][j]=sum; + } + } + free_vector(wti,1,ma); +} +#undef NRANSI diff --git a/doc/Programs/SVD/Fortran/README.txt b/doc/Programs/SVD/Fortran/README.txt new file mode 100644 index 000000000..7e58f9d49 --- /dev/null +++ b/doc/Programs/SVD/Fortran/README.txt @@ -0,0 +1,23 @@ +This folder contains a simple benchmark case on how to run the program and how the +input file should look like. + +The best test case is to run the simplefit.dat file using the simplefit.f90 code. +This input file is set up as follows: + + 8 2 # number of data entries and order of polynomial + 0.001 -2.89017 0.00073621 + 0.002 -2.88946 0.00052732 + 0.005 -2.89067 0.00055038 + 0.010 -2.89091 0.00040973 + 0.015 -2.89084 0.00034278 + 0.02 -2.89086 0.00029315 + 0.025 -2.89059 0.00034278 + 0.03 -2.89077 0.00025017 + The present input file stems from a variational Monte Carlo calculation of the ground + state energy (second column). It contains also the standard deviation (3rd column) + The first column is the time step used in importance sampling. Extrapolating to zero + should give the best estimate for the VMC energy. + + +Run the code as executable < inputfile > outputfile after having compiled and linked the source +files. diff --git a/doc/Programs/SVD/Fortran/poly.dat b/doc/Programs/SVD/Fortran/poly.dat new file mode 100755 index 000000000..038960d48 --- /dev/null +++ b/doc/Programs/SVD/Fortran/poly.dat @@ -0,0 +1,6 @@ +5 3 +1.0 3.0 0.0001 +2.0 17.0 0.0001 +3.0 49.0 0.0001 +4.0 105.0 0.0001 +5.0 191.0 0.0001 \ No newline at end of file diff --git a/doc/Programs/SVD/Fortran/simplefit.dat b/doc/Programs/SVD/Fortran/simplefit.dat new file mode 100755 index 000000000..2ed9a6136 --- /dev/null +++ b/doc/Programs/SVD/Fortran/simplefit.dat @@ -0,0 +1,9 @@ +8 3 +0.001 -2.89017 0.00073621 +0.002 -2.88946 0.00052732 +0.005 -2.89067 0.00055038 +0.010 -2.89091 0.00040973 +0.015 -2.89084 0.00034278 +0.02 -2.89086 0.00029315 +0.025 -2.89059 0.00034278 +0.03 -2.89077 0.00025017 \ No newline at end of file diff --git a/doc/Programs/SVD/Fortran/simplefit.f90 b/doc/Programs/SVD/Fortran/simplefit.f90 new file mode 100755 index 000000000..81b6d1684 --- /dev/null +++ b/doc/Programs/SVD/Fortran/simplefit.f90 @@ -0,0 +1,444 @@ +! Program to fit a given set of data in terms of a given polynomial using +! Singular Value Decomposition, see Golub and Van Loan, Matrix Computations +! Numerical Recipes, chapter 15. +! This explict example fits a polynomial to a given degree based. The example +! input should look like +! 8 2 # number of data entries and order of polynomial +! 0.001 -2.89017 0.00073621 +! 0.002 -2.88946 0.00052732 +! 0.005 -2.89067 0.00055038 +! 0.010 -2.89091 0.00040973 +! 0.015 -2.89084 0.00034278 +! 0.02 -2.89086 0.00029315 +! 0.025 -2.89059 0.00034278 +! 0.03 -2.89077 0.00025017 +! The present input file stems from a variational Monte Carlo calculation of the ground +! state energy (second column). It contains also the standard deviation (3rd column) +! The first column is the time step used in importance sampling. Extrapolating to zero +! should give the best estimate for the VMC energy. +! +! Author: Morten Hjorth-Jensen, Department of Physics, University of Oslo, +! POB 1048 Blindern, N-0316 Oslo, Norway. email: morten.hjorth-jensen@fys.uio.no +! Latest update: October 1999. +! Run in old fashioned way, executable < inputfile > outputfile + +! +! Main program starts here +! + PROGRAM fitting + IMPLICIT NONE + CALL vmcfit + + END PROGRAM fitting +! +! This subroutine fits a data set with in terms +! of a polynomial expansion. The number of terms in the polynomial +! expansion is given by the variable number_terms +! The number of data ( delta and corresponding energy +! and the standard error) are defined by the variable number_of_data +! + SUBROUTINE vmcfit + IMPLICIT NONE + INTEGER :: i, number_of_data, number_terms + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: w, n, e,sig,polynom_terms, variance + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: cvm, v, u + DOUBLE PRECISION :: chisq, kappa, v_over_c, adiabatic + EXTERNAL eos_param + + READ (5,*) number_of_data, number_terms + ALLOCATE (cvm(number_terms,number_terms), & + v(number_terms,number_terms) ) + ALLOCATE ( u(number_of_data,number_terms) ) + ALLOCATE ( w(number_terms) ) + ALLOCATE ( variance(number_terms) ) + ALLOCATE ( polynom_terms(number_terms) ) + ALLOCATE (n(number_of_data), e(number_of_data), sig(number_of_data) ) + DO i=1,number_of_data + READ (5,*)n(i) , e(i), sig(i) + ENDDO + CALL svdfit(n,e,sig,number_of_data,polynom_terms, & + number_terms, & + u,v,w,number_of_data,number_terms,chisq) + CALL svdvar(v,number_terms,number_terms,w,cvm,number_terms) +! Compute now the variance for each coefficient using Eq. 15.4.19 +! of numerical recipes + variance = 0.0; w = 1.0/(w*w) + DO i = 1, number_terms + variance(i) = SUM(v(i,:)*v(i,:)*w(:)) + ENDDO + WRITE(*,*) ' Quality of fit:' + WRITE(*,*) ' Number of terms in polynomial expansion:', number_terms + WRITE (*,'(''CHISQ = '',F12.5)') chisq + WRITE(*,*) 'Extrapolated energy with error' + WRITE(*,*) polynom_terms(1), SQRT(variance(1)) + DEALLOCATE (cvm, v) + DEALLOCATE ( u ) + DEALLOCATE ( w,polynom_terms, variance ) + DEALLOCATE (n, e, sig) + + END SUBROUTINE vmcfit + + +! This function encodes the actual functional form of the polynomial. +! If you need to change the functional form, this is the function to adjust. +! This specific version is a polynomial in x^i, with i >=0 and being an integer + + SUBROUTINE eos_param(x,afunc,ma) + IMPLICIT NONE + INTEGER :: ma, i + DOUBLE PRECISION :: afunc, x + DIMENSION afunc(ma) + afunc(1)=1.D0 + DO i=2,ma + afunc(i)=afunc(i-1)*x + ENDDO + + END SUBROUTINE eos_param + + + + SUBROUTINE svdvar(v,ma,np,w,cvm,ncvm) + IMPLICIT NONE + INTEGER :: ma,ncvm,np,MMAX + DOUBLE PRECISION :: cvm(ncvm,ncvm),v(np,np),w(np) + PARAMETER (MMAX=100) + INTEGER :: i,j,k + DOUBLE PRECISION :: sum,wti(MMAX) + DO i=1,ma + wti(i)=0. + IF (w(i) /= 0.) wti(i)=1./(w(i)*w(i)) + ENDDO + DO i=1,ma + DO j=1,i + sum=0. + DO k=1,ma + sum=sum+v(i,k)*v(j,k)*wti(k) + ENDDO + cvm(i,j)=sum + cvm(j,i)=sum + ENDDO + ENDDO + + END SUBROUTINE svdvar + + + SUBROUTINE svbksb(u,w,v,m,n,mp,np,b,x) + IMPLICIT NONE + INTEGER :: m,mp,n,np,NMAX + DOUBLE PRECISION :: b(mp),u(mp,np),v(np,np),w(np),x(np) + PARAMETER (NMAX=5000) + INTEGER :: i,j,jj + DOUBLE PRECISION :: s,tmp(NMAX) + DO j=1,n + s=0. + IF (w(j) /= 0.) THEN + DO i=1,m + s=s+u(i,j)*b(i) + ENDDO + s=s/w(j) + ENDIF + tmp(j)=s + ENDDO + DO j=1,n + s=0. + DO jj=1,n + s=s+v(j,jj)*tmp(jj) + ENDDO + x(j)=s + ENDDO + + END SUBROUTINE svbksb + + + SUBROUTINE svdfit(x,y,sig,ndata,a,ma,u,v,w,mp,np,chisq) + IMPLICIT NONE + INTEGER :: ma,mp,ndata,np,NMAX,MMAX + DOUBLE PRECISION :: chisq,a(ma),sig(ndata),u(mp,np),v(np,np),w(np), & + x(ndata), y(ndata),TOL + PARAMETER (NMAX=5000,MMAX=100,TOL=1.e-5) + INTEGER :: i,j + DOUBLE PRECISION :: sum,thresh,tmp,wmax,afunc(MMAX), & + b(NMAX) + + DO i=1,ndata + CALL eos_param(x(i),afunc,ma) + tmp=1./sig(i) + DO j=1,ma + u(i,j)=afunc(j)*tmp + ENDDO + b(i)=y(i)*tmp + ENDDO + CALL svdcmp(u,ndata,ma,mp,np,w,v) + wmax=0. + DO j=1,ma + IF (w(j) > wmax)wmax=w(j) + ENDDO + thresh=TOL*wmax + DO j=1,ma + IF (w(j) < thresh) w(j)=0. + ENDDO + CALL svbksb(u,w,v,ndata,ma,mp,np,b,a) + chisq=0. + DO i=1,ndata + CALL eos_param(x(i),afunc,ma) + sum=0. + DO j=1,ma + sum=sum+a(j)*afunc(j) + ENDDO + chisq=chisq+((y(i)-sum)/sig(i))**2 + ENDDO + + END SUBROUTINE svdfit + + + FUNCTION pythag(a,b) + IMPLICIT NONE + DOUBLE PRECISION :: a, b, pythag + DOUBLE PRECISION :: absa, absb + absa=ABS(a) + absb=ABS(b) + IF (absa > absb) THEN + pythag=absa*SQRT(1.+(absb/absa)**2) + ELSE + IF(absb == 0.) THEN + pythag=0. + ELSE + pythag=absb*SQRT(1.+(absa/absb)**2) + ENDIF + ENDIF + + END FUNCTION pythag + + SUBROUTINE svdcmp(a,m,n,mp,np,w,v) + IMPLICIT NONE + INTEGER :: m,mp,n,np,NMAX + DOUBLE PRECISION :: a(mp,np),v(np,np),w(np) + PARAMETER (NMAX=5000) + INTEGER :: i,its,j,jj,k,l,nm + DOUBLE PRECISION :: anorm,c,f,g,h,s,scale,x,y,z, & + rv1(NMAX),pythag + g=0.0 + scale=0.0 + anorm=0.0 + DO i=1,n + l=i+1 + rv1(i)=scale*g + g=0.0 + s=0.0 + scale=0.0 + IF(i <= m)THEN + DO k=i,m + scale=scale+ABS(a(k,i)) + ENDDO + IF(scale /= 0.0)THEN + DO k=i,m + a(k,i)=a(k,i)/scale + s=s+a(k,i)*a(k,i) + ENDDO + f=a(i,i) + g=-sign(SQRT(s),f) + h=f*g-s + a(i,i)=f-g + DO j=l,n + s=0.0 + DO k=i,m + s=s+a(k,i)*a(k,j) + ENDDO + f=s/h + DO k=i,m + a(k,j)=a(k,j)+f*a(k,i) + ENDDO + ENDDO + DO k=i,m + a(k,i)=scale*a(k,i) + ENDDO + ENDIF + ENDIF + w(i)=scale *g + g=0.0 + s=0.0 + scale=0.0 + IF((i <= m).and.(i /= n))THEN + DO k=l,n + scale=scale+ABS(a(i,k)) + ENDDO + IF(scale /= 0.0)THEN + DO k=l,n + a(i,k)=a(i,k)/scale + s=s+a(i,k)*a(i,k) + ENDDO + f=a(i,l) + g=-sign(SQRT(s),f) + h=f*g-s + a(i,l)=f-g + DO k=l,n + rv1(k)=a(i,k)/h + ENDDO + DO j=l,m + s=0.0 + DO k=l,n + s=s+a(j,k)*a(i,k) + ENDDO + DO k=l,n + a(j,k)=a(j,k)+s*rv1(k) + ENDDO + ENDDO + DO k=l,n + a(i,k)=scale*a(i,k) + ENDDO + ENDIF + ENDIF + anorm=max(anorm,(ABS(w(i))+ABS(rv1(i)))) + ENDDO + DO i=n,1,-1 + IF(i < n)THEN + IF(g /= 0.0)THEN + DO j=l,n + v(j,i)=(a(i,j)/a(i,l))/g + ENDDO + DO j=l,n + s=0.0 + DO k=l,n + s=s+a(i,k)*v(k,j) + ENDDO + DO k=l,n + v(k,j)=v(k,j)+s*v(k,i) + ENDDO + ENDDO + ENDIF + DO j=l,n + v(i,j)=0.0 + v(j,i)=0.0 + ENDDO + ENDIF + v(i,i)=1.0 + g=rv1(i) + l=i + ENDDO + DO i=min(m,n),1,-1 + l=i+1 + g=w(i) + DO j=l,n + a(i,j)=0.0 + ENDDO + IF(g /= 0.0)THEN + g=1.0/g + DO j=l,n + s=0.0 + DO k=l,m + s=s+a(k,i)*a(k,j) + ENDDO + f=(s/a(i,i))*g + DO k=i,m + a(k,j)=a(k,j)+f*a(k,i) + ENDDO + ENDDO + DO j=i,m + a(j,i)=a(j,i)*g + ENDDO + ELSE + DO j= i,m + a(j,i)=0.0 + ENDDO + ENDIF + a(i,i)=a(i,i)+1.0 + ENDDO + DO k=n,1,-1 + DO its=1,30 + DO l=k,1,-1 + nm=l-1 + IF((ABS(rv1(l))+anorm) == anorm) goto 2 + IF((ABS(w(nm))+anorm) == anorm) goto 1 + ENDDO +1 c=0.0 + s=1.0 + DO i=l,k + f=s*rv1(i) + rv1(i)=c*rv1(i) + IF((ABS(f)+anorm) == anorm) goto 2 + g=w(i) + h=pythag(f,g) + w(i)=h + h=1.0/h + c= (g*h) + s=-(f*h) + DO j=1,m + y=a(j,nm) + z=a(j,i) + a(j,nm)=(y*c)+(z*s) + a(j,i)=-(y*s)+(z*c) + ENDDO + ENDDO +2 z=w(k) + IF(l == k)THEN + IF(z < 0.0)THEN + w(k)=-z + DO j=1,n + v(j,k)=-v(j,k) + ENDDO + ENDIF + goto 3 + ENDIF + IF(its == 30) THEN + WRITE(*,*) 'no convergence in svdcmp'; STOP + ENDIF + 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.D0) + f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x + c=1.0 + s=1.0 + DO j=l,nm + 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=-(x*s)+(g*c) + h=y*s + y=y*c + DO jj=1,n + x=v(jj,j) + z=v(jj,i) + v(jj,j)= (x*c)+(z*s) + v(jj,i)=-(x*s)+(z*c) + ENDDO + z=pythag(f,h) + w(j)=z + IF(z /= 0.0)THEN + z=1.0/z + c=f*z + s=h*z + ENDIF + f= (c*g)+(s*y) + x=-(s*g)+(c*y) + DO jj=1,m + y=a(jj,j) + z=a(jj,i) + a(jj,j)= (y*c)+(z*s) + a(jj,i)=-(y*s)+(z*c) + ENDDO + ENDDO + rv1(l)=0.0 + rv1(k)=f + w(k)=x + ENDDO +3 CONTINUE + ENDDO + + END SUBROUTINE svdcmp + + + + + +