c
c  note for The IBM parallel machines, 
c  it appears that having lower case for the 
c  function names and/or declaring EXTERNAL's 
c  carefully is necessary for successful running.
c
c
C  The coefficient of the exponential
C  is gotten from a study of a model
C  with simplified longitudinal structure.
C
C  This subroutine is used by the programs
C  monster.f and extrapolate.c
c
c   f77 -c weights.for
c
c  The gamma function is defined below 
c  to increase portability.
c
#define NOSU 0 /* for OSU alpha's, this should be zero */ 

#ifdef _CRAY
#define DGESVD SGESVD
#endif

#ifndef EISSVD
#define EISSVD 0
#endif

C
C  "interface" function for the gamma function
C
      double precision function mygamma(x)
c      integer i
c      data i/0/
#if 0 
      double precision s14aae,x
      integer ifail
      external s14aae
      mygamma=S14AAE(x,ifail)
#elif 0
      double precision s14aaf,x
      integer ifail
      external s14aaf
      mygamma=S14AAF(x,ifail)
c
c
#else
      double precision dgamma,x
      external dgamma
      mygamma=dgamma(x)
#endif
c      i=i+1
c      if(i.lt.10)then
c      write(*,*)'   mygamma(',x,')=',mygamma
c      end if
      end
C
C  "interface" function for the Bessel function J_n(x)
C  for negative order, must use recursion.
C
      double precision function mybesselj(n,x)
      integer ifail
c      integer i
c      data i/0/
C
C       There is a NAG routine, but out of laziness:
C
      double precision n,x,y(1)
c      save i
      external dbesj
      call dbesj(x,n,1,y,ifail)
      if(ifail.ne.0)then
         write(*,*)'Error in dbesj(n,x), order n=',n,' x=',x
      end if
      mybesselj=y(1)
c      i=i+1
c      if(i.lt.100)then
c      write(*,*)'   besselj(',n,',',x,')=',mybesselj
c      end if
      end
c
c
c  Find constant c used in extrapolating in particle 
c  truncation.  The form is the one discussed in the 2+1 paper.
c
c
      double precision function findc(beta,la1,la2)
c
      implicit none
      double precision pi,la1,la2,beta,linc
c                                          ,t,ln2,y,aijmax(maxsects)
      double precision mygamma,gumby
      external mygamma
      pi=4.0d0*datan(1.0d0)
      gumby=0.5d0
c
c   The following is d_p/d_{p-2} for a p-particle
c   basis state to leading order in p and 
c   exact in beta.
c
      linc=(2.0d0*(1.0d0 - gumby)*la1 + gumby*la2)*
     -     mygamma(0.5d0 + beta)**3*
     -     mygamma(1.0d0 + 2*beta)**2*
     -     mygamma(1.0d0 + 4.0d0*beta)/
     -     (mygamma(1.5d0 + 3.0d0*beta)*
c  The extra factor of 0.5 for the g^2 term is
c  from the new (standard) definition of g^2 subsequent to
c  The first paper 
     -     ((0.5d0 *(1.0d0 + 4.0d0*beta) + 2.0d0*(1.0d0 - gumby)*la1 
     -     + gumby*la2)*
     -     mygamma(0.5d0 + beta)**4*
     -     mygamma(1.0d0 + 4.0d0*beta)
c            turned 4 into 2 to match new coupling definition.
     -     + 2.0d0*pi*
     -     mygamma(1.0d0 + 2.0d0*beta)**4*Tan(beta*pi)))
c
c    This is the coefficient in the exponential
c    for the convergence in p.  The factor of 2
C    is from the number of states.
c
      linc= Log(2.0d0*Abs(linc))
c
c      write(*,*)'linc=',linc
c
      findc = MIN(linc,-0.1d0)
      return 
      end
c
c
c
c
c  Calculate the pseudoinverse of matrix a.
c  It assumes that the size of w is pretty small
c  since it has fixed memory set aside.
c
c
      SUBROUTINE PSEUDOINVERSE(A,TH,NSECTS)
c
      implicit none
      integer maxsects,nsects,maxth,th
      parameter(maxsects=40,maxth=10)
      integer i,j,k,ifail
      double precision s(maxsects),work(maxsects,maxth)
      double precision a(nsects,th)
#if EISSVD 
      logical matu,matv
      parameter(matu=.true.,matv=.true.)
      double precision v(maxsects*maxth)
      external SVD
#else
      double precision vt(maxth,maxsects)
      external DGESVD
