added codes
This commit is contained in:
mhjensen
Executable
+183
@@ -0,0 +1,183 @@
|
||||
#include <math.h>
|
||||
#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
|
||||
Executable
+41
@@ -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
|
||||
Executable
+22
@@ -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
|
||||
@@ -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.
|
||||
Executable
+6
@@ -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
|
||||
Executable
+9
@@ -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
|
||||
Executable
+444
@@ -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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user