      integer function continued_fraction(a,b,fraction,mask,size,iax)
      implicit none
      include '/usr/include/cm/paris-configuration-fort.h'
c
c number of equations per system
      integer size
c axis number along which to solve
      integer iax
c id's of real arrays
      integer a,b,fraction
c id of logical mask
      integer mask
c id's of shifted copies of a,b
      integer ashift,bshift
c id's of local real arrays
      integer alpha,beta,gamma
c and shifted copies
      integer alphashift,betashift,gammashift
c id's of local real arrays
      integer alpha1,beta1,gamma1,temp1,temp2
c id's of local logical arrays
      integer mask1,mask2
c
      integer i,j,k
      integer ilog2
c
c prime the pump with the first loop
c
      ashift=CM_allocate_stack_field(32)
      bshift=CM_allocate_stack_field(32)
      call CM_get_from_news_always_1L(ashift,a,iax,
     1    CM_downward,32)
      call CM_get_from_news_always_1L(bshift,b,iax,
     1    CM_downward,32)
c
      alpha=CM_allocate_stack_field(32)
      beta=CM_allocate_stack_field(32)
      gamma=CM_allocate_stack_field(32)
      call CM_set_context
      call CM_f_divinto_const_always_3_1L(gamma,a,1.0,23,8)
      call CM_f_mult_add_1L(beta,b,gamma,ashift,23,8)
      call CM_f_news_mult_add_4_1L(alpha,b,gamma,a,iax,
     1   CM_downward,23,8)
      call CM_f_move_always_1L(gamma,ashift,23,8)
      call CM_load_context(mask)
      call CM_f_move_1L(alpha,a,23,8)
      call CM_set_context
c
      call CM_get_from_news_always_1L(b,b,iax,
     1   CM_upward,32)
c
      mask2=CM_allocate_stack_field(1)
      call CM_get_from_power_two_always_1L(mask2,mask,iax,1,
     1   CM_downward,1)
      alphashift=CM_allocate_stack_field(32)
      betashift=CM_allocate_stack_field(32)
      gammashift=CM_allocate_stack_field(32)
      call CM_my_news_coordinate_1L(alphashift,iax,32)
      call CM_u_le_constant_1L(alphashift,1,32)
      call CM_logand_context_with_test
      call CM_clear_bit(mask2)
c
      call CM_get_from_power_two_always_1L(alphashift,alpha,iax,1,
     1   CM_downward,32)
      call CM_get_from_power_two_always_1L(bshift,b,iax,1,
     1   CM_downward,32)
c
      alpha1=CM_allocate_stack_field(32)
      beta1=CM_allocate_stack_field(32)
      gamma1=CM_allocate_stack_field(32)
      call CM_load_context(mask2)
      call CM_f_divide_3_1L(beta,bshift,beta,23,8)
      call CM_f_divide_3_1L(gamma,bshift,gamma,23,8)
      call CM_f_add_3_1L(beta1,alphashift,beta,23,8)
      call CM_f_add_3_1L(gamma1,alphashift,gamma,23,8)
      call CM_f_divide_3_1L(alpha1,beta1,gamma1,23,8)
      call CM_f_multiply_2_1L(alpha1,alpha,23,8)
      call CM_f_move_1L(alpha,alpha1,23,8)
      call CM_f_move_1L(beta,beta1,23,8)
      call CM_f_move_1L(gamma,gamma1,23,8)
c
      call CM_load_context(mask)
      call CM_logior_context(mask2)
      call CM_get_from_news_always_1L(mask,mask,iax,CM_downward,1)
      call CM_logior_context(mask)
      call CM_store_context(mask)
