diff --git a/src/gmat.jl b/src/gmat.jl
index 62df0880..4e2786ce 100644
--- a/src/gmat.jl
+++ b/src/gmat.jl
@@ -57,9 +57,10 @@ function gmat!(mx, θ, ::AR_)
         mx[i, i] = de
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = de * θ[2] ^ (n - m)
+        @inbounds θ2 = θ[2]
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = de * θ2 ^ (n - m)
             end
         end
     end
@@ -72,14 +73,16 @@ function gmat!(mx, θ, ::ARH_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = mx[m, m] * mx[n, n] * θ[end] ^ (n - m)
+        θe = last(θ)
+        for n = 2:s
+            mxnn = mx[n, n]
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = mx[m, m] * mxnn * θe ^ (n - m)
             end
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m] 
     end
     mx
 end
@@ -105,14 +108,14 @@ function gmat!(mx, θ, ::CSH_)
     end
     if s > 1
         for n = 2:s
-            @inbounds mxnθe = mx[n, n] * θ[end]
+            @inbounds mxnθe = mx[n, n] * last(θ)
             @inbounds @simd for m = 1:n-1
                 mx[m, n] = mx[m, m] * mxnθe
             end
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
@@ -125,9 +128,10 @@ function gmat!(mx, θ, ::ARMA_)
         mx[i, i] = de
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = de * θ[2] * θ[3] ^ (n - m - 1)
+        deθ2 = de * θ[2]
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = deθ2 * θ[3] ^ (n - m - 1)
             end
         end
     end
@@ -141,8 +145,8 @@ function gmat!(mx, θ, ::TOEP_)
         mx[i, i] = de
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] = de * θ[n-m+1]
             end
         end
@@ -171,14 +175,14 @@ function gmat!(mx, θ, ::TOEPH_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] = mx[m, m] * mx[n, n] * θ[n-m+s]
             end
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
@@ -196,7 +200,7 @@ function gmat!(mx, θ, ct::TOEPHP_)
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
@@ -207,15 +211,14 @@ function gmat!(mx, θ, ::UN_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] = mx[m, m] * mx[n, n] * θ[s+tpnum(m, n, s)]
             end
         end
     end
     @inbounds @simd for m = 1:s
-        v = mx[m, m]
-        mx[m, m] = v * v
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
diff --git a/src/rmat.jl b/src/rmat.jl
index 852aba34..1966305a 100644
--- a/src/rmat.jl
+++ b/src/rmat.jl
@@ -3,22 +3,6 @@
 ################################################################################
 
 ################################################################################
-#=
-function rmat_base_inc_b!(mx, θ, zrv, covstr)
-    rmat!(mx, θ, zrv, covstr.repeated.covtype.s)
-end
-=#
-################################################################################
-################################################################################
-#=
-function rmat_base_inc!(mx, θ, covstr, block, sblock)
-    zblock    = view(covstr.rz, block, :)
-    @simd for i ∈ axes(sblock[end], 1)
-        @inbounds rmat_base_inc_b!(view(mx, sblock[end][i],  sblock[end][i]), θ, view(zblock,  sblock[end][i], :), covstr)
-    end
-    mx
-end
-=#
 @noinline function rmat_base_inc!(mx, θ, covstr, bi)
     en        = covstr.rn + 1
     block     = covstr.vcovblock[bi]
@@ -43,23 +27,23 @@ Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_)
 end
 #DIAG
 function rmat!(mx, θ, rz, ::DIAG_)
-    #=@turbo=#  for i ∈ axes(mx, 1)
+    for i ∈ axes(mx, 1)
         @inbounds @simd for c ∈ axes(θ, 1)
-            mx[i, i] += rz[i, c] * θ[c] * θ[c]
+            mx[i, i] += rz[i, c] * θ[c] ^ 2
         end
     end
     mx
 end
 #AR
 function rmat!(mx, θ, ::AbstractMatrix, ::AR_)
-    rn  = size(mx, 1)
+    s  = size(mx, 1)
     de  = θ[1] ^ 2
-    @inbounds @simd for m = 1:rn
+    @inbounds @simd for m = 1:s
         mx[m, m] += de
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @inbounds @simd for n = m + 1:rn
+    if s > 1
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] += de * θ[2] ^ (n - m)
             end
         end
@@ -69,10 +53,10 @@ end
 #ARH
 function rmat!(mx, θ, rz, ::ARH_)
     vec = tmul_unsafe(rz, θ)
-    rn    = size(mx, 1)
-    if rn > 1
-        for m = 1:rn - 1
-            @inbounds @simd for n = m + 1:rn
+    s    = size(mx, 1)
+    if s > 1
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                  mx[m, n] += vec[m] * vec[n] * θ[end] ^ (n - m)
             end
         end
