This page was generated from unit-2.1.2-blockjacobi/blockjacobi.ipynb.

2.1.2 Building blocks for programming preconditioners

In \(\S\)2.1.1, we saw the Jacobi method given as a local preconditioner. We now delve deeper into such local preconditioners, which are often useful ingredients/smoothers in more complicated preconditioning strategies. In this tutorial we will see

• Jacobi and Gauss-Seidel smoothers/preconditioners,

• their block versions,

• how to select your own smoothing blocks, and

• how to additively combine preconditioners.

from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.la import EigenValues_Preconditioner

def Setup(h=0.1, p=3):
    mesh = Mesh(unit_square.GenerateMesh(maxh=h))
    fes = H1(mesh, order=p, dirichlet="left|bottom")
    u, v = fes.TnT()
    a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx)
    f = LinearForm(v*dx)
    gfu = GridFunction(fes)
    return mesh, fes, a, f, gfu

mesh, fes, a, f, gfu = Setup(h=0.1, p=3)
a.Assemble()
f.Assemble();

Create a Jacobi-preconditioner

Let \(A=\) a.mat be the assembled matrix, which can be decomposed based on FreeDofs (\(F\)) the remainder (\(D\)), as in \(\S\)1.3,

\[\begin{split}A = \left( \begin{array}{cc} A_{FF} & A_{FD} \\ A_{DF} & A_{DD} \end{array} \right).\end{split}\]

Then the matrix form of the point Jacobi preconditioner is

\[\begin{split}J = \left( \begin{array}{cc} \text{diag}(A_{FF})^{-1} & 0 \\ 0 & 0 \end{array} \right),\end{split}\]

which can be obtained in NGSolve using CreateSmoother:

preJpoint = a.mat.CreateSmoother(fes.FreeDofs())

NGSolve also gives us a facility to quickly check an estimate of the condition number of the preconditioned matrix by applying the Lanczos algorithm on the preconditioned system.

lams = EigenValues_Preconditioner(mat=a.mat, pre=preJpoint)
lams

 0.0145464
 0.0568008
 0.0922746
 0.111764
 0.141313
 0.192416
 0.243329
 0.344161
 0.455739
 0.522676
 0.639076
 0.755497
 0.831726
 0.917639
 1.02986
 1.15959
 1.24795
 1.34881
 1.45217
 1.56496
 1.69223
 1.80351
  1.9083
 2.01703
 2.10665
 2.20345
 2.28227
 2.36363
  2.4083
 2.46259
 2.47731
 2.51528
 2.80058

An estimate of the condition number \(\kappa = \lambda_{\text{max}} / \lambda_{\text{min}}\) is therefore given as follows:

max(lams)/min(lams)

192.52692231522082

One might wonder if we have gained anything by this point Jacobi preconditioning. What if we did not precondition at all?

Not preconditioning is the same as preconditioning by the identity operator on \(F\)-dofs. One way to realize this identity operator in NGSolve is through the projection into the space of free dofs (i.e., the space spanned by the shape functions corresponding to free dofs). NGSolve provides

Projector(mask, range)   # mask: bit array; range: bool

which projects into the space spanned by the shape functions of the degrees of freedom marked as range in mask.

preI = Projector(mask=fes.FreeDofs(), range=True)

Note that another way to obtain the identity matrix in NGSolve is

IdentityMatrix(fes.ndof, complex=False).

lams = EigenValues_Preconditioner(mat=a.mat, pre=preI)
max(lams)/min(lams)

1562.655376007982

Clearly the point Jacobi preconditioner has reduced the condition number.

We can use preconditioners within iterative solvers provided by NGSolve’s solvers module (which has MinRes, QMR etc.) Here is an illustration of its use within CG, or the conjugate gradient solver:

solvers.CG(mat=a.mat, pre=preJpoint, rhs=f.vec, sol=gfu.vec)
Draw(gfu)

CG iteration 1, residual = 0.05538055852845106
CG iteration 2, residual = 0.09315935397110503
CG iteration 3, residual = 0.10351189020825113
CG iteration 4, residual = 0.09055344263695937
CG iteration 5, residual = 0.09690655381638001
CG iteration 6, residual = 0.08685499565444357
CG iteration 7, residual = 0.09375145900717156
CG iteration 8, residual = 0.08102545067766499
CG iteration 9, residual = 0.05141273300638216
CG iteration 10, residual = 0.0431870326518925
CG iteration 11, residual = 0.03478485636733268
CG iteration 12, residual = 0.027955532180578252
CG iteration 13, residual = 0.021786413075808392
CG iteration 14, residual = 0.01956256250251828
CG iteration 15, residual = 0.013925054803092866
CG iteration 16, residual = 0.010717691595519545
CG iteration 17, residual = 0.009072032767387475
CG iteration 18, residual = 0.006853407280989169
CG iteration 19, residual = 0.004560661678349433
CG iteration 20, residual = 0.003188722005885534
CG iteration 21, residual = 0.001827684585575292
CG iteration 22, residual = 0.0012509916643804668
CG iteration 23, residual = 0.001025491414944295
CG iteration 24, residual = 0.0007220602502288801
CG iteration 25, residual = 0.0004749807632415533
CG iteration 26, residual = 0.0003512488512668743
CG iteration 27, residual = 0.0002911953880497042
CG iteration 28, residual = 0.0002182695359589094
CG iteration 29, residual = 0.00015584865841437626
CG iteration 30, residual = 0.00011137197511770437
CG iteration 31, residual = 7.923058395297498e-05
CG iteration 32, residual = 6.678061443918647e-05
CG iteration 33, residual = 5.4140486060236344e-05
CG iteration 34, residual = 4.3002947423831404e-05
CG iteration 35, residual = 3.217682856039787e-05
CG iteration 36, residual = 2.2537832752118193e-05
CG iteration 37, residual = 1.731251235250529e-05
CG iteration 38, residual = 1.4474897710791329e-05
CG iteration 39, residual = 1.1166971332630283e-05
CG iteration 40, residual = 7.78946322045134e-06
CG iteration 41, residual = 5.4300175642087405e-06
CG iteration 42, residual = 3.955667003062854e-06
CG iteration 43, residual = 2.744469036729557e-06
CG iteration 44, residual = 2.049745474309032e-06
CG iteration 45, residual = 1.3824160075542e-06
CG iteration 46, residual = 9.404899136462762e-07
CG iteration 47, residual = 6.193600569040406e-07
CG iteration 48, residual = 4.6248282363918463e-07
CG iteration 49, residual = 3.3644149196602156e-07
CG iteration 50, residual = 2.1125565850464987e-07
CG iteration 51, residual = 1.40776959336719e-07
CG iteration 52, residual = 1.0212975776385104e-07
CG iteration 53, residual = 7.829956407494647e-08
CG iteration 54, residual = 5.491378131315398e-08
CG iteration 55, residual = 3.911575811331202e-08
CG iteration 56, residual = 2.7095063137655022e-08
CG iteration 57, residual = 1.8715447617397395e-08
CG iteration 58, residual = 1.1445619024755785e-08
CG iteration 59, residual = 7.494637549995804e-09
CG iteration 60, residual = 5.2694613053402356e-09
CG iteration 61, residual = 3.944326777051754e-09
CG iteration 62, residual = 2.7013286089511052e-09
CG iteration 63, residual = 1.6319138253283868e-09
CG iteration 64, residual = 1.0500554674678948e-09
CG iteration 65, residual = 7.384473405359392e-10
CG iteration 66, residual = 5.838223572826365e-10
CG iteration 67, residual = 4.3745559364989083e-10
CG iteration 68, residual = 3.006836760858395e-10
CG iteration 69, residual = 1.9118334699242257e-10
CG iteration 70, residual = 1.088794487217495e-10
CG iteration 71, residual = 7.671367833741455e-11
CG iteration 72, residual = 5.314368906574211e-11
CG iteration 73, residual = 3.560578465222463e-11
CG iteration 74, residual = 2.3612549339269696e-11
CG iteration 75, residual = 1.5841966756569016e-11
CG iteration 76, residual = 1.0798034307016644e-11
CG iteration 77, residual = 7.143662364767825e-12
CG iteration 78, residual = 4.933643130631897e-12
CG iteration 79, residual = 3.447652275556076e-12
CG iteration 80, residual = 2.2686743245400916e-12
CG iteration 81, residual = 1.3026588591991905e-12
CG iteration 82, residual = 8.401833878076998e-13
CG iteration 83, residual = 5.917311690681452e-13
CG iteration 84, residual = 4.373819936341534e-13
CG iteration 85, residual = 2.9291689508133064e-13
CG iteration 86, residual = 1.7380642827410246e-13
CG iteration 87, residual = 1.0062790058575541e-13
CG iteration 88, residual = 7.175414172010894e-14
CG iteration 89, residual = 5.353935895756829e-14

BaseWebGuiScene

Gauss-Seidel smoothing

The same point Jacobi smoother object can also used to perform point Gauss-Seidel smoothing. One step of the classical Gauss-Seidel iteration is realized by the method preJpoint.Smooth(). It is well known that this iteration converges for matrices like \(A\). Below we show how to use it as a linear iterative solver.

gfu.vec[:] = 0
res = f.vec.CreateVector()              # residual
projres = f.vec.CreateVector()          # residual projected to freedofs
proj = Projector(fes.FreeDofs(), True)

