Функция Fortran для перегрузки умножения между производными типами с выделяемыми компонентами

Предисловие

Для хранения полосатых матриц, полные аналоги которых могут иметь как строки, так и столбцы, проиндексированные из индексов, отличных от 1, я определил производный тип данных как

TYPE CDS
  REAL, DIMENSION(:,:), ALLOCATABLE :: matrix
  INTEGER, DIMENSION(2) :: lb, ub
  INTEGER :: ld, ud
END TYPE CDS

где CDS означает сжатое диагональное хранилище. Учитывая объявление TYPE(CDS) :: A,

1. Компонент ранга 2 matrix должен содержать в качестве столбцов диагонали фактической полной матрицы (например, здесь, за исключением того, что я сохраняю диагонали как столбцы, а не как строки).
2. Предполагается, что компоненты ld и ud содержат количество нижних и верхних диагоналей соответственно, то есть -lbound(A%matrix,2) и +ubound(A%matrix,2).
3. Предполагается, что двухэлементные компоненты lb и ub содержат нижние и верхние границы фактической полной матрицы по двум измерениям. В частности, lb(1) и ub(1) должны совпадать с lbound(A%matrix,1) и lbound(A%matrix,2).

Как вы можете видеть в пунктах 2. и 3., производный тип содержит некоторую избыточную информацию, но меня это не волнует, так как это всего лишь 3 пары целых чисел. Кроме того, в коде, который я пишу, информация о границах и диапазоне фактической полной матрицы известна до заполнения матрицы. Поэтому я сначала присвоил значения компонентам ld, ud, lb и ub, а затем использовал эти компоненты для ALLOCATE компонента matrix (тогда я смогу правильно его заполнить).

Проблема

Мне нужно выполнить умножение матриц между такими разреженными матрицами, поэтому я написал FUNCTION для выполнения такого произведения и использовал его для перегрузки оператора *.

На данный момент функция выглядит следующим образом:

FUNCTION CDS_mat_x_CDS_mat(A, B)
IMPLICIT NONE
TYPE(CDS), INTENT(IN) :: A, B
TYPE(CDS) :: cds_mat_x_cds_mat

! determine the lower and upper bounds and band of the result based on those of the operands
CDS_mat_x_CDS_mat%lb(1) = A%lb(1)
CDS_mat_x_CDS_mat%ub(1) = A%ub(1)
CDS_mat_x_CDS_mat%lb(2) = B%lb(2)
CDS_mat_x_CDS_mat%ub(2) = B%ub(2)
CDS_mat_x_CDS_mat%ld = A%ld + B%ld
CDS_mat_x_CDS_mat%ud = A%ud + B%ud

! allocate the matrix component
ALLOCATE(CDS_mat_x_CDS_mat%matrix(CDS_mat_x_CDS_mat%lb(1):CDS_mat_x_CDS_mat%ub(1),&
                               & -CDS_mat_x_CDS_mat%ld:+CDS_mat_x_CDS_mat%ud))

! perform the product
:
:

END FUNCTION

Это означает, что если мне нужно выполнить продукт несколько раз, распределение будет выполнено несколько раз внутри функции. Я думаю, что это не хорошо с точки зрения производительности.

Я прошу предложений о том, как выполнить задачу ленточной разреженной матрицы, умноженной на лентовидную разреженную матрицу. Я хотел бы использовать тип, который я определил, потому что мне нужно, чтобы он был таким же общим, с точки зрения границ, как на данный момент. Но я мог бы изменить процедуру выполнения продукта (с FUNCTION на SUBROUTINE, если нужно).

Идеи

Я мог бы переписать процедуру как SUBROUTINE, чтобы объявить CDS_mat_x_CDS_mat с помощью INTENT(INOUT), выполнить назначение компонентов, отличных от matrix, а также распределение вне SUBROUTINE. Недостатком было бы то, что я не мог перегрузить оператор *.

Я заметил, что встроенная функция MATMUL может работать с любыми операндами ранга 2, какими бы ни были верхняя и нижняя границы двух измерений. Это означает, что выделение выполняется внутри функции. Я полагаю, что это эффективно (поскольку это встроенная функция). Отличие моей функции в том, что она принимает массивы ранга 2 любой формы, тогда как моя принимает объекты производных типов данных с компонентом массива ранга 2 любой формы.