c
c       REQUIRED IMPROVEMENTS:
c       (1)     Move the parameters which must be here for this
c               version of fortran to include files
c       (2)     Move the common array definitions to include files
c       (3)     Remove any extraneous arrays which are not used
c               in the subroutines.
c       (4)     Calculate the maximum number of loops on the fly.
c       (5)     
c
c       This routine factors a tridiagonal system which is then solved 
c       in the subroutine factapp (rename this function).  The tridiagonal
c       system is assumed to have the following format:
c
c       a1  b1  c1                                       k1
c           a2  b2  c2                                   k2
c               a3  b3  c3                               k3
c                   .                                    .
c                   .                         [X]     =  .
c                   .                                    .
c                          an-2  bn-2  cn-2              kn-2
c                                an-1  bn-1  cn-1        kn-1
c                                      an    bn    cn    kn
c
c       The following restrictions apply to this subroutine
c       (1)     n must be the size declared for the arrays a, b, c, and k
c       (2)     The value for the parameters nequat, l2nequat, float_size,
c               and nsys must be defined by equivalent values in this routine
c               the factapp routine, and the program where the arrays
c               a, b, c, and k are defined.
c       (3)     The arrays a, b, and c are destroyed during this processing
c               so if they are going to be used again, then you need to
c               pass in temporary arrays
c       (4)     The storage locations pointed to be factor1 and factor2
c               must be declared externally to this subroutine.
c       (5)     The arrays a, b, c, and k must all be in the same vp set
c               and the locations pointed to by factor1 and factor2
c               must be in the same vp set as these four arrays or 
c               the results will be unpredictable and almost certainly 
c               wrong.
c
c       The following features exist for this subroutine
c       (1)     The subroutine checks the values of the diagonal dominance
c               to stop before the complete set of log2(l) iterations are
c               complete.
c
c       The inputs to this subroutine are:
c       (1)     a, b, c, and k as defined in the introduction.  These
c               can be two dimensional arrays with the rows representing
c               multiple tridiagonal systems to solve.  To repeat
c               each row is assumed to be an independent tridiagonal
c               system to solve.
c       (2)     factor1 is a front-end integer array of size log2(l) + 1
c               which holds addresses to put the values of the factors which
c               will be applied in factapp.  The values which are stored in
c               the location defined by factor1(k) k<= loop_cnt are the
c               factors used to reduce the lower tridiagonal system.  
c               The values
c               stored at factor1(loop_cnt + 1) is the normalizer at the
c               end of the solution.
c       (3)     factor2 is the same as factor1 except it is used to 
c               eliminate the upper portion of the tridiagonal system.
c       (4)     l is log2(number of variables in this system)  (I am
c               passing it here, but later I may calculate it on the fly.
c               7/18/89)
c       (5)     epsilon is the error bound which is used to check if the
c               process can stop before the complete log2 number of 
c               iterations.
c       Outputs
c       (1)     At the locations pointed to by factor1 and factor2 are
c               the factors which are applied during the factapp subroutine
c       (2)     loop_cnt is the number of iterations which are needed
c               to find the solution in factapp.
c
c       See (reference) for a description of the algorithm which
c       is used in this derivation.  Generally this would be called
c       a parallel cyclic reduction or a generalized even-odd 
c       cyclic reduction.c
c
        subroutine trifact(a,b,c,nsys,nequat,
     1             factor1,factor2,log2_nequat,
     1             epsilon,loop_cnt)
        implicit none
        include '/usr/include/cm/paris-configuration-fort.h'
        integer nsys, nequat, log2_nequat
c       real a(nsys,nequat), b(nsys,nequat)
        integer a,b
c       real c(nsys,nequat)
        integer c
cmf$    layout a(:news,:news), b(:news,:news)
cmf$    layout c(:news,:news)
        integer factor1(log2_nequat+1)
        integer factor2(log2_nequat+1)
cmf$    layout factor1(:serial)
cmf$    layout factor2(:serial)
        integer j,i
        integer vp
c       real t0(nsys,nequat)
        integer t0
        integer t2
c       real t4(nsys,nequat), t5(nsys,nequat)
        integer t4,t5
c       real t6(nsys,nequat)
        integer t6
cmf$    layout t0(:news,:news)
cmf$    layout t4(:news,:news)
cmf$    layout t5(:news,:news), t6(:news,:news)

        integer lbound
        real epsilon
	double precision tepsilon
        integer loop_cnt
        logical rpt
	call CM_set_context
        t0=CM_allocate_stack_field(32)
        t2=CM_allocate_stack_field(32)
        t4=CM_allocate_stack_field(32)
        t5=CM_allocate_stack_field(32)
        t6=CM_allocate_stack_field(32)
        rpt = .true.
        i=1
c
c       This starts the main processing loop
c
        lbound=1
	call CM_my_news_coordinate_1L(t4,0,32)
        do while ( (rpt .eq. .true.) .and. (i .le. log2_nequat))
        factor1(i)=cm_allocate_heap_field(64)
        factor2(i)=cm_add_offset_to_field_id(factor1(i),32)
c
c       NOTICE:  t0 is defined with a negative sign so that all of the
c       solutions can use addition instead of subtraction.  This must
c       also be carried through to the factapp routine.  Therefore if
c       this changes you must also change the factapp routine.
c       
c          t0=-1.0/b
c          t6=a
c          t5=c
c          t1=a*eoshift(t0,dim=2,shift=-lbound)
c          tepsilon=maxval(abs(t1))
       call cm_f_divinto_const_always_3_1L(t0,b,-1.0,23,8)
       call cm_f_move_always_1L(t6,a,23,8)
       call cm_f_move_always_1L(t5,c,23,8)
       call cm_get_from_power_two_always_1L(factor1(i),t0,0,i-1,
     1      cm_downward,32)
       call cm_u_lt_constant_1L(t4,lbound,32)
       call cm_logand_context_with_test
       call cm_f_move_zero_1L(factor1(i),23,8)
       call cm_set_context
       call cm_f_multiply_always_2_1L(factor1(i),a,23,8)
       call cm_f_abs_2_1L(t2,factor1(i),23,8)
       tepsilon=cm_global_f_max_1L(t2,23,8)
c
c       the factor application actuatlly performs the following step
c       which is commented out.  This is for the lower diagonal.
c       
c       
c
c       Note the addition in the next two steps, in the normal derivation
c       this is subtraction but I negated t0 first to eliminate the need
c       for subtraction.
c
c          b = b + t1*cshift(t5,dim=2,shift=-lbound)
c          a = t1*cshift(t6,dim=2,shift=-lbound)
       call cm_get_from_power_two_always_1L(t2,t5,0,i-1,
     1      cm_downward,32)
       call cm_f_mult_add_1L(b,factor1(i),t2,b,23,8)
       call cm_get_from_power_two_always_1L(a,t6,0,i-1,
     1      cm_downward,32)
       call cm_f_multiply_always_2_1L(a,factor1(i),23,8)
c       
c          t1=c*eoshift(t0,dim=2,shift=lbound)
c          tepsilon=max(tepsilon,maxval(abs(t1)))
       call cm_get_from_power_two_always_1L(factor2(i),t0,0,i-1,
     1      cm_upward,32)
       call cm_u_ge_constant_1L(t4,nequat-lbound,32)
       call cm_logand_context_with_test
       call cm_f_move_zero_1L(factor2(i),23,8)
       call cm_set_context
       call cm_f_multiply_always_2_1L(factor2(i),c,23,8)
       call cm_f_abs_2_1L(t2,factor2(i),23,8)
       tepsilon=max(cm_global_f_max_1L(t2,23,8),
     2              tepsilon)
c
c          b = b + t1*cshift(t6,dim=2,shift=lbound)
c          c = +t1*cshift(t5,dim=2,shift=lbound)
c
       call cm_get_from_power_two_always_1L(t2,t6,0,i-1,
     1      cm_upward,32)
       call cm_f_mult_add_1L(b,factor2(i),t2,b,23,8)
       call cm_get_from_power_two_always_1L(c,t5,0,i-1,
     1      cm_upward,32)
       call cm_f_multiply_always_2_1L(c,factor2(i),23,8)

        lbound=lbound*2
           if(tepsilon .le. epsilon) rpt = .false.
           i=i+1
        end do
        loop_cnt=i-1
C       t1=1.0/b
c
c       Move the normalizer into to last address pointed to
c       by factor1.  This is the reason that factor1 is dimensioned
c       one more than "necessary"
c       
        factor1(i)=cm_allocate_heap_field(32)
	call cm_f_divinto_const_always_3_1L(factor1(i),b,1.0,23,8)
        call cm_deallocate_stack_through(t0)
        return
        end
