I'm implementing an algebraic multigrid solver based on "A Multigrid Tutorial, 2ed", and in chapter 3, the formation of the prolongation operator is described as in the attached text below. My confusion arises from the equations that describe how the coefficients in the prolongation matrix P = I_{2h}^{h}
should be chosen, and it's not clear to me how one would iterate through j = 0:(n_div2_minus1)
to formulate P
. Any tips would be great, since the formulation of operations as matrices comes up a lot in multigrid literature (and in numerical methods in general) and I always find conversion of these operations into matrices moderately unintuitive (though for cases like finite differences/finite elements, the logic can be found in many texts).
julia> n_div2_min1 = 3
julia> n_min1 = 7
julia> P = zeros(n_min1, n_div2_min1)
julia> for i 1:n_min1 # rows in P, obviously this code won't work
for j in 0:n_div2_min1
P[i, j] = 1 # this would fail for P[1, 0]
P[i, j] = 1/2 # would overwrite the previous assignment
P[i, j+1] = 1/2
end
end
Please explain your logic and how it relates to the formulas in the description. I know that if I were simply shown a picture of the matrix, I could translate that to code, so that's not the part I'm as interested in, and I intentionally wrote the sample code to mimic (though obviously incorrectly) the equations in question.
Obviously the other answer didn't satisfy. Perhaps the desired answer is to fix the code in the question. This could be done in the following way:
for j in 1:n_div2_min1
P[2j, j] = 1
P[2j-1, j] = 1/2
P[2j+1, j] = 1/2
end
The equation description of the matrix in the question needs to be translated delicately:
j
runs over the short vectors and thus over the 1:n_div2_min1
range.j
and two j
s appear in RHS, one is for the contribution of the equation of the previous j
s value. Hence the use of 2j-1
.In any case, the mental model is a conversion of one vector into another, and therefore one for
loop needed. The for
loop helps create a matrix which is later used (this is why this is a "diagonal" matrix, as explained in the other answer).