for i in range(500):
    preJpoint.Smooth(gfu.vec, f.vec)    # one step of point Gauss-Seidel
    res.data = f.vec - a.mat*gfu.vec
    projres.data = proj * res
    print ("it#", i, ", res =", Norm(projres))
Draw (gfu)

it# 0 , res = 0.08158454491029796
it# 1 , res = 0.07701420259744884
it# 2 , res = 0.07344552937844082
it# 3 , res = 0.07022797813293108
it# 4 , res = 0.067335311402751
it# 5 , res = 0.06470423687226663
it# 6 , res = 0.06228529657252073
it# 7 , res = 0.060041797551721836
it# 8 , res = 0.05794597823795032
it# 9 , res = 0.05597643413016141
it# 10 , res = 0.05411638255690982
it# 11 , res = 0.05235246169022937
it# 12 , res = 0.05067388538239118
it# 13 , res = 0.04907183792095946
it# 14 , res = 0.04753903282470946
it# 15 , res = 0.0460693855286243
it# 16 , res = 0.044657766510979774
it# 17 , res = 0.04329981230624786
it# 18 , res = 0.04199177895713715
it# 19 , res = 0.04073042713003674
it# 20 , res = 0.039512931225209505
it# 21 , res = 0.03833680691588803
it# 22 , res = 0.03719985300247283
it# 23 , res = 0.03610010449189233
it# 24 , res = 0.03503579454882256
it# 25 , res = 0.0340053235053189
it# 26 , res = 0.03300723351756637
it# 27 , res = 0.03204018776229419
it# 28 , res = 0.0311029532977578
it# 29 , res = 0.030194386893718894
it# 30 , res = 0.029313423274774392
it# 31 , res = 0.02845906533125371
it# 32 , res = 0.027630375938709244
it# 33 , res = 0.026826471095984954
it# 34 , res = 0.02604651414687739
it# 35 , res = 0.02528971089449334
it# 36 , res = 0.024555305452852363
it# 37 , res = 0.02384257670886459
it# 38 , res = 0.023150835290920452
it# 39 , res = 0.022479420959054936
it# 40 , res = 0.021827700346859897
it# 41 , res = 0.02119506499769431
it# 42 , res = 0.020580929647832393
it# 43 , res = 0.019984730717435357
it# 44 , res = 0.01940592497697394
it# 45 , res = 0.018843988362255343
it# 46 , res = 0.01829841491574992
it# 47 , res = 0.017768715835639097
it# 48 , res = 0.017254418617082018
it# 49 , res = 0.01675506627273352
it# 50 , res = 0.016270216621642367
it# 51 , res = 0.01579944163739242
it# 52 , res = 0.015342326847788895
it# 53 , res = 0.01489847077958863
it# 54 , res = 0.01446748444276868
it# 55 , res = 0.014048990849655548
it# 56 , res = 0.013642624564934733
it# 57 , res = 0.013248031283138761
it# 58 , res = 0.01286486743070027
it# 59 , res = 0.01249279979006747
it# 60 , res = 0.012131505143723527
it# 61 , res = 0.011780669936243044
it# 62 , res = 0.011439989952767765
it# 63 , res = 0.011109170012491036
it# 64 , res = 0.010787923675918175
it# 65 , res = 0.010475972964825048
it# 66 , res = 0.010173048093962071
it# 67 , res = 0.00987888721366541
it# 68 , res = 0.009593236162631505
it# 69 , res = 0.009315848230190882
it# 70 , res = 0.009046483927492372
it# 71 , res = 0.008784910767064285
it# 72 , res = 0.008530903050277679
it# 73 , res = 0.00828424166227764
it# 74 , res = 0.008044713873993787
it# 75 , res = 0.007812113150872186
it# 76 , res = 0.007586238968003393
it# 77 , res = 0.00736689663134775
it# 78 , res = 0.007153897104784373
it# 79 , res = 0.006947056842728936
it# 80 , res = 0.006746197628084737
it# 81 , res = 0.006551146415310406
it# 82 , res = 0.0063617351783984
it# 83 , res = 0.006177800763575091
it# 84 , res = 0.00599918474654228
it# 85 , res = 0.005825733294094131
it# 86 , res = 0.005657297029949266
it# 87 , res = 0.005493730904648926
it# 88 , res = 0.005334894069379847
it# 89 , res = 0.005180649753586581
it# 90 , res = 0.005030865146245809
it# 91 , res = 0.004885411280680435
it# 92 , res = 0.0047441629227982975
it# 93 , res = 0.0046069984626418625
it# 94 , res = 0.004473799809145551
it# 95 , res = 0.00434445228799705
it# 96 , res = 0.004218844542504399
it# 97 , res = 0.00409686843737689
it# 98 , res = 0.003978418965327425
it# 99 , res = 0.0038633941564108123
it# 100 , res = 0.0037516949900136495
it# 101 , res = 0.003643225309414549
it# 102 , res = 0.0035378917388391405
it# 103 , res = 0.003435603602931567
it# 104 , res = 0.0033362728485721075
it# 105 , res = 0.0032398139689707923
it# 106 , res = 0.0031461439299674066
it# 107 , res = 0.003055182098474388
it# 108 , res = 0.002966850172997948
it# 109 , res = 0.002881072116177055
it# 110 , res = 0.0027977740892785813
it# 111 , res = 0.0027168843885944077
it# 112 , res = 0.002638333383681428
it# 113 , res = 0.0025620534573925764
it# 114 , res = 0.002487978947645425
it# 115 , res = 0.0024160460908778693
it# 116 , res = 0.0023461929671415595
it# 117 , res = 0.0022783594467852067
it# 118 , res = 0.002212487138681771
it# 119 , res = 0.0021485193399545647
it# 120 , res = 0.0020864009871577136
it# 121 , res = 0.002026078608870218
it# 122 , res = 0.001967500279661614
it# 123 , res = 0.00191061557538926
it# 124 , res = 0.0018553755297888953
it# 125 , res = 0.0018017325923216233
it# 126 , res = 0.0017496405872394634
it# 127 , res = 0.0016990546738353764
it# 128 , res = 0.0016499313078423227
it# 129 , res = 0.0016022282039501678
it# 130 , res = 0.001555904299404273
it# 131 , res = 0.0015109197186591984
it# 132 , res = 0.0014672357390530606
it# 133 , res = 0.0014248147574748356
it# 134 , res = 0.0013836202579954754
it# 135 , res = 0.0013436167804349976
it# 136 , res = 0.0013047698898390869
it# 137 , res = 0.0012670461468371225
it# 138 , res = 0.0012304130788577119
it# 139 , res = 0.001194839152177292
it# 140 , res = 0.001160293744775447
it# 141 , res = 0.0011267471199764146
it# 142 , res = 0.0010941704008518064
it# 143 , res = 0.0010625355453635344
it# 144 , res = 0.0010318153222260088
it# 145 , res = 0.0010019832874657193
it# 146 , res = 0.000973013761657863
it# 147 , res = 0.0009448818078227786
it# 148 , res = 0.0009175632099601754
it# 149 , res = 0.0008910344522036725
it# 150 , res = 0.0008652726985796375
it# 151 , res = 0.0008402557733492908
it# 152 , res = 0.0008159621419208075
it# 153 , res = 0.0007923708923120122
it# 154 , res = 0.0007694617171498034
it# 155 , res = 0.0007472148961896871
it# 156 , res = 0.0007256112793408188
it# 157 , res = 0.0007046322701816679
it# 158 , res = 0.0006842598099526338
it# 159 , res = 0.0006644763620104103
it# 160 , res = 0.0006452648967339671
it# 161 , res = 0.0006266088768638825
it# 162 , res = 0.0006084922432685012
it# 163 , res = 0.0005908994011198445
it# 164 , res = 0.0005738152064693205
it# 165 , res = 0.0005572249532125081
it# 166 , res = 0.0005411143604297035
it# 167 , res = 0.0005254695600933729
it# 168 , res = 0.0005102770851301749
it# 169 , res = 0.0004955238578287787
it# 170 , res = 0.00048119717858276917
it# 171 , res = 0.00046728471495857194
it# 172 , res = 0.00045377449107991815
it# 173 , res = 0.0004406548773187243
it# 174 , res = 0.00042791458028482014
it# 175 , res = 0.000415542633104162
it# 176 , res = 0.0004035283859787939
it# 177 , res = 0.0003918614970195488
it# 178 , res = 0.00038053192334361156
it# 179 , res = 0.00036952991242965267
it# 180 , res = 0.0003588459937228913
it# 181 , res = 0.00034847097048297987
it# 182 , res = 0.00033839591186708365
it# 183 , res = 0.00032861214524251245
it# 184 , res = 0.0003191112487213219
it# 185 , res = 0.0003098850439106512
it# 186 , res = 0.0003009255888726686
it# 187 , res = 0.00029222517128845305
it# 188 , res = 0.00028377630181871066
it# 189 , res = 0.00027557170765746605
it# 190 , res = 0.00026760432627143544
it# 191 , res = 0.0002598672993206221
it# 192 , res = 0.00025235396675431525
it# 193 , res = 0.0002450578610785725
it# 194 , res = 0.00023797270178899424
it# 195 , res = 0.0002310923899639278
it# 196 , res = 0.00022441100301515767
it# 197 , res = 0.0002179227895891631
it# 198 , res = 0.0002116221646168121
it# 199 , res = 0.00020550370450540654
it# 200 , res = 0.00019956214247014324
it# 201 , res = 0.00019379236400073367
it# 202 , res = 0.00018818940245816907
it# 203 , res = 0.0001827484348000783
it# 204 , res = 0.0001774647774287972
it# 205 , res = 0.00017233388215949966
it# 206 , res = 0.00016735133230563656
it# 207 , res = 0.0001625128388765407
it# 208 , res = 0.000157814236885419
it# 209 , res = 0.00015325148176491962
it# 210 , res = 0.00014882064588451184
it# 211 , res = 0.00014451791517014722
it# 212 , res = 0.0001403395858210718
it# 213 , res = 0.00013628206112149068
it# 214 , res = 0.00013234184834463088
it# 215 , res = 0.00012851555574587692
it# 216 , res = 0.0001247998896437193
it# 217 , res = 0.00012119165158360017
it# 218 , res = 0.00011768773558597085
it# 219 , res = 0.00011428512547146749
it# 220 , res = 0.00011098089226527009
it# 221 , res = 0.00010777219167547494
it# 222 , res = 0.00010465626164515458
it# 223 , res = 0.00010163041997422051
it# 224 , res = 9.869206201096987e-05
it# 225 , res = 9.583865840998952e-05
it# 226 , res = 9.3067752954595e-05
it# 227 , res = 9.037696044276308e-05
it# 228 , res = 8.776396463388748e-05
it# 229 , res = 8.522651625499745e-05
it# 230 , res = 8.276243106458327e-05
it# 231 , res = 8.036958797203013e-05
it# 232 , res = 7.804592721245146e-05
it# 233 , res = 7.578944857315928e-05
it# 234 , res = 7.359820967186025e-05
it# 235 , res = 7.147032428494175e-05
it# 236 , res = 6.940396072362522e-05
it# 237 , res = 6.739734025713993e-05
it# 238 , res = 6.544873558196441e-05
it# 239 , res = 6.355646933430308e-05
it# 240 , res = 6.171891264706871e-05
it# 241 , res = 5.9934483747105876e-05
it# 242 , res = 5.820164659365177e-05
it# 243 , res = 5.651890955629979e-05
it# 244 , res = 5.488482413117977e-05
it# 245 , res = 5.329798369346851e-05
it# 246 , res = 5.175702228727788e-05
it# 247 , res = 5.026061344943178e-05
it# 248 , res = 4.880746906739233e-05
it# 249 , res = 4.739633827119828e-05
it# 250 , res = 4.602600635589202e-05
it# 251 , res = 4.469529373676427e-05
it# 252 , res = 4.3403054933081085e-05
it# 253 , res = 4.214817758266315e-05
it# 254 , res = 4.092958148389811e-05
it# 255 , res = 3.9746217666637415e-05
it# 256 , res = 3.8597067488216773e-05
it# 257 , res = 3.748114175753226e-05
it# 258 , res = 3.639747988322248e-05
it# 259 , res = 3.5345149046246083e-05
it# 260 , res = 3.432324339801092e-05
it# 261 , res = 3.333088327949151e-05
it# 262 , res = 3.236721446476762e-05
it# 263 , res = 3.143140742553918e-05
it# 264 , res = 3.0522656616784636e-05
it# 265 , res = 2.9640179783605735e-05
it# 266 , res = 2.878321728785916e-05
it# 267 , res = 2.7951031454392612e-05
it# 268 , res = 2.7142905935279094e-05
it# 269 , res = 2.6358145094471567e-05
it# 270 , res = 2.5596073407678055e-05
it# 271 , res = 2.4856034882048463e-05
it# 272 , res = 2.4137392490481393e-05
it# 273 , res = 2.3439527623947663e-05
it# 274 , res = 2.276183955874233e-05
it# 275 , res = 2.2103744939333284e-05
it# 276 , res = 2.146467727623399e-05
it# 277 , res = 2.084408645848844e-05
it# 278 , res = 2.0241438280177926e-05
it# 279 , res = 1.965621398017024e-05
it# 280 , res = 1.90879097961431e-05
it# 281 , res = 1.85360365303218e-05
it# 282 , res = 1.800011912895102e-05
it# 283 , res = 1.7479696273088705e-05
it# 284 , res = 1.697431998150314e-05
it# 285 , res = 1.6483555225008542e-05
it# 286 , res = 1.6006979552068665e-05
it# 287 , res = 1.5544182725429892e-05
it# 288 , res = 1.5094766368262612e-05
it# 289 , res = 1.4658343622137905e-05
it# 290 , res = 1.423453881305167e-05
it# 291 , res = 1.3822987128936247e-05
it# 292 , res = 1.342333430518749e-05
it# 293 , res = 1.3035236319652174e-05
it# 294 , res = 1.265835909651356e-05
it# 295 , res = 1.2292378219103142e-05
it# 296 , res = 1.1936978650294078e-05
it# 297 , res = 1.1591854461068472e-05
it# 298 , res = 1.1256708567819e-05
it# 299 , res = 1.0931252476223177e-05
it# 300 , res = 1.0615206032788631e-05
it# 301 , res = 1.0308297183967661e-05
it# 302 , res = 1.0010261742061279e-05
it# 303 , res = 9.720843157272578e-06
it# 304 , res = 9.439792297552055e-06
it# 305 , res = 9.166867233401467e-06
it# 306 , res = 8.90183303054385e-06
it# 307 , res = 8.644461546710797e-06
it# 308 , res = 8.394531236113314e-06
it# 309 , res = 8.151826958042768e-06
it# 310 , res = 7.91613979187625e-06
it# 311 , res = 7.687266857923437e-06
it# 312 , res = 7.4650111414826734e-06
it# 313 , res = 7.2491813245328434e-06
it# 314 , res = 7.0395916201008e-06
it# 315 , res = 6.83606161293897e-06
it# 316 , res = 6.638416103887351e-06
it# 317 , res = 6.446484959347123e-06
it# 318 , res = 6.260102964380015e-06
it# 319 , res = 6.079109680926496e-06
it# 320 , res = 5.903349309692706e-06
it# 321 , res = 5.732670555553196e-06
it# 322 , res = 5.566926497892747e-06
it# 323 , res = 5.405974463960224e-06
it# 324 , res = 5.24967590561079e-06
it# 325 , res = 5.0978962808264416e-06
it# 326 , res = 4.950504937252275e-06
it# 327 , res = 4.807375000306979e-06
it# 328 , res = 4.668383263117248e-06
it# 329 , res = 4.533410081375718e-06
it# 330 , res = 4.402339269854127e-06
it# 331 , res = 4.275058002566916e-06
it# 332 , res = 4.151456715372918e-06
it# 333 , res = 4.031429011861027e-06
it# 334 , res = 3.914871572243141e-06
it# 335 , res = 3.801684063362447e-06
it# 336 , res = 3.6917690530804396e-06
it# 337 , res = 3.585031926463136e-06
it# 338 , res = 3.4813808039547113e-06
it# 339 , res = 3.3807264623897484e-06
it# 340 , res = 3.2829822581423883e-06
it# 341 , res = 3.1880640528380667e-06
it# 342 , res = 3.095890140693216e-06
it# 343 , res = 3.006381178191567e-06
it# 344 , res = 2.919460115749521e-06
it# 345 , res = 2.8350521315287564e-06
it# 346 , res = 2.753084567002455e-06
it# 347 , res = 2.6734868641370716e-06
it# 348 , res = 2.5961905052304734e-06
it# 349 , res = 2.5211289532966214e-06
it# 350 , res = 2.4482375951586986e-06
it# 351 , res = 2.3774536857566343e-06
it# 352 , res = 2.3087162941635875e-06
it# 353 , res = 2.241966251172684e-06
it# 354 , res = 2.177146097960148e-06
it# 355 , res = 2.1142000373417927e-06
it# 356 , res = 2.0530738851072737e-06
it# 357 , res = 1.9937150237696942e-06
it# 358 , res = 1.936072357181074e-06
it# 359 , res = 1.8800962662495933e-06
it# 360 , res = 1.8257385667085815e-06
it# 361 , res = 1.7729524671792297e-06
it# 362 , res = 1.7216925294828362e-06
it# 363 , res = 1.6719146286956745e-06
it# 364 , res = 1.6235759160513422e-06
it# 365 , res = 1.5766347814217245e-06
it# 366 , res = 1.531050817848509e-06
it# 367 , res = 1.4867847864129542e-06
it# 368 , res = 1.4437985828425085e-06
it# 369 , res = 1.4020552046501457e-06
it# 370 , res = 1.3615187187343345e-06
it# 371 , res = 1.3221542315604774e-06
it# 372 , res = 1.2839278577939486e-06
it# 373 , res = 1.246806692270947e-06
it# 374 , res = 1.210758780938851e-06
it# 375 , res = 1.1757530934774583e-06
it# 376 , res = 1.1417594971409268e-06
it# 377 , res = 1.1087487303012274e-06
it# 378 , res = 1.0766923768355604e-06
it# 379 , res = 1.045562842647791e-06
it# 380 , res = 1.0153333315852946e-06
it# 381 , res = 9.859778219439997e-07
it# 382 , res = 9.574710442775193e-07
it# 383 , res = 9.297884599455395e-07
it# 384 , res = 9.029062397952039e-07
it# 385 , res = 8.768012434373425e-07
it# 386 , res = 8.514509997174116e-07
it# 387 , res = 8.268336869995626e-07
it# 388 , res = 8.029281147706311e-07
it# 389 , res = 7.797137048987637e-07
it# 390 , res = 7.571704744975274e-07
it# 391 , res = 7.35279018228663e-07
it# 392 , res = 7.140204918689466e-07
it# 393 , res = 6.93376596058694e-07
it# 394 , res = 6.733295604632773e-07
it# 395 , res = 6.53862128650529e-07
it# 396 , res = 6.349575429155277e-07
it# 397 , res = 6.165995300271977e-07
it# 398 , res = 5.98772287532169e-07
it# 399 , res = 5.814604697391092e-07
it# 400 , res = 5.646491744036758e-07
it# 401 , res = 5.483239304642691e-07
it# 402 , res = 5.32470684947866e-07
it# 403 , res = 5.170757916970794e-07
it# 404 , res = 5.021259983958332e-07
it# 405 , res = 4.876084362859546e-07
it# 406 , res = 4.735106088191733e-07
it# 407 , res = 4.5982038026472787e-07
it# 408 , res = 4.46525966382994e-07
it# 409 , res = 4.336159230082819e-07
it# 410 , res = 4.2107913725664225e-07
it# 411 , res = 4.089048174050019e-07
it# 412 , res = 3.970824838286707e-07
it# 413 , res = 3.8560195978840776e-07
it# 414 , res = 3.744533629083909e-07
it# 415 , res = 3.636270962427332e-07
it# 416 , res = 3.5311384078062575e-07
it# 417 , res = 3.4290454641639865e-07
it# 418 , res = 3.329904252073484e-07
it# 419 , res = 3.2336294296469057e-07
it# 420 , res = 3.1401381226269816e-07
it# 421 , res = 3.0493498533732626e-07
it# 422 , res = 2.961186472202231e-07
it# 423 , res = 2.8755720881975843e-07
it# 424 , res = 2.792433002902755e-07
it# 425 , res = 2.7116976503668966e-07
it# 426 , res = 2.6332965340875473e-07
it# 427 , res = 2.5571621659478066e-07
it# 428 , res = 2.483229008691228e-07
it# 429 , res = 2.41143342016248e-07
it# 430 , res = 2.3417136005705055e-07
it# 431 , res = 2.2740095330561052e-07
it# 432 , res = 2.2082629389733077e-07
it# 433 , res = 2.1444172224571714e-07
it# 434 , res = 2.082417423980099e-07
it# 435 , res = 2.0222101777280838e-07
it# 436 , res = 1.9637436540555348e-07
it# 437 , res = 1.9069675252856704e-07
it# 438 , res = 1.8518329187760713e-07
it# 439 , res = 1.7982923739096188e-07
it# 440 , res = 1.7462998044525073e-07
it# 441 , res = 1.6958104531865657e-07
it# 442 , res = 1.6467808595555063e-07
it# 443 , res = 1.5991688200268945e-07
it# 444 , res = 1.552933347339495e-07
it# 445 , res = 1.5080346438472974e-07
it# 446 , res = 1.4644340606054435e-07
it# 447 , res = 1.4220940649303992e-07
it# 448 , res = 1.3809782129503953e-07
it# 449 , res = 1.3410511089589528e-07
it# 450 , res = 1.3022783844615977e-07
it# 451 , res = 1.264626664611372e-07
it# 452 , res = 1.2280635396912627e-07
it# 453 , res = 1.192557534197149e-07
it# 454 , res = 1.1580780843191293e-07
it# 455 , res = 1.1245955113747845e-07
it# 456 , res = 1.0920809921238367e-07
it# 457 , res = 1.060506539613815e-07
it# 458 , res = 1.0298449740153936e-07
it# 459 , res = 1.0000699016813586e-07
it# 460 , res = 9.711556908401632e-08
it# 461 , res = 9.430774531725115e-08
it# 462 , res = 9.158110191286881e-08
it# 463 , res = 8.893329167061978e-08
it# 464 , res = 8.636203554906392e-08
it# 465 , res = 8.386511997738113e-08
it# 466 , res = 8.144039579703933e-08
it# 467 , res = 7.908577556887003e-08
it# 468 , res = 7.679923262836436e-08
it# 469 , res = 7.457879873058425e-08
it# 470 , res = 7.242256224946768e-08
it# 471 , res = 7.032866741207166e-08
it# 472 , res = 6.829531169007436e-08
it# 473 , res = 6.63207447933448e-08
it# 474 , res = 6.440326680450478e-08
it# 475 , res = 6.254122729288852e-08
it# 476 , res = 6.073302353275807e-08
it# 477 , res = 5.8977098741561216e-08
it# 478 , res = 5.7271941736539335e-08
it# 479 , res = 5.561608445499275e-08
it# 480 , res = 5.400810172947907e-08
it# 481 , res = 5.244660922702858e-08
it# 482 , res = 5.0930262945152884e-08
it# 483 , res = 4.945775745474459e-08
it# 484 , res = 4.802782550385657e-08
it# 485 , res = 4.663923593991956e-08
it# 486 , res = 4.5290793429019615e-08
it# 487 , res = 4.398133744510376e-08
it# 488 , res = 4.270974068475564e-08
it# 489 , res = 4.147490854888346e-08
it# 490 , res = 4.027577816718289e-08
it# 491 , res = 3.91113172788408e-08
it# 492 , res = 3.7980523413718034e-08
it# 493 , res = 3.6882423346904056e-08
it# 494 , res = 3.581607173244975e-08
it# 495 , res = 3.478055071392342e-08
it# 496 , res = 3.377496886965999e-08
it# 497 , res = 3.2798460497119535e-08
it# 498 , res = 3.185018523784621e-08
it# 499 , res = 3.092932658343176e-08

