The coefficient matrices in M3D are reformed here to have symmetric structures. They are further categorized into 3 types: weak diagonally-dominant matrix, moderate diagonally-dominant matrix, and strong diagonally-dominant matrix. The weak diagonally-dominant matrix corresponds to the solution of auxiliary quantity F of the perturbed toroidal flux Ĩ with Neumann boundary conditions. The moderate diagonally-dominant matrix corresponds to the solution of the toroidal current C and the scalar potential Φ with Dirichlet boundary conditions. The strong diagonally-dominant matrix corresponds to the solution of the perturbed toroidal flux Ĩ, the poloidal flux ψ, the pressure p, and the 3 components of velocity: χ, U, and v ϕ with Dirichlet boundary conditions, respectively. We compare LU, GMRES, and ICCG algorithms for linear equation with these 3 types of matrices. It is observed that ICCG greatly accelerates the solution process as compared to GMRES, especially for the weak diagonally-dominant matrix. In this case we achieved 4 to 44 times speedup when the matrix order ranges from ∼10 1 to ∼10 4. For the moderate diagonally-dominant case, there is a 4 to 24 times speedup. For the strong diagonally-dominant case, an average 2 times speedup is observed. It is also shown that the GMRES algorithm is much slower for the weak diagonally-dominant type than for the other 2 types: 4.4 times slower than the moderate diagonally-dominant case, 333 times slower than the strong diagonally-dominant case. The LU algorithm is faster than GMRES for the weak diagonally-dominant matrix since the matrix is strong ill-conditioned, but GMRES outperforms LU for the strong diagonally-dominant matrix since the matrix is well-conditioned. However, the ICCG algorithm outperforms both ILU and GMRES for all matrix types.