c
      ilog2=0
      temp1=CM_allocate_stack_field(32)
      temp2=CM_allocate_stack_field(32)
 10   continue
      if(1.eq.CM_global_logand_bit_always(mask))goto 20
      ilog2=ilog2+1
      call CM_get_from_power_two_always_1L(alphashift,alpha,
     1   iax,ilog2,CM_downward,32)
      call CM_get_from_power_two_always_1L(betashift,beta,
     1   iax,ilog2,CM_downward,32)
      call CM_get_from_power_two_always_1L(gammashift,gamma,
     1   iax,ilog2,CM_downward,32)
      call CM_get_from_power_two_always_1L(bshift,b,
     1   iax,ilog2,CM_downward,32)
      call CM_load_context(mask)
      call CM_invert_context
      call CM_f_divide_3_1L(temp1,bshift,alphashift,23,8)
      call CM_f_divide_3_1L(temp2,betashift,gammashift,23,8)
      call CM_f_multiply_2_1L(temp2,temp1,23,8)
c
      call CM_f_divide_3_1L(beta1,bshift,beta,23,8)
      call CM_f_divide_3_1L(gamma1,bshift,gamma,23,8)
      call CM_f_add_mult_1L(alpha1,alphashift,beta1,alpha,23,8)
      call CM_f_add_3_1L(alpha,alphashift,gamma1,23,8)
      call CM_f_divide_2_1L(alpha1,alpha,23,8)
      call CM_f_add_mult_1L(beta1,beta,temp1,betashift,23,8)
      call CM_f_add_2_1L(beta,temp2,23,8)
      call CM_f_divide_2_1L(beta1,beta,23,8)
      call CM_f_add_mult_1L(gamma1,gamma,temp1,betashift,23,8)
      call CM_f_add_2_1L(gamma,temp2,23,8)
      call CM_f_divide_2_1L(gamma1,gamma,23,8)
      call CM_f_move_1L(alpha,alpha1,23,8)
      call CM_f_move_1L(beta,beta1,23,8)
      call CM_f_move_1L(gamma,gamma1,23,8)
c
      call CM_load_context(mask)
      call CM_get_from_power_two_always_1L(mask,mask,iax,ilog2,
     1   CM_downward,1)
      call CM_logior_context(mask)
      call CM_store_context(mask)
      goto 10
 20   continue
c note context is always 1 on exit from the above loop
      call CM_f_move_always_1L(fraction,alpha,23,8)
      call CM_get_from_news_always_1L(b,b,iax,CM_downward,32)
      call CM_deallocate_stack_through(ashift)
      continued_fraction=1
      return
      end
c
c  This subroutine factors a tridiagonal system of equations
c
      integer function tridiag_factor(a,b,c,size,iax)
      implicit none
c number of equations per system
      integer size
c news axis along which each system runs
      integer iax
c id's of three input/output diagonals
      integer a,b,c
c
c id's of local real arrays
      integer a1,b1,fraction
c id of local logical array
      integer mask1
      integer error,i
      integer continued_fraction
      include '/usr/include/cm/paris-configuration-fort.h'
c
      call CM_set_context
      a1=CM_allocate_stack_field(32)
      b1=CM_allocate_stack_field(32)
      fraction=CM_allocate_stack_field(32)
      mask1=CM_allocate_stack_field(1)
c
      call CM_my_news_coordinate_1L(a1,iax,32)
      call CM_u_eq_zero_1L(a1,32)
      call CM_store_test_always(mask1)
c
      call CM_f_news_mult_always_3_1L(b1,a,c,iax,CM_downward,23,8)
      call CM_f_negate_1_1L(b1,23,8)
      call CM_load_context(mask1)
      call CM_f_move_zero_1L(b1,23,8)
      call CM_set_context
      call CM_f_eq_zero_1L(b1,23,8)
      call CM_store_test_always(mask1)
c
      call CM_f_move_always_1L(a1,b,23,8)
c
      error=continued_fraction(a1,b1,fraction,mask1,size,iax)
      call CM_f_divinto_const_always_3_1L(b,fraction,1.0,23,8)
      call CM_f_news_mult_always_2_1L(a,b,iax,CM_downward,23,8)
      call CM_f_multiply_always_2_1L(c,b,23,8)
      call CM_deallocate_stack_through(a1)
c
      tridiag_factor=1
      return
      end
