face_splitting.F90

Go to the documentation of this file.

!! Copyright (C) 2024 A. Buccheri
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
 
#include "global.h"
 
module face_splitting_oct_m
  use, intrinsic :: iso_fortran_env
 
  use debug_oct_m
  use global_oct_m
  use io_function_oct_m
  use io_oct_m,              only: io_open, io_close
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
  private
  public :: column_wise_khatri_rao_product
 
contains
 
  subroutine column_wise_khatri_rao_product(y, x, z)
    real(real64), intent(in ), contiguous :: y(:, :)
    real(real64), intent(in ), contiguous :: x(:, :)
    real(real64), intent(out), contiguous :: z(:, :)
 
    real(real64), parameter :: alpha = 1.0_real64  
    integer                 :: nrows_x, nrows_y, ncols, k, ij
 
    nrows_x = size(x, 1)
    nrows_y = size(y, 1)
    ncols = size(x, 2)
 
    ! x and y arrays should have the same number of columns
    assert(size(y, 2) == ncols)
    ! The number of rows of z should equal the size(x, 1) * size(y, 1)
    assert(size(z, 1) == nrows_x * nrows_y)
    ! z should have the same number of columns as both x and y arrays
    assert(size(z, 2) == ncols)
 
    ! dger performs A := alpha*x*y**T + A, requiring A to be zeroed
    !$omp parallel do collapse(2) default(shared)
    do k = 1, ncols
      do ij = 1, nrows_x * nrows_y
        z(ij, k) = 0.0_real64
      enddo
    enddo
    !$omp end parallel do
 
    do k = 1, ncols
      ! Abuse the fact that I can pass a 1D array, z(:, k), of size (nrows_x * nrows_y) to
      ! lapack''s interface as long as it is contiguous in memory
      call dger(nrows_x, nrows_y, alpha, x(:, k), 1, y(:, k), 1, z(:, k), nrows_y)
    enddo
 
  end subroutine column_wise_khatri_rao_product
 
end module face_splitting_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: