!*******************************************************
!*  Program to demonstrate multidimensional operation  *
!*  of the multi-nonlinear regression subroutine with  *
!*  iterative error reduction                          *
!* --------------------------------------------------- *
!*  Ref.:  BASIC Scientific Subroutines Vol. II, by    *
!*         F.R. Ruckdeschel, Byte/McGRAW-HILL, 1981    *
!*         [BIBLI 01].                                 *
!*                                                     *
!*                  F90 version by J-P Moreau, Paris   *
!*                         (www.jpmoreau.fr)           *
!* --------------------------------------------------- *
!* SAMPLE RUN:                                         *
!*                                                     *
!* MULTI-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION    *
!*                                                     *
!* How many data points ? 10                           *
!*                                                     *
!* How many dimensions ? 2                             *
!*                                                     *
!* What is the fit for dimension 1 ? 2                 *
!* What is the fit for dimension 2 ? 1                 *
!*                                                     *
!* Input the data as prompted:                         *
!*                                                     *
!* Y( 1) = ? 7                                         *
!* X( 1, 1) = ? 1                                      *
!* X( 1, 2) = ? 6                                      *
!*                                                     *
!* Y( 2) = ? 7                                         *
!* X( 2, 1) = ? 6                                      *
!* X( 2, 2) = ? 1                                      *
!*                                                     *
!* Y( 3) = ? 6                                         *
!* X( 3, 1) = ? 3                                      *
!* X( 3, 2) = ? 3                                      *
!*                                                     *
!* Y( 4) = ? 8                                         *
!* X( 4, 1) = ? 2                                      *
!* X( 4, 2) = ? 6                                      *
!*                                                     *
!* Y( 5) = ? 9                                         *
!* X( 5, 1) = ? 1                                      *
!* X( 5, 2) = ? 8                                      *
!*                                                     *
!* Y( 6) = ? 9                                         *
!* X( 6, 1) = ? 7                                      *
!* X( 6, 2) = ? 2                                      *
!*                                                     *
!* Y( 7) = ? 6                                         *
!* X( 7, 1) = ? 3                                      *
!* X( 7, 2) = ? 3                                      *
!*                                                     *
!* Y( 8) = ? 7                                         *
!* X( 8, 1) = ? 3                                      *
!* X( 8, 2) = ? 4                                      *
!*                                                     *
!* Y( 9) = ? 7                                         *
!* X( 9, 1) = ? 4                                      *
!* X( 9, 2) = ? 3                                      *
!*                                                     *
!* Y( 10) = ? 2                                        *
!* X( 10, 1) = ? 0                                     *
!* X( 10, 2) = ? 2                                     *
!*                                                     *
!*                                                     *
!* The calculated coefficients are:                    *
!*                                                     *
!*  1    0.0000                                        *
!*  2    1.0000                                        *
!*  3    0.0000                                        *
!*  4    1.0000                                        *
!*  5    0.0000                                        *
!*  6    0.0000                                        *
!*                                                     *
!* Standard deviation:  0.000000                       *
!*                                                     *
!* Number of iterations: 2                             *
!*                                                     *
!*******************************************************
PROGRAM REGITER

parameter (LMAX=9,MMAX=25,NMAX=25)

real*8 D(NMAX)
real*8 Y(MMAX)

integer i,j,l,l1,m,n
real*8 dd

real*8 X(MMAX,LMAX)              ! maxi l:=9  number of dimensions
real*8 Z(MMAX,NMAX+1)            ! maxi m:=25 number of data points
real*8 A(MMAX,MMAX)              ! maxi n:=25 order of regression
real*8 B(MMAX,2*MMAX)
real*8 C(MMAX,MMAX)

integer  M1(LMAX)


  print *,'MULTI-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION'
  
  write (*,'(" How many data points ? ")',advance='no')
  read *, m
  
  write (*,'(" How many dimensions  ? ")',advance='no') 
  read *, l

  do i = 1, l
    write(*,25,advance='no') i
    read *, M1(i)
  end do

  n = 1
  do i = 1, l  
    n = n * (M1(i) + 1)
  end do
  if (m < n+2) m = n+2;   ! m must be >= n+2 for LstSqrN()

  print *,'Input the data as prompted:'
  
  do i = 1, m
    write(*,30,advance='no') i
    read *, Y(i)
    
    do j = 1, l
      write(*,35,advance='no') i,j
      read *, X(i,j)
    end do
    print *, ' '
  end do

  !Call iterations supervisor subroutine
  call Iter_Supervisor(l,l1,m,n,dd,D,M1,X,Y)

  print *,'The calculated coefficients are:'
  
  do i = 1, n
    write(*,50)  i, D(i) 
  end do
  
  write(*,60) dd
  write(*,70) l1
  stop

