Hi,

is there a way to print the matrices of an MPO/MPS? What I mean by that is, an MPO is a formed by multiplying matrices on each site. However, ITensors, stores these matrices as 4 indexed objects in case of MPOs and 3 indexed objects in case of MPS. Now, what I want is to print these 4 indexed and 3 indexed objects as matrices. In other words, I want to print the constituent matrices of MPO/MPS.

Hello,

You can easily print the constituent matrices using Julia array access notation. For example

```
julia> siteinds("S=1/2", 10)
julia> psi = randomMPS(s)
julia> @show psi[5]
psi[5] = ITensor ord=3
Dim 1: (dim=1|id=416|"Link,l=4")
Dim 2: (dim=2|id=192|"S=1/2,Site,n=5")
Dim 3: (dim=1|id=967|"Link,l=5")
Dense{Float64, Vector{Float64}}
1Ă—2Ă—1
[:, :, 1] =
-0.8970638949108982 0.4419008581653681
ITensor ord=3 (dim=1|id=416|"Link,l=4") (dim=2|id=192|"S=1/2,Site,n=5") (dim=1|id=967|"Link,l=5")
Dense{Float64, Vector{Float64}}
```

Hi Bhargava,

Karlâ€™s answer is right, but you may want to see the output organized so that the printing always shows the â€śmatricesâ€ť in the sense of how a â€śmatrix product stateâ€ť is a collection of d matrices on each site, one for every setting of the â€śphysical indexâ€ť of dimension d. And similar for an MPO thought of as a collection of d^2 indices on each site.

So I made a small function that permutes the index ordering of each tensor to put the physical or Site indices last. This way the printout shows each matrix grouped together. Please give it a try:

```
function print_matrices(A::AbstractMPS)
s = siteinds(A)
for j=1:length(s)
Aj = A[j]
links = uniqueinds(Aj,s[j])
Aj = permute(Aj,links...,collect(s[j])...)
println(Aj)
end
end
```

For example, when I use it to print a Hamiltonian MPO (XY model for S=1/2) I get:

```
Dim 1: (dim=4|id=324|"Link,l=1")
Dim 2: (dim=2|id=640|"Qubit,Site,n=1")'
Dim 3: (dim=2|id=640|"Qubit,Site,n=1")
NDTensors.Dense{Float64, Vector{Float64}}
4Ă—2Ă—2
[:, :, 1] =
0.0 0.0
0.0 0.0
0.0 0.5
1.0 0.0
[:, :, 2] =
0.0 0.0
0.5 0.0
0.0 0.0
0.0 1.0
ITensor ord=4
Dim 1: (dim=4|id=324|"Link,l=1")
Dim 2: (dim=4|id=576|"Link,l=2")
Dim 3: (dim=2|id=762|"Qubit,Site,n=2")'
Dim 4: (dim=2|id=762|"Qubit,Site,n=2")
NDTensors.Dense{Float64, Vector{Float64}}
4Ă—4Ă—2Ă—2
[:, :, 1, 1] =
1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0
[:, :, 2, 1] =
0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.5 0.0
[:, :, 1, 2] =
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0
0.0 0.5 0.0 0.0
[:, :, 2, 2] =
1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0
...
```

Hi Karl and Miles,

Thank you very much for the detailed reply. Your answers were very useful.