Функция на 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. Предполага се, че 2-елементните компоненти 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 с всякаква форма.