25 format(' What is the fit for dimension ',I1,' ? ')
30 format(' Y(',I2,') = ? ')
35 format(' X(',I2,',',I2,') = ? ')
50 format('  ',I2,'  ',F10.4)
60 format(' Standard deviation: ',F10.6)
70 format(' Number of iterations: ',I2)

end ! main program


! Matrix transpose B = Tr(A)
! A and B are n by m matrices
Subroutine Transpose(n,m,A,B)
  parameter(MMAX=25,NMAX=25)
  integer n,m
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX)
  integer i,j
  do i = 1, n
    do j = 1, m
      B(i,j) = A(j,i)
    end do
  end do
  return
end

! Matrix save (A in C) 
Subroutine A_IN_C(n1,n2,A,C) 
  parameter(MMAX=25,NMAX=25)
  integer n1,n2
  real*8  A(MMAX,MMAX), C(MMAX,MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        C(i1,i2) = A(i1,i2)
      end do
    end do
  end if
  return
end

! Matrix save (B in A) 
Subroutine B_IN_A(n1,n2,A,B) 
  parameter(MMAX=25,NMAX=25)
  integer n1,n2
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        A(i1,i2) = B(i1,i2)
      end do
    end do
  end if
  return
end

! Matrix save (C in B)
Subroutine C_IN_B(n1,n2,B,C) 
  parameter(MMAX=25,NMAX=25)
  integer n1,n2
  real*8  C(MMAX,MMAX), B(MMAX,2*MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        B(i1,i2) = C(i1,i2)
      end do
    end do
  end if
  return
end 

! Matrix save (C in A)
Subroutine C_IN_A(n1,n2,A,C) 
  parameter(MMAX=25,NMAX=25)
  integer n1,n2
  real*8  A(MMAX,MMAX), C(MMAX,MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        A(i1,i2) = C(i1,i2)
      end do
    end do
  end if
  return
end 

! Matrix multiplication C = A * B
! A is m1 by n1, B is m2 by n2 ( m2 = n1 )
! Result C is m1 by n2.
Subroutine Matmult(m1,n1,n2,A,B,C) 
  parameter(MMAX=25,NMAX=25)
  integer m1,n1,n2
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX), C(MMAX,MMAX)  
  integer i,j,k
  do i = 1, m1
    do j = 1, n2	 
      C(i,j) = 0
      do k = 1, n1
        C(i,j) = C(i,j) + A(i,k) * B(k,j)
      end do
    end do
  end do
  return
end   

! Matrix inversion: B = Inv(A) by Gauss-Jordan method
! A and B are n by n matrices
Subroutine Matinv(n,A,B)  
  ! Labels: 10, 20, 30
  parameter(MMAX=25,NMAX=25)
  integer n
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX)  
  integer i,j,k
  real*8 bb

  do i = 1, n
    do j = 1, n
      B(i,j + n) = 0.d0
      B(i,j) = A(i,j)
    end do
    B(i,i + n) = 1.d0
  end do  

  do k = 1, n
    if (k.eq.n) goto 10
    m = k
    do i = k+1, n
      if (abs(B(i,k)) > abs(B(m,k)))  m = i
    end do
    if (m.eq.k) goto 10
    do j = k, 2*n
      bb = B(k,j)
      B(k,j) = B(m,j)
      B(m,j) = bb
    end do
10  do j = k+1, 2*n 
      B(k,j) = B(k,j) / B(k,k)
    end do
    if (k.eq.1) goto 20
    do i = 1, k-1
      do j = k+1, 2*n
        B(i,j) = B(i,j) - B(i,k) * B(k,j)
      end do
    end do
    if (k.eq.n) goto 30
20  do i = k+1, n
    do j = k+1, 2*n
      B(i,j) = B(i,j) - B(i,k) * B(k,j)
      end do
    end do
  end do    ! k loop

30  do i = 1, n
      do j = 1, n   
        B(i,j) = B(i,j + n)
      end do
    end do
  return
end ! Matinv()


!**********************************************************
!*   Least squares fitting subroutine, general purpose    *
!* subroutine for multidimensional, nonlinear regression. *
!* ------------------------------------------------------ *
!* The equation fitted has the form:                      *
!*     Y = D(1)X1 + D(2)X2 + ... + D(n)Xn                 *
!* The coefficients are returned by the program in D(i).  *  
!* The X(i) can be simple powers of x, or functions.      *      
!* Note that the X(i) are assumed to be independent.      *   
!* The measured responses are Y(i), there are m of them.  * 
!* Y is a m row column vector, Z(i,j) is a m by n matrix. *
!* m must be >= n+2. The subroutine inputs are m, n, Y(i) *
!* and Z(i,j) previously calculated. The subroutine calls *
!* several other matrix routines during the calculation.  *
!**********************************************************
Subroutine LstSqrN(m,n,A,B,C,D,Y,Z)
  parameter(MMAX=25,NMAX=25)
  integer m,n
  real*8 A(MMAX,MMAX),B(MMAX,2*MMAX),C(MMAX,MMAX)
  real*8 D(NMAX), Y(MMAX), Z(MMAX,NMAX+1) 
  integer i,j,m4,n4
  m4 = m
  n4 = n
  do i = 1, m
    do j = 1, n 
      A(i,j) = Z(i,j)
    end do
  end do
  call Transpose(n,m,A,B)                ! B=Transpose(A) 
  call A_IN_C(m,n,A,C)                   ! move A to C
  call B_IN_A(n,m,A,B)                   ! move B to A
  call C_IN_B(m,n,B,C)                   ! move C to B
  call Matmult(n,m,n,A,B,C)              ! multiply A and B, result in C
  call C_IN_A(n,n,A,C)                   ! move C to A
  call Matinv(n,A,B)                     ! B=Inverse(A) 
  m = m4                                 ! restore m
  call B_IN_A(n,n,A,B)                   ! move B to A
  do i = 1, m
    do j = 1, n
      B(j,i) = Z(i,j)
    end do
  end do
  call Matmult(n,n,m,A,B,C)              ! multiply A and B, result in C 
  call C_IN_A(n,m,A,C)                   ! move C to A
  do i = 1, m
    B(i,1) = Y(i)
  end do
  call Matmult(n,m,1,A,B,C)              ! multiply A and B

  ! Product C is n by 1 - Regression coefficients are in C(i,1)
  ! and put for convenience into vector D(i).
  do i = 1, n
    D(i) = C(i,1)
  end do
  return
end ! LstSqrN()

Subroutine S40(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1,j1;
  cc = bb
  j1 = j
  do i1 = 0, M1(1)
    j1 = j1 + 1  
    Z(i,j1) = bb
    bb = bb * X(i,1)
  end do
  bb = cc
  j = j1 
  return
end

Subroutine S50(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(2)
    call S40(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,2)
  end do
  bb = cc
  return
end

Subroutine S60(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(3)
    call S50(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,3)
  end do
  bb = cc
  return
end
  
Subroutine S70(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(4)
    call S60(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,4)
  end do
  bb = cc
  return
end

Subroutine S80(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(5)
    call S70(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,5)
  end do
  bb = cc
  return
end

Subroutine S90(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(6)
    call S80(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,6)
  end do
  bb = cc
  return
end

Subroutine S100(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(7)
    call S90(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,7)
  end do
  bb = cc
  return
end

Subroutine S110(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(8)
    call S100(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,8)
  end do
  bb = cc
  return
end

Subroutine S120(i,j,bb,cc,M1,X,Z) 
parameter (LMAX=9,MMAX=25,NMAX=25)
integer i,j
real*8  bb,cc
real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
integer M1(LMAX)   
integer i1;
  cc = bb
  do i1 = 0, M1(9)
    call S110(i,j,bb,cc,M1,X,Z)
    bb = bb * X(i,9)
  end do
  bb = cc
  return
end


!**********************************************************
!*   Least squares fitting subroutine, general purpose    *
!* subroutine for multidimensional, nonlinear regression  *
!* ------------------------------------------------------ *
!* The equation fitted has the form:                      *
!*     Y = D(1)X1 + D(2)X2 + ... + D(n)Xn                 *
!* The coefficients are returned by the program in D(i).  *
!* The X(i) can be simple powers of x, or functions.      *
!* Note that the X(i) are assumed to be independent.      *
!* The measured responses are Y(i), there are m of them.  *
!* Y is a m row column vector, Z(i,j) is a m by n matrix. *
!* m must be >= n+2. The subroutine inputs are m, n, Y(i) *
!* and Z(i,j) previously calculated. The subroutine calls *
!* several other matrix routines during the calculation.  *
!**********************************************************
Subroutine Gene_Coeffs_SD(m,l,bb,cc,dd,D,M1,X,Y,Z)  
  parameter (LMAX=9,MMAX=25,NMAX=25)
  integer i,ii,j,k,m,n,l
  real*8  bb,cc,dd,yy
  real*8  D(NMAX), Y(MMAX)
  real*8  X(MMAX,LMAX), Z(MMAX,NMAX+1)
  integer M1(LMAX)   
  
  ! Calculate the total number of dimensions
  n = 1
  do i = 1, l 
    n = n * (M1(i) + 1)
  end do

  Z =0.d0

  dd = 0.d0

  do i = 1, m
    !Branch according to dimension l (return if l > 9) 
    if (l < 1) then
      l = 0
      return
    end if
    if (l > 9) then
      l = 0
      return
    end if
    j = 0; ii = i
    if (l.eq.1) then
      bb=1
      call S40(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.2) then
      bb=1
      call S50(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.3) then
      bb=1
      call S60(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.4) then
      bb=1
      call S70(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.5) then
      bb=1
      call S80(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.6) then
      bb=1
      call S90(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.7) then
      bb=1
      call S100(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.8) then
      bb=1
      call S110(ii,j,bb,cc,M1,X,Z)
    end if
    if (l.eq.9) then
      bb=1
      call S120(ii,j,bb,cc,M1,X,Z)
    end if

    yy = 0
    do k = 1, n
      yy = yy + D(k) * Z(i,k)
    end do
    dd = dd + (Y(i) - yy) * (Y(i) - yy)
  end do  !i loop

  !Calculate standard deviation (if m > n) 
  if (m - n < 1) then 
    dd = 0
  else 
    dd = dd / (m - n)
    dd = dsqrt(dd)
  end if

  return

end ! Gene_Coeffs_SD()


!********************************************************************
!*   Multi-dimensional polynomial regression iteration subroutine   *
!* ---------------------------------------------------------------- *
!* This routine supervises the calling of several other subroutines *
!* in order to iteratively fit least squares polynomials in more    * 
!* one dimension.                                                   *
!* The routine repeatedly calculates improved coefficients until    *
!* the standard deviation is no longer reduced. The inputs to the   *
!* subroutine are the number of dimensions l, the degree of fit for *
!* each dimension m(i), and the input data, x(i) and y(i).          *
!* The coefficients are returned in d(i), with the standard devia-  *
!* tion in d. Also returned is the number of iterations tried, l1.  *
!* y1(i), d1(i) and d1 are used respectively to save the original   *
!* values of y(i) and the current values of d(i) and d.             *
!********************************************************************
Subroutine Iter_Supervisor(l,l1,m,n,dd,D,M1,X,Y)
!Labels: 50, 100
parameter (LMAX=9,MMAX=25,NMAX=25)

integer i,l,l1,m,n
real*8  A(MMAX,MMAX)              
real*8  B(MMAX,2*MMAX)
real*8  C(MMAX,MMAX)
real*8  D(NMAX), D1(NMAX)
real*8  X(MMAX,LMAX)  
real*8  Y(MMAX), Y1(MMAX)
real*8  Z(MMAX,NMAX+1)    
integer M1(LMAX)
real*8  bb,cc,dd,dd1

  l1 = 0
  !Save the y(i)
  do i = 1, m  
    Y1(i) = Y(i)
  end do
  !Zero d1(i) 
  do i = 1, n  
    D1(i) = 0.d0
  end do
  !Set the initial standard deviation high
  dd1 = 1.d7
  !Call coefficients subroutine
50 call Gene_Coeffs_SD(m,l,bb,cc,dd,D,M1,X,Y,Z)

  !Call regression subroutine
  call LstSqrN(m,n,A,B,C,D,Y,Z)

  !Get standard deviation
  call Gene_Coeffs_SD(m,l,bb,cc,dd,D,M1,X,Y,Z)

  !if standard deviation is decreasing, continue
  if (dd1 > dd)  goto 100
  !Terminate iteration
  do i = 1, n
    D(i) = D1(i)
  end do
  !Restore the y(i)
  do i = 1, m
    Y(i) = Y1(i)
  end do
  !Get the final standard deviation
  call Gene_Coeffs_SD(m,l,bb,cc,dd,D,M1,X,Y,Z)
  return
  !Save the standard deviation
100 dd1 = dd
  l1 = l1 + 1
  !Augment coefficient matrix
  do i = 1, n
    D(i)  = D1(i) + D(i);
    D1(i) = D(i)
  end do
  !Restore y(i)
  do i = 1, m
    Y(i) = Y1(i)
  end do
  !Reduce y(i) according to the d(i) 
  call S250(l,m,n,bb,cc,D,M1,X,Y,Z)
  !We now have a set of error values
  goto 50
end ! Iter_Supervisor()

  !begin internal subroutines used by Iter_Supervisor()
  Subroutine S160(i,j,bb,cc,M1,X,Z) 
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)  
  real*8 Z(MMAX,NMAX+1)    
  integer M1(LMAX)
    cc = bb
    j1 = j
    do i1 = 0, M1(1)
      j1 = j1 + 1
      Z(i,j1) = bb
      bb = bb * X(i,1)
    end do
    bb = cc
    j = j1
    return
  end

  Subroutine S170(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i2 = 0, M1(2)
      call S160(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,2)
    end do
    bb = cc
    return
  end

  Subroutine S180(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i3 = 0, M1(3)
      call S170(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,3)
    end do
    bb = cc
    return
  end

  Subroutine S190(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i4 = 0, M1(4)
      call S180(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,4)
    end do
    bb = cc
    return
  end

  Subroutine S200(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i5 = 0, M1(5)
      call S190(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,5)
    end do
    bb = cc
    return
  end

  Subroutine S210(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i6 = 0, M1(6)
      call S200(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,6)
    end do
    bb = cc
    return
  end

  Subroutine S220(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i7 = 0, M1(7)
      call S210(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,7)
    end do
    bb = cc
    return
  end

  Subroutine S230(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    cc = bb
    do i8 = 0, M1(8)
      call S220(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,8)
    end do
    bb = cc
    return
  end

  Subroutine S240(i,j,bb,cc,M1,X,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)   
    do i9 = 0, M1(9)
      call S230(i,j,bb,cc,M1,X,Z)
      bb = bb * X(i,9)
    end do
    bb = cc
    return
  end

  Subroutine S250(l,m,n,bb,cc,D,M1,X,Y,Z)
  parameter (LMAX=9,MMAX=25,NMAX=25)  
  real*8 bb,cc
  real*8 X(MMAX,LMAX)
  real*8 D(NMAX),Y(MMAX)
  real*8 Z(MMAX,NMAX+1)   
  integer M1(LMAX)      
  real*8 yy
	do i = 1, m
      j = 0
      if (l.eq.1) then
	bb=1
	call S160(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.2) then
	bb=1
	call S170(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.3) then
	bb=1
	call S180(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.4) then
	bb=1
	call S190(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.5) then
	bb=1
	call S200(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.6) then
	bb=1
	call S210(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.7) then
	bb=1
	call S220(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.8) then
	bb=1
	call S230(i,j,bb,cc,M1,X,Z)
      end if
      if (l.eq.9) then
	bb=1
	call S240(i,j,bb,cc,M1,X,Z)
      end if
      !Array generated for row i
      yy = 0.d0
      do k = 1, n
        yy = yy + D(k) * Z(i,k)
      end do
      Y(i) = Y(i) - yy
    end do !i loop
    return
  end !S250()

! End of file regiter.f90