BaseWebGuiScene

Implement a forward-backward GS preconditioner

The same point Jacobi smoother object is also able to perform a Gauss-Seidel iteration after reversing the ordering of the points, i.e., a backward Gauss-Seidel sweep. One can combine the forward and backward sweeps to construct a symmetric preconditioner, often called the symmetrized Gauss-Seidel preconditioner. This offers a good simple illustration of how to construct NGSolve preconditioners from within python.

class SymmetricGS(BaseMatrix):
    def __init__ (self, smoother):
        super(SymmetricGS, self).__init__()
        self.smoother = smoother
    def Mult (self, x, y):
        y[:] = 0.0
        self.smoother.Smooth(y, x)
        self.smoother.SmoothBack(y,x)
    def Height (self):
        return self.smoother.height
    def Width (self):
        return self.smoother.height

preGS = SymmetricGS(preJpoint)
solvers.CG(mat=a.mat, pre=preGS, rhs=f.vec, sol=gfu.vec)
Draw(gfu)

CG iteration 1, residual = 0.09431736772142078
CG iteration 2, residual = 0.12202935492847881
CG iteration 3, residual = 0.08635033881376865
CG iteration 4, residual = 0.054724963439555664
CG iteration 5, residual = 0.02653117244015519
CG iteration 6, residual = 0.014817837346444674
CG iteration 7, residual = 0.006568692977315905
CG iteration 8, residual = 0.0028179384297480617
CG iteration 9, residual = 0.000975490712281912
CG iteration 10, residual = 0.00038788789222606035
CG iteration 11, residual = 0.00012801282895042707
CG iteration 12, residual = 6.138184244284248e-05
CG iteration 13, residual = 2.9576601747088142e-05
CG iteration 14, residual = 1.352015176068953e-05
CG iteration 15, residual = 5.668185475388735e-06
CG iteration 16, residual = 2.4123963208265913e-06
CG iteration 17, residual = 8.67459981293994e-07
CG iteration 18, residual = 4.017571016763011e-07
CG iteration 19, residual = 1.3703872761080615e-07
CG iteration 20, residual = 5.4896914357743264e-08
CG iteration 21, residual = 2.0516648155353695e-08
CG iteration 22, residual = 6.144069986364002e-09
CG iteration 23, residual = 2.7864832241461385e-09
CG iteration 24, residual = 8.307004904433915e-10
CG iteration 25, residual = 3.3016506758272495e-10
CG iteration 26, residual = 1.2444445435727118e-10
CG iteration 27, residual = 4.1685221594516925e-11
CG iteration 28, residual = 1.855470440241164e-11
CG iteration 29, residual = 5.3869614737237855e-12
CG iteration 30, residual = 1.834661509401809e-12
CG iteration 31, residual = 6.73909315211797e-13
CG iteration 32, residual = 2.004237186203173e-13
CG iteration 33, residual = 6.690135963709771e-14

