module chain1_mod {
  logical, parameter, private :: T = .true., F = .false.
  interface chain0
     module procedure chain20
     module procedure chain30
  end interface
contains
	
  subroutine row0(op1,op2, adj,add, m1, m2, d, eps) { #        ROW  d = Am1+epsBm2
    interface {
       integer function op1(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
       integer function op2(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
    }
    logical, intent (in) :: adj, add
    real, dimension (:)  :: m1,m2,d
    real                 :: eps
    integer              :: st

    if (adj) { st = op1 (T, add, m1, d)        # m1 = A'd
               st = op2 (T, add, m2, d)        # m2 = B'd
               m2 = eps*m2                     # m2 = eps(B'd)
    }
    else {     st = op2 (F, add, eps*m2, d)    # d = epsBm2
               st = op1 (F,   T,     m1, d)    # d = Am1+epsBm2
    }
  }	

 
  subroutine chain20(op1,op2, adj,add, m,d,t1) { #        CHAIN 2
    interface {
       integer function op1(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
       integer function op2(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
    }
    logical, intent(in) :: adj, add
    real, dimension(:)  :: m,d, t1
    integer             :: st
    if(adj) {  st = op1( T,   F, t1,d )		#        t =    A' d
               st = op2( T, add, m, t1)		# m = B' t = B' A' d
    } 
    else    {  st = op2( F,   F, m, t1)		#        t =    B  m
               st = op1( F, add, t1,d )		# d = A  t = A  B  m
    }
  }

  subroutine chain30(op1,op2,op3, adj,add, m,d,t1,t2) { # CHAIN 3
    interface {
       integer function op1(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
       integer function op2(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
       integer function op3(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
    }
    logical, intent(in) :: adj, add
    real, dimension(:)  :: m,d, t1,t2
    integer             :: st
    if(adj) {  st = op1( T,   F, t2, d )     #                 t1 =       A' d
               st = op2( T,   F, t1, t2)     #         t2 = B' t1 =    B' A' d
               st = op3( T, add, m , t1)     # m  = C' t2 =         C' B' A' d
    } 
    else    {  st = op3( F,   F, m , t1)     #                 t1 =       C  m
               st = op2( F,   F, t1, t2)     #         t2 = B  t1 =    B  C  m
               st = op1( F, add, t2, d )     # d  = A  t2 =         A  B  C  m
    }
  }

}