#endif
      if(nsects.gt.maxsects) write(*,*)
     c                  'too many sectors in pseudoinverse'
      if(th.gt.maxth) write(*,*)'too many terms in function'
c
c  The pseudoinverse is only calculated when nsects is large
c  enough.  If nsects is one, then a weighting of one is assigned
c  to each element of w.
c
      if(nsects.ge.th)then
c
c
C      Use LAPACK subroutine DGESVD.
c
#if EISSVD 
c     Use EISPACK routine svd
         call  svd(nsects,nsects,th,a,s,matu,a,matv,v,ifail,work)
         if(ifail.ne.0)then
            write(*,*)'Error using SVD in weights.f, ifail=',ifail
         end if
#else 
         call DGESVD('O','A',nsects,th,a,nsects,s,a,nsects,vt,maxth,
     1        work,maxsects*maxth,ifail)
         if(ifail.ne.0)then
            write(*,*)'Error using DGESVD in weights.f, ifail=',ifail
            write(*,*)'nsects=',nsects,' th=',th
            write(*,*)'a=',a(1,1),a(2,1),a(1,2),a(2,2)
         end if
#endif
c      write(*,*)'optimal length=',rv1(1)
c
c
         do k=1,th
            do i=1,nsects
               work(i,k)=0.0
               do j=1,th
#if EISSVD
                  work(i,k)=work(i,k)+v((j-1)*nsects+k)*a(i,j)/s(j)
#else 
                  work(i,k)=work(i,k)+vt(j,k)*a(i,j)/s(j)
#endif
c
c      write(*,*)' a=',a(i,j),' s=',s(j),' vt=',vt(j,k)
c
               end do
c
c      write(*,*)'weights.F:  work(',i,',',k,')=',work(i,k)
c
            end do
         end do
      do i=1,nsects
         do j=1,th
            a(i,j)=work(i,j)
         end do
      end do
      else
         write(*,*)' Too few sectors in weights.F'
      end if
      return
      end
c
c
c
#if 0  
c        The following is just for testing:
c        At least on GNU it runs substantially slower 
c        than gcc  
C
C  These are FORTRAN versions of the routines called by
C  LASO.  They have not been tested.
C
C
      SUBROUTINE FORELOP(N,M,P,Q,ELIST,ILIST,JLIST,NLIST)
C
C THIS IS THE MATRIX MULTIPLY FOR THE MATRIX A.
c The sparse format is meant for C so that we
c have to shift indices in ILIST & JLIST by 1
c
      INTEGER I,J,JTEMP,K,M,N,NLIST
      INTEGER ILIST(*),JLIST(*)
      DOUBLE PRECISION P(N,M),Q(N,M),ELIST(*)
      IF(N.NE.NLIST)THEN
         WRITE(*,*)'ERROR IN FORTRAN ELOP'
         RETURN
      END IF
      DO K=1,M         
         DO I=1,N
            Q(I,K)=0.0D+00
         END DO
         DO I=1,N
            JTEMP=JLIST(I)+1
            IF(ILIST(JTEMP)+1.EQ.I)THEN
               Q(I,K)=Q(I,K)+ELIST(JTEMP)*P(I,K)
               JTEMP=JTEMP+1
            END IF
            DO J=JTEMP,JLIST(I+1)
               Q(ILIST(J)+1,K)=Q(ILIST(J)+1,K)+ELIST(J)*P(I,K)
               Q(I,K)=Q(I,K)+ELIST(J)*P(ILIST(J)+1,K)
            END DO
         END DO
      END DO
      RETURN
      END
C 
C     The following was simply copied from the 
C     documentation for DNLASO with QS added to
C     the argument list    
C
#endif
#if 0 /* Turned off */
      SUBROUTINE IOVECT(N,M,Q,J,K,QS)
C     
C     THIS SUBROUTINE STORES AND RECALLS VECTORS.  IT STORES THE VECTORS
C     IN CORE ALTHOUGH MOST REAL PROBLEMS WOULD USE A DISK.
      INTEGER M,N,K,L1,L,JJ
      DOUBLE PRECISION Q(N,M),QS(N,M)
      IF(K.EQ.1)GO TO 30
      DO 20 L=1,M
         L1=J-M+L
         DO 10 I=1,N
            QS(I,L1)=Q(I,L)
 10      CONTINUE
 20   CONTINUE
      RETURN
C     
 30   DO 50 L=1,M
         L1=J-M+L
         DO 40 I=1,N
            Q(I,L)=QS(I,L1)
 40      CONTINUE
 50   CONTINUE
      RETURN
      END
      
#endif