PathToValues="C:/users/s_pauli/Desktop/";

ReadFull1ArrayVle[handle_InputStream,systemname_:"",namefromfile_:""]:=Module[{k,whereto},
line=Read[handle,String];Print[line];
If[namefromfile=!=""&&namefromfile=!=line,
Print["Wrong name ",line," expected ", namefromfile];
Abort[]];
name=If[systemname==="",line,systemname];
prec=Read[handle,Number];
start=Read[handle,Number];
end=Read[handle,Number];
step=Read[handle,Number];
For[k=start,k≤end,k=k+step,
val=Read[handle,String];
pos=StringPosition[val,"e"];
If[pos=!={},
val=StringReplace[val,"e"->"*10^"]];
ToExpression[name<>"["<>ToString[k,InputForm]<>"]="<>val]
];
]

Loading precomputed values of Log[n] as NLogk[n]

Map[Close,Drop[Streams[],2]];
handle=OpenRead[PathToValues<>"nlogk.vle",BinaryFormat->True]
ReadFull1ArrayVle[handle]
Close[handle];

Loading Imaginary parts of zeta zeroes as NImZetaZero[n]

Map[Close,Drop[Streams[],2]];
handle=OpenRead[PathToValues<>"imLzero_001_001.vle",BinaryFormat->True]
ReadFull1ArrayVle[handle]
Close[handle];

Loading Imaginary parts of L[3,2] as imLzero32[n]

Map[Close,Drop[Streams[],2]];
handle=OpenRead[PathToValues<>"imLzero_003_002.vle",BinaryFormat->True]
ReadFull1ArrayVle[handle,"imLzero32"]
Close[handle];

TrueZeroes[ ], Multiplicities[ ], and ns[ ] contain input data

NImZetaZero[2]//N

GeneralFillDeltaMatrixComplex[TrueZeroes_,Multiplicities_,ns_,Nsize_]:=
Module[{answer,nplace,n,zeroplace,truezero,m},
answer=Table[0,{nplace,1,Nsize},{zeroplace,1,Nsize}];
For[nplace=1,nplace<=Nsize,nplace++,
Prin["nplace=",nplace];
n=ns[[nplace]];
zeroplace=1;zeronumber=1;
While[zeroplace≤Nsize,Prin["zeroplace=",zeroplace," mult=",Multiplicities[[zeronumber]]];
  truezero=TrueZeroes[[zeronumber]];
m=1;
While[zeroplace≤Nsize&&m≤Multiplicities[[zeronumber]],
Prin["m=",m];
answer[[nplace,zeroplace]]=n^(-truezero)*If[n===1&&m===1,1,(NLogk[n])^(m-1)];
zeroplace++;m++];
zeronumber++]];
answer]

Setting parameters

MatrixSize=17; Nprecision=100;

ActualCaliculations

tz=Table[N[{1/2+I*NImZetaZero[k],1/2-I*NImZetaZero[k]},Nprecision],{k,1,MatrixSize/2+1}]//Flatten;
a=GeneralFillDeltaMatrixComplex[tz,Table[1,{MatrixSize}],Table[k,{k,1,MatrixSize}],MatrixSize];
DIAGONAL=diag=Table[0,{MatrixSize+2}]; L=Table[0,{MatrixSize+2},{MatrixSize+2}];
For[i=1,i≤MatrixSize,i++,
diag [[i]]=   a [[i,i]];
    DIAGONAL [[i]]=   diag [[i ]]*If[i===1,1,   DIAGONAL [[i-1  ]]];
For[j=i+1,j≤MatrixSize,j++,
z=-   a [[j,i  ]]/   diag [[i  ]];
    L [[j,i]]=z;    a [[j,i]]=z;
For[k=i+1,k≤MatrixSize,k++,
  a [[j,k]]=   a [[j,k   ]]+   a [[i,k  ]]*   a [[j,i ]];L[[j,k]]=a[[j,k]];
];
]];
For[k=1,k≤MatrixSize,k++,
For[i=k+1,i≤MatrixSize,i++,For[j=i+1,j≤MatrixSize,j++,
L [[j,k]]=   L [[j,k  ]]+   L [[i,k  ]]*   a [[j,i]];];
];
];
For[j=1,j≤MatrixSize,j++,
For[k=1,k<j,k++,
LinVectElemNew[j,k]=L[[j,k]]*If[j===1,1,DIAGONAL[[j-1]]]];
LinVectElemNew[j,j]=If[j===1,1,DIAGONAL[[j-1]]];
InvNTildedelta[j]=1/LinVectElemNew[j,1];
For[k=1,k≤j,k++,
Ndelta[j,k]=LinVectElemNew[j,k]*InvNTildedelta[j]];
];

Viewing the results

(tt=Table[Ndelta[MatrixSize,n]//Re,{n,1,MatrixSize}])//N

tt//ListPlot