derived data type with allocatable element

[tservas]

hi,

i would like to implement a derived data type for square matrices. the dimension of the square matrices are not known at compilation time ... so i used allocatable arrays and implemented some basic new and delete functions ... see below ... my problem is that i would like to implement an operator .x. corresponding to the direct matrix product (kronecker product), put my implementation for this function (called "direct_product" in the code below) does not work properly ... the result of the function should be of type square_matrix this means that i have to allocate space for it ... after the calculation is done it should give back the resulting square_matrix and after that free the allocated memory ... i dont know how i can achieve such a behavior, i.e. how i can free the memory after the function direct_product returned the result ... i get the following message if i run a program which uses this module:

Remaining memory: xxx bytes allocated at line xx of program.f95

MODULE square_matrices
!
! Purpose:
! Define a derived data type square_matrix and the operations that can
! performed on it. The module defines the operations:
!
! Operation Operator
! ==================== ===========================
! Memory allocation for a square matrix new
! Memory deallocation for a square matrix delete
! Creation of a square matrix from an array =
! Conversion of a square matrix to an array =
! Assigning one matrix to another =
! Matrix scalar product (4 cases) *
! Matrix scalar division (2 cases) /
! Addition of matrices +
! Subtraction of matrices -
! Matrix multiplication *
! Direct matrix product (tensor product) .x.
!
!
IMPLICIT NONE
PRIVATE
PUBLIC :: square_matrix, new, delete, dimension_of, unit_matrix, ASSIGNMENT(=), OPERATOR(.x.)

TYPE :: square_matrix
COMPLEX, ALLOCATABLE, DIMENSION, :: array
INTEGER :: dim
END TYPE


INTERFACE ASSIGNMENT(=)
MODULE PROCEDURE matrix_assignment
END INTERFACE

INTERFACE OPERATOR(.x.)
MODULE PROCEDURE direct_product
END INTERFACE

CONTAINS

SUBROUTINE matrix_assignment(sq_matrix_1, sq_matrix_2)
TYPE (square_matrix), INTENT(INOUT) :: sq_matrix_1
TYPE (square_matrix), INTENT(IN) :: sq_matrix_2
INTEGER :: i, j

IF (sq_matrix_1%dim == sq_matrix_2%dim) THEN
FORALL (i=1:sq_matrix_1%dim, j=1:sq_matrix_1%dim)
sq_matrix_1%array(i, j) = sq_matrix_2%array(i, j)
END FORALL
sq_matrix_1%dim = sq_matrix_2%dim
ELSE
WRITE(*,*) "Matrix_assignment: matrix dimensions do not match"
RETURN
ENDIF
END SUBROUTINE matrix_assignment

FUNCTION new(sq_matrix, matrix_size)
INTEGER :: new
TYPE (square_matrix), INTENT(INOUT) :: sq_matrix
INTEGER :: matrix_size
INTEGER :: alloc_status

ALLOCATE ( sq_matrix%array(matrix_size,matrix_size), STAT = alloc_status)
sq_matrix%array = 0
IF (alloc_status /= 0) THEN
WRITE (*,*) "Memory allocation error!"
RETURN
ELSE
new = alloc_status
sq_matrix%dim = matrix_size
ENDIF
END FUNCTION new

FUNCTION delete(sq_matrix)
INTEGER :: delete
TYPE (square_matrix), INTENT(INOUT) :: sq_matrix
INTEGER :: dealloc_status
DEALLOCATE(sq_matrix%array, STAT=dealloc_status)
IF (dealloc_status /= 0) WRITE(*,*) "Matrix deletion: matrix could not be deleted"
delete = dealloc_status
END FUNCTION delete

FUNCTION dimension_of (sq_matrix)
INTEGER :: dimension_of
TYPE (square_matrix), INTENT(IN) :: sq_matrix
dimension_of = sq_matrix%dim
END FUNCTION dimension_of

FUNCTION direct_product (sq_matrix_1, sq_matrix_2)
TYPE (square_matrix) :: direct_product
TYPE (square_matrix), INTENT(IN) :: sq_matrix_1, sq_matrix_2
INTEGER :: dim_direct_product, dim_1, dim_2
INTEGER :: alloc_status
INTEGER :: n, m, i, j, i_prime, j_prime, row_min, row_max, column_min, column_max

dim_1 = sq_matrix_1%dim
dim_2 = sq_matrix_2%dim

dim_direct_product = dim_1 * dim_2

alloc_status = new(direct_product, dim_direct_product)

! The following nested do loops are used to
! calculate the direct product of two square
! matrices. This calculation is based on the following
! formula: C(i,j) = A(n,m)*B(i_prime, j_prime), with
! (n-1)*num_rows_B + 1 <= i <= n*num_rows_B,
! (m-1)*num_columns_B + 1 <= j <= m*num_columns_B,
! i_prime = i - (n-1)*num_rows_B + 1,
! j_prime = j - (m-1)*num_columns_B + 1,
! The resulting matrix C is build up block by block, according to the
! definition of the direct matrix product, i.e.
! the (num_rows_B x num_columns_B) block A(1,1)*B gets in the upper left corner, etc.
! The two outer do loops ( iterating n, m) thus go from one block to the next one.
! The two inner do loops ( iterating i, j) thus operate only within one block
!
DO n=1,dim_1, 1
DO m=1, dim_1, 1
row_min = (n-1)*dim_2 + 1
row_max = n*dim_2
column_min = (m-1)*dim_2 + 1
column_max = m*dim_2
DO i=row_min, row_max, 1
DO j=column_min, column_max, 1
i_prime = i - row_min + 1
j_prime = j - column_min +1
direct_product%array(i,j) = sq_matrix_1%array(n, m) * sq_matrix_2%array(i_prime, j_prime)
ENDDO
ENDDO
ENDDO
ENDDO


!
!
! the above do loop construct is not equivalent to the following forall construct ... i dont know why ...
!
!
! FORALL ( n=1:dim_1, m=1:dim_1 )
! row_min = (n-1)*dim_2 + 1
! row_max = n*dim_2
! column_min = (m-1)*dim_2 + 1
! column_max = m*dim_2
! FORALL ( i=(n-1)*dim_2:n*dim_2, j=(m-1)*dim_2:m*dim_2 )
! i_prime = i - row_min + 1
! j_prime = j - column_min +1
! direct_product%array(i,j) = sq_matrix_1%array(n, m) * sq_matrix_2%array(i_prime, j_prime)
! END FORALL
! END FORALL

END FUNCTION direct_product

FUNCTION unit_matrix (dim_matrix)
TYPE (square_matrix) :: unit_matrix
INTEGER, INTENT(IN) :: dim_matrix
INTEGER :: alloc_status, i

alloc_status = new(unit_matrix, dim_matrix)

FORALL ( i=1:dim_matrix )
unit_matrix%array(i,i) = (1.0, 0.0)
END FORALL
END FUNCTION unit_matrix

END MODULE square_matrices

... thx for your help ... it would also be very helpful if someone could give me a reference to the source code of an already existing implementation of such a data type with the corresponding operations ...

tservas