BaseWebGuiScene

lams = EigenValues_Preconditioner(mat=a.mat, pre=preGS)
max(lams)/min(lams)

20.242155787942128

Note that the condition number now is better than that of the system preconditioned by point Jacobi.

A Block Jacobi preconditioner

The point Jacobi preconditioner is based on inverses of 1 x 1 diagonal blocks. Condition numbers can be improved by using larger blocks. It is possible to group dofs into blocks within python and construct an NGSolve preconditioner based on the blocks.

Here is an example that constructs vertex-based blocks, using the mesh querying techniques we learnt in \(\S\)1.8.

def VertexPatchBlocks(mesh, fes):
    blocks = []
    freedofs = fes.FreeDofs()
    for v in mesh.vertices:
        vdofs = set()
        for el in mesh[v].elements:
            vdofs |= set(d for d in fes.GetDofNrs(el)
                         if freedofs[d])
        blocks.append(vdofs)
    return blocks

blocks = VertexPatchBlocks(mesh, fes)
print(blocks)

[{160, 873, 159}, {13, 207, 48, 145, 146, 208, 143, 211, 212, 144, 889, 891}, {257, 2, 258, 908, 147, 148, 21, 22, 149, 150}, {65, 312, 924, 315, 926, 316, 311, 152, 151, 155, 156, 30}, {160, 161, 162, 165, 166, 40, 873, 874, 361, 362, 942, 159}, {161, 162, 165, 166, 167, 40, 168, 874, 41, 171, 172, 875, 363, 364, 876}, {167, 168, 41, 42, 875, 171, 172, 173, 174, 877, 177, 178, 878, 369, 370}, {375, 376, 42, 43, 877, 173, 174, 879, 177, 178, 179, 180, 880, 183, 184}, {382, 43, 44, 879, 881, 882, 179, 180, 381, 183, 184, 185, 186, 189, 190}, {192, 195, 196, 387, 388, 44, 45, 881, 883, 884, 185, 186, 189, 190, 191}, {192, 195, 196, 197, 198, 201, 202, 393, 394, 45, 46, 883, 885, 886, 191}, {197, 198, 201, 202, 203, 204, 205, 206, 399, 400, 46, 47, 885, 887, 888}, {203, 204, 205, 206, 207, 208, 145, 146, 405, 406, 47, 48, 887, 889, 890}, {13, 14, 143, 144, 145, 146, 211, 212, 209, 210, 213, 214, 217, 218, 411, 412, 48, 49, 891, 892, 893}, {13, 14, 15, 209, 210, 213, 214, 215, 216, 217, 218, 219, 220, 223, 224, 416, 415, 49, 50, 892, 894, 895}, {896, 897, 14, 15, 16, 215, 216, 219, 220, 221, 222, 223, 224, 225, 226, 229, 230, 421, 422, 50, 51, 894}, {896, 898, 899, 15, 16, 17, 221, 222, 225, 226, 227, 228, 229, 230, 231, 232, 235, 236, 427, 428, 51, 52}, {898, 900, 901, 16, 17, 18, 227, 228, 231, 232, 233, 234, 235, 236, 237, 238, 241, 242, 433, 52, 53, 434}, {900, 902, 903, 17, 18, 19, 439, 440, 233, 234, 237, 238, 239, 240, 241, 242, 243, 244, 53, 54, 247, 248}, {446, 902, 904, 905, 18, 19, 20, 55, 249, 239, 240, 243, 244, 245, 54, 247, 248, 246, 250, 253, 254, 445}, {256, 259, 260, 451, 452, 904, 906, 907, 19, 20, 21, 245, 246, 55, 56, 249, 250, 251, 252, 253, 254, 255}, {256, 257, 2, 259, 260, 258, 263, 264, 906, 908, 147, 20, 21, 148, 22, 149, 150, 943, 56, 251, 252, 255}, {257, 2, 258, 259, 261, 262, 260, 263, 265, 266, 264, 908, 269, 270, 909, 457, 458, 147, 148, 21, 22, 149, 150, 23, 943, 944, 56, 57}, {261, 262, 265, 266, 267, 268, 269, 270, 909, 910, 271, 272, 275, 276, 461, 22, 23, 24, 462, 911, 57, 58}, {267, 268, 910, 271, 272, 912, 273, 275, 276, 274, 277, 23, 24, 281, 278, 25, 282, 465, 466, 913, 58, 59}, {471, 472, 912, 273, 274, 914, 915, 277, 278, 279, 24, 281, 25, 282, 284, 280, 26, 283, 288, 287, 59, 60}, {914, 916, 917, 279, 280, 25, 26, 283, 284, 285, 27, 287, 288, 289, 290, 286, 477, 293, 294, 478, 60, 61}, {916, 918, 919, 26, 27, 28, 285, 286, 289, 290, 291, 292, 293, 294, 295, 296, 483, 484, 299, 300, 61, 62}, {918, 920, 921, 27, 28, 29, 291, 292, 295, 296, 297, 298, 299, 300, 301, 302, 489, 490, 305, 306, 62, 63}, {64, 920, 922, 923, 28, 29, 30, 297, 298, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 495, 496, 63}, {64, 65, 151, 152, 922, 155, 924, 29, 30, 156, 925, 303, 304, 307, 308, 309, 310, 311, 312, 501, 502}, {65, 321, 66, 322, 155, 156, 926, 927, 928, 508, 507, 315, 316, 317, 318}, {512, 321, 66, 322, 323, 67, 324, 327, 328, 927, 929, 930, 317, 318, 511}, {323, 67, 324, 68, 327, 328, 329, 330, 517, 518, 333, 334, 929, 931, 932}, {68, 69, 329, 330, 523, 524, 333, 334, 335, 336, 339, 340, 931, 933, 934}, {69, 70, 335, 336, 529, 530, 339, 340, 341, 342, 345, 346, 933, 935, 936}, {70, 71, 341, 342, 535, 536, 345, 346, 347, 348, 351, 352, 935, 937, 938}, {71, 72, 347, 348, 541, 542, 351, 352, 353, 354, 359, 360, 937, 939, 940}, {72, 353, 354, 357, 358, 359, 360, 40, 361, 939, 362, 941, 365, 366, 945}, {160, 161, 162, 357, 358, 40, 873, 361, 362, 941, 942, 159}, {159, 160, 161, 162, 547, 548, 165, 166, 167, 40, 41, 168, 941, 942, 945, 947, 72, 74, 357, 358, 359, 360, 361, 874, 363, 364, 876, 362, 365, 366, 367, 368, 1006, 373, 374}, {165, 166, 167, 168, 41, 40, 171, 172, 42, 173, 174, 551, 552, 946, 947, 948, 73, 74, 875, 363, 364, 876, 367, 878, 369, 370, 371, 372, 368, 373, 374, 377, 378}, {385, 386, 41, 42, 43, 171, 172, 173, 174, 553, 177, 178, 179, 180, 946, 73, 100, 554, 1001, 1002, 877, 878, 880, 369, 370, 371, 372, 375, 376, 377, 378, 379, 380}, {384, 385, 386, 389, 390, 42, 43, 44, 177, 178, 179, 180, 949, 565, 183, 184, 185, 186, 566, 75, 100, 1001, 1003, 879, 880, 882, 375, 376, 379, 380, 381, 382, 383}, {384, 387, 388, 389, 390, 391, 392, 395, 396, 43, 44, 45, 563, 564, 949, 950, 183, 184, 185, 186, 951, 189, 190, 191, 192, 75, 76, 881, 882, 884, 381, 382, 383}, {387, 388, 391, 392, 393, 394, 395, 396, 397, 398, 401, 402, 44, 45, 46, 950, 952, 953, 571, 572, 189, 190, 191, 192, 195, 196, 197, 198, 76, 77, 883, 884, 886}, {393, 394, 397, 398, 399, 400, 401, 402, 403, 404, 407, 408, 45, 46, 47, 952, 954, 955, 577, 578, 195, 196, 197, 198, 201, 202, 203, 204, 77, 78, 885, 886, 888}, {399, 400, 403, 404, 405, 406, 407, 408, 409, 410, 413, 414, 46, 47, 48, 954, 956, 957, 583, 584, 201, 202, 203, 204, 205, 206, 79, 207, 208, 78, 887, 888, 890}, {13, 143, 144, 145, 146, 405, 406, 409, 410, 411, 412, 413, 414, 417, 418, 47, 48, 49, 956, 958, 205, 206, 79, 207, 208, 211, 212, 213, 214, 889, 890, 891, 893}, {13, 14, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 425, 426, 48, 49, 50, 958, 960, 961, 589, 590, 79, 209, 210, 211, 212, 213, 214, 81, 217, 218, 219, 220, 892, 893, 895}, {897, 14, 15, 415, 416, 419, 420, 421, 422, 423, 424, 425, 426, 429, 430, 49, 50, 51, 959, 960, 962, 591, 80, 81, 592, 215, 216, 217, 218, 219, 220, 223, 224, 225, 226, 894, 895}, {896, 897, 899, 15, 16, 421, 422, 423, 424, 427, 428, 429, 430, 431, 432, 50, 51, 52, 435, 436, 959, 963, 964, 80, 593, 82, 594, 221, 222, 223, 224, 225, 226, 229, 230, 231, 232}, {898, 899, 901, 16, 17, 427, 428, 431, 432, 433, 434, 51, 52, 53, 435, 436, 438, 437, 441, 442, 963, 965, 966, 82, 83, 601, 602, 227, 228, 229, 230, 231, 232, 235, 236, 237, 238}, {900, 901, 903, 17, 18, 433, 434, 52, 53, 54, 439, 440, 441, 442, 438, 437, 443, 444, 447, 448, 965, 967, 968, 83, 84, 607, 608, 233, 234, 235, 236, 237, 238, 241, 242, 243, 244}, {902, 903, 905, 18, 19, 53, 54, 55, 439, 440, 443, 444, 445, 446, 447, 448, 449, 450, 453, 454, 967, 969, 970, 84, 85, 613, 614, 239, 240, 241, 242, 243, 244, 247, 248, 249, 250}, {256, 904, 905, 907, 19, 20, 54, 55, 56, 445, 446, 449, 450, 451, 452, 453, 454, 455, 456, 969, 459, 460, 971, 972, 85, 86, 619, 620, 245, 246, 247, 248, 249, 250, 253, 254, 255}, {256, 257, 258, 259, 260, 263, 264, 265, 906, 907, 266, 20, 21, 22, 943, 944, 55, 56, 57, 451, 452, 455, 456, 457, 458, 459, 460, 971, 973, 463, 464, 86, 251, 252, 253, 254, 255}, {261, 262, 263, 264, 265, 266, 269, 270, 909, 271, 272, 911, 22, 23, 944, 56, 57, 58, 457, 458, 459, 460, 461, 462, 463, 464, 973, 977, 467, 468, 86}, {267, 268, 269, 910, 271, 272, 270, 911, 275, 276, 277, 278, 23, 24, 913, 57, 58, 59, 461, 462, 975, 463, 465, 466, 464, 977, 469, 470, 467, 88, 468, 86, 475, 476, 1004, 625, 626}, {912, 273, 274, 275, 276, 277, 278, 913, 24, 281, 25, 282, 283, 284, 58, 59, 60, 974, 975, 976, 465, 466, 469, 470, 471, 472, 473, 474, 87, 475, 88, 476, 479, 480, 627, 628, 915}, {914, 915, 917, 279, 280, 25, 26, 283, 284, 281, 282, 287, 288, 289, 290, 59, 60, 61, 974, 978, 979, 471, 472, 473, 474, 87, 89, 477, 478, 479, 480, 481, 482, 485, 486, 629, 630}, {640, 916, 917, 919, 26, 27, 285, 286, 287, 288, 289, 290, 293, 294, 295, 296, 60, 61, 62, 978, 980, 981, 89, 90, 477, 478, 481, 482, 483, 484, 485, 486, 487, 488, 491, 492, 639}, {645, 646, 918, 919, 921, 27, 28, 291, 292, 293, 294, 295, 296, 299, 300, 301, 302, 61, 62, 63, 980, 982, 983, 90, 91, 483, 484, 487, 488, 489, 490, 491, 492, 493, 494, 497, 498}, {651, 652, 920, 921, 923, 28, 29, 297, 298, 299, 300, 301, 302, 305, 306, 307, 308, 62, 63, 64, 982, 984, 985, 91, 92, 489, 490, 493, 494, 495, 496, 497, 498, 499, 500, 503, 504}, {657, 658, 922, 923, 29, 30, 925, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 63, 64, 65, 984, 986, 987, 92, 93, 495, 496, 499, 500, 501, 502, 503, 504, 505, 506, 509, 510}, {513, 514, 151, 152, 155, 924, 156, 30, 925, 926, 928, 309, 310, 311, 312, 315, 316, 317, 318, 64, 65, 66, 986, 988, 93, 501, 502, 505, 506, 507, 508, 509, 510}, {512, 513, 514, 515, 516, 521, 522, 663, 664, 927, 928, 930, 315, 316, 317, 318, 321, 66, 322, 65, 67, 323, 324, 988, 93, 990, 95, 991, 507, 508, 509, 510, 511}, {512, 515, 516, 517, 518, 519, 520, 521, 522, 525, 526, 665, 666, 929, 930, 932, 321, 66, 323, 67, 324, 322, 327, 328, 68, 329, 330, 989, 94, 990, 95, 992, 511}, {517, 518, 519, 520, 523, 524, 525, 526, 527, 528, 531, 532, 667, 668, 931, 932, 934, 67, 68, 69, 327, 328, 329, 330, 333, 334, 335, 336, 989, 94, 96, 993, 994}, {523, 524, 527, 528, 529, 530, 531, 532, 533, 534, 537, 538, 675, 676, 933, 934, 936, 68, 69, 70, 333, 334, 335, 336, 339, 340, 341, 342, 96, 993, 97, 995, 996}, {529, 530, 533, 534, 535, 536, 537, 538, 539, 540, 543, 544, 935, 936, 681, 938, 682, 69, 70, 71, 339, 340, 341, 342, 345, 346, 347, 348, 97, 98, 995, 997, 998}, {535, 536, 539, 540, 541, 542, 543, 544, 545, 546, 549, 550, 937, 938, 940, 687, 688, 70, 71, 72, 345, 346, 347, 348, 351, 352, 353, 354, 98, 99, 997, 999, 1000}, {541, 542, 545, 546, 547, 548, 549, 550, 40, 939, 940, 559, 560, 945, 71, 72, 74, 351, 352, 353, 354, 99, 357, 358, 359, 360, 999, 365, 366, 1005, 367, 368, 1006}, {1027, 1046, 1048, 551, 552, 41, 42, 553, 554, 555, 556, 557, 558, 561, 562, 946, 948, 699, 700, 73, 74, 119, 100, 1002, 111, 369, 370, 371, 372, 373, 374, 759, 760, 377, 378, 379, 380}, {1026, 1027, 547, 548, 549, 550, 551, 40, 41, 552, 555, 556, 559, 560, 561, 562, 947, 948, 693, 694, 72, 73, 74, 99, 363, 364, 365, 1006, 1005, 367, 368, 366, 371, 372, 373, 374, 111}, {384, 385, 386, 389, 390, 391, 392, 1043, 1044, 43, 44, 563, 564, 949, 565, 567, 568, 951, 566, 569, 570, 697, 698, 575, 576, 709, 710, 75, 76, 100, 102, 1003, 1008, 118, 381, 382, 383}, {387, 388, 389, 390, 391, 392, 395, 396, 397, 398, 44, 45, 563, 564, 950, 951, 567, 953, 568, 571, 572, 573, 574, 575, 576, 701, 702, 579, 580, 75, 76, 77, 101, 102, 1007, 1008, 1009}, {1025, 393, 394, 395, 396, 397, 398, 401, 402, 403, 404, 1042, 45, 46, 952, 953, 955, 571, 572, 573, 574, 704, 577, 578, 579, 580, 581, 582, 703, 587, 76, 77, 78, 588, 101, 103, 1007}, {1025, 399, 400, 401, 402, 403, 404, 407, 408, 409, 410, 46, 47, 954, 955, 957, 577, 578, 581, 582, 583, 584, 585, 586, 587, 588, 77, 78, 79, 589, 590, 81, 599, 600, 103, 1010, 1011}, {405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 417, 418, 419, 420, 47, 48, 49, 956, 957, 958, 961, 583, 584, 585, 586, 589, 78, 79, 590, 81, 1010}, {1045, 1052, 421, 422, 423, 424, 425, 426, 429, 430, 431, 432, 50, 51, 959, 962, 964, 715, 716, 591, 80, 81, 592, 593, 82, 594, 597, 598, 595, 596, 599, 600, 605, 606, 103, 106, 1015}, {1045, 415, 416, 417, 418, 419, 420, 423, 424, 425, 426, 49, 50, 960, 961, 962, 583, 584, 585, 586, 587, 588, 589, 590, 79, 591, 81, 80, 592, 78, 595, 596, 599, 600, 103, 1010, 1011}, {427, 428, 429, 430, 431, 432, 51, 52, 435, 436, 437, 438, 963, 964, 966, 719, 80, 593, 82, 594, 83, 597, 598, 720, 601, 602, 603, 604, 605, 606, 609, 610, 104, 106, 1012, 1015, 1016}, {433, 434, 435, 52, 53, 438, 437, 436, 441, 442, 443, 444, 965, 966, 968, 717, 718, 82, 83, 84, 601, 602, 603, 604, 607, 608, 609, 610, 611, 612, 615, 104, 616, 105, 1012, 1013, 1014}, {53, 54, 439, 440, 441, 442, 443, 444, 447, 448, 449, 450, 967, 968, 970, 83, 84, 85, 725, 726, 607, 608, 611, 612, 613, 614, 615, 616, 105, 617, 618, 107, 623, 624, 1013, 1017, 1018}, {1041, 54, 55, 445, 446, 447, 448, 449, 450, 453, 454, 455, 456, 969, 970, 972, 84, 85, 86, 88, 613, 614, 635, 617, 618, 619, 620, 107, 621, 623, 624, 622, 625, 626, 1017, 1019, 636}, {55, 56, 57, 58, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 971, 972, 973, 463, 464, 461, 462, 977, 85, 86, 467, 468, 469, 470, 88, 619, 620, 1004, 621, 622, 625, 626, 1019}, {1024, 643, 644, 637, 59, 60, 974, 976, 979, 471, 87, 89, 472, 473, 474, 475, 88, 476, 479, 480, 481, 482, 747, 748, 109, 110, 627, 628, 629, 630, 631, 632, 633, 634, 1021, 638, 1023}, {1041, 1049, 107, 58, 59, 975, 976, 465, 466, 467, 468, 469, 470, 87, 88, 473, 474, 475, 476, 86, 85, 735, 736, 635, 619, 1004, 620, 621, 622, 110, 625, 626, 627, 628, 623, 624, 633, 634, 1019, 636, 637, 638, 1023}, {640, 641, 642, 643, 644, 647, 648, 60, 61, 978, 979, 981, 87, 89, 90, 477, 478, 479, 480, 481, 482, 739, 740, 485, 486, 487, 488, 108, 109, 629, 630, 631, 632, 1020, 1021, 1022, 639}, {640, 641, 642, 1028, 645, 646, 647, 648, 649, 650, 1029, 653, 654, 61, 62, 980, 981, 983, 89, 90, 91, 483, 484, 485, 486, 487, 488, 741, 742, 491, 492, 493, 494, 108, 112, 1020, 639}, {1028, 645, 646, 1030, 1031, 649, 650, 651, 652, 653, 654, 655, 656, 661, 662, 62, 63, 982, 983, 985, 90, 91, 92, 489, 490, 491, 492, 493, 494, 112, 497, 498, 499, 500, 113, 761, 762}, {1030, 1032, 651, 652, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 1050, 673, 674, 63, 64, 984, 985, 91, 92, 987, 93, 95, 495, 496, 497, 498, 499, 500, 113, 503, 504, 505, 506}, {513, 514, 515, 516, 1032, 657, 658, 659, 660, 663, 664, 64, 65, 66, 986, 987, 92, 93, 988, 991, 95, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510}, {768, 517, 518, 519, 520, 521, 522, 1036, 525, 526, 527, 528, 665, 666, 667, 668, 1051, 670, 671, 672, 673, 674, 669, 1053, 679, 680, 67, 68, 989, 94, 95, 992, 96, 994, 113, 116, 767}, {512, 513, 514, 515, 516, 519, 520, 521, 522, 1032, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 1050, 1051, 669, 670, 673, 674, 66, 67, 92, 93, 990, 95, 991, 992, 94, 113, 511}, {773, 774, 1033, 523, 524, 525, 526, 527, 528, 1036, 1037, 531, 532, 533, 534, 667, 668, 671, 672, 675, 676, 677, 678, 679, 680, 683, 684, 68, 69, 94, 96, 993, 994, 97, 996, 114, 116}, {771, 772, 1033, 1034, 1035, 529, 530, 531, 532, 533, 534, 537, 538, 539, 540, 675, 676, 677, 678, 681, 682, 683, 684, 685, 686, 689, 690, 69, 70, 96, 97, 98, 995, 996, 998, 114, 115}, {1034, 779, 780, 1038, 1039, 535, 536, 537, 538, 539, 540, 543, 544, 545, 546, 681, 682, 685, 686, 687, 688, 689, 690, 691, 692, 695, 696, 70, 71, 97, 98, 99, 997, 998, 1000, 115, 117}, {1026, 1038, 1040, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 687, 688, 559, 560, 561, 562, 693, 694, 691, 692, 695, 696, 71, 72, 74, 98, 99, 999, 1000, 1005, 111, 117, 757, 758}, {384, 385, 386, 1044, 1048, 793, 794, 553, 42, 43, 554, 557, 558, 565, 566, 1079, 569, 698, 697, 570, 699, 700, 73, 75, 119, 100, 1001, 1002, 1003, 118, 375, 376, 377, 378, 379, 380, 383}, {129, 1042, 1062, 1074, 1075, 571, 572, 701, 702, 703, 576, 704, 705, 579, 580, 581, 582, 706, 573, 574, 707, 715, 76, 77, 716, 575, 711, 712, 733, 734, 101, 102, 103, 106, 708, 1007, 1009}, {129, 130, 1043, 795, 796, 563, 564, 1075, 1077, 567, 568, 569, 570, 1078, 701, 702, 573, 576, 574, 575, 707, 708, 709, 710, 711, 712, 713, 714, 75, 76, 849, 850, 101, 102, 1008, 1009, 118}, {1025, 715, 1042, 1045, 1052, 1062, 703, 704, 577, 578, 579, 580, 581, 582, 705, 706, 585, 586, 587, 588, 77, 78, 591, 80, 81, 592, 595, 596, 716, 597, 599, 600, 598, 101, 103, 106, 1011}, {1056, 1057, 1058, 807, 808, 717, 718, 719, 720, 721, 82, 83, 722, 723, 724, 727, 728, 601, 602, 603, 604, 605, 606, 731, 732, 609, 610, 611, 612, 104, 105, 106, 1012, 1014, 1016, 121, 122}, {1056, 1059, 1060, 809, 810, 717, 718, 721, 722, 83, 84, 725, 726, 727, 728, 729, 730, 607, 608, 609, 610, 611, 612, 737, 738, 615, 616, 105, 104, 617, 107, 618, 1013, 1014, 121, 1018, 123}, {129, 1052, 1057, 1062, 813, 814, 1074, 1076, 703, 704, 705, 706, 707, 708, 715, 716, 719, 80, 593, 82, 594, 720, 597, 598, 595, 596, 723, 724, 603, 604, 605, 606, 731, 732, 733, 734, 101, 103, 104, 106, 1015, 1016, 122}, {1041, 1049, 1059, 1061, 84, 85, 725, 726, 88, 729, 730, 735, 736, 737, 738, 613, 614, 615, 616, 617, 618, 107, 105, 621, 622, 623, 624, 110, 123, 753, 754, 1017, 1018, 635, 636, 637, 638}, {640, 641, 642, 643, 644, 1029, 647, 648, 649, 650, 1020, 766, 1063, 1067, 1068, 825, 826, 89, 90, 739, 740, 741, 742, 743, 744, 745, 746, 108, 109, 749, 750, 112, 126, 124, 765, 1022, 639}, {1024, 641, 642, 643, 644, 1063, 1064, 1066, 823, 824, 125, 87, 89, 739, 740, 743, 744, 747, 108, 109, 748, 110, 749, 750, 751, 752, 755, 629, 630, 631, 632, 633, 634, 756, 124, 1021, 1022}, {1024, 1049, 1061, 1064, 1065, 819, 820, 125, 87, 88, 756, 735, 736, 737, 738, 635, 747, 748, 109, 110, 107, 751, 753, 754, 627, 628, 752, 755, 631, 632, 633, 634, 123, 636, 637, 638, 1023}, {1026, 1027, 1040, 789, 1046, 1047, 790, 551, 552, 555, 556, 557, 558, 559, 560, 561, 562, 693, 694, 695, 696, 73, 74, 99, 759, 111, 760, 117, 757, 119, 758}, {769, 770, 1028, 645, 646, 647, 648, 649, 650, 1029, 1031, 653, 654, 655, 656, 1055, 801, 802, 1067, 1069, 90, 91, 126, 741, 742, 745, 746, 108, 112, 113, 120, 761, 762, 763, 764, 765, 766}, {768, 769, 770, 1030, 1031, 651, 652, 653, 654, 655, 656, 785, 786, 659, 660, 661, 662, 665, 1050, 666, 1051, 669, 670, 671, 672, 673, 674, 1053, 1054, 1055, 91, 92, 94, 95, 112, 113, 116, 120, 761, 762, 763, 764, 767}, {128, 771, 772, 773, 774, 775, 776, 1033, 777, 1035, 778, 1037, 781, 782, 787, 788, 675, 676, 677, 678, 679, 680, 683, 684, 685, 686, 1070, 1071, 1073, 837, 838, 96, 97, 114, 115, 116, 127}, {771, 772, 131, 775, 776, 1034, 1035, 779, 780, 781, 1039, 782, 783, 784, 791, 792, 681, 682, 683, 684, 685, 686, 1070, 689, 690, 691, 692, 1080, 1081, 841, 842, 97, 98, 114, 115, 117, 127}, {768, 769, 770, 128, 773, 774, 777, 778, 1036, 1037, 785, 786, 787, 788, 667, 668, 1053, 669, 671, 672, 670, 1054, 803, 804, 677, 678, 679, 680, 1071, 1072, 94, 96, 113, 114, 116, 120, 767}, {131, 779, 780, 1038, 1039, 1040, 783, 784, 789, 790, 1047, 791, 792, 799, 800, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 1081, 1082, 98, 99, 759, 111, 760, 115, 117, 757, 119, 758}, {130, 131, 1043, 1044, 793, 794, 795, 796, 797, 798, 799, 800, 1077, 565, 566, 567, 568, 697, 698, 699, 700, 569, 570, 1079, 1087, 1091, 709, 710, 713, 714, 75, 853, 854, 100, 102, 118, 119}, {131, 789, 1046, 1047, 790, 1048, 793, 794, 791, 792, 797, 799, 800, 798, 553, 554, 555, 556, 557, 558, 1079, 697, 698, 699, 700, 1082, 1087, 73, 100, 759, 111, 757, 117, 119, 760, 758, 118}, {768, 769, 770, 128, 135, 766, 785, 786, 787, 788, 1054, 1055, 801, 802, 803, 804, 805, 806, 1069, 1072, 835, 836, 1092, 1093, 845, 846, 112, 113, 116, 120, 761, 762, 763, 764, 765, 126, 767}, {132, 1056, 1058, 1060, 807, 808, 809, 810, 811, 812, 815, 816, 821, 822, 1083, 1088, 717, 718, 721, 722, 723, 724, 727, 728, 729, 730, 104, 105, 121, 122, 123}, {129, 132, 133, 1057, 1058, 807, 808, 811, 812, 813, 814, 815, 816, 817, 818, 1076, 1085, 1088, 1094, 719, 720, 721, 722, 723, 724, 851, 852, 731, 732, 733, 734, 859, 860, 104, 106, 121, 122}, {132, 1059, 1060, 1061, 809, 810, 1065, 811, 812, 819, 820, 821, 822, 1083, 1084, 831, 832, 725, 726, 727, 728, 729, 730, 735, 736, 737, 738, 105, 107, 110, 753, 754, 755, 756, 121, 123, 125}, {134, 135, 1063, 1066, 1068, 823, 824, 825, 826, 827, 828, 829, 830, 833, 834, 1089, 835, 836, 1098, 1099, 739, 740, 867, 868, 743, 744, 745, 746, 108, 109, 749, 750, 751, 752, 124, 125, 126}, {132, 134, 1064, 1065, 1066, 819, 820, 821, 822, 823, 824, 827, 1084, 828, 831, 832, 833, 834, 1089, 1090, 861, 862, 747, 748, 109, 110, 751, 752, 753, 754, 755, 756, 750, 749, 123, 124, 125}, {135, 766, 801, 802, 805, 806, 1067, 1068, 1069, 825, 826, 829, 830, 835, 836, 1092, 1098, 741, 742, 743, 744, 745, 746, 108, 112, 124, 120, 763, 764, 765, 126}, {128, 130, 771, 772, 131, 775, 776, 777, 778, 136, 781, 782, 783, 784, 1070, 1073, 1080, 837, 838, 1095, 840, 841, 842, 839, 1096, 843, 844, 847, 848, 1097, 853, 854, 857, 858, 114, 115, 127}, {128, 773, 774, 775, 776, 777, 778, 135, 136, 785, 786, 787, 788, 803, 804, 805, 806, 1071, 1072, 1073, 837, 838, 1093, 1096, 843, 844, 845, 846, 847, 848, 1103, 871, 872, 114, 116, 120, 127}, {129, 130, 133, 813, 814, 817, 1074, 1075, 1076, 818, 1078, 701, 702, 1085, 1086, 705, 706, 707, 708, 711, 712, 713, 714, 849, 850, 851, 852, 855, 856, 731, 732, 733, 734, 101, 102, 106, 122}, {129, 130, 131, 133, 136, 843, 844, 795, 796, 797, 798, 1100, 1077, 1078, 1086, 1091, 709, 710, 711, 712, 713, 714, 839, 840, 841, 842, 1095, 1097, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 865, 866, 102, 118, 127}, {130, 131, 779, 780, 781, 782, 783, 784, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 1080, 1081, 1082, 1087, 1091, 839, 840, 841, 842, 1095, 853, 854, 115, 117, 118, 119, 127}, {132, 133, 134, 807, 808, 809, 810, 811, 812, 815, 816, 817, 818, 819, 820, 821, 822, 1083, 1084, 831, 832, 1088, 833, 834, 1090, 1094, 1101, 859, 860, 861, 862, 863, 864, 121, 122, 123, 125}, {129, 130, 132, 133, 134, 136, 813, 814, 815, 816, 817, 818, 1085, 1086, 1094, 1100, 1101, 1102, 849, 850, 851, 852, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 866, 869, 870, 122}, {132, 133, 134, 135, 136, 823, 824, 827, 828, 829, 830, 831, 832, 833, 834, 1089, 1090, 1099, 1101, 1102, 1104, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 124, 125}, {128, 134, 135, 136, 801, 802, 803, 804, 805, 806, 825, 826, 827, 828, 829, 830, 835, 836, 1092, 1093, 1098, 1099, 845, 846, 847, 848, 1103, 1104, 867, 868, 869, 870, 871, 872, 120, 124, 126}, {128, 130, 133, 134, 135, 136, 837, 838, 839, 1096, 840, 1097, 843, 844, 1100, 1102, 847, 848, 845, 846, 1103, 1104, 855, 856, 857, 858, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 127}]

