Replace the diagonal of a list of matrices

Let's generate sample data:

array = ConstantArray[
   RandomInteger[{0, 10}, {100, 100}], {100, 100}];
replace = ConstantArray[RandomInteger[{0, 10}, {100, 100}], 100];

array1 = array2 = array3 = array4 = array;

and apply different methods:

RepeatedTiming[
 Table[array1[[i, i]] = replace[[i]], {i, Length@replace}];]

{0.0064, Null}

RepeatedTiming[(array2[[#, #]] = replace[[#]]) & /@ 
   Range[Length@replace];]

{0.0071, Null}

RepeatedTiming[array3=ReplacePart[array3, {i_, i_} :> replace[[i]]];]

{12.9, Null}

RepeatedTiming[
 ParallelTable[array4[[i, i]] = replace[[i]], {i, Length@replace}];]

{0.016, Null}

array1 == array2 == array3 == array4

True

Suggested by Alx (the fastest so far):

RepeatedTiming[
 Do[array1[[i, i]] = replace[[i]], {i, Length@replace}];]

{0.0018, Null}

We can also use SparseArrays.

m1 = {{{{1, 2}, {3, 4}}, {{5, 6}, {7, 8}}}, {{{9, 10}, {11, 12}}, {{13, 14}, {15, 16}}}};
v1 = {{{11, 3}, {7, 5}}, {{13, 97}, {3, 16}}};
res = {{{{11, 3}, {7, 5}}, {{5, 6}, {7, 8}}}, {{{9, 10}, {11, 12}}, {{13, 97}, {3, 16}}}};

blockmatrixreplace[mat_, rep_] := Block[
 {
  diagonalmatrixindices = Flatten[Table[{i, i, j, k}, ##] & @@ ({#, Dimensions[mat][[1]]} & /@ {i, j, k}), 2],
  emptydiag,
  repten
 },
 emptydiag = SparseArray[mat] - SparseArray[diagonalmatrixindices -> Flatten@Diagonal@mat, Dimensions[mat]];
 repten = SparseArray[diagonalmatrixindices -> Flatten@rep, Dimensions[mat]];
 emptydiag + repten // Normal
]

{#1, #2 == res} & @@ RepeatedTiming@blockmatrixreplace[m1, v1]

{0.000055, True}

Or we could use Band in the SparseArray since

With an array a of the same rank as the whole sparse array, Band[start]->a by default inserts a at the position specified by start.

bandblockmatrixreplace[mat_, rep_] := Block[
 {diagreps, offdiagreps, reps},
 diagreps = MapIndexed[Band[{#2[[1]], #2[[1]], 1, 1}] -> {{#1}} &, rep];
 offdiagreps = Table[Band[{j, k, 1, 1}] -> {{mat[[j, k]]}}, {j, #1}, {k, Range[#2]~Complement~{j}}] & @@ Dimensions[mat][[;; 2]] // Flatten;
 reps = Join[diagreps, offdiagreps];
 Normal@SparseArray[reps, Dimensions@mat]
]
{#1, #2 == res} & @@ RepeatedTiming@bandblockmatrixreplace[m1, v1]

{0.000070, True}