program specompnl
c     minimize (S-Sici)^2+alpha/sum(ci^4)
      implicit none
      integer*4 ns,nsn,np1,np,N,j,js,maxit,nt,m,ix,ib
      real*8 alpha,xs1,xs0,xs2,wmin,wmax,ITOL,tol,damp,
     1norm,n1,na,utol,g,an,dev,db
      integer*4,allocatable::u(:)
      real*8,allocatable::sx(:),sy(:),sally(:,:),c(:),
     1cc(:),sfit(:),syt(:),sf(:,:),sat(:,:),satt(:,:)
      character*170,allocatable::sn(:)
      character*10 st
      logical obyc,lwr,lstable,unscaled,usem,usec,ortho,lmax,lcon,
     1ladj,lran,lcg
c      
      write(6,6000)
6000  format(/,' Specompnl compares spectra listed in SP.LST',/,
     1        ' to the first spectrum in the list.',/)
      call iopar(alpha,xs0,xs1,xs2,wmin,wmax,obyc,lwr,tol,maxit,
     1ITOL,damp,lstable,unscaled,usem,usec,1,ortho,utol,lmax,lcon,g,
     1ladj,lran,lcg)
c
c     getn number of spectra:
      ns=nsn('SP.LST')
      write(6,*)ns,' spectra in SP.LST'
      N=ns-1

      allocate(sn(ns))

c     load spectral names:
      call lsn(ns,sn)

c     get number of points within wmin wmax in first
      np=np1(sn(1),wmin,wmax)
      write(6,*)np,' points in the first spectrum'
      if(np.lt.ns)write(6,*)' Warning: np < ns!'

c     load detailed scaling table:
      if(lstable)then
       nt=nsn('STABLE.TXT')
       allocate(sf(2,nt))
       call lstab(nt,sf)
      else
       nt=1
       allocate(sf(2,nt))
      endif

      allocate(sx(np),sy(np),sally(N,np),c(N),sfit(np),syt(np))

c     read decomposed spectrum
      call ls(sx,sy,sn(1),wmin,wmax)
      n1=norm(np,sy)
      write(6,599)n1
 599  format(' Decomposed spectrum norm:',e12.4)
      if(lwr)call wrs('ref_.prn',np,sx,sy,lwr)

c     read rest, adjust number of points and scale:
      na=0.0d0
      ib=2
      db=1.0d9
      do 7 js=2,ns
      call lsr(ns,sx,syt,np,js,sn,xs0,xs1,xs2,lstable,nt,sf,lcon,g)
      write(st,700)js
700   format(i10)
      if(lwr)then
       do 73 j=1,len(st)
