columns_in_partition_group Subroutine

    subroutine columns_in_partition_group(me,igroup,n_cols,cols,nonzero_rows,indices,status_ok)

    implicit none

    class(sparsity_pattern),intent(in)           :: me
    integer,intent(in)                           :: igroup       !! group number. Should be `>0` and `<=me%mxgrp`
    integer,intent(out)                          :: n_cols       !! number of columns in the `igroup` group.
    integer,dimension(:),allocatable,intent(out) :: cols         !! the column numbers in the `igroup` group.
                                                                 !! (if none, then it is not allocated)
    integer,dimension(:),allocatable,intent(out) :: nonzero_rows !! the row numbers of all the nonzero
                                                                 !! Jacobian elements in this group
    integer,dimension(:),allocatable,intent(out) :: indices      !! nonzero indices in `jac` for a group
    logical,intent(out)                          :: status_ok    !! true if the partition is valid

    integer :: i  !! counter
    integer :: num_nonzero_elements_in_col          !! number of nonzero elements in a column
    integer :: num_nonzero_elements_in_group        !! number of nonzero elements in a group
    integer,dimension(:),allocatable :: col_indices !! nonzero indices in `jac` for a column

    if (me%maxgrp>0 .and. allocated(me%ngrp)) then

        status_ok = .true.

        n_cols = count(me%ngrp==igroup)
        if (n_cols>0) then
            allocate(cols(n_cols))
            cols = pack([(i,i=1,size(me%ngrp))],mask=me%ngrp==igroup)
        end if

        ! get all the non-zero elements in each column:
        num_nonzero_elements_in_group = 0  ! initialize
        do i = 1, n_cols
            num_nonzero_elements_in_col = count(me%icol==cols(i))
            if (num_nonzero_elements_in_col/=0) then ! there are functions to
                                                     ! compute in this column
                num_nonzero_elements_in_group = num_nonzero_elements_in_group + &
                                                num_nonzero_elements_in_col
                if (allocated(col_indices)) deallocate(col_indices)
                allocate(col_indices(num_nonzero_elements_in_col))
                !col_indices = pack(me%indices,mask=me%icol==cols(i))
                block
                    integer :: j,n
                    n = 0
                    do j = 1, size(me%icol)
                        if (me%icol(j)==cols(i)) then
                            n = n + 1
                            col_indices(n) = j
                        end if
                    end do
                end block
                if (allocated(nonzero_rows)) then
                    nonzero_rows = [nonzero_rows,me%irow(col_indices)]
                    indices = [indices,col_indices]
                else
                    nonzero_rows = me%irow(col_indices)
                    indices = col_indices
                end if
            end if
        end do

    else
        status_ok = .false.
    end if

    end subroutine columns_in_partition_group