I am trying to define the (+) operator between Fortran derived types that describe matrices (linear operators). My goal is to implicitly define a matrix M = M1 + M2 + M3 such that, given a vector x, Mx = M1x + M2x + M3x.
First, I defined an abstract type (abs_linop) with the abstract interface for a matrix vector multiplication (y = M *x). Then, I built an derived type (add_linop), extending the abstract type (abs_linop). The operator (+) is defined for the type (add_linop). I then create an example of concrete type (eye) extending the abstract type (abs_linop) that describes the identity matrix. This type is used in the main program. This is the source code
module LinearOperator
implicit none
private
public :: abs_linop,multiplication
type, abstract :: abs_linop
integer :: nrow=0
integer :: ncol=0
character(len=20) :: name='empty'
contains
!> Procedure for computation of (matrix) times (vector)
procedure(multiplication), deferred :: Mxv
end type abs_linop
abstract interface
!>-------------------------------------------------------------
!> Abstract procedure defining the interface for a general
!<-------------------------------------------------------------
subroutine multiplication(this,vec_in,vec_out,info,lun_err)
import abs_linop
implicit none
class(abs_linop), intent(inout) :: this
real(kind=8), intent(in ) :: vec_in(this%ncol)
real(kind=8), intent(inout) :: vec_out(this%nrow)
integer, optional, intent(inout) :: info
integer, optional, intent(in ) :: lun_err
end subroutine multiplication
end interface
!>---------------------------------------------------------
!> Structure variable for Identity matrix
!> (rectangular case included)
!>---------------------------------------------------------
type, extends(abs_linop), public :: eye
contains
!> Static constructor
procedure, public, pass :: init => init_eye
!> Compute matrix times vector operatoration
procedure, public, pass :: Mxv => apply_eye
end type eye
!>----------------------------------------------------------------
!> Structure variable to build implicit matrix defined
!> as composition and sum of linear operator
!>----------------------------------------------------------------
public :: add_linop, operator(+)
type, extends(abs_linop) :: add_linop
class(abs_linop) , pointer :: matrix_1
class(abs_linop) , pointer :: matrix_2
real(kind=8), allocatable :: scr(:)
contains
procedure, public , pass:: Mxv => add_Mxv
end type add_linop
INTERFACE OPERATOR (+)
module PROCEDURE mmsum
END INTERFACE OPERATOR (+)
contains
!>------------------------------------------------------
!> Function that give two linear operator A1 and A2
!> defines, implicitely, the linear operator
!> A=A1+A2
!> (public procedure for class add_linop)
!>
!> usage:
!> 'var' = A1 + A2
!<-------------------------------------------------------------
function mmsum(matrix_1,matrix_2) result(this)
implicit none
class(abs_linop), target, intent(in) :: matrix_1
class(abs_linop), target, intent(in) :: matrix_2
type(add_linop) :: this
! local
integer :: res
character(len=20) :: n1,n2
if (matrix_1%nrow .ne. matrix_2%nrow) &
write(*,*) 'Error mmproc dimension must agree '
if (matrix_1%ncol .ne. matrix_2%ncol) &
write(*,*) 'Error mmproc dimension must agree '
this%matrix_1 => matrix_1
this%matrix_2 => matrix_2
this%nrow = matrix_1%nrow
this%ncol = matrix_2%ncol
this%name=etb(matrix_1%name)//'+'//etb(matrix_2%name)
write(*,*) 'Sum Matrix initialization '
write(*,*) 'M1 : ',this%matrix_1%name
write(*,*) 'M2 : ',this%matrix_2%name
write(*,*) 'sum : ',this%name
allocate(this%scr(this%nrow),stat=res)
contains
function etb(strIn) result(strOut)
implicit none
! vars
character(len=*), intent(in) :: strIn
character(len=len_trim(adjustl(strIn))) :: strOut
strOut=trim(adjustl(strIn))
end function etb
end function mmsum
recursive subroutine add_Mxv(this,vec_in,vec_out,info,lun_err)
implicit none
class(add_linop), intent(inout) :: this
real(kind=8), intent(in ) :: vec_in(this%ncol)
real(kind=8), intent(inout) :: vec_out(this%nrow)
integer, optional, intent(inout) :: info
integer, optional, intent(in ) :: lun_err
write(*,*) 'Matrix vector multipliction',&
'matrix:',this%name,&
'M1: ',this%matrix_1%name,&
'M2: ',this%matrix_2%name
select type (mat=>this%matrix_1)
type is (add_linop)
write(*,*) 'is allocated(mat%scr) ?', allocated(mat%scr)
end select
call this%matrix_1%Mxv(vec_in,this%scr,info=info,lun_err=lun_err)
call this%matrix_2%Mxv(vec_in,vec_out,info=info,lun_err=lun_err)
vec_out = this%scr + vec_out
end subroutine add_Mxv
subroutine init_eye(this,nrow)
implicit none
class(eye), intent(inout) :: this
integer, intent(in ) :: nrow
this%nrow = nrow
this%ncol = nrow
end subroutine init_eye
subroutine apply_eye(this,vec_in,vec_out,info,lun_err)
class(eye), intent(inout) :: this
real(kind=8), intent(in ) :: vec_in(this%ncol)
real(kind=8), intent(inout) :: vec_out(this%nrow)
integer, optional, intent(inout) :: info
integer, optional, intent(in ) :: lun_err
! local
integer :: mindim
vec_out = vec_in
if (present(info)) info=0
end subroutine apply_eye
end module LinearOperator
program main
use LinearOperator
implicit none
real(kind=8) :: x(2),y(2),z(2),t(2)
type(eye) :: id1,id2,id3
type(add_linop) :: sum12,sum23,sum123_ok,sum123_ko
integer :: i
call id1%init(2)
id1%name='I1'
call id2%init(2)
id2%name='I2'
call id3%init(2)
id3%name='I3'
x=1.0d0
y=1.0d0
z=1.0d0
write(*,*) ' Vector x =', x
call id1%Mxv(x,t)
write(*,*) ' Vector t = I1 *x', t
write(*,*) ' '
sum12 = id1 + id2
call sum12%Mxv(x,t)
write(*,*) ' Vector t = (I1 +I2) *x', t
write(*,*) ' '
sum23 = id2 + id3
sum123_ok = id1 + sum23
call sum123_ok%Mxv(x,t)
write(*,*) ' Vector t = ( I1 + (I2 + I3) )*x', t
write(*,*) ' '
sum123_ko = id1 + id2 + id3
call sum123_ko%Mxv(x,t)
write(*,*) ' Vector t = ( I1 +I2 + I3) *x', t
end program main
I compile this code with gfortran version 7.5.0 and flags "-g -C -Wall -fcheck=all -O -ffree-line-length-none -mcmodel=medium " and this is what I get
Vector x = 1.0000000000000000 1.0000000000000000 Vector t = I1 *x 1.0000000000000000 1.0000000000000000 Sum Matrix initialization M1 : I1 M2 : I2 sum : I1+I2 Matrix vector multiplictionmatrix:I1+I2 M1: I1 M2: I2 Vector t = (I1 +I2) *x 2.0000000000000000 2.0000000000000000 Sum Matrix initialization M1 : I2 M2 : I3 sum : I2+I3 Sum Matrix initialization M1 : I1 M2 : I2+I3 sum : I1+I2+I3 Matrix vector multiplictionmatrix:I1+I2+I3 M1: I1 M2: I2+I3 Matrix vector multiplictionmatrix:I2+I3 M1: I2 M2: I3 Vector t = ( I1 + (I2 + I3) )*x 3.0000000000000000 3.0000000000000000 Sum Matrix initialization M1 : I1 M2 : I2 sum : I1+I2 Sum Matrix initialization M1 : I1+I2 M2 : I3 sum : I1+I2+I3 Matrix vector multiplictionmatrix:I1+I2+I3 M1: I1+I2 M2: I3 is allocated(mat%scr) ? F Matrix vector multiplictionmatrix:I1+I2 M1: I1 M2: I2 At line 126 of file LinearOperator.f90 Fortran runtime error: Allocatable actual argument 'this' is not allocated
Everthing works fine when I use the (+) operator with 2 terms. But when 3 terms are used there is an issue with the allocatable array scr, member of type (add_linop), that is not allocated.
Does anybody knows the reason of this issue and how to solve it? I include the Makefile used for compiling the code.
#Gfortran compiler
FC = gfortran
OPENMP = -fopenmp
MODEL = -mcmodel=medium
OFLAGS = -O5 -ffree-line-length-none
DFLAGS = -g -C -Wall -fcheck=all -O -ffree-line-length-none
#DFLAGS = -g -C -Wall -ffree-line-length-none -fcheck=all
PFLAGS = -pg
CPPFLAGS = -D_GFORTRAN_COMP
ARFLAGS =
ODIR = objs
MDIR = mods
LDIR = libs
INCLUDE = -J$(MODDIR)
OBJDIR = $(CURDIR)/$(ODIR)
MODDIR = $(CURDIR)/$(MDIR)
LIBDIR = $(CURDIR)/$(LDIR)
INCLUDE += -I$(MODDIR)
FFLAGS = $(OFLAGS) $(MODEL) $(INCLUDE)
LIBSRCS =
DEST = .
EXTHDRS =
HDRS =
LIBS = -llapack -lblas
LIBMODS =
LDFLAGS = $(MODEL) $(INCLUDE) -L. -L/usr/lib -L/usr/local/lib -L$(LIBDIR)
LINKER = $(FC)
MAKEFILE = Makefile
PRINT = pr
CAT = cat
PROGRAM = main.out
SRCS = LinearOperator.f90
OBJS = LinearOperator.f90
PRJS= $(SRCS:jo=.prj)
OBJECTS = $(SRCS:%.f90=$(OBJDIR)/%.o)
MODULES = $(addprefix $(MODDIR)/,$(MODS))
.SUFFIXES: .prj .f90
print-% :
@echo $* = $($*)
.f.prj:
ftnchek -project -declare -noverbose $<
.f90.o:
$(FC) $(FFLAGS) $(INCLUDE) -c $<
all::
@make dirs
@make $(PROGRAM)
$(PROGRAM): $(LIBS) $(MODULES) $(OBJECTS)
$(LINKER) -o $(PROGRAM) $(LDFLAGS) $(OBJECTS) $(LIBS)
$(LIBS):
@set -e; for i in $(LIBSRCS); do cd $$i; $(MAKE) --no-print-directory -e CURDIR=$(CURDIR); cd $(CURDIR); done
$(OBJECTS): $(OBJDIR)/%.o: %.f90
$(FC) $(CPPFLAGS) $(FFLAGS) -o $@ -c $<
dirs:
@-mkdir -p $(OBJDIR) $(MODDIR) $(LIBDIR)
clean-emacs:
@-rm -f $(CURDIR)/*.*~
@-rm -f $(CURDIR)/*\#*
check: $(PRJS)
ftnchek -noverbose -declare $(PRJS) -project -noextern -library > $(PROGRAM).ftn
profile:; @make "FFLAGS=$(PFLAGS) $(MODEL) " "CFLAGS=$(PFLAGS) $(MODEL)" "LDFLAGS=$(PFLAGS) $(LDFLAGS)" $(PROGRAM)
debug:; @make "FFLAGS=$(DFLAGS) $(MODEL) $(INCLUDE)" "LDFLAGS=$(DFLAGS) $(LDFLAGS)" $(PROGRAM)
openmp:; @make "FFLAGS=$(OFLAGS) $(OPENMP) $(MODEL) $(INCLUDE)" "LDFLAGS=$(LDFLAGS) $(OPENMP)" $(PROGRAM)
clean:; @rm -f $(OBJECTS) $(MODULES) $(PROGRAM).cat $(PROGRAM).ftn
@set -e; for i in $(LIBSRCS); do cd $$i; $(MAKE) --no-print-directory clean; cd $(CURDIR); done
clobber:; @rm -f $(OBJECTS) $(MODULES) $(PROGRAM).cat $(PROGRAM).ftn $(PROGRAM)
@-rm -rf $(OBJDIR) $(MODDIR) $(LIBDIR)
@-rm -f $(CURDIR)/*.*~
@-rm -f $(CURDIR)/*\#*
.PHONY: mods
index:; ctags -wx $(HDRS) $(SRCS)
install: $(PROGRAM)
install -s $(PROGRAM) $(DEST)
print:; $(PRINT) $(HDRS) $(SRCS)
cat:; $(CAT) $(HDRS) $(SRCS) > $(PROGRAM).cat
program: $(PROGRAM)
profile: $(PROFILE)
tags: $(HDRS) $(SRCS); ctags $(HDRS) $(SRCS)
update: $(DEST)/$(PROGRAM)
main.o: linearoperator.mod
# DO NOT EDIT --- auto-generated file
linearoperator.mod : LinearOperator.f90
$(FC) $(FCFLAGS) -c $<