CreateBlockSmoother can now take these blocks and construct a block Jacobi preconditioner.

blockjac = a.mat.CreateBlockSmoother(blocks)
lams = EigenValues_Preconditioner(mat=a.mat, pre=blockjac)
max(lams)/min(lams)

34.880387978454515

Multiplicative smoothers and its symmetrized version often yield better condition numbers in practice. We can apply the same code we wrote above for symmetrization (SymmetricGS) to the block smoother:

blockgs = SymmetricGS(blockjac)
lams = EigenValues_Preconditioner(mat=a.mat, pre=blockgs)
max(lams)/min(lams)

2.9623083716671723

Add a coarse grid correction

Dependence of the condition number on degrees of freedom can often be reduced by preconditioners that appropriately use a coarse grid correction. We can experiment with coarse grid corrections using NGSolve’s python interface. Here is an example on how to precondition with a coarse grid correction made using the lowest order subspace of fes.

In the example below, note that we use fes.GetDofNrs again. Previously we used it with argument el of type ElementId, while now we use it with an argument v of type MeshNode.

def VertexDofs(mesh, fes):
    vertexdofs = BitArray(fes.ndof)
    vertexdofs[:] = False
    for v in mesh.vertices:
        for d in fes.GetDofNrs(v):
            vertexdofs[d] = True
    vertexdofs &= fes.FreeDofs()
    return vertexdofs

