时间戳

首页 > 解决方案 > 使用派生类型在 Fortran 中定义具有多个加数的运算符（+）。可分配数组的问题

oop - 使用派生类型在 Fortran 中定义具有多个加数的运算符（+）。可分配数组的问题

问题描述

我正在尝试在描述矩阵（线性运算符）的 Fortran 派生类型之间定义 (+) 运算符。我的目标是隐式定义一个矩阵 M = M1 + M2 + M3，这样，给定一个向量 x，M x = M1 x + M2 x + M3 x。

首先，我用矩阵向量乘法 (y = M *x) 的抽象接口定义了一个抽象类型 (abs_linop)。然后，我构建了一个派生类型 (add_linop)，扩展了抽象类型 (abs_linop)。运算符 (+) 是为类型 (add_linop) 定义的。然后，我创建了一个具体类型（眼睛）的示例，它扩展了描述单位矩阵的抽象类型（abs_linop）。这种类型在主程序中使用。这是源代码

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

我用 gfortran 7.5.0 版和标志“-g -C -Wall -fcheck=all -O -ffree-line-length-none -mcmodel=medium”编译这段代码，这就是我得到的

  向量 x = 1.0000000000000000 1.0000000000000000     
  向量 t = I1 *x 1.0000000000000000 1.0000000000000000     
  
 和矩阵初始化
 M1 : I1                  
 M2 : I2                  
 总和：I1+I2               
 矩阵向量乘法matrix:I1+I2 M1: I1 M2: I2                  
  向量 t = (I1 +I2) *x 2.0000000000000000 2.0000000000000000     
  
 和矩阵初始化
 M1 : I2                  
 M2 : I3                  
 总和：I2+I3               
 和矩阵初始化
 M1 : I1                  
 M2 : I2+I3               
 总和：I1+I2+I3            
 矩阵向量乘法矩阵：I1+I2+I3 M1：I1 M2：I2+I3               
 矩阵向量乘法matrix:I2+I3 M1: I2 M2: I3                  
  向量 t = ( I1 + (I2 + I3) )*x 3.0000000000000000 3.0000000000000000     
  
 和矩阵初始化
 M1 : I1                  
 M2 : I2                  
 总和：I1+I2               
 和矩阵初始化
 M1 : I1+I2               
 M2 : I3                  
 总和：I1+I2+I3            
 矩阵向量乘法矩阵：I1+I2+I3 M1：I1+I2 M2：I3                  
 已分配（mat%scr）？F
 矩阵向量乘法matrix:I1+I2 M1: I1 M2: I2                  
在文件 LinearOperator.f90 的第 126 行
Fortran 运行时错误：未分配可分配的实际参数“this”

当我使用带有 2 个术语的 (+) 运算符时，一切正常。但是，当使用 3 个术语时，未分配的可分配数组 scr（类型为 (add_linop) 的成员）存在问题。

有谁知道这个问题的原因以及如何解决它？我包括用于编译代码的 Makefile。

#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 $<

标签： oopfortranoperator-overloadinggfortranderived-types

解决方案

您的程序不是有效的 Fortran。

的函数结果mmsum有一个指针组件，在函数执行期间，它是与虚拟参数相关联的指针。此虚拟参数（正确用于此用途）具有目标属性。但是，实际参数没有目标属性：当函数执行完成时，指针组件变为未定义的指针关联状态。

在子例程add_Mxv中尝试取消引用该指针。这是不允许的。

有必要重新审视在您的数据类型中如何处理操作数。请特别注意，表达式不能具有目标属性：在表达式的情况下，id1+id2+id3将id1+id2不会有用地保留为稍后引用的内容。

推荐阅读