c
c  the subroutines for working with diagonally dominant systems
c
c This subroutine factors a tridiagonal system using a step
c to take care of diagonal dominance.  This algorithm will
c prevent the scans from going to zero
c
      subroutine tridiag_diag_dom_factor_step(a,b,c,
     1   fact1,ovfact1,fact2,ovfact2,
     2   mask1,mask2,mask3,mask4,
     3   onep,sizep,
     3   size,step,iax)
      implicit none
      include '/usr/include/cm/paris-configuration-fort.h'
      integer step, size,iax
c id's of real arrays
      integer a,b,c,fact1,ovfact1,fact2,ovfact2
c id's of logical arrays
      integer mask1,mask2,mask3,mask4,onep,sizep
c
c id's of local real arrays
      integer k,k1,a1,b1,c1,f
c id's of local integer arrays
      integer coord,coord2
c id's of local logical arrays
      integer segbit,tmask1,tmask2
c
      integer i, error
      integer tridiag_factor
c
      k=CM_allocate_stack_field(32)
      k1=CM_allocate_stack_field(32)
      a1=CM_allocate_stack_field(32)
      b1=CM_allocate_stack_field(32)
      c1=CM_allocate_stack_field(32)
      f=CM_allocate_stack_field(32)
      segbit=CM_allocate_stack_field(1)
      tmask1=CM_allocate_stack_field(1)
      tmask2=CM_allocate_stack_field(1)
      coord=k
      coord2=k1
      call CM_set_context
      call CM_f_eq_zero_1L(a,23,8)
      call CM_store_test_always(tmask1)
      call CM_f_eq_zero_1L(c,23,8)
      call CM_store_test_always(tmask2)
c
      call CM_get_from_news_always_1L(mask1,tmask1,iax,CM_upward,1)
      call CM_get_from_news_always_1L(mask3,tmask2,iax,CM_downward,1)
c
      call CM_enumerate_1L(coord,iax,32,CM_upward,CM_inclusive,
     1   CM_none,0)
      call CM_u_eq_constant_1L(coord,size,32)
      call CM_store_test_always(sizep)
      call CM_u_eq_constant_1L(coord,1,32)
      call CM_store_test_always(onep)
c
      call CM_s_subfrom_constant_3_1L(coord2,coord,size+1,32)
      call CM_s_mod_constant_2_1L(coord,step,32)
      call CM_s_mod_constant_2_1L(coord2,step,32)
c
      call CM_s_eq_zero_1L(coord,32)
      call CM_logior_test_always(mask1)
      call CM_logior_test_always(sizep)
      call CM_store_test_always(mask1)
c
      call CM_s_eq_constant_1L(coord,1,32)
      call CM_logior_test_always(tmask1)
      call CM_store_test_always(mask2)
c
      call CM_s_eq_zero_1L(coord2,32)
      call CM_logior_test_always(mask3)
      call CM_logior_test_always(onep)
      call CM_store_test_always(mask3)
c
      call CM_s_eq_constant_1L(coord2,1,32)
      call CM_logior_test_always(tmask2)
      call CM_store_test_always(mask4)
c
      error=tridiag_factor(a,b,c,size,iax)
c
      call CM_f_move_zero_always_1L(k,23,8)
      call CM_f_news_sub_always_2_1L(k,a,iax,CM_upward,23,8)
      call CM_load_context(mask1)
      call CM_f_move_constant_1L(k,1.0,23,8)
      call CM_f_move_always_1L(f,k,23,8)
      call CM_set_context
c
      call CM_scan_with_f_multiply_1L(fact1,f,iax,23,8,
     1   CM_downward,CM_inclusive,CM_start_bit,mask1)
      call CM_logior_2_1L(mask1,onep,1)
c
c  get the exra factor for the second factor
c
      call CM_f_move_zero_always_1L(k,23,8)
      call CM_f_news_sub_always_2_1L(k,c,iax,CM_downward,23,8)
      call CM_load_context(mask3)
      call CM_logior_context(onep)
      call CM_f_move_constant_1L(k,1.0,23,8)
      call CM_f_move_always_1L(f,k,23,8)