vertexdofs = VertexDofs(mesh, fes)
print(vertexdofs)   # bit array, printed 50 chars/line

0: 00100000000001111111111111111110000000001111111111
50: 11111111111111111111111111111111111111111111111111
100: 11111111111111111111111111111111111110000000000000
150: 00000000000000000000000000000000000000000000000000
200: 00000000000000000000000000000000000000000000000000
250: 00000000000000000000000000000000000000000000000000
300: 00000000000000000000000000000000000000000000000000
350: 00000000000000000000000000000000000000000000000000
400: 00000000000000000000000000000000000000000000000000
450: 00000000000000000000000000000000000000000000000000
500: 00000000000000000000000000000000000000000000000000
550: 00000000000000000000000000000000000000000000000000
600: 00000000000000000000000000000000000000000000000000
650: 00000000000000000000000000000000000000000000000000
700: 00000000000000000000000000000000000000000000000000
750: 00000000000000000000000000000000000000000000000000
800: 00000000000000000000000000000000000000000000000000
850: 00000000000000000000000000000000000000000000000000
900: 00000000000000000000000000000000000000000000000000
950: 00000000000000000000000000000000000000000000000000
1000: 00000000000000000000000000000000000000000000000000
1050: 00000000000000000000000000000000000000000000000000
1100: 00000

