Search code examples
matrixjulia

Forming a matrix (i.e., linear operator) based on a description of coefficients


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.

enter image description here


Solution

  • 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:

    1. The equations are linear, and therefore translated into a matrix relating the elements in the LHS long vectors, to the RHS short vectors.
    2. j runs over the short vectors and thus over the 1:n_div2_min1 range.
    3. The tricky bit is to translate the indices in the appearance of two elements in RHS of the second equation. Since we are looping over j and two js appear in RHS, one is for the contribution of the equation of the previous js 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).