关于最小二乘法的那个程序

zxjcool

最小二乘法的subroutine.



SUBROUTINE mrqmin(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,&
           COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA)
INTEGER ma,nca,ndata,lista(ma),MMAX
REAL alamda,chisq,funcs,a(ma),alpha(nca,nca),&
     covar(nca,nca),sig(ndata),x(ndata),y(ndata)
PARAMETER (MMAX=20)
!USES covsrt,gaussj,mrqcof
INTEGER j,k,l,mfit
REAL ochisq,atry(MMAX),beta(MMAX),da(MMAX)
if(alamda<0.) then
  kk=mfit+1
  do j=1,ma
    ihit=0
    do k=1,mfit
      if(lista(k)==j) ihit=ihit+1
    end do
    if(ihit==0) then
      lista(kk)=j
      kk=kk+1
    else if(ihit>1) then
      pause 'improper permutation in lista'
    endif
  end do
  if(kk/=(ma+1)) pause 'improper permutation in lista'
  alamda=0.001
  call mrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,&
                     beta,nca,chisq,funcs)
  ochisq=chisq
  do j=1,ma
    atry(j)=a(j)
  end do
endif
do j=1,mfit
  do k=1,mfit
    covar(j,k)=alpha(j,k)
  end do
  covar(j,j)=alpha(j,j)*(1.+alamda)
  da(j)=beta(j)
end do
call gaussj(covar,mfit,da)
if(alamda==0.) then
  call covsrt(covar,nca,ma,lista,mfit)
  return
endif
do j=1,mfit
  atry(lista(j))=a(lista(j))+da(j)
end do
call mrqcof(x,y,sig,ndata,atry,ma,lista,mfit,covar,da,&
                nca,chisq,funcs)
if(chisq<ochisq) then
  alamda=0.1*alamda
  ochisq=chisq
  do j=1,mfit
    do k=1,mfit
      alpha(j,k)=covar(j,k)
    end do
    beta(j)=da(j)
    a(lista(j))=atry(lista(j))
  end do
else
  alamda=10.*alamda
  chisq=ochisq
endif
END SUBROUTINE mrqmin

SUBROUTINE mrqcof(x,y,sig,ndata,a,ma,lista,mfit,&
             alpha,beta,nalp,chisq,funcs)
PARAMETER (mmax=20)
DIMENSION x(ndata),y(ndata),alpha(nalp,nalp),&
   beta(ma),dyda(mmax),lista(mfit),sig(ndata),a(ma)
do j=1,mfit
  do k=1,j
    alpha(j,k)=0.
  end do
  beta(j)=0.
end do
chisq=0.
do i=1,ndata
  call funcs(x(i),a,ymod,dyda,ma)
  sig2i=1./(sig(i)*sig(i))
  dy=y(i)-ymod
  do j=1,mfit
    wt=dyda(lista(j))*sig2i
    do k=1,j
      alpha(j,k)=alpha(j,k)+wt*dyda(lista(k))
    end do
    beta(j)=beta(j)+dy*wt
  end do
  chisq=chisq+dy*dy*sig2i
end do
do j=2,mfit
  do k=1,j-1
    alpha(k,j)=alpha(j,k)
  end do
end do
END SUBROUTINE mrqcof

SUBROUTINE covsrt(covar,ncvm,ma,lista,mfit)
INTEGER ma,mfit,ncvm,lista(mfit)
REAL covar(ncvm,ncvm)
INTEGER i,j,k
REAL swap
do j=1,ma-1
  do i=j+1,ma
    covar(i,j)=0.
  end do
end do
do i=1,mfit-1
  do j=i+1,mfit
    if(lista(j)>lista(i)) then
      covar(lista(j),lista(i))=covar(i,j)
    else
      covar(lista(i),lista(j))=covar(i,j)
    endif
  end do
end do
swap=covar(1,1)
do j=1,ma
  covar(1,j)=covar(j,j)
  covar(j,j)=0.
end do
covar(lista(1),lista(1))=swap
do j=2,mfit
  covar(lista(j),lista(j))=covar(1,j)
end do
do j=2,ma
  do i=1,j-1
    covar(i,j)=covar(j,i)
  end do
end do
END SUBROUTINE covsrt

SUBROUTINE gaussj(a,n,b)
INTEGER n,NMAX
REAL a(n,n),b(n)
PARAMETER (NMAX=50)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),&
        indxr(NMAX),ipiv(NMAX)
REAL big,dum,pivinv
do j=1,n
  ipiv(j)=0
end do
do i=1,n
  big=0.
  do j=1,n
    if(ipiv(j)/=1) then
      do k=1,n
        if (ipiv(k)==0) then
          if (abs(a(j,k))>=big) then
            big=abs(a(j,k))
            irow=j
            icol=k
          endif
        else if (ipiv(k)>1) then
          pause 'singular matrix in gaussj'
        endif
      end do
    endif
  end do
  ipiv(icol)=ipiv(icol)+1
  if (irow/=icol) then
    do l=1,n
      dum=a(irow,l)
      a(irow,l)=a(icol,l)
      a(icol,l)=dum
    end do
    dum=b(irow)
        b(irow)=b(icol)
        b(icol)=dum
  endif
  indxr(i)=irow
  indxc(i)=icol
  if (a(icol,icol)==0.) pause 'singular matrix in gaussj'
  pivinv=1./a(icol,icol)
  a(icol,icol)=1.
  do l=1,n
    a(icol,l)=a(icol,l)*pivinv
  end do
  b(icol)=b(icol)*pivinv
  do ll=1,n
    if(ll/=icol) then
      dum=a(ll,icol)
      a(ll,icol)=0.
      do l=1,n
        a(ll,l)=a(ll,l)-a(icol,l)*dum
      end do
      b(ll)=b(ll)-b(icol)*dum
    endif
  end do
end do
do l=n,1,-1
  if(indxr(l)/=indxc(l)) then
    do k=1,n
      dum=a(k,indxr(l))
      a(k,indxr(l))=a(k,indxc(l))
      a(k,indxc(l))=dum
    end do
  endif
end do
END SUBROUTINE gaussj


这是之前有人贴上去的最小二乘法的代码,现在我想单独用这个程序,但是编译时老是出现这样的问题,
--------------------Configuration: gaussj - Win32 Debug--------------------
Linking...
main.obj : error LNK2001: unresolved external symbol [email=_GAUSSJ@0]_GAUSSJ@0[/email]
Debug/gaussj.exe : fatal error LNK1120: 1 unresolved externals
Error executing link.exe.

gaussj.exe - 2 error(s), 0 warning(s)

有谁知道这是怎么回事么??