Init_ColumnOperators3D_t Subroutine

  subroutine Init_ColumnOperators3D_t(this,interp,mesh,geometry)
    !! Binds the column operators to an extruded mesh and its geometry and, on
    !! Gauss-Lobatto nodes, precomputes the inverses of the two sweep matrices
    !!
    !!   A = D + M^{-1} e_0 e_0^T ,  At = D - M^{-1} e_N e_N^T .
    !!
    !! Both are (N+1) x (N+1). They are built from interp%dMatrix and inverted by
    !! Gauss-Jordan elimination with partial pivoting in double precision, then rounded
    !! to prec. The geometry must stay alive while the operators are used.
    !!
    !!  Input
    !!    - interp : the interpolant of the fields and of the geometry, any N >= 1
    !!    - mesh : a mesh built by Mesh3D%ExtrudedMesh (nLayers >= 1)
    !!    - geometry : SEMHex generated from mesh with interp
    implicit none
    class(ColumnOperators3D_t),intent(out) :: this
    type(Lagrange),target,intent(in) :: interp
    type(Mesh3D),intent(in) :: mesh
    type(SEMHex),target,intent(in) :: geometry
    ! Local
    integer :: i,m,np
    real(real64),allocatable :: a(:,:),aInv(:,:)

    if(mesh%nLayers <= 0) then
      error stop "ColumnOperators3D : the mesh is not extruded (nLayers = 0); build it with Mesh3D%ExtrudedMesh"
    endif
    if(mod(mesh%nElem,mesh%nLayers) /= 0) then
      error stop "ColumnOperators3D : the element count of the mesh is not a multiple of its layer count"
    endif
    if(geometry%nElem /= mesh%nElem) then
      error stop "ColumnOperators3D : the geometry does not have the element count of the mesh"
    endif
    if(geometry%dxds%N /= interp%N) then
      error stop "ColumnOperators3D : the geometry was not generated with this interpolant degree"
    endif

    this%interp => interp
    this%geometry => geometry
    this%N = interp%N
    this%nLayers = mesh%nLayers
    this%nElem = mesh%nElem
    this%nColumns = mesh%nElem/mesh%nLayers
    this%lobatto = interp%controlNodeType == GAUSS_LOBATTO

    if(this%lobatto) then
      np = interp%N+1
      allocate(this%upInverse(1:np,1:np),this%downInverse(1:np,1:np))
      allocate(a(1:np,1:np),aInv(1:np,1:np))

      ! A(i,m) = D(i,m) + delta_{i,0} delta_{m,0}/w_0, in the usual row, column order
      do m = 1,np
        do i = 1,np
          a(i,m) = real(interp%dMatrix(m,i),real64)
        enddo
      enddo
      a(1,1) = a(1,1)+1.0_real64/real(interp%qWeights(1),real64)
      call InvertMatrix(a,aInv,np)
      do m = 1,np
        do i = 1,np
          this%upInverse(m,i) = real(aInv(i,m),prec)
        enddo
      enddo

      ! At(i,m) = D(i,m) - delta_{i,N} delta_{m,N}/w_N
      do m = 1,np
        do i = 1,np
          a(i,m) = real(interp%dMatrix(m,i),real64)
        enddo
      enddo
      a(np,np) = a(np,np)-1.0_real64/real(interp%qWeights(np),real64)
      call InvertMatrix(a,aInv,np)
      do m = 1,np
        do i = 1,np
          this%downInverse(m,i) = real(aInv(i,m),prec)
        enddo
      enddo

      deallocate(a,aInv)
    endif

  endsubroutine Init_ColumnOperators3D_t