Is this the correct procedure for calculating matrix spectrum?

Question

I am not sure if my question is on topic but I have a piece of Fortran code that is used to perform successive over relaxation. Prior to performing successive over relaxation the author is calculating the spectrum of the matrix.There is very little documentation so I am just trying to figure out all the details for documentation. It appears from the code comments that the following code calculates matrix spectrum. Would it be possible to get a mathematical background for this code ? Here is the equation for the successive over relaxation as requested.

The linearized second order partial differential equation(Laplacian) that I am trying to solve is this

$$\Delta_x(A.\Delta_x\psi) + \Delta_y(B.\Delta_y\psi) + \Delta_z(C.\Delta_z\psi) + S = 0$$

where $$\Delta_x(A.\Delta_x\psi) = A_{i+1/2,j,k} .(\psi_{i+1,j,k}-\psi_{i,j,k}) - A_{i-1/2,j,k} .(\psi_{i,j,k}-\psi_{i-1,j,k})$$

$$\Delta_y(B.\Delta_y\psi) = B_{i,j+1/2,k} .(\psi_{i,j+1,k}-\psi_{i,j,k}) - B_{i,j-1/2,k} .(\psi_{i,j,k}-\psi_{i,j-1,k})$$

$$\Delta_z(C.\Delta_z\psi) = C_{i,j,k+1/2} .(\psi_{i,j,k+1}-\psi_{i,j,k}) - A_{i,j,k-1/2} .(\psi_{i,j,k}-\psi_{i,j,k-1})$$

The additive constant is

$$ S_{i,j,k} = \Delta_x .\Delta_y . \Delta_z * \sigma_{i,j,k} $$

$$ \psi_i^{(n+1)} = \omega\left( b_i - \sum_{j=1}^{i-1} A_{ij} \psi_j^{(n+1)} - \sum_{j=i}^{m}A_{ij}\psi_j^{(n)}\right) + (1-\omega)\psi_i^{(n)}. $$

  i=nx/2
  specx=4.*a(i,0,0)/
 >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
  specy=2.*(b(i,0,0)+b(i,1,0))/
 >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
  specz=2.*(c(i,0,0)+c(i,0,1))/
 >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
  do k=1,nz-2
     do j=1,ny-2

        helpx=4.*a(i,j,k)/
 >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
        if (helpx.gt.specx) specx=helpx

        helpy=2.*(b(i,j,k)+b(i,j+1,k))/
 >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
        if (helpy.gt.specy) specy=helpy

        helpz=2.*(c(i,j,k)+c(i,j,k+1))/
 >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
        if (helpz.gt.specz) specz=helpz

     enddo
  enddo

The following code fragments calculate the a,b,c matrices

  do i=0,nx
     do j=0,ny
        do k=0,nz
           a(i,j,k)=dy*dz*rhoref(2*k)/(dx*coriol(i,j))
           if (j.lt.ny) then
              b(i,j,k)=dx*dz*rhoref(2*k)/
 >                     (dy*0.5*(coriol(i,j)+coriol(i,j+1)))
           else
              b(i,j,k)=0.
           endif
           if (k.lt.nz) then
              c(i,j,k)=dx*dy*rhoref(2*k+1)*coriol(i,j)/
 >                     (dz*nsq(2*k+1))
           else
              c(i,j,k)=0.
           endif
        enddo
     enddo
  enddo

What are $a$, $b$, and $c$, and what linear system is it solving? — Kirill, Jun 29 '15 at 05:24
No, I meant what are the matrices precisely? Give details, expressions, anything! I don't know that equation. How do the matrices enter into it? What is the equation exactly, and how does the code try to solve it? — Kirill, Jun 29 '15 at 05:41
I converted your formula from JPEG to LaTeX (BTW, compare it with the similar one at https://en.wikipedia.org/wiki/Successive_over-relaxation), but what has it got to do with the code at hand? See my first comment. — Kirill, Jun 29 '15 at 11:26
@Kirill - Better late than never. I answered your comment. I hope :-) — , Sep 07 '16 at 15:18

score 3 · Accepted Answer · edited Apr 13 '17 at 12:53

The spectrum of a matrix is a set of numbers, so since there isn't a set of numbers being computed here, it's almost certainly not the spectrum. Also relevant is that you have three tensors ($a$, $b$, and $c$), so which spectrum do you think it would calculate? Calculating the spectrum of a matrix numerically also cannot be done in a finite number of rational operations and requires an iterative algorithm, unless something very special is known about the matrix.

Why do you think it calculates a spectrum of any kind?

You mention it does this before performing successive over-relaxation. SOR requires a guess as to the best value of the relaxation parameter, which is related to the spectral norm of a matrix (see, e.g., Are there any heuristics for optimizing the successive over-relaxation (SOR) method?). So I guess it's possible it's trying to compute a heuristic estimate of something related to the spectral norm of the matrix, and hence to the optimal relaxation parameter.

But since there are three different matrices present here, and since we don't know the form of the linear system it's trying to solve, or how it uses the computed values in the solution, it's just impossible to say what it's doing.

Is this the correct procedure for calculating matrix spectrum?

1 Answers1