@@ -1186,39 +1186,27 @@ function _dibimul!(C, A, B, _add)
1186
1186
end
1187
1187
function _dibimul_nonzeroalpha! (C, A, B, _add)
1188
1188
n = size (A,1 )
1189
- if n <= 3
1190
- # For simplicity, use a naive multiplication for small matrices
1191
- # that loops over all elements.
1192
- for I in CartesianIndices (C)
1193
- C[I] += _add (A. diag[I[1 ]] * B[I[1 ], I[2 ]])
1194
- end
1195
- return C
1196
- end
1197
1189
Ad = A. diag
1198
1190
Bl = _diag (B, - 1 )
1199
1191
Bd = _diag (B, 0 )
1200
1192
Bu = _diag (B, 1 )
1201
1193
@inbounds begin
1202
1194
# first row of C
1203
- C[1 ,1 ] += _add (A[1 ,1 ]* B[1 ,1 ])
1204
- C[1 ,2 ] += _add (A[1 ,1 ]* B[1 ,2 ])
1205
- # second row of C
1206
- C[2 ,1 ] += _add (A[2 ,2 ]* B[2 ,1 ])
1207
- C[2 ,2 ] += _add (A[2 ,2 ]* B[2 ,2 ])
1208
- C[2 ,3 ] += _add (A[2 ,2 ]* B[2 ,3 ])
1209
- for j in 3 : n- 2
1195
+ C[1 ,1 ] += _add (Ad[1 ]* Bd[1 ])
1196
+ if n >= 2
1197
+ C[1 ,2 ] += _add (Ad[1 ]* Bu[1 ])
1198
+ end
1199
+ for j in 2 : n- 1
1210
1200
Ajj = Ad[j]
1211
1201
C[j, j- 1 ] += _add (Ajj* Bl[j- 1 ])
1212
1202
C[j, j ] += _add (Ajj* Bd[j])
1213
1203
C[j, j+ 1 ] += _add (Ajj* Bu[j])
1214
1204
end
1215
- # row before last of C
1216
- C[n- 1 ,n- 2 ] += _add (A[n- 1 ,n- 1 ]* B[n- 1 ,n- 2 ])
1217
- C[n- 1 ,n- 1 ] += _add (A[n- 1 ,n- 1 ]* B[n- 1 ,n- 1 ])
1218
- C[n- 1 ,n ] += _add (A[n- 1 ,n- 1 ]* B[n- 1 ,n ])
1219
- # last row of C
1220
- C[n,n- 1 ] += _add (A[n,n]* B[n,n- 1 ])
1221
- C[n,n ] += _add (A[n,n]* B[n,n ])
1205
+ if n >= 2
1206
+ # last row of C
1207
+ C[n,n- 1 ] += _add (Ad[n]* Bl[n- 1 ])
1208
+ C[n,n ] += _add (Ad[n]* Bd[n])
1209
+ end
1222
1210
end # inbounds
1223
1211
C
1224
1212
end
0 commit comments