scale_mat Subroutine

   pure subroutine scale_mat(n, m, wt, ldwt, ld2wt, t, wtt)
   !! Scale matrix `t` using `wt`, i.e., compute `wtt = wt*t`.

      use odrpack_kinds, only: zero

      integer, intent(in) :: n
         !! Number of rows of data in `t`.
      integer, intent(in) :: m
         !! Number of columns of data in `t`.
      real(wp), intent(in), target :: wt(..)
         !! Array of shape conformable to `(ldwt,ld2wt,m)` holding the weights.
      integer, intent(in) :: ldwt
         !! Leading dimension of array `wt`.
      integer, intent(in) :: ld2wt
         !! Second dimension of array `wt`.
      real(wp), intent(in), target :: t(..)
         !! Array of shape conformable to `(n,m)` being scaled by `wt`.
      real(wp), intent(out), target :: wtt(..)
         !! Array of shape conformable to `(n,m)` holding the result of weighting array `t` by
         !! array `wt`. Array `wtt` can be the same as `t` only if the arrays in `wt` are upper
         !! triangular with zeros below the diagonal.

      ! Local scalars
      integer :: i, j
      real(wp), pointer :: wt_(:, :, :), t_(:, :), wtt_(:, :)

      ! Variable Definitions (alphabetically)
      !  I:       An indexing variable.
      !  J:       An indexing variable.

      if (n == 0 .or. m == 0) return

      select rank (wt)
      rank (1); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
      rank (2); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
      rank (3); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
      rank default; error stop "Invalid rank of `wt`."
      end select

      select rank (t)
      rank (1); t_(1:n, 1:m) => t
      rank (2); t_(1:n, 1:m) => t
      rank default; error stop "Invalid rank of `t`."
      end select

      select rank (wtt)
      rank (1); wtt_(1:n, 1:m) => wtt
      rank (2); wtt_(1:n, 1:m) => wtt
      rank default; error stop "Invalid rank of `wtt`."
      end select

      if (wt_(1, 1, 1) >= zero) then
         if (ldwt >= n) then
            if (ld2wt >= m) then
               ! WT is an N-array of M by M matrices
               do j = 1, m
                  do i = 1, n
                     wtt_(i, j) = dot_product(wt_(i, j, :), t_(i, :))
                  end do
               end do
            else
               ! WT is an N-array of diagonal matrices
               do j = 1, m
                  wtt_(:, j) = wt_(:, 1, j)*t_(:, j)
               end do
            end if
         else
            if (ld2wt >= m) then
               ! WT is an M by M matrix
               do j = 1, m
                  do i = 1, n
                     wtt_(i, j) = dot_product(wt_(1, j, :), t_(i, :))
                  end do
               end do
            else
               ! WT is a diagonal matrix
               do j = 1, m
                  wtt_(:, j) = wt_(1, 1, j)*t_(:, j)
               end do
            end if
         end if
      else
         ! WT is a scalar
         wtt_ = abs(wt_(1, 1, 1))*t_
      end if

   end subroutine scale_mat