# Issues using linsolve (cont'd discussion)

Continuing the discussion from Constant shift in MPS shift-and-invert algo:

I finally have time to come back to this problem. In short, I’m trying to variationally solve the following linear system:
(H - \lambda) |\psi\rangle = |\phi\rangle

Or that

o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o-      -o--o--o--o--o--o |psi>

=
o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o |phi>


(For now let’s assume \lambda=0, since I believe you can always add back the cross terms ( \lambda \langle \psi | H | \psi\rangle , \lambda^2\langle \psi|\psi\rangle, \lambda\langle \psi|\phi\rangle ) which doesn’t change the problem)
(Or there’s a way to add constant terms?)

The way I approach this, is that I defined a new position! function:

function position!(P::AbstractProjMPO, psi::MPS, phi::MPS, pos::Int)
L = makeL!(P, psi, phi, pos - 1)
R = makeR!(P, psi, phi, pos + nsite(P))

# start of the left end
if L == nothing
return phi[pos] * phi[pos + 1] * R

# start of the right end
elseif R == nothing
return L * phi[pos] * phi[pos + 1]

# in the middle
else
return L * phi[pos] * phi[pos + 1] * R

end
end


with appropriate modification to makeL! and makeR!

function _makeR!(P::AbstractProjMPO, psi::MPS, phi::MPS, k::Int)
"""everything were kept the same except"""
R = R * phi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))
"""instead of R = R * psi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))"""
end


Eventually I defined a hamiltonian ‘squared’:

  H2 = contract(H', H; cutoff=1e-12)
H2 = replaceprime(H2, 2 => 1)


And I calculate the LHS and RHS and a guess

RHS = position!(PH, psi, phi, b)
position!(PH2, psi, b)
guess = psi[b] * psi[b + 1]


Put it through linsolve

          vecs, conv = linsolve(
PH2,
RHS,
guess,  0, 1;
ishermitian=ishermitian,
tol=linsolve_tol,
krylovdim=linsolve_krylovdim,
maxiter=linsolve_maxiter,
)


This error initially occurred
LoadError: DimensionMismatch("In scalar(T) or T[], ITensor T is not a scalar

After I removed ‘guess’ from the linsolve arguments, another error

LoadError: The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, this is probably due to an index mismatch.
Common reasons for this error:
(1) You are trying to multiply the ProjMPO with the 2-site wave-function at the wrong position.
(2) orthogonalize! was called, changing the MPS without updating the ProjMPO.`

Essentially, the LHS (The PH2 here) uses all the unmodified methods of ProjMPS, it’s just that I used a different Hamiltonian, which I believe is done correctly. I can run the LHS Hamiltonian through eigsolve as well.

I guess there’s something wrong in the way I modified the RHS methods, from contraction between the same state \langle \psi | H|\psi\rangle to \langle \psi | H |\phi\rangle, but it’s not immediately clear to me what was the issue.