1

我正在尝试在 Fortran90（使用 ifort 编译）中编写一个代码，其中我将两个矩阵相乘。我正在为此编写代码，因为其中一个矩阵是稀疏的，因此您可以在不为整个矩阵分配内存的情况下进行乘法运算。

我有两个子程序。第一个，将稀疏矩阵（k在对角线上，在对角线l两侧）和一个向量（b）相乘。结果通过指针传递给主函数r。我选择使用子例程而不使用函数，因为这样，我不需要在子例程中再次分配内存。

    subroutine matSparXVec(k, l, b, r)

    implicit none
    real, intent(in)        ::k, l
    real, intent(in), dimension(:)  ::b
    real, intent(out), dimension(:) ::r
    integer             ::ierr, i

    r = (/&
        k*b(1)+l*b(2), &
        (l*b(i-1)+k*b(i)+l*b(i+1),i=2,size(b)-1),& 
        l*b(size(b)-1)+k*b(size(b))&
        /)

end subroutine matSparXVec

第二个子程序使用第一个子程序将稀疏矩阵 ( k, l) 与另一个矩阵 ( B) 相乘：

subroutine matSparXMat(k, l, B, R)
    implicit none
    real, intent(in)            ::k, l
    real, intent(in), dimension(:,:)    ::B
    real, intent(out), dimension(:,:)   ::R

    call  matsparXVec(k, l, B(1, :), R(1, :))
    R(1, :) = R(1, :) + l * B(1, :)

    do i = 2, (size(R)-2)
        call matsparXVec(k,l,B(i,:),R(i,:))
        R(i,:)=R(i-1,:)+R(i,:)
    enddo 

    call matsparXVec(k,l,B(size(B),:),R(size(R),:))
    R(size(R),:) = R(size(R),:) + l * B(size(B)-1,:)


end subroutine matSparXmat

现在的问题是，在子例程的某个地方matSparXmat，我正在修改B(:,:). 例如，使用代码：

    implicit none
    real, dimension(:,:), allocatable   ::BB, RR
    integer                 ::i, j, ierr, n, m
    real, parameter             ::k = 4.0, l = 1.0

    n=3 !Dimension vector
    m=3 !Dimension del segundo orden

    allocate(RR(n,m), stat=ierr)
    allocate(BB(n,m), stat=ierr)

    forall(i = 1:size(BB(:,1)), j=1:size(BB(1,:)))
        BB(i,j)=i+j
        RR(i,j) = 0 
    endforall

    do i=1,size(BB(:,1))
        print *, BB(:,i)
    enddo

    call matSparXMat(k, l, BB, RR)

    do i=1,size(BB(:,1))
        print *, BB(:,i)
    enddo

我得到输出：

 2.000000       3.000000       4.000000    
   3.000000       4.000000       5.000000    
   4.000000       5.000000       6.000000  

  4.6526092E+33   3.000000      1.9366391E+31
   3.000000       4.000000       5.000000    
   4.000000       5.000000       6.000000

您可以在其中看到 的值BB已被修改。

4

1 回答 1

1

看起来这里出错了：

do i = 2, (size(R)-2)                      !<---- here
   call matsparXVec(k,l,B(i,:),R(i,:))
   R(i,:)=R(i-1,:)+R(i,:)
enddo

size(R)=9但是size(R(:,1))=3。但是，使用 gfortran 我确实得到了错误的正确输出

*** glibc detected *** ./sparse_mults: free(): invalid next size (fast): 0x00000000014e4010 ***

但是当我使用 ifort 2013 时，我得到了1.93...E+31右上角的值，否则它是正确的。不知道发生了什么，但如果我想出一些东西，我会告诉你的。

于 2013-08-29T20:11:27.233 回答