另外，按照你的代码，子程序必须是内部子程序。

 parameter(SIZE=2) !方阵的阶数，待定，可修改。
      parameter(EPSILON=0.0001) !小量，待定，可修改。（矩阵元小于EPSILON时认为等于0）
      parameter(ETA=0.0001) !小量，待定，可修改。作为w的虚部。

 complex*16 a(SIZE,SIZE),b(SIZE,SIZE),c(SIZE,SIZE),w(SIZE,SIZE),H_00(SIZE,SIZE),H_01(SIZE,SIZE)
 complex*16 G_00(SIZE,SIZE),G00bar(SIZE,SIZE),G_b(SIZE,SIZE),w_ep_s(SIZE,SIZE),wepsbar(SIZE,SIZE)
 complex*16 trG00(SIZE,SIZE),trG00b(SIZE,SIZE),trGb(SIZE,SIZE)    
 complex m(SIZE,SIZE)


 data a/(1.0,0.0),(1.0,0.0),(1.0,0.0),(0.0,0.0)/
 data b/(1.0,0.0),(0.0,0.0),(1.0,0.0),(1.0,0.0)/
 c=matmul(a,b)
 b=transpose(b)
 call prtmtx(a)
 pause
 call prtmtx(b)
 pause
 call prtmtx(c)

 contains

 
 subroutine prtmtx(m)
 complex*16 m(:,:)
 integer i,j
 do i=1,SIZE
 do j=1,SIZE
 if(abs(real(m(i,j))).GT.EPSILON.and.aimag(m(i,j)).GT.EPSILON)then
write(*,10)i,j,real(m(i,j)),aimag(m(i,j))
10 format(1x,I3,'行',I3,'列',e12.4,'+',e12.4,'i     ')
else if(abs(real(m(i,j))).GT.EPSILON.and.aimag(m(i,j)).LT.-EPSILON)then
write(*,20)i,j,real(m(i,j)),aimag(m(i,j))
20 format(1x,I3,'行',I3,'列',e12.4,e13.4,'i     ')
else if(abs(real(m(i,j))).GT.EPSILON)then
write(*,30)i,j,real(m(i,j))
30 format(1x,I3,'行',I3,'列',e12.4,'                   ')
else if(abs(aimag(m(i,j))).GT.EPSILON)then
write(*,40)i,j,aimag(m(i,j))
40 format(1x,I3,'行',I3,'列','            ',e13.4,'i     ')
else
write(*,50)
50 format(1x,I3,'行',I3,'列','           0                   ')
endif
 end do
 end do

 end subroutine
 end program