c
      call CM_set_context
      call CM_scan_with_f_multiply_1L(fact2,f,iax,23,8,
     1   CM_upward,CM_inclusive,CM_start_bit,mask3)
      call CM_logior_2_1L(mask3,sizep,1)
c
c
      call CM_f_negate_1_1L(a,23,8)
      call CM_load_context(mask1)
      call CM_logior_context(onep)
      call CM_f_move_constant_1L(a,1.0,23,8)
c
      call CM_set_context
      call CM_f_negate_1_1L(c,23,8)
      call CM_scan_with_f_multiply_1L(k,a,iax,23,8,
     1   CM_upward,CM_inclusive,CM_start_bit,mask1)
      call CM_f_move_always_1L(a,k,23,8)
      call CM_load_context(mask1)
      call CM_f_move_zero_1L(a,23,8)
c
c going up second factor
c
      call CM_load_context(sizep)
      call CM_logior_context(mask3)
      call CM_f_move_constant_1L(c,1.0,23,8)
c
      call CM_set_context
      call CM_scan_with_f_multiply_1L(k,c,iax,23,8,
     1   CM_downward,CM_inclusive,CM_start_bit,mask3)
c
      call CM_f_move_always_1L(c,k,23,8)
      call CM_load_context(mask3)
      call CM_logior_context(sizep)
      call CM_f_move_zero_1L(c,23,8)
c
c store reciprocal of factors
c
      call CM_set_context
      call CM_f_divinto_const_always_3_1L(ovfact1,fact1,1.0,23,8)
      call CM_f_divinto_const_always_3_1L(ovfact2,fact2,1.0,23,8)
      call CM_deallocate_stack_through(k)
      return
      end
      subroutine tridiag_diag_dom_apply_step(a,b,c,k,
     1   fact1,ovfact1,fact2,ovfact2,
     2   mask1,mask2,mask3,mask4,
     3   onep,sizep,
     4   size,step, iax)
c
      implicit none
      include '/usr/include/cm/paris-configuration-fort.h'
      integer size,step
c axis along which to apply factors
      integer iax
c id's of real arrays
      integer a,b,c,k
c id's of real arrays
      integer fact1, ovfact1, fact2, ovfact2
c id's of logical arrays
      integer mask1,mask2,mask3,mask4
c id's of logical arrays
      integer onep,sizep
c
c id of local real array
      integer d,e
      integer i
c
      call CM_set_context
      d=CM_allocate_stack_field(32)
      e=CM_allocate_stack_field(32)
c
      call CM_f_multiply_always_3_1L(d,fact1,k,23,8)
      call CM_scan_with_f_add_1L(e,d,iax,23,8,CM_upward,
     1   CM_inclusive,CM_start_bit,mask2)
c
      call CM_load_context(mask1)
      call CM_f_move_1L(k,e,23,8)
      call CM_load_context(onep)
      call CM_f_move_zero_1L(k,23,8)
c
      call CM_set_context
      call CM_f_multiply_always_3_1L(d,e,ovfact1,23,8)
      call CM_scan_with_copy_1L(e,k,iax,32,CM_upward,
     1   CM_inclusive,CM_start_bit,mask1)
c
      call CM_f_mult_add_1L(k,a,e,d,23,8)
      call CM_f_multiply_always_2_1L(k,b,23,8)
      call CM_f_multiply_always_3_1L(d,k,fact2,23,8)
c
      call CM_scan_with_f_add_1L(e,d,iax,23,8,CM_downward,
     1   CM_inclusive,CM_start_bit,mask4)
c
      call CM_load_context(mask3)
      call CM_f_move_1L(k,e,23,8)
      call CM_load_context(sizep)
      call CM_f_move_zero_1L(k,23,8)
      call CM_set_context
c
      call CM_f_multiply_always_3_1L(d,e,ovfact2,23,8)
      call CM_scan_with_copy_1L(e,k,iax,32,CM_downward,
     1   CM_inclusive,CM_start_bit,mask3)
      call CM_f_mult_add_1L(k,c,e,d,23,8)
c
      call CM_deallocate_stack_through(d)
      return
      end