Thus we have made a mask vertexdofs which reveals all free dofs associated to vertices. If these are labeled \(c\) (and the remainder is labeled \(f\)), then the matrix \(A\) can partitioned into

\[\begin{split}A = \left( \begin{array}{cc} A_{cc} & A_{cf} \\ A_{fc} & A_{ff} \end{array} \right).\end{split}\]

The matrix coarsepre below represents

\[\begin{split}\left( \begin{array}{cc} A_{cc}^{-1} & 0 \\ 0 & 0 \end{array} \right).\end{split}\]

coarsepre = a.mat.Inverse(vertexdofs)

This matrix can be used for coarse grid correction.

Pitfall!

Note that coarsepre is not appropriate as a preconditioner by itself as it has a large null space. You might get the wrong idea from the results of a Lanczos eigenvalue estimation:

EigenValues_Preconditioner(mat=a.mat, pre=coarsepre)

       1

But this result only gives the Laczos eigenvalue estimates on the range of the preconditioner. The preconditioned operator in this case is simply

\[\begin{split}\left( \begin{array}{cc} A_{cc}^{-1} & 0 \\ 0 & 0 \end{array} \right) \left( \begin{array}{cc} A_{cc} & A_{cf} \\ A_{fc} & A_{ff} \end{array} \right) = \left( \begin{array}{cc} I & A_{cc}^{-1} A_{cf} \\ 0 & 0 \end{array} \right),\end{split}\]

