get_sparsity_pattern Subroutine

    subroutine get_sparsity_pattern(me,irow,icol,&
                                    linear_irow,linear_icol,linear_vals,&
                                    maxgrp,ngrp)

    implicit none

    class(mission_type),intent(inout) :: me
    integer,dimension(:),allocatable,intent(out)  :: irow  !! sparsity pattern nonzero elements row indices
    integer,dimension(:),allocatable,intent(out)  :: icol  !! sparsity pattern nonzero elements column indices
    integer,dimension(:),allocatable,intent(out),optional  :: linear_irow !! linear sparsity pattern
                                                                          !! nonzero elements row indices
    integer,dimension(:),allocatable,intent(out),optional  :: linear_icol !! linear sparsity pattern nonzero
                                                                          !! elements column indices
    real(wp),dimension(:),allocatable,intent(out),optional :: linear_vals !! linear sparsity values (constant
                                                                          !! elements of the Jacobian)
    integer,intent(out),optional                           :: maxgrp      !! DSM sparsity partition
    integer,dimension(:),allocatable,intent(out),optional  :: ngrp        !! DSM sparsity partition

    integer :: k,ii,jj,icol_start,irow_start
    integer :: n_nonzero  !! number of nonzero elements in the jacobian
    integer :: n  !! number of optimization variables
    integer :: m  !! number of functions
    integer :: iblock
    integer,dimension(:),allocatable :: cols_to_remove !! the columns from the full pattern to remove
    logical,dimension(:),allocatable :: mask
    integer,dimension(:),allocatable :: icol_tmp

    ! get the size of the full problem, we will first construct the full pattern:
    call me%define_problem_size(n=n, m=m, n_nonzero=n_nonzero, full_problem=.true.)

    !-------------------------------------
    ! test... dense pattern
    !-------------------------------------
    ! allocate(irow(n*m))
    ! allocate(icol(n*m))
    ! k = 0
    ! do ii = 1, n
    !     do jj = 1, m
    !         k = k + 1
    !         icol(k) = ii
    !         irow(k) = jj
    !     end do
    ! end do
    ! return
    !-------------------------------------

    allocate(irow(n_nonzero))
    allocate(icol(n_nonzero))

    k = 0
    do iblock = 1, me%n_revs*4 ! block loop

        ! see Figure 22b in Modern Fortran paper
        icol_start = (iblock-1)*7 + 1  ! [1-14, 8-21, 15-29, ...]
        irow_start = (iblock-1)*6 + 1  ! [1-6,  7-12, 13-18, ...]

        do ii = icol_start,icol_start+13
            do jj = irow_start,irow_start+5
                k = k + 1
                irow(k) = jj
                icol(k) = ii
            end do
        end do

    end do

    ! the full pattern was generated above.
    ! Here we just remove the columns (optimization variables) we don't need.

    allocate(cols_to_remove(0))
    if (me%fix_initial_time) &
        cols_to_remove = [cols_to_remove, 1]
    if (me%fix_initial_r) &
        cols_to_remove = [cols_to_remove, [2,3,4]]

    if (me%fix_ry_at_end_of_rev > 0) then
        ! remove Ry at end of specified rev
        ! e.g. rev 1 : SEG8 Ry (km) : column 31
        !
        ! patch point: [t rx ry rz vx vy vz]
        !                    ^^
        cols_to_remove = [cols_to_remove, me%fix_ry_at_end_of_rev*7*5 - 4]
    end if

    if (me%fix_final_ry_and_vx) &
        cols_to_remove = [cols_to_remove, [n-4,n-2]] ! last segment Ry & Vx

    if (size(cols_to_remove) > 0) then

        allocate(mask(size(icol)))
        mask = .true. ! keep these
        do k = 1, size(cols_to_remove)
            where(icol == cols_to_remove(k)) mask = .false. ! remove these
        end do

        ! remove the cols:
        irow = pack(irow, mask)
        icol = pack(icol, mask)

        ! decrement ones > the ones removed
        icol_tmp = icol ! so we have the original indices
        do k = 1, size(cols_to_remove)
            where (icol_tmp > cols_to_remove(k)) icol = icol - 1
        end do

        call me%define_problem_size(n_nonzero=n_nonzero)

    end if

    end subroutine get_sparsity_pattern