c-----------------------------------------------------------------------
c Orthonormal functions. Numerical approach based on LU decomposition.
c Link with either Numerical Recipes (QTRAP, CHOLDC, and GAUSSJ) or
c                  IMSL              (QDAGS, LCHRG,  and LINRG)
c      for integration, matrix decomposition, and matrix inverse.
c Link with either Numerical Recipes (BESSJ0) or
c                  IMSL              (BSJ0)
c      for Bessel function of first kind of order 0
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

      parameter (n=5)
c     dimension P(n,n), diag(n), A(n,n), b(n)  !for Numerical Recipes
      dimension P(n,n), U(n,n),  A(n,n)        !for IMSL
      external prod
      common i, j

c Program Header -------------------------------------------------------
      print *,'Construct orthonormal functions e from the following ',
     *         'given functions'
      print *,'  f1(x)=J0(x) '
      print *,'  f2(x)=J0(2*x)'
      print *,'  f3(x)=J0(3*x)'
      print *,'  f4(x)=J0(4*x)'
      print *,'  f5(x)=J0(5*x)'
      print *,' '
      print *,'The scalar product (f,g) is defined as'
      print *,'          ô1             '
      print *,'  (f,g) = ³ f(x)*g(x)*dx '
      print *,'          õ0             '
      print *,' '

c Find the scalar product between each pair of f
      do i=1, n
        do j=1, n
c         call a routine to find the integral
c         call QTRAP(prod, 0., 1., P(i,j))  !from Numerical Recipes, wrong answers!
          call QDAGS(prod, 0., 1., 1.e-2, 1.e-2, P(i,j), error) !from IMSL
        enddo
      enddo

      print *, 'Scalar product matrix P'
      do i=1, n
        print '(1p5e15.6)', (P(i,j),j=1,n)
      enddo

c IMSL's LCHRG decomposes P into U, an upper triangular matrix, such that P=U'*U

c Eigenvalue lookup when Cholesky LU decomposition does not work.
c Examine the eigenvalues to make sure that the matrix is positive semidefinite.
c     print *, (U(i,1),i=1,n)
c     call EVLRG(n, P, n, U)

      call LCHRG(n, P, n, .false., ipvt, U, n)
c A=U'  (extract only the upper triangular elements)
      do i=1, n
        do j=1, i
          A(i,j)=U(j,i)
        enddo
      enddo
c Invert A
      call LINRG(n, A, n, A, n)

c-----------------------------------------------------------------------
c Numerical Recipes' choldc() returns the lower triangular matrix L such that L*L'=P
c     call CHOLDC(P, n, n, diag)
c     The lower off-diagonal elements are returned in P, and the diagonal elements are returned in diag
c     A=P
c     do i=1, n
c       A(i,i)=diag(i)
c       do j=i+1, n
c         A(i,j)=0.
c       enddo
c     enddo
c Invert A
c     call GAUSSJ(A, n, n, b, 1, 1)
c-----------------------------------------------------------------------

c Output results
      print *, 'Orthonormal basis functions e=A*f'
      print *, 'where A is...'
      do i=1, n
        print 650, (A(i,j), j=1, n)
      enddo

650   format (5f9.3)
      end



c-----------------------------------------------------------------------
      function prod(x)
c-----------------------------------------------------------------------
c Return the product of a set of given basis functions.
c-----------------------------------------------------------------------
      dimension f(5)
      common i, j

c For Numerical Recipes
c     f(1)=BESSJ0(x)
c     f(2)=BESSJ0(2.*x)
c     f(3)=BESSJ0(3.*x)
c     f(4)=BESSJ0(4.*x)
c     f(5)=BESSJ0(5.*x)
c For IMSL
      f(1)=BSJ0(x)
      f(2)=BSJ0(2.*x)
      f(3)=BSJ0(3.*x)
      f(4)=BSJ0(4.*x)
      f(5)=BSJ0(5.*x)

      prod = f(i)*f(j)

      end