which is a projection into the \(c\)-dofs. Hence its no surprise that Lanczos estimated the eigenvalues of this operator (on its range) to be just one. But this does not suggest that coarsepre has any utility as a preconditioner by itself.

Additive two-grid preconditioner

One well-known correct way to combine the coarse grid correction with one of the previous smoothers is to combine them additively, to get an additive two-grid preconditioner as follows.

twogrid = coarsepre + blockgs

This addition of two operators (of type BaseMatrix) results in another operator, which is stored as an expression, to be evaluated only when needed. The 2-grid preconditioner has a very good condition number:

lam = EigenValues_Preconditioner(mat=a.mat, pre=twogrid)
lam

 0.993407
 0.998068
 0.999962
 1.34402
 1.81233
 1.82051
 1.95709
 1.98356
 1.99985

Combining multigrid with block smoothers

The twogrid preconditioner becomes more expensive on finer meshes due to the large matrix inversion required for the computation of coarsepre. This can be avoided by replacing coarsepre by the multigrid preconditioner we saw in 2.1.1. It is easy to combine your own block smoothers on the finest grid with the built-in multigrid on coarser levels, as the next example shows.

mesh, fes, a, f, gfu = Setup(h=0.5, p=3)
Draw(gfu)
mg = Preconditioner(a, 'multigrid')  # register mg to a
a.Assemble()                         # assemble on coarsest mesh

<ngsolve.comp.BilinearForm at 0x7f026c13fcb0>

Let us refine a few times to make the problem size larger.

for i in range(6):   # finer mesh updates & assembly
    mesh.Refine()
    fes.Update()
    a.Assemble()
print('NDOF = ', fes.ndof)

NDOF =  111361

Next, reset the block smoother to use the blocks on the finest grid.

blocks = VertexPatchBlocks(mesh, fes)
blockjac = a.mat.CreateBlockSmoother(blocks)
blockgs = SymmetricGS(blockjac)

Finally, add the multigrid and symmetrized block Gauss-Seidel preconditioners together to form a new preconditioner.

mgblock = mg.mat + blockgs

lam = EigenValues_Preconditioner(mat=a.mat, pre=mgblock)
lam

 0.806035
 0.844919
 0.895001
 0.931765
 0.996094
 1.11259
 1.36937
 1.42792
 1.54132
 1.65624
 1.76296
 1.84948
  1.9073
 1.97197
 1.99987

Although mgblock is similarly conditioned as twogrid, it is much more efficient.