73     if(st(j:j).ne.' ')goto 72
72     call wrs(st(j:len(st))//'_.prn',np,sx,syt,lwr)
      endif
      an=norm(np,syt)
      na=na+an
      dev=0.0d0
      do 77 ix=1,np
77    dev=dev+dabs(sy(ix)/n1-syt(ix)/an)
      if(js.eq.2)then
       write(6,76)
76     format(' Loading subspectra',/,
     1 ' Spectrum    norm    deviation from the first:')
       db=dev
      else
       if(dev.lt.db)then
        db=dev
        ib=js
       endif
      endif
      write(6,78)js,an,dev
78    format(i6,e12.4,f12.6)
c     save it to the buffer:
7     call sav(js-1,np,syt,sally,N)
      write(6,597)ib,sn(ib)(1:40)
597   format(' Best is number',i4,1x,a40)
      na=na/dble(N)
      write(6,598)N,na
598   format(i6,' subspectra loaded, average norm:',E12.4)

      if(ortho)then
c      orthogonalize subspectra:
       allocate(u(N),satt(N,np))
       call dort(u,sally,satt,N,np,utol,na,m)
       allocate(sat(m,np),cc(m))
c      fill m non-redundant into sat:
       call fills(satt,sat,N,m,np,na)
       deallocate(satt)
c      decompose sy into m orthogonal subspectra in sat:
       if(lran)then
        call e2(m,np,sat,sy,cc,n1,na,sx,ladj,alpha,maxit,lwr)
       else
        if(lcg)then
         call dcg(m,np,cc,alpha,sy,sfit,sat,n1,na,lwr,maxit,usem,sx)
        else
         call decomp(m,np,sat,sy,cc,n1,na,sx,sn,usem,obyc,usec,lmax)
        endif
       endif
c      reconstruct coefficients:
       c=0.0d0
       do 11 j=1,m
   11  c(u(j))=cc(j)
      else
c      decompose sy into N subspectra in sally:
       if(lran)then
        call e2(N,np,sally,sy,c,n1,na,sx,ladj,alpha,maxit,lwr)
       else
        if(lcg)then
         call dcg(N,np,c,alpha,sy,sfit,sally,n1,na,lwr,maxit,usem,sx)
        else
         call decomp(N,np,sally,sy,c,n1,na,sx,sn,usem,obyc,usec,lmax)
        endif
       endif
      endif

      call sort(c,N,sn,obyc)

      call wrc(N,c,sn)
c
      if(unscaled)then
       write(6,*)'Loading unscaled subspectra ...'
       do 71 js=2,ns
       call lsr(ns,sx,syt,np,js,sn,0.0d0,1.0d0,0.0d0,.false.,0,sf,
     1 lcon,g)
71     call sav(js-1,np,syt,sally,N)
       write(6,*)ns-1,' .. loaded'
       call mkfit(N,np,sfit,c,sally,obyc)
       call wrs('fitu.prn',np,sx,sfit,lwr)
      endif

      end

      function ef(sy,s,np)
      implicit none
      integer*4 i,np
      real*8 s(np),sy(np),ef
      ef=0.0d0
      do 1 i=1,np
  1   ef=ef+(s(i)-sy(i))**2
      return
      end

      subroutine adj(N,np,c,sally,obyc,n1)
c     adjust coefficients to norm:
      implicit none
      integer*4 N,np,i
      real*8 c(N),sally(N,np),norm,nf,n1
      real*8 ,allocatable::sfit(:)
      logical obyc
      allocate(sfit(np))
      call mkfit(N,np,sfit,c,sally,obyc)
      nf=norm(np,sfit)
      do 800 i=1,N
 800  c(i)=c(i)*n1/nf
      return
      end

      subroutine nrc(a,c,N)
      implicit none
      integer*4 N
      real*8 c(N)
      character*2 a
      write(6,600)a,c(1),c(2),c(3),c(N-1),c(N)
 600  format(a2,3E12.4,' ... ',2E12.4)
      return
      end

      subroutine sete(e,es,ec,sy,sfit,np,N,cold,as)
      implicit none
      integer*4 N,np
      real*8 e,es,ec,sy(np),sfit(np),cold(N),psu,p1,p2,as,ef
      es=ef(sy,sfit,np)
      p1=psu(1,cold,N)
      p2=psu(2,cold,N)
      ec=as*p1**2/p2
      e=ec+es
      return
      end

      subroutine setg(g,ssi,sij,N,c,as,gnorm)
      implicit none
      integer*4 N,i,j
      real*8 g(N),c(N),psu,p1,p2,as,p22,gnorm,s,sij(N,N),ssi(N)
      p1=psu(1,c,N)
      p2=psu(2,c,N)
      p22=p2**2
      gnorm=0.0d0
      do 1 j=1,N
      s=0.0d0
      do 2 i=1,N
   2  s=s+c(i)*sij(i,j)
      s=2.0d0*(s-ssi(j)+as*p1*(p2-p1*c(j))/p22)
      gnorm=gnorm+s**2
  1   g(j)=s
      gnorm=dsqrt(gnorm)
      return
      end

      subroutine dcg(N,np,c,alpha,sy,sfit,sally,n1,na,lwr,maxit,usem,
     1sx)
      implicit none
      integer*4 N,nstep,lstep,i,maxit,np
      real*8 c(N),as,maxstep,sy(np),sfit(np),sx(np),
     1e0,gnorm0,es,ec,e00,dh,secder,eold,h,hh0,sally(N,np),hh1,e2,hh2,
     1e1,gradh,alpha,n1,na,etot,enew,gnorm1
      logical lwr,usem
      real*8,allocatable::c0(:),g0(:),h0(:),ssi(:),sij(:,:),g(:)
      as=alpha*n1**2
      maxstep=0.1d0*n1/na
      nstep=0
      allocate(c0(N),g0(N),h0(N),sij(N,N),ssi(N),g(N))
      call getm(N,np,sij,ssi,sally,sy,usem)
      call inic(N,c0,.true.,n1,na,.false.)
      call adj(N,np,c0,sally,.true.,n1)
      call mkfit(N,np,sfit,c0,sally,.true.)
      call sete(etot,es,ec,sy,sfit,np,N,c0,as)
      call setg(g0,ssi,sij,N,c0,as,gnorm0)
      e0=etot
      h0=g0

222   e00=etot
      dh=maxstep/gnorm0
      secder=-1.0d0
      nstep=nstep+1
      if(lwr)write(6,600)nstep,etot
 600  format(' Linear search at point ',i6,e12.4)
      eold=etot
      e0=etot
      h=0.0d0
      hh0=0.0d0
      lstep=0

111   lstep=lstep+1
      h=h+dh
      do 1 i=1,N
1     c(i)=max(0.0d0,c0(i)-h*h0(i))
      call mkfit(N,np,sfit,c,sally,.true.)
      call sete(etot,es,ec,sy,sfit,np,N,c,as)
      enew=etot
c     possibly end of linear search?:
      if(dabs(enew-eold)<1.0d-5*eold)goto 444
      if( lstep.eq.1) then
       e1=etot
       hh1=h
      else
       if(lstep.eq.2) then
        e2=etot
        hh2=h
       else 
        e0=e1
        hh0=hh1
        e1=e2
        hh1=hh2
        e2=etot
        hh2=h
       endif
       if(hh1.eq.hh0.or.hh2.eq.hh0.or.hh2.eq.hh1)then
        secder=-1.0d0
       else
        secder=-2.0d0*((hh2-hh1)*(e1-e0)-(e2-e1)*(hh1-hh0))
     1  /((hh2-hh0)*(hh1-hh0)*(hh2-hh1))
        gradh=(e0-e1-secder/2.0d0*(hh0**2-hh1**2))/(hh0-hh1)
       endif
      endif
      if(secder.gt.0.0d0)then
       dh=-gradh/secder-h
       if(dabs(dh).gt.maxstep/gnorm0)then
        if(dh.gt.0.0d0)then
         dh=maxstep/gnorm0
        else
         dh=-maxstep/gnorm0
        endif
       endif
       eold=enew
       goto 111
      endif
      if(enew.gt.eold)then
       eold=enew
       dh=-dh/2
       goto 111
      endif
      eold=enew
      goto 111

444   if(lwr)write(6,601)lstep,etot
 601  format(i6,' linear steps ',e12.4)
      call mkfit(N,np,sfit,c,sally,.true.)
      call sete(etot,es,ec,sy,sfit,np,N,c,as)
      call setg(g,ssi,sij,N,c,as,gnorm1)
      c0=c
      do 2 i=1,n
      h0(i)=g(i)+(gnorm1/gnorm0)**2*h0(i)
   2  g0(i)=g(i)
      if(lwr)write(6,602)gnorm1,e00-etot
 602  format(' grad:',e12.4,', dE:',e12.4)
      call wrs('fit.prn',np,sx,sfit,.false.)
      if(dabs(e00-etot).lt.1.0d-5*e00)then
       write(6,*)'Optimization converged'
       goto 555
      endif
      if(nstep.gt.maxit)then
       write(6,*)'too many steps'
       goto 555
      endif
      e0=etot
      gnorm0=gnorm1
      goto 222
555   write(6,605)es,ec,etot
 605  format(' Errors, spectra:',e12.4,', coef:',e12.4,', total:',e12.4)
      return
      end

      subroutine e2(N,np,sally,sy,c,n1,na,sx,ladj,alpha,maxit,lwr)
      implicit none
      integer*4 N,np,iter,maxit,i,im
      real*8 sally(N,np),sy(np),c(N),n1,na,sx(np),e,alpha,as,en,dh,
     1p1,es,ec,psu,frandom,eold,h,error
      logical ladj,lwr
      real*8,allocatable::cold(:),sfit(:),cr(:),cn(:),del(:)
      as=alpha*n1**2
      iter=0
      allocate(cold(N),sfit(np),cr(N),cn(N),del(N))
      call inic(N,cold,.true.,n1,na,.false.)

      eold=0.0d0
9999  iter=iter+1
c     now fitted spectrum = sum(i) c(i)^xx S(i)
      if(ladj)call adj(N,np,cold,sally,.true.,n1)
      call mkfit(N,np,sfit,cold,sally,.true.)
      if(lwr)call nrc(' c',cold,N)
      call wrs('fit.prn',np,sx,sfit,.false.)
      call sete(e,es,ec,sy,sfit,np,N,cold,as)
      if(lwr)then
       write(6,6792)iter,es,ec
6792   format(i6,2e11.4)
      else
       write(6,6791)e
6791   format(e11.4,$)
       if(mod(iter,7).eq.0)write(6,*)
      endif

c     random change:
      p1=psu(1,cold,N)
      do 11 i=1,N
 11   cr(i)=max(0.0d0,cold(i)+frandom(i)*p1*0.01d0)
      if(ladj)call adj(N,np,cr,sally,.true.,n1)
      call mkfit(N,np,sfit,cr,sally,.true.)
      call sete(en,es,ec,sy,sfit,np,N,cr,as)
      write(6,601)en
 601  format(' After random change:',e12.4)

c     line search along cr-cold
      do 3 i=1,N
 3    del(i)=cr(i)-cold(i)
      if(en.lt.e)then
       dh=1.0d0
      else
       dh=-1.0d0
      endif
      h=0.0d0
      im=0
88    im=im+1
      h=h+dh
      if(im.gt.maxit)goto 882
      do 2 i=1,N
 2    cn(i)=max(0.0d0,cold(i)+h*del(i))
      if(ladj)call adj(N,np,cn,sally,.true.,n1)
c     if(lwr)call nrc('cn',cn,N)
      call mkfit(N,np,sfit,cn,sally,.true.)
      call sete(en,es,ec,sy,sfit,np,N,cn,as)
      write(6,603)im,h,dh,en,es,ec
603   format(i5,5E12.4)
      if(en.lt.e)then
       error=dabs((en-e)/e)
       e=en
       if(error.lt.1.0d-5)goto 882
       write(6,*)'match'
      else
       h=h-dh
       dh=0.5d0*dh
       if(dabs(dh).lt.1.0d-4)goto 882
      endif
      goto 88
882   write(6,602)im,es,ec,e
602   format(i6,' miniiterations ',3e12.4)
      e=en
      error=dabs((eold-e)/e)*100.0d0
      eold=e
      cold=cn
      if(error.gt.1.0d-4.and.iter.lt.maxit)goto 9999
      c=cn
      write(6,*)
      write(6,679)iter,error,es,ec
679   format(' Iter:',i4,' er:',f6.3,'% ser:',G12.4,' aer:',G12.4)
      if(lwr)call nrc(' c',cold,N)
      return
      end

      subroutine decomp(N,np,sally,sy,c,n1,na,sx,sn,usem,
     1obyc,usec,lmax)
      implicit none
      integer*4 i,j,IERR,N,np,iter,maxit
      real*8 sally(N,np),sy(np),c(N),n1,alpha,xs0,xs1,xs2,wmin,wmax,nf,
     1tol,utol,as,p,t,cn,error,sx(np),e1,e2,ser,psu,cer,na,ITOL,damp,g,
     1norm
      logical usem,obyc,lwr,lstable,unscaled,usec,ortho,lmax,lcon,
     1ladj,lran,lcg
      real*8,allocatable::sij(:,:),ssi(:),cold(:),e(:,:),si(:,:),
     1sfit(:),so(:,:),vo(:)
      character*170 sn(*)
      allocate(sij(N,N),ssi(N),cold(N),so(N,N),vo(N),e(N,N+N),si(N,N),
     1sfit(np))
c     correlation matrix:
      call getm(N,np,sij,ssi,sally,sy,usem)
c     initial values of coefficients and normalization factor:
      call inic(N,cold,obyc,n1,na,usec)

      iter=0
c
c     IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

2     iter=iter+1
      call iopar(alpha,xs0,xs1,xs2,wmin,wmax,obyc,lwr,tol,maxit,
     1ITOL,damp,lstable,unscaled,usem,usec,0,ortho,utol,lmax,lcon,g,
     1ladj,lran,lcg)
      as=alpha*n1**2
      so=sij
      vo=ssi
      p=psu(2,cold,N)
      t=psu(4,cold,N)

c     make correlation matrix and vector:
      if(.not.obyc)then
       do 200 i=1,N
       vo(i)=vo(i)*cold(i)
       do 200 j=1,N
200    so(i,j)=cold(i)*so(i,j)*cold(j)
       if(lmax)then
        do 202 i=1,N
202     so(i,i)=so(i,i)+as*cold(i)**2
       else
        do 201 i=1,N
201     so(i,i)=so(i,i)-as*p*(p*cold(i)**2-t)/t**2
       endif
      endif

c     invert and make new coefficients:
      IERR=0
      if(lwr)then
       write(6,*)'Matrix:'
       call wrm1(N,so)
      endif
      call INV(N,so,si,N,ITOL,e,IERR)
      if(lwr)then
       write(6,*)'Inverted:'
       call wrm1(N,si)
      endif
      if(IERR.ne.0) call report('unsuccessful inversion')

      do 11 i=1,N
      cn=0.0d0
      do 111 j=1,N
 111  cn=cn+si(j,i)*vo(j)
c     if(cn*cold(i).lt.0.0d0)cn=-cn
 11   c(i)=(1.0d0-damp)*cn + damp*cold(i)
      error=cer(cold,c,N)
      cold=c
      open(67,file='C.SCR',form='unformatted')
      write(67)c
      close(67)

c     now fitted spectrum = sum(i) c(i)^xx S(i)
      call mkfit(N,np,sfit,c,sally,obyc)

      if(ladj)then
c      norm of fitted spectrum:
       nf=norm(np,sfit)
c      adjust coefficients to norm:
       if(obyc)then
        do 800 i=1,N
 800    c(i)=c(i)*n1/nf
       else
        do 900 i=1,N
 900    c(i)=c(i)*dsqrt(n1/nf)
       endif
       call mkfit(N,np,sfit,c,sally,obyc)
      endif

      call wrs('fit.prn',np,sx,sfit,lwr)

      e1=ser(np,sy,sfit)
      if(lmax)then
       e2=as*t
      else
       e2=as*p**2/t
      endif
      if(lwr)then
       if(.not.ortho)call sort(c,N,sn,obyc)
       write(6,679)iter,error,e1,e2
679    format(' Iter:',i4,' cer(%):',G10.4,' ser:',G12.4,' aer:',G12.4)
      else
       write(6,6791)e1+e2
6791   format(e11.4,$)
       if(mod(iter,7).eq.0)write(6,*)
      endif
      if(dabs(error).gt.tol.and.iter.lt.maxit)goto 2
c     IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
      write(6,*)
      write(6,679)iter,error,e1,e2
      return
      end

      subroutine wrs(s,np,x,y,lwr)
      implicit none
      character*(*) s
      integer*4 i,np
      real*8 x(*),y(*)
      logical lwr
      open(67,file=s)
      do 1 i=1,np
1     write(67,900)x(i),y(i)
900   format(f9.2,g13.5)
      close(67)
      if(lwr)write(6,*)s
      return
      end

      function ser(np,sy,sfit)
      implicit none
      integer*4 np,i
      real*8 sy(*),sfit(*),s,ser
      s=0.0d0
      do 1 i=1,np
 1    s=s+dabs(sfit(i)-sy(i))
      ser=s
      return
      end

      subroutine retract(i,np,sy,sally,N)
      implicit none
      integer*4 i,N,np,j
      real*8 sy(*),sally(N,np)
      do 3 j=1,np
3     sy(j)=sally(i,j)
      return
      end

      function sp(np,si,sj)
c     scalar product of spectra i j
      implicit none
      integer*4 i,np
      real*8 si(*),sj(*),sp,s
      s=0.0d0
      do 1 i=1,np
  1   s=s+si(i)*sj(i)
      sp=s
      return
      end



      function cer(cold,c,n)
      implicit none
      integer*4 i,n
      real*8 cer,s,cold(*),c(*),sn
      s=0.0d0
      sn=0.0d0
      do 1 i=1,n
      sn=sn+cold(i)**2
 1    s=s+(cold(i)-c(i))**2
      cer=s/sn*100.0d0
      return
      end

      subroutine sav(is,np,sy,sally,N)
      integer*4 is,i,np,N
      real*8 sy(*),sally(N,np)
      do 1 i=1,np
 1    sally(is,i)=sy(i)
      return
      end

      subroutine ls(xs,ys,name,wmin,wmax)
c     loads in the spectrum  within wmin wmax    
      implicit none
      integer*4 i
      real*8 xs(*),ys(*),wmin,wmax,x,y
      character*(*) name
      i=0
      open(9,file=name)
 11   read(9,*,end=99,err=99)x,y
      if(wmin.le.x.and.x.le.wmax)then
       i=i+1       
       xs(i)=x
       ys(i)=y
      endif
      goto 11
  99  close(9)
      return
      end

      function xjf(xj,nt,sf)
      implicit none
      integer*4 nt,i
      real*8 xjf,xj,sf(2,nt),xe1,xe2,xc1,xc2
      do 1 i=1,nt-1
      xe1=sf(1,i )
      xe2=sf(1,i+1)
      xc1=sf(2,i )
      xc2=sf(2,i+1)
c     xe1,xe2 ... experimental values, xc1,xc2 ... corresponding calc.
      if((xj.ge.xc1.and.xj.le.xc2).or.
     1   (xj.le.xc1.and.xj.ge.xc2))then
       xjf=xe1+(xj-xc1)*(xe2-xe1)/(xc2-xc1)
       return
      endif
 1    continue
      if(dabs(xj-sf(2,nt)).lt.dabs(xj-sf(2,1)))then
       xjf=sf(1,nt)+xj-sf(2,nt)
      else
       xjf=sf(1,1)+xj-sf(2,1)
      endif
      return
      end

      function npall(s)
      implicit none
      integer*4 npall,n
      character*(*) s
      n=0
      open(9,file=s)
101   read(9,*,end=99,err=99)
      n=n+1
      goto 101
99    close(9)
      npall=n
      return
      end

      subroutine lsr(ns,xs,ys,np,js,sn,xs0,xs1,xs2,lstable,nt,sf,
     1lcon,g)
c     loads in the full spectrum, adjust number of points to np
c     with the scaled xs
      implicit none
      integer*4 np,i,js,npn,j,nt,npall,ns
      real*8 xs(*),ys(*),xi,xj,xjp,xs0,xs1,xs2,sf(2,nt),xjf,g,
     1ef,pi
      real*8 ,allocatable::sx(:),sy(:),sp(:)
      character*170 sn(ns)
      logical lstable,lcon
      pi=4.0d0*datan(1.0d0)
      npn=npall(sn(js))
      allocate(sx(npn),sy(npn),sp(npn))
c     loads in the full spectrum:
      call lsf(npn,sx,sy,sn(js))
c     convolute
      if(lcon)then
       do 3 i=1,npn
       sp(i)=0.0d0
       do 31 j=1,npn-1
       ef=((sx(i)-sx(j))/g)**2
31     if(ef.lt.20.0d0)sp(i)=sp(i)+exp(-ef)*sy(j)*dabs(sx(j)-sx(j+1))
3      sp(i)=sp(i)/dsqrt(pi)/g
       sy=sp
      endif
      do 1 i=1,np
      xi=xs(i)
      ys(i)=0.0d0
      do 2 j=1,npn-1
c     scaled values:
      if(lstable)then
       xj =xjf(sx(j  ),nt,sf)
       xjp=xjf(sx(j+1),nt,sf)
      else
       xj =xs0+xs1*sx(j  )+xs2*sx(j  )*2
       xjp=xs0+xs1*sx(j+1)+xs2*sx(j+1)*2
      endif
      if(xj.le.xi.and.xjp.ge.xi)then
       ys(i)=sy(j)+(sy(j+1)-sy(j))*(xi-xj)/(xjp-xj)
       goto 1
      endif
 2    continue
 1    continue
      return
      end

      function np1(s,wmin,wmax)
      implicit none
      integer*4 np1,n
      real*8 wmin,wmax,x
      character*(*) s
      open(9,file=s)
      n=0
101   read(9,*,end=99,err=99)x
      if(x.ge.wmin.and.x.le.wmax)n=n+1
      goto 101
99    close(9)
      np1=n
      return
      end

      function nsn(fn)
      implicit none
      integer*4 nsn,ns
      character*(*) fn
      open(9,file=fn,status='old')
      ns=0
101   read(9,*,end=99,err=99)
      ns=ns+1
      goto 101
99    close(9)
      nsn=ns
      return
      end

      subroutine lsn(ns,sn)
      implicit none
      integer*4 j,ns
      character*170 sn(ns)
      open(9,file='SP.LST')
      do 1 j=1,ns
1     read(9,9000,end=99,err=99)sn(j)
9000  format(a170)
99    close(9)
      write(6,*)'names loaded'
      return
      end

      subroutine lstab(ns,sf)
      implicit none
      integer*4 j,ns
      real*8 sf(2,ns)
      open(9,file='STABLE.TXT')
      do 1 j=1,ns
1     read(9,*,end=99,err=99)sf(1,j),sf(2,j)
99    close(9)
      write(6,*)ns,' assignments loaded'
      return
      end

      subroutine lsf(n,xs,ys,name)
c     loads in the full spectrum      
      implicit none
      integer*4 i,n
      real*8 xs(*),ys(*)
      character*(*) name
      open(9,file=name)
      do 11 i=1,n
 11   read(9,*)xs(i),ys(i)
      close(9)
      return
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine iopar(alpha,xs0,xs1,xs2,wmin,wmax,obyc,lwr,tol,maxit,
     1ITOL,damp,lstable,unscaled,usem,usec,ic,ortho,utol,lmax,lcon,g,
     1ladj,lran,lcg)
      implicit none
      integer*4 maxit,ic,iran,i
      real*8 alpha,xs0,xs1,xs2,wmin,wmax,ITOL,tol,damp,utol,g,frandom
      logical lex,obyc,lwr,lstable,unscaled,usem,usec,ortho,lmax,lcon,
     1ladj,lran,lcg
      character*10 s10
      integer*4 ix, iy, iz
      common /randc/ ix, iy, iz
      ix=1
      iy=1000
      iz=99
      alpha=1.0d-6
      damp=0.5d0
      wmin=0.0d0
      wmax=4000.0d0
      xs0=0.0d0
      xs1=1.0d0
      xs2=0.0d0
      obyc=.false.
      lwr=.false.
      tol=0.001d0
      utol=0.001d0
      ITOL=1.0d-20
      maxit=1000
      lstable=.false.
      unscaled=.true.
      usem=.false.
      usec=.false.
      ortho=.false.
      lmax=.false.
      lcon=.false.
      ladj=.false.
      lran=.false.
      g=20.0d0
      iran=0
      lcg=.false.
      inquire(file='SP.OPT',exist=lex)
      if(lex)then
       open(9,file='SP.OPT')
1      read(9,90,end=99,err=99)s10
90     format(a10)
       if(s10(1:4).eq.'alph')read(9,*)alpha
       if(s10(1:4).eq.'wmax')read(9,*)wmax
       if(s10(1:4).eq.'wmin')read(9,*)wmin
       if(s10(1:4).eq.'damp')read(9,*)damp
       if(s10(1:3).eq.'xs0')read(9,*)xs0
       if(s10(1:3).eq.'xs1')read(9,*)xs1
       if(s10(1:3).eq.'xs2')read(9,*)xs2
       if(s10(1:3).eq.'lwr')read(9,*)lwr
       if(s10(1:4).eq.'obyc')read(9,*)obyc
       if(s10(1:3).eq.'tol')read(9,*)tol
       if(s10(1:4).eq.'lsta')read(9,*)lstable
       if(s10(1:4).eq.'itol')read(9,*)ITOL
       if(s10(1:4).eq.'utol')read(9,*)utol
       if(s10(1:4).eq.'maxi')read(9,*)maxit
       if(s10(1:4).eq.'unsc')read(9,*)unscaled
       if(s10(1:4).eq.'usem')read(9,*)usem
       if(s10(1:4).eq.'usec')read(9,*)usec
       if(s10(1:4).eq.'orth')read(9,*)ortho
       if(s10(1:4).eq.'lmax')read(9,*)lmax
       if(s10(1:4).eq.'lcon')read(9,*)lcon
       if(s10(1:4).eq.'ladj')read(9,*)ladj
       if(s10(1:4).eq.'lran')read(9,*)lran
       if(s10(1:4).eq.'iran')read(9,*)iran
       if(s10(1:4).eq.'conj')read(9,*)lcg
       if(s10(1:3).eq.'ggg')read(9,*)g
       goto 1
99     close(9)
      endif
      if(ic.eq.1)write(6,6000) wmin,wmax,alpha,xs0,xs1,xs2,tol,ITOL,
     1utol,damp,g,lwr,lstable,unscaled,obyc,lmax,usem,usec,ortho,
     1ladj,lran,iran,lcg,lcon
6000  format(
     1' wmin:                         : ',f8.1,/,
     1' wmax:                         : ',f8.1,/,
     2' alpha                         : ',g8.2,/,
     1' xs0                           : ',f8.3,/,
     1' xs1                           : ',f8.3,/,
     1' xs2                           : ',f8.2,/,
     1' tolerance                     : ',g9.2,/,
     1' itolerance                    : ',g9.2,/,
     1' utolerance                    : ',g9.2,/,
     1' damping                       : ',g9.2,/,
     1' g                             : ',g9.2,/,
     1' lwr                           : ',l8,/,
     1' scaling table                 : ',l8,/,
     1' unscaling at the end          : ',l8,/,
     1' obyc                          : ',l8,/,
     1' lmax                          : ',l8,/,
     1' Use old M.SCR                 : ',l8,/,
     1' Use old C.SCR                 : ',l8,/,
     1' Pre-orthonormalize subspectra : ',l8,/,
     1' Adjust to norm                : ',l8,/,
     1' Random search                 : ',l8,/,
     1' Random seed                   : ',i8,/,
     1' Congujate gradient            : ',l8,/,
     1' Convoluted subspectra         : ',l8,/)
       do 2 i=1,iran
  2    if(frandom(i).gt.10.0d0)return
       return
      end


      SUBROUTINE INV(nr,a,ai,n,TOL,e,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension ai(nr,nr),a(nr,nr),e(nr,2*nr)
C
10000 IERR=0
      DO 1  ii=1,n
      DO 1  jj=1,n
      e(ii,jj)=a(ii,jj)
1     e(ii,jj+n)=0.0D0
      do 13 ii=1,n
13    e(ii,ii+n)=1.0D0
c
      DO 2 ii=1,n-1
      if (ABS(e(ii,ii)).LE.TOL) then
       DO 3 io=ii+1,n
3      if (ABS(e(io,ii)).GT.TOL) goto 11
       IERR=1
       write(6,*)ii
       write(6,*)'tol = ',tol
       tol=tol*0.50d0
       if(tol.gt.1.0d-20)goto 10000
       RETURN
c
11     CONTINUE
       DO 4 kk=1,2*n
       w=e(ii,kk)
       e(ii,kk)=e(io,kk)
4      e(io,kk)=w
      ENDIF
      eii=e(ii,ii)
      DO 5 jj=ii+1,n
      e1=e(jj,ii)/eii
      DO 6 kk=ii+1, 2*n
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e1
5     e(jj,ii)=0.0D0
2     CONTINUE
c
      DO 7  i2=n,2,-1
      eii=e(i2,i2)
      DO 7 j2=i2-1,1,-1
      e1=e(j2,i2)/eii
      DO 9 kk=1, n
9     e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1
7     e(j2,i2)=0.0d0
c
      DO 10 ii=1,n
      ei=1.0d0/e(ii,ii)
      DO 12  jj=1,n
12    ai(ii,jj)=e(ii,jj+n)*ei
10    CONTINUE
c
      RETURN
      END

      function psu(e,c,N)
      implicit none
      integer*4 e,N,i
      real*8 c(*),psu,p
      p=0.0d0
      do 4 i=1,N
4     p=p+c(i)**e
      psu=p
      return
      end

      subroutine wrc(N,c,sn)
      implicit none
      integer*4 N,i
      real*8 c(*)
      character*170 sn(*)
      open(61,file='COEF.TXT')
      do 3 i=1,N
3     write(61,678)i,c(i),sn(i+1)(1:40)
678   format(i4,f19.4,1x,a40) 
      close(61)
      return
      end

      subroutine sort(co,N,sn,obyc)
      implicit none
      real*8 co(*),ci,cmax,s
      integer*4 n,i,imax,j,it,is,ie
      integer*4 ,allocatable::ind(:)
      character*170 sn(*)
      real*8,allocatable::c(:)
      logical obyc
      allocate(c(N))
      if(obyc)then
       ie=1
      else
       ie=2
      endif
      s=0.0d0
      do 5 i=1,N
      s=s+co(i)**ie
5     c(i)=co(i)
      allocate(ind(N))
      do 10 i=1,N
  10  ind(i)=i
      do 1 i=1,N-1
      ci=c(i)
      imax=i
      cmax=dabs(ci)
      it=ind(i)
      do 2 j=i+1,N
      if(dabs(c(j)).gt.cmax)then
       imax=j
       cmax=dabs(c(j))
      endif
2     continue
      if(i.ne.imax)then
       c(i)=c(imax)
       c(imax)=ci
       ind(i)=ind(imax)
       ind(imax)=it
      endif
1     continue
      write(6,*)'Largest cofficients absulute/normalized: '
      do 3 i=1,min(10,N)
      is=ind(i)+1
3     write(6,679)i,ind(i),c(i)**ie,c(i)**ie/s,sn(is)(1:50)
679   format(2i4,E12.4,' /',f12.6,1x,a50) 
      return
      end      

      subroutine mkfit(N,np,sfit,c,sally,obyc)
      implicit none
      integer*4 ii,np,i,N
      real*8 sfit(*),c(*),sally(N,np)
      logical obyc
      if(obyc)then
       do 3 ii=1,np
       sfit(ii)=0.0d0
       do 3 i=1,N
3      sfit(ii)=sfit(ii)+c(i)*sally(i,ii)
      else
       do 1 ii=1,np
       sfit(ii)=0.0d0
       do 1 i=1,N
1      sfit(ii)=sfit(ii)+c(i)**2*sally(i,ii)
      endif
      return
      end

      subroutine getm(N,np,sij,ssi,sally,sy,usem)
      implicit none
      integer*4 N,np,N_u,i,j
      real*8 sij(N,N),sally(N,np),
     1sy(np),ssi(N),sp
      logical usem,lex
      real*8,allocatable::syi(:),syj(:)

      inquire(file='M.SCR',exist=lex)
77    if(lex.and.usem)then
       write(6,*)'Correlation matrix from M.SCR'
       open(9,file='M.SCR',form='unformatted')
       read(9)N_u
       if(N_u.ne.N)then
        write(6,*)'M.SCR not usable'
        usem=.false.
        close(9)
        goto 77
       endif
       read(9)((sij(i,j),i=1,N),j=1,N)
       read(9)(ssi(i),i=1,N)
       close(9)
      else
       write(6,*)'Making correlation matrix '
       allocate(syi(np),syj(np))
       do 5 i=1,N
       call retract(i,np,syi,sally,N)
       ssi(i)=sp(np,sy,syi)
       do 5 j=1,N
       call retract(j,np,syj,sally,N)
5      sij(i,j)=sp(np,syi,syj)
       open(9,file='M.SCR',form='unformatted')
       write(9)N
       write(9)((sij(i,j),i=1,N),j=1,N)
       write(9)(ssi(i),i=1,N)
       close(9)
      endif
      return
      end

      subroutine inic(N,c,obyc,n1,na,usec)
      implicit none
      integer*4 N,i
      real*8 c(N),n1,na
      logical obyc,usec
      do 1 i=1,N
      if(obyc)then
       c(i)=n1/na/dble(N)
      else
       c(i)=dsqrt(n1/na)/dsqrt(dble(N))
      endif
1     continue
      if(usec)then
       open(67,file='C.SCR',form='unformatted',status='old')
       read(67)c
       close(67)
      endif
      return
      end

      subroutine wrm(N,so)
      implicit none
      integer*4 N,i,j
      real*8 so(N,N)
      do 1 i=1,N
    1 write(6,600)(so(i,j),j=1,N)
  600 format(8E10.4)
      return
      end

      subroutine wrm1(N,so)
      implicit none
      integer*4 N
      real*8 so(N,N)
      if(N.lt.5)then
       call wrm(N,so)
      else
       write(6,600)so(1  ,1),so(1  ,2),so(1  ,N-1),so(1  ,N),
     1             so(2  ,1),so(2  ,2),so(2  ,N-1),so(2  ,N),
     1             so(N-1,1),so(N-1,2),so(N-1,N-1),so(N-1,N),
     1             so(N  ,1),so(N  ,2),so(N  ,N-1),so(N  ,N)
  600  format('  1,1   1,2 ..  1,N-1   1,N',2E11.3,' ..',2e11.3,/,
     1        '  2,1   2,2 ..  2,N-1   2,N',2E11.3,' ..',2e11.3,/,
     1        '   ..    .. ..     ..    ..',2(9x,'..'),' ..',
     1                                                 2(9x,'..'),/,
     1        'N-1,1 N-1,2 ..N-1,N-1 N-1,N',2E11.3,' ..',2e11.3,/,
     1        '  N,1   N,2 ..  N,N-1   N,N',2E11.3,' ..',2e11.3,/)
      endif
      return
      end

      function norm(np,s)
      implicit none
      integer*4 i,np
      real*8 norm,s(*)
      norm=0.0d0
      do 1 i=1,np
   1  norm=norm+dabs(s(i))
      return
      end

      function normq(np,s)
      implicit none
      integer*4 i,np
      real*8 normq,s(*)
      normq=0.0d0
      do 1 i=1,np
   1  normq=normq+s(i)**2
      normq=dsqrt(normq)
      return
      end

      subroutine dort(u,sally,satt,N,np,utol,na,m)
      implicit none
      integer*4 N,np,m,i,ii,j,u(N)
      real*8 sally(N,np),no,na,sp,sij,utol,normq,satt(N,np)
      real*8,allocatable::si(:),sj(:),sat(:,:)
      write(6,*)'orthogonalizing spectra'
      allocate(si(np),sj(np),sat(N,np))
      m=0
      sat=sally
      do 1 i=1,N
      call retract(i,np,si,sat,N)
      do 2 j=1,i-1
      call retract(j,np,sj,sat,N)
      sij=sp(np,si,sj)
      do 2 ii=1,np
2     si(ii)=si(ii)-sij*sj(ii)
      no=normq(np,si)
      if(no/na.gt.utol)then
       do 3 ii=1,np
3      si(ii)=si(ii)/no
c      take the original spectrum if something is left:
       m=m+1
       u(m)=i
c      save i-th spectrum to satt for further use:
       call retract(i,np,sj,sally,N)
       call sav(m,np,sj,satt,N)
      else
       si=0.0d0
      endif
c     save i-th spectrum to sat just for further orthogonalization:
1     call sav(i,np,si,sat,N)
      write(6,600)m,n
600   format(i6,' non-redundant spectra from ',i6)
      return
      end

      subroutine fills(satt,sat,N,m,np,na)
      implicit none
      integer*4 N,m,np,i,ii
      real*8 satt(N,np),sat(m,np),na,norm
      real*8,allocatable::si(:)
      allocate(si(np))
      na=0.0d0
      do 1 i =1, m
      call retract(i,np,si,satt,N)
      na=na+norm(np,si)
      do 1 ii=1,np
  1   sat(i,ii)=si(ii)
      na=na/dble(m)
      write(6,598)na
 598  format(' average norm:',E12.4)
      return
      end

      function frandom(i)
      implicit none
      real*8 frandom
      integer*4 i
c
c     Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2
c
c     Returns a pseudo-random number rectangularly distributed
c     between 0 and 1.   The cycle length is 6.95E+12 (See page 123
c     of Applied Statistics (1984) vol.33), not as claimed in the
c     original article.
c
c     IX, IY and IZ should be set to integer values between 1 and
c     30000 before the first entry.
c
      integer*4 ix, iy, iz
      common /randc/ ix, iy, iz
c
      ix = mod(171 * ix, 30269)
      iy = mod(172 * iy, 30307)
      iz = mod(170 * iz, 30323)
      frandom = mod(dble(ix) / 30269.0d0 + dble(iy) / 30307.0d0 +
     1                        dble(iz) / 30323.0d0, 1.0d0)
      return
      end