csr_transpose Subroutine

   subroutine csr_transpose(self, at, err)
      !! Build the transpose of this matrix into `at`.
      !!
      !! `at` is an `n_cols` by `n_rows` matrix whose rows are, as always,
      !! sorted ascending by column index. Transposing twice reproduces the
      !! original structure and values exactly. The source is checked with
      !! `is_valid` first; a malformed source reports `ERROR_VALIDATION` and
      !! leaves `at` empty.
      class(csr_t), intent(in) :: self
      type(csr_t), intent(out) :: at
         !! Receives the transpose
      type(error_t), intent(out), optional :: err

      integer(default_int) :: i, k, j, n_entries
      integer(default_int), allocatable :: next_free(:)

      call at%destroy()
      if (.not. self%is_valid(err)) return

      n_entries = self%nnz()
      at%nrow = self%ncol
      at%ncol = self%nrow

      allocate (at%row_ptr(at%nrow + 1))
      allocate (at%col_idx(n_entries))
      allocate (at%values(n_entries))
      allocate (next_free(at%nrow))

      at%row_ptr = 0
      do k = 1, n_entries
         j = self%col_idx(k) + 1
         at%row_ptr(j) = at%row_ptr(j) + 1
      end do
      at%row_ptr(1) = 1
      do i = 1, at%nrow
         at%row_ptr(i + 1) = at%row_ptr(i + 1) + at%row_ptr(i)
         next_free(i) = at%row_ptr(i)
      end do

      ! walking the source in row order leaves every target row ascending
      do i = 1, self%nrow
         do k = self%row_ptr(i), self%row_ptr(i + 1) - 1
            j = self%col_idx(k)
            at%col_idx(next_free(j)) = i
            at%values(next_free(j)) = self%values(k)
            next_free(j) = next_free(j) + 1
         end do
      end do
   end subroutine csr_transpose