Finding matrix in Krylov subspace (Lanczos method)

The Lanczos method for finding the smallest eigenvalue of a hermiteian matrix $$H$$ is based on the construction of a vector subspace (Krylov space) where one can build a matrix $$H_{Krylov}$$ which is tridiagonal and which shares the smallest eigenvalue and the corresponding eigenvector with $$H$$. However through the functions Eigenvalues[] and Eigensystem[] one can only get the eigenvalue/eigenvector, but not the matrix $$H_{Krylov}$$. My question is: is there a simple way to get $$H_{Krylov}$$ from any Mathematica function or package? Of course I can write the Lanczos algorithm myself - and I did it - but I'm looking for some built in function to check out my code. I also know that "Krylov" is a Method for the function LinearSolve[], but still I can't figure out how to get $$H_{Krylov}$$... Thanks for helping!

Update

Following a comment by CA Trevillian, I paste below my Lanczos code. Although comments about this code are most welcome, I would really appreciate answers that are related to the question above.

    Swap[x_,y_]:=Module[{},Return[{y,x}]];

Lanczos[H_,threshold_,miniter_,maxiter_,startingvector_]:=Module[
{
a=ConstantArray[0,maxiter+1],(*initialize array of Subscript[a, n]*)
(* a[1]=Subscript[a, 0] , a[1+n]=Subscript[a, n]*)
b=ConstantArray[0,maxiter],(*initialize array of Subscript[b, n]*)
(* b[n]=Subscript[b, n] *)
dim=Length@H,
(*dimension of H*)
a0,b1,v,w,HKrylov,E0old,E0new,nfinal
},
(*initialize the starting vector: if startingvector=Null drop a default v*)
If[startingvector===Null,
v=SparseArray[{1->1.0},{dim}],
(*else*)v=startingvector
];

v=v/Norm[v];(*normalize it*)
w=H.v;
a0=(Conjugate@v).w;
w=w-a0*v;
b1=Norm[w];
a[[1]]=a0;b[[1]]=b1;

E0old=a0;
nfinal=maxiter;

Do[
If[b[[n]]<threshold&&n>=miniter,
nfinal=n;
Break[];
];
w=w/b[[n]];
v=-b[[n]]*v;
{v,w}=Swap[v,w];
w=w+H.v;
a[[n+1]]=(Conjugate@v).w;
w=w-a[[n+1]]*v;
If[n<maxiter,
b[[n+1]]=Norm[w];
];
(*Build hamiltonian in the Krylov subspace*)
HKrylov=DiagonalMatrix[
Table[b[[k]],{k,1,n}],
1];
HKrylov=HKrylov+HKrylov\[ConjugateTranspose];
HKrylov+=DiagonalMatrix[
Table[a[[k]],{k,1,1+n}]
];
E0new=Min@Eigenvalues@HKrylov;
If[Abs[E0new-E0old]<threshold&&n>miniter,
nfinal=n;
Break[];
];
E0old=E0new;
,{n,1,maxiter}];
a=Take[a,nfinal+1];
b=Take[b,nfinal];
Return[{E0new,a,b}]
];


The function requires a hermiteian matrix H, a convergence threshold (I suggest threshold=10^(-8) or less), a minimum and maximum number of Lanczos iterations and a starting vector which must have the same dimension of H. If you set startingvector=Null, then a default vector is dropped. The output is a list with the lowest eigenvalue of H (E0new), the list of parameters in the Krylov matrix (a and b).

• Can you, please, provide the code for your self-written Lanzcos algorithm? This will help others to be able to compare their attempts for accuracy, as well as helping future users to see the value of this QA by being able to reproduce the comparisons & so-on on their own. Thanks! Commented Jul 14, 2021 at 14:44
• I have added my code to the question following your suggestion. Please be careful since this code might have several problems... In fact I am still working on that and this is just a rough draft Commented Jul 14, 2021 at 15:12
• To test you code, check if your HKrylof matrix has the same Eigen-values/vectors as H Commented Jul 14, 2021 at 15:53
• Yes, the lowest eigenvalue of $H$ and $H_{Krylov}$ coincide, but I need a fast and well optimized algorithm to compute $H_{Krylov}$, that's why I'm looking for some built in function Commented Jul 14, 2021 at 19:54