@@ -108,12 +108,13 @@ contains
108108 p = .true.
109109 do k = 2, order_
110110 c = c * (order_ - k + 1) / (k * (2*order_ - k + 1))
111- X_tmp = X
112- #:if rt.startswith('complex')
113- call gemm("N", "N", n, n, n, one_c${rk}$, A2, n, X_tmp, n, zero_c${rk}$, X, n)
114- #:else
115- call gemm("N", "N", n, n, n, one_${rk}$, A2, n, X_tmp, n, zero_${rk}$, X, n)
116- #:endif
111+ X = matmul(A2, X)
112+ ! X_tmp = X
113+ ! #:if rt.startswith('complex')
114+ ! call gemm("N", "N", n, n, n, one_c${rk}$, A2, n, X_tmp, n, zero_c${rk}$, X, n)
115+ ! #:else
116+ ! call gemm("N", "N", n, n, n, one_${rk}$, A2, n, X_tmp, n, zero_${rk}$, X, n)
117+ ! #:endif
117118 do concurrent(i=1:n, j=1:n)
118119 A(i, j) = A(i, j) + c*X(i, j) ! E = E + c*X
119120 enddo
@@ -140,12 +141,13 @@ contains
140141 block
141142 ${rt}$ :: E_tmp(n, n)
142143 do k = 1, s
143- E_tmp = A
144- #:if rt.startswith('complex')
145- call gemm("N", "N", n, n, n, one_c${rk}$, E_tmp, n, E_tmp, n, zero_c${rk}$, A, n)
146- #:else
147- call gemm("N", "N", n, n, n, one_${rk}$, E_tmp, n, E_tmp, n, zero_${rk}$, A, n)
148- #:endif
144+ A = matmul(A, A)
145+ ! E_tmp = A
146+ ! #:if rt.startswith('complex')
147+ ! call gemm("N", "N", n, n, n, one_c${rk}$, E_tmp, n, E_tmp, n, zero_c${rk}$, A, n)
148+ ! #:else
149+ ! call gemm("N", "N", n, n, n, one_${rk}$, E_tmp, n, E_tmp, n, zero_${rk}$, A, n)
150+ ! #:endif
149151 enddo
150152 end block
151153 endif
0 commit comments