@@ -84,15 +68,15 @@ function rmat!(mx, θ, rz, ::ARH_)
 end
 #CS
 function rmat!(mx, θ, ::AbstractMatrix,  ::CS_)
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     θsq   =  θ[1]*θ[1]
     θsqp  =  θsq*θ[2]
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] += θsq
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @inbounds @simd for n = m + 1:rn
+    if s > 1
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] += θsqp
             end
         end
@@ -120,15 +104,16 @@ end
 ################################################################################
 #ARMA
 function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_)
-    rn  = size(mx, 1)
+    s  = size(mx, 1)
     de  = θ[1] ^ 2
-    @inbounds @simd for m = 1:rn
+    @inbounds @simd for m = 1:s
         mx[m, m] += de
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @inbounds @simd for n = m + 1:rn
-                mx[m, n] += de * θ[2] * θ[3] ^ (n - m - 1)
+    if s > 1
+        deθ2 = de * θ[2]
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] += deθ2 * θ[3] ^ (n - m - 1)
             end
         end
     end
@@ -195,13 +180,13 @@ function rmat!(mx, θ, rz,  ::SPEXP_)
     θe    = θ[2]
     θe    = abs(θe) < eps() ? sqrt(eps()) : abs(θe)
     #θe    = abs(θ[2])
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     @simd for i = 1:size(mx, 1)
         @inbounds mx[i, i] += σ²
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @simd for n = m + 1:rn
+    if s > 1
+        for m = 1:s - 1
+            @simd for n = m + 1:s
                 mx[m, n] += σ² * exp(-edistance(rz, m, n) / θe)
             end
         end
@@ -213,13 +198,13 @@ end
 function rmat!(mx, θ, rz,  ::SPPOW_)
     σ²    = θ[1]^2
     ρ     = θ[2]
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     @simd for i = 1:size(mx, 1)
         @inbounds mx[i, i] += σ²
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @simd for n = m + 1:rn
+    if s > 1
+        for m = 1:s - 1
+            @simd for n = m + 1:s
                 ed = edistance(rz, m, n)
                 mx[m, n] += ρ > 0 ? σ² * ρ^ed : σ² * ρ^ed * cos(π * ed)
                 #mx[m, n] += ρ > 0 ? σ² * ρ^ed : σ² * abs(ρ)^ed * sign(ρ)
@@ -237,13 +222,13 @@ function rmat!(mx, θ, rz,  ::SPGAU_)
     θe    = abs(θe) < eps() ? sqrt(eps()) : θe^2
     #θe    = θ[2]
     #θe    = abs(θ[2])
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     @simd for i = 1:size(mx, 1)
         @inbounds mx[i, i] += σ²
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @simd for n = m + 1:rn
+    if s > 1
+        for m = 1:s - 1
+            @simd for n = m + 1:s
                 mx[m, n] += σ² * exp(- (edistance(rz, m, n)^2) / θe)
             end
         end
@@ -257,13 +242,13 @@ function rmat!(mx, θ, rz,  ::SPEXPD_)
     σ²d   = θ[1]^2 + σ²
     θe    = θ[3]
     θe    = iszero(θe) ? 1e-16 : abs(θe)
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     @simd for i = 1:size(mx, 1)
         @inbounds mx[i, i] += σ²d
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @simd for n = m + 1:rn
+    if s > 1
+        for m = 1:s - 1
+            @simd for n = m + 1:s
                 mx[m, n] += σ² * exp(-edistance(rz, m, n) / θe)
             end
         end
@@ -275,13 +260,13 @@ function rmat!(mx, θ, rz,  ::SPPOWD_)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     ρ     = θ[3]
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     @simd for i = 1:size(mx, 1)
         @inbounds mx[i, i] += σ²d
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @simd for n = m + 1:rn
+    if s > 1
+        for m = 1:s - 1
+            @simd for n = m + 1:s
                 mx[m, n] += σ² * ρ^edistance(rz, m, n)
             end
         end
@@ -294,13 +279,13 @@ function rmat!(mx, θ, rz,  ::SPGAUD_)
     σ²d   = θ[1]^2 + σ²
     θe    = θ[3]
     θe    = iszero(θe) ? 1e-16 : θe^2
-    rn    = size(mx, 1)
+    s    = size(mx, 1)
     @simd for i = 1:size(mx, 1)
         @inbounds mx[i, i] += σ²d
     end
-    if rn > 1
-        for m = 1:rn - 1
-            @simd for n = m + 1:rn
+    if s > 1
+        for m = 1:s - 1
+            @simd for n = m + 1:s
                 mx[m, n] += σ² * exp(- (edistance(rz, m, n)^2) / θe)
             end
         end