I have finally succeeded in doing the step 2 of the first method. It wasn’t that hard, but it is a solution that only works for me, because I only fill diagonals in the matrix:
To insert values of an array of 2*n*(n-1) in the matrix, I do it this way :
>>> n = 4
>>> W = torch.rand(2,n*(n-1), requires_grad=True)
tensor([[ 0.4963,  0.7682,  0.0885,  0.1320,  0.3074,  0.6341,  0.4901, 0.8964,  0.4556,  0.6323,  0.3489,  0.4017],
        [ 0.0223,  0.1689,  0.2939,  0.5185,  0.6977,  0.8000,  0.1610, 0.2823,  0.6816,  0.9152,  0.3971,  0.8742]])
>>> D = W
>>> D1 = torch.nn.ConstantPad1d((0, 1), 0)(D[0].view(n, n - 1)).view(-1)[:-1]
>>> Dn = D[1]
>>> data = (Dn, D1, D1, Dn)
>>> offsets = torch.LongTensor([-n, -1, 1, n])
>>> M = [torch.diag(datum, offset) for datum, offset in zip(data, offsets)]
>>> A = torch.stack(M).sum(0)
# Calling .backward() on for example Z = torch.sum(A) works !
>>> A
tensor([[ 0.0000,  0.4963,  0.0000,  0.0000,  0.0223,  0.0000,  0.0000, 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.4963,  0.0000,  0.7682,  0.0000,  0.0000,  0.1689,  0.0000, 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.7682,  0.0000,  0.0885,  0.0000,  0.0000,  0.2939, 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0885,  0.0000,  0.0000,  0.0000,  0.0000, 0.5185,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0223,  0.0000,  0.0000,  0.0000,  0.0000,  0.1320,  0.0000, 0.0000,  0.6977,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.1689,  0.0000,  0.0000,  0.1320,  0.0000,  0.3074, 0.0000,  0.0000,  0.8000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.2939,  0.0000,  0.0000,  0.3074,  0.0000, 0.6341,  0.0000,  0.0000,  0.1610,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.5185,  0.0000,  0.0000,  0.6341, 0.0000,  0.0000,  0.0000,  0.0000,  0.2823,  0.0000,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.6977,  0.0000,  0.0000, 0.0000,  0.0000,  0.4901,  0.0000,  0.0000,  0.6816,  0.0000, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.8000,  0.0000, 0.0000,  0.4901,  0.0000,  0.8964,  0.0000,  0.0000,  0.9152, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.1610, 0.0000,  0.0000,  0.8964,  0.0000,  0.4556,  0.0000,  0.0000, 0.3971,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.2823,  0.0000,  0.0000,  0.4556,  0.0000,  0.0000,  0.0000, 0.0000,  0.8742],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.6816,  0.0000,  0.0000,  0.0000,  0.0000,  0.6323, 0.0000,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000,  0.9152,  0.0000,  0.0000,  0.6323,  0.0000, 0.3489,  0.0000],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000,  0.0000,  0.3971,  0.0000,  0.0000,  0.3489, 0.0000,  0.4017],
        [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000, 0.0000,  0.0000,  0.0000,  0.0000,  0.8742,  0.0000,  0.0000, 0.4017,  0.0000]])

Further explanations:
# Pad 1 zero every n-1 value, and remove the last zero padded.
>>> D1 = torch.nn.ConstantPad1d((0, 1), 0)(D[0].view(n, n-1)).view(-1)[:-1]
# There is the right number of values in D to fill the nth diagonal and have a resulting N*N matrix (N=N**2)
>>> Dn = D[1]
# List of diagonal data
>>> data = (Dn, D1, D1, Dn)
# List of offsets
>>> offsets = torch.LongTensor([-n, -1, 1, n])
# Creates a list of diagonal matrices for each diagonal in data.
>>> M = [torch.diag(datum, offset) for datum, offset in zip(data, offsets)]
# Sum all these matrices to get the final matrix
>>> A = torch.stack(M).sum(0)

I see a lot of functions in Pytorch to index in big arrays to get smaller arrays (torch.index_select(), torch.masked_select(), torch.take(), torch.where()), but very few (i.e. torch.diag()) to do the opposite, which is, placing values of a small arrays in a bigger one.
Does anyone know a reason for that ?
Thanks.