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.0144918
 0.0565868
 0.0913563
 0.112048
 0.141222
 0.195761
 0.247646
 0.351665
 0.446235
 0.524795
 0.632006
 0.756244
 0.835742
 0.927184
 1.03588
 1.16687
 1.25179
 1.36664
 1.47109
 1.57336
 1.68741
 1.81492
 1.93739
 2.01627
 2.13865
 2.22417
 2.31655
 2.37388
  2.4388
 2.49019
 2.50489
 2.53176
 2.80787

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

max(lams)/min(lams)

193.75502520680777

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)

1561.8273803217198

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.05514173105997197
CG iteration 2, residual = 0.09362431014556014
CG iteration 3, residual = 0.10183798655456612
CG iteration 4, residual = 0.0894640151539913
CG iteration 5, residual = 0.09864154662595719
CG iteration 6, residual = 0.08752084132820961
CG iteration 7, residual = 0.09196141307218679
CG iteration 8, residual = 0.07931974777210136
CG iteration 9, residual = 0.05262078540101388
CG iteration 10, residual = 0.04323706125676107
CG iteration 11, residual = 0.035206451492069794
CG iteration 12, residual = 0.027795575409704187
CG iteration 13, residual = 0.022311587802186686
CG iteration 14, residual = 0.019516985933755487
CG iteration 15, residual = 0.013880031216222444
CG iteration 16, residual = 0.0108194142234518
CG iteration 17, residual = 0.009432534548762429
CG iteration 18, residual = 0.00701882909110787
CG iteration 19, residual = 0.004666218187002317
CG iteration 20, residual = 0.003164824338234744
CG iteration 21, residual = 0.0017730431581364221
CG iteration 22, residual = 0.0011842359412064284
CG iteration 23, residual = 0.0010076069006142196
CG iteration 24, residual = 0.0007466659458352517
CG iteration 25, residual = 0.00046973504901749746
CG iteration 26, residual = 0.0003348113460042977
CG iteration 27, residual = 0.00028297480579273854
CG iteration 28, residual = 0.00021828093905188484
CG iteration 29, residual = 0.0001544991773773806
CG iteration 30, residual = 0.00011300774305439136
CG iteration 31, residual = 8.206625685090396e-05
CG iteration 32, residual = 6.899528154493067e-05
CG iteration 33, residual = 5.795236525379527e-05
CG iteration 34, residual = 4.3057180145713665e-05
CG iteration 35, residual = 3.1203668635216604e-05
CG iteration 36, residual = 2.345226320419534e-05
CG iteration 37, residual = 1.786138853014384e-05
CG iteration 38, residual = 1.416656738033428e-05
CG iteration 39, residual = 1.0900081367054276e-05
CG iteration 40, residual = 7.91390948920572e-06
CG iteration 41, residual = 5.1320784199512385e-06
CG iteration 42, residual = 3.953008836662163e-06
CG iteration 43, residual = 2.7810541448692e-06
CG iteration 44, residual = 1.9779747184497486e-06
CG iteration 45, residual = 1.3918628904742241e-06
CG iteration 46, residual = 9.503960877490242e-07
CG iteration 47, residual = 6.381955833858475e-07
CG iteration 48, residual = 4.5071157835581943e-07
CG iteration 49, residual = 3.239790436646943e-07
CG iteration 50, residual = 1.9978682354894165e-07
CG iteration 51, residual = 1.378203320762726e-07
CG iteration 52, residual = 1.0164602197566733e-07
CG iteration 53, residual = 7.211310475331023e-08
CG iteration 54, residual = 5.1854466909551945e-08
CG iteration 55, residual = 3.7315868525952276e-08
CG iteration 56, residual = 2.7811764945272955e-08
CG iteration 57, residual = 1.8544161271782038e-08
CG iteration 58, residual = 1.1904926948915224e-08
CG iteration 59, residual = 7.2770794329562116e-09
CG iteration 60, residual = 4.878300158371808e-09
CG iteration 61, residual = 3.6605711575865344e-09
CG iteration 62, residual = 2.6228740834933517e-09
CG iteration 63, residual = 1.7311231693992392e-09
CG iteration 64, residual = 1.047684632674657e-09
CG iteration 65, residual = 7.030664661358563e-10
CG iteration 66, residual = 5.363882276071544e-10
CG iteration 67, residual = 4.1979431527714086e-10
CG iteration 68, residual = 2.9107645901670605e-10
CG iteration 69, residual = 1.7349682850543667e-10
CG iteration 70, residual = 1.0534698503698504e-10
CG iteration 71, residual = 7.302316135366631e-11
CG iteration 72, residual = 5.6033041247913775e-11
CG iteration 73, residual = 4.089963021344319e-11
CG iteration 74, residual = 2.4897862383281646e-11
CG iteration 75, residual = 1.450886464273219e-11
CG iteration 76, residual = 1.0177902500853327e-11
CG iteration 77, residual = 7.705336097582845e-12
CG iteration 78, residual = 5.342321883059619e-12
CG iteration 79, residual = 3.495020490721236e-12
CG iteration 80, residual = 2.0168020920495257e-12
CG iteration 81, residual = 1.202799859371847e-12
CG iteration 82, residual = 8.409546049057451e-13
CG iteration 83, residual = 6.237776378508858e-13
CG iteration 84, residual = 4.2057902489793705e-13
CG iteration 85, residual = 2.698437985464794e-13
CG iteration 86, residual = 1.7604807654260919e-13
CG iteration 87, residual = 1.1316954415130668e-13
CG iteration 88, residual = 7.890045683921838e-14
CG iteration 89, residual = 5.3139344034525494e-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.08154902419641553
it# 1 , res = 0.07698831978567876
it# 2 , res = 0.07343571088587378
it# 3 , res = 0.07022987452159993
it# 4 , res = 0.06734784593987088
it# 5 , res = 0.06472557072763736
it# 6 , res = 0.0623133582367283
it# 7 , res = 0.060074879731911196
it# 8 , res = 0.057982859669015126
it# 9 , res = 0.0560163166745898
it# 10 , res = 0.05415877995370275
it# 11 , res = 0.052397096854838995
it# 12 , res = 0.050720611703502234
it# 13 , res = 0.04912058411072688
it# 14 , res = 0.04758976638505557
it# 15 , res = 0.04612208963244292
it# 16 , res = 0.04471242607808056
it# 17 , res = 0.04335640617752858
it# 18 , res = 0.04205027600466388
it# 19 , res = 0.040790784835190906
it# 20 , res = 0.039575095747334
it# 21 , res = 0.03840071401068153
it# 22 , res = 0.03726542937710803
it# 23 , res = 0.03616726933615241
it# 24 , res = 0.035104461082522534
it# 25 , res = 0.0340754004485919
it# 26 , res = 0.033078626433638486
it# 27 , res = 0.032112800249860504
it# 28 , res = 0.031176688027209687
it# 29 , res = 0.030269146491763287
it# 30 , res = 0.029389111067791962
it# 31 , res = 0.028535585960640916
it# 32 , res = 0.027707635862510554
it# 33 , res = 0.026904378991044618
it# 34 , res = 0.026124981225003133
it# 35 , res = 0.02536865114502759
it# 36 , res = 0.024634635822794855
it# 37 , res = 0.023922217230401895
it# 38 , res = 0.023230709164974075
it# 39 , res = 0.022559454602299387
it# 40 , res = 0.021907823408608817
it# 41 , res = 0.021275210352110618
it# 42 , res = 0.020661033366095424
it# 43 , res = 0.020064732023773794
it# 44 , res = 0.019485766191853885
it# 45 , res = 0.018923614835480957
it# 46 , res = 0.018377774951776744
it# 47 , res = 0.017847760613016548
it# 48 , res = 0.017333102103614963
it# 49 , res = 0.0168333451376772
it# 50 , res = 0.01634805014601414
it# 51 , res = 0.01587679162328915
it# 52 , res = 0.015419157527437737
it# 53 , res = 0.01497474872472232
it# 54 , res = 0.01454317847480195
it# 55 , res = 0.01412407195104894
it# 56 , res = 0.013717065792048405
it# 57 , res = 0.01332180768081772
it# 58 , res = 0.012937955948776208
it# 59 , res = 0.012565179201917694
it# 60 , res = 0.012203155966991194
it# 61 , res = 0.011851574355791086
it# 62 , res = 0.011510131745914694
it# 63 , res = 0.011178534476553498
it# 64 , res = 0.010856497558072505
it# 65 , res = 0.01054374439428077
it# 66 , res = 0.010240006516432807
it# 67 , res = 0.009945023328111691
it# 68 , res = 0.009658541860242453
it# 69 , res = 0.009380316535567276
it# 70 , res = 0.009110108941986206
it# 71 , res = 0.008847687614229388
it# 72 , res = 0.008592827823380365
it# 73 , res = 0.008345311373817449
it# 74 , res = 0.008104926407179146
it# 75 , res = 0.007871467212997961
it# 76 , res = 0.007644734045675322
it# 77 , res = 0.007424532947499581
it# 78 , res = 0.007210675577431655
it# 79 , res = 0.007002979045406332
it# 80 , res = 0.006801265751912558
it# 81 , res = 0.006605363232636882
it# 82 , res = 0.006415104007966036
it# 83 , res = 0.0062303254371586046
it# 84 , res = 0.006050869577009595
it# 85 , res = 0.0058765830448392516
it# 86 , res = 0.005707316885649529
it# 87 , res = 0.005542926443300176
it# 88 , res = 0.00538327123556271
it# 89 , res = 0.005228214832920554
it# 90 , res = 0.005077624740985917
it# 91 , res = 0.004931372286416372
it# 92 , res = 0.004789332506213412
it# 93 , res = 0.0046513840402931245
it# 94 , res = 0.00451740902722616
it# 95 , res = 0.004387293003042702
it# 96 , res = 0.004260924803010059
it# 97 , res = 0.004138196466285686
it# 98 , res = 0.004019003143360063
it# 99 , res = 0.0039032430062018084
it# 100 , res = 0.003790817161022731
it# 101 , res = 0.003681629563582871
it# 102 , res = 0.003575586936960342
it# 103 , res = 0.0034725986917093134
it# 104 , res = 0.0033725768483363907
it# 105 , res = 0.003275435962026372
it# 106 , res = 0.003181093049548195
it# 107 , res = 0.003089467518280411
it# 108 , res = 0.0030004810972891036
it# 109 , res = 0.002914057770401241
it# 110 , res = 0.002830123711212999
it# 111 , res = 0.0027486072199768563
it# 112 , res = 0.002669438662312354
it# 113 , res = 0.0025925504096869926
it# 114 , res = 0.0025178767816157696
it# 115 , res = 0.0024453539895301023
it# 116 , res = 0.002374920082265284
it# 117 , res = 0.002306514893122229
it# 118 , res = 0.0022400799884551973
it# 119 , res = 0.0021755586177435644
it# 120 , res = 0.0021128956651026657
it# 121 , res = 0.0020520376021928724
it# 122 , res = 0.0019929324424863507
it# 123 , res = 0.0019355296968528882
it# 124 , res = 0.001879780330424566
it# 125 , res = 0.0018256367207050569
it# 126 , res = 0.0017730526168851283
it# 127 , res = 0.0017219831003313863
it# 128 , res = 0.0016723845462129105
it# 129 , res = 0.0016242145862339033
it# 130 , res = 0.0015774320724403345
it# 131 , res = 0.0015319970420685444
it# 132 , res = 0.001487870683407928
it# 133 , res = 0.0014450153026455206
it# 134 , res = 0.0014033942916667076
it# 135 , res = 0.0013629720967834164
it# 136 , res = 0.0013237141883634587
it# 137 , res = 0.0012855870313332957
it# 138 , res = 0.0012485580565323362
it# 139 , res = 0.001212595632891253
it# 140 , res = 0.0011776690404121562
it# 141 , res = 0.001143748443926802
it# 142 , res = 0.0011108048676108555
it# 143 , res = 0.001078810170232208
it# 144 , res = 0.0010477370211120218
it# 145 , res = 0.0010175588767790855
it# 146 , res = 0.0009882499582944116
it# 147 , res = 0.0009597852292321164
it# 148 , res = 0.0009321403742918262
it# 149 , res = 0.0009052917785286286
it# 150 , res = 0.0008792165071802845
it# 151 , res = 0.0008538922860765319
it# 152 , res = 0.000829297482611615
it# 153 , res = 0.0008054110872656348
it# 154 , res = 0.0007822126956577797
it# 155 , res = 0.0007596824911161527
it# 156 , res = 0.0007378012277509984
it# 157 , res = 0.0007165502140129352
it# 158 , res = 0.000695911296728008
it# 159 , res = 0.000675866845589638
it# 160 , res = 0.0006563997380989742
it# 161 , res = 0.0006374933449383231
it# 162 , res = 0.00061913151576634
it# 163 , res = 0.0006012985654219236
it# 164 , res = 0.0005839792605254148
it# 165 , res = 0.0005671588064661007
it# 166 , res = 0.0005508228347642415
it# 167 , res = 0.0005349573907977536
it# 168 , res = 0.0005195489218806896
it# 169 , res = 0.000504584265687486
it# 170 , res = 0.0004900506390085957
it# 171 , res = 0.00047593562683111574
it# 172 , res = 0.0004622271717336364
it# 173 , res = 0.0004489135635858828
it# 174 , res = 0.0004359834295469073
it# 175 , res = 0.00042342572434875384
it# 176 , res = 0.0004112297208621661
it# 177 , res = 0.00039938500093308747
it# 178 , res = 0.0003878814464829784
it# 179 , res = 0.0003767092308659766
it# 180 , res = 0.0003658588104745064
it# 181 , res = 0.00035532091658703187
it# 182 , res = 0.00034508654745059244
it# 183 , res = 0.00033514696059091343
it# 184 , res = 0.00032549366534458715
it# 185 , res = 0.00031611841560625676
it# 186 , res = 0.00030701320278448763
it# 187 , res = 0.00029817024896047523
it# 188 , res = 0.00028958200024409064
it# 189 , res = 0.00028124112032107716
it# 190 , res = 0.00027314048418642273
it# 191 , res = 0.0002652731720577016
it# 192 , res = 0.00025763246346381224
it# 193 , res = 0.0002502118315049001
it# 194 , res = 0.00024300493727650376
it# 195 , res = 0.00023600562445428026
it# 196 , res = 0.000229207914035672
it# 197 , res = 0.00022260599923279697
it# 198 , res = 0.000216194240510804
it# 199 , res = 0.00020996716077177903
it# 200 , res = 0.00020391944067588703
it# 201 , res = 0.00019804591409642238
it# 202 , res = 0.00019234156370810147
it# 203 , res = 0.0001868015167001111
it# 204 , res = 0.0001814210406143494
it# 205 , res = 0.00017619553930200784
it# 206 , res = 0.00017112054899849432
it# 207 , res = 0.00016619173450954522
it# 208 , res = 0.00016140488550807952
it# 209 , res = 0.0001567559129379746
it# 210 , res = 0.00015224084552077757
it# 211 , res = 0.00014785582636353791
it# 212 , res = 0.00014359710966396427
it# 213 , res = 0.00013946105751092876
it# 214 , res = 0.00013544413677638774
it# 215 , res = 0.00013154291609756934
it# 216 , res = 0.00012775406294634033
it# 217 , res = 0.0001240743407813531
it# 218 , res = 0.00012050060628418139
it# 219 , res = 0.00011702980667384358
it# 220 , res = 0.00011365897709933486
it# 221 , res = 0.00011038523810642975
it# 222 , res = 0.00010720579317854559
it# 223 , res = 0.00010411792634767921
it# 224 , res = 0.00010111899987404832
it# 225 , res = 9.820645199346894e-05
it# 226 , res = 9.53777947284792e-05
it# 227 , res = 9.263061176340115e-05
it# 228 , res = 8.99625563799107e-05
it# 229 , res = 8.737134945269809e-05
it# 230 , res = 8.485477750261787e-05
it# 231 , res = 8.241069080562388e-05
it# 232 , res = 8.003700155646018e-05
it# 233 , res = 7.773168208551489e-05
it# 234 , res = 7.549276312627871e-05
it# 235 , res = 7.331833213351093e-05
it# 236 , res = 7.120653164919626e-05
it# 237 , res = 6.915555771600876e-05
it# 238 , res = 6.716365833642106e-05
it# 239 , res = 6.522913197611547e-05
it# 240 , res = 6.335032610999069e-05
it# 241 , res = 6.152563581116582e-05
it# 242 , res = 5.975350237961197e-05
it# 243 , res = 5.803241201113983e-05
it# 244 , res = 5.636089450328863e-05
it# 245 , res = 5.473752200074414e-05
it# 246 , res = 5.3160907774555937e-05
it# 247 , res = 5.1629705038003324e-05
it# 248 , res = 5.014260579634901e-05
it# 249 , res = 4.86983397289116e-05
it# 250 , res = 4.729567310456153e-05
it# 251 , res = 4.5933407727361e-05
it# 252 , res = 4.461037991323045e-05
it# 253 , res = 4.332545949601001e-05
it# 254 , res = 4.2077548861641145e-05
it# 255 , res = 4.086558201112279e-05
it# 256 , res = 3.9688523649569016e-05
it# 257 , res = 3.85453683016029e-05
it# 258 , res = 3.743513945318645e-05
it# 259 , res = 3.6356888716439906e-05
it# 260 , res = 3.530969502045334e-05
it# 261 , res = 3.42926638238473e-05
it# 262 , res = 3.330492635112865e-05
it# 263 , res = 3.234563885029479e-05
it# 264 , res = 3.141398187151549e-05
it# 265 , res = 3.050915956841697e-05
it# 266 , res = 2.9630399017100668e-05
it# 267 , res = 2.8776949556564516e-05
it# 268 , res = 2.7948082147171685e-05
it# 269 , res = 2.7143088747750425e-05
it# 270 , res = 2.6361281711144098e-05
it# 271 , res = 2.5601993196574352e-05
it# 272 , res = 2.486457459926337e-05
it# 273 , res = 2.414839599616004e-05
it# 274 , res = 2.3452845608048666e-05
it# 275 , res = 2.2777329276979837e-05
it# 276 , res = 2.212126995854093e-05
it# 277 , res = 2.1484107229179095e-05
it# 278 , res = 2.0865296807028026e-05
it# 279 , res = 2.0264310087510022e-05
it# 280 , res = 1.9680633691464185e-05
it# 281 , res = 1.9113769026573634e-05
it# 282 , res = 1.8563231861888455e-05
it# 283 , res = 1.8028551913545032e-05
it# 284 , res = 1.7509272443238964e-05
it# 285 , res = 1.7004949868596263e-05
it# 286 , res = 1.651515338349012e-05
it# 287 , res = 1.6039464590383427e-05
it# 288 , res = 1.5577477143326494e-05
it# 289 , res = 1.5128796399655212e-05
it# 290 , res = 1.4693039084370336e-05
it# 291 , res = 1.4269832961721445e-05
it# 292 , res = 1.385881651754313e-05
it# 293 , res = 1.3459638650514473e-05
it# 294 , res = 1.3071958371928187e-05
it# 295 , res = 1.2695444514873325e-05
it# 296 , res = 1.2329775450892118e-05
it# 297 , res = 1.1974638815679472e-05
it# 298 , res = 1.1629731241894019e-05
it# 299 , res = 1.1294758099956226e-05
it# 300 , res = 1.0969433246766563e-05
it# 301 , res = 1.0653478781223555e-05
it# 302 , res = 1.0346624806177666e-05
it# 303 , res = 1.0048609198968241e-05
it# 304 , res = 9.759177386223277e-06
it# 305 , res = 9.478082127585961e-06
it# 306 , res = 9.205083303781558e-06
it# 307 , res = 8.939947711947309e-06
it# 308 , res = 8.682448866029702e-06
it# 309 , res = 8.432366803522394e-06
it# 310 , res = 8.189487897421104e-06
it# 311 , res = 7.953604674076511e-06
it# 312 , res = 7.724515635701855e-06
it# 313 , res = 7.502025088089216e-06
it# 314 , res = 7.285942973921041e-06
it# 315 , res = 7.076084709880889e-06
it# 316 , res = 6.872271029337415e-06
it# 317 , res = 6.67432782923665e-06
it# 318 , res = 6.482086020986079e-06
it# 319 , res = 6.29538138659759e-06
it# 320 , res = 6.114054437685901e-06
it# 321 , res = 5.937950279934657e-06
it# 322 , res = 5.766918480527648e-06
it# 323 , res = 5.600812939295695e-06
it# 324 , res = 5.439491764367114e-06
it# 325 , res = 5.282817150945344e-06
it# 326 , res = 5.130655263201839e-06
it# 327 , res = 4.9828761204626765e-06
it# 328 , res = 4.839353485683533e-06
it# 329 , res = 4.699964757990287e-06
it# 330 , res = 4.564590867627099e-06
it# 331 , res = 4.433116174650422e-06
it# 332 , res = 4.3054283698132924e-06
it# 333 , res = 4.181418378650014e-06
it# 334 , res = 4.060980268496734e-06
it# 335 , res = 3.94401115789135e-06
it# 336 , res = 3.830411128705417e-06
it# 337 , res = 3.720083140572494e-06
it# 338 , res = 3.6129329483195256e-06
it# 339 , res = 3.508869021532425e-06
it# 340 , res = 3.407802465867934e-06
it# 341 , res = 3.3096469474337217e-06
it# 342 , res = 3.2143186192184127e-06
it# 343 , res = 3.1217360491737653e-06
it# 344 , res = 3.0318201507558456e-06
it# 345 , res = 2.9444941154434957e-06
it# 346 , res = 2.8596833468371776e-06
it# 347 , res = 2.7773153974419775e-06
it# 348 , res = 2.697319906004142e-06
it# 349 , res = 2.619628538545977e-06
it# 350 , res = 2.544174928821312e-06
it# 351 , res = 2.4708946223660645e-06
it# 352 , res = 2.3997250211257222e-06
it# 353 , res = 2.3306053300706394e-06
it# 354 , res = 2.2634765053544437e-06
it# 355 , res = 2.198281203807783e-06
it# 356 , res = 2.134963733591044e-06
it# 357 , res = 2.073470007246095e-06
it# 358 , res = 2.0137474952864376e-06
it# 359 , res = 1.955745181080915e-06
it# 360 , res = 1.899413517516605e-06
it# 361 , res = 1.8447043844596275e-06
it# 362 , res = 1.7915710479762233e-06
it# 363 , res = 1.7399681199862072e-06
it# 364 , res = 1.6898515201829538e-06
it# 365 , res = 1.6411784371726523e-06
it# 366 , res = 1.5939072934037181e-06
it# 367 , res = 1.5479977084870128e-06
it# 368 , res = 1.5034104653253024e-06
it# 369 , res = 1.4601074759647723e-06
it# 370 , res = 1.4180517500992914e-06
it# 371 , res = 1.3772073625916025e-06
it# 372 , res = 1.3375394228747665e-06
it# 373 , res = 1.2990140457361015e-06
it# 374 , res = 1.2615983215359587e-06
it# 375 , res = 1.225260288890629e-06
it# 376 , res = 1.1899689068565927e-06
it# 377 , res = 1.1556940284555925e-06
it# 378 , res = 1.1224063753903157e-06
it# 379 , res = 1.0900775121497274e-06
it# 380 , res = 1.05867982276121e-06
it# 381 , res = 1.0281864864177523e-06
it# 382 , res = 9.98571454864159e-07
it# 383 , res = 9.698094301132213e-07
it# 384 , res = 9.418758428884068e-07
it# 385 , res = 9.147468316431556e-07
it# 386 , res = 8.883992219039149e-07
it# 387 , res = 8.62810506855781e-07
it# 388 , res = 8.379588278791083e-07
it# 389 , res = 8.138229562411623e-07
it# 390 , res = 7.903822741245798e-07
it# 391 , res = 7.676167581023888e-07
it# 392 , res = 7.455069609863247e-07
it# 393 , res = 7.240339963135654e-07
it# 394 , res = 7.031795210344196e-07
it# 395 , res = 6.829257207424407e-07
it# 396 , res = 6.632552942546613e-07
it# 397 , res = 6.441514382868382e-07
it# 398 , res = 6.255978340744685e-07
it# 399 , res = 6.075786323223278e-07
it# 400 , res = 5.90078440914256e-07
it# 401 , res = 5.730823105435421e-07
it# 402 , res = 5.565757225816124e-07
it# 403 , res = 5.405445766423993e-07
it# 404 , res = 5.24975178564602e-07
it# 405 , res = 5.098542285493482e-07
it# 406 , res = 4.951688097564827e-07
it# 407 , res = 4.809063777066817e-07
it# 408 , res = 4.6705474886287745e-07
it# 409 , res = 4.5360209092639554e-07
it# 410 , res = 4.405369121590279e-07
it# 411 , res = 4.2784805200412186e-07
it# 412 , res = 4.155246712565314e-07
it# 413 , res = 4.0355624292501993e-07
it# 414 , res = 3.9193254330675177e-07
it# 415 , res = 3.806436430414755e-07
it# 416 , res = 3.696798990759136e-07
it# 417 , res = 3.590319455474361e-07
it# 418 , res = 3.4869068686795014e-07
it# 419 , res = 3.386472892191472e-07
it# 420 , res = 3.2889317337383455e-07
it# 421 , res = 3.1942000696388067e-07
it# 422 , res = 3.102196977233625e-07
it# 423 , res = 3.0128438668781334e-07
it# 424 , res = 2.926064408322276e-07
it# 425 , res = 2.8417844740400097e-07
it# 426 , res = 2.759932068834431e-07
it# 427 , res = 2.68043727196597e-07
it# 428 , res = 2.603232178645734e-07
it# 429 , res = 2.5282508353274463e-07
it# 430 , res = 2.455429193502768e-07
it# 431 , res = 2.3847050452922816e-07
it# 432 , res = 2.3160179767358632e-07
it# 433 , res = 2.24930931461046e-07
it# 434 , res = 2.1845220732370615e-07
it# 435 , res = 2.1216009092983153e-07
it# 436 , res = 2.060492075002507e-07
it# 437 , res = 2.0011433686188444e-07
it# 438 , res = 1.9435040928430636e-07
it# 439 , res = 1.8875250122155405e-07
it# 440 , res = 1.8331583060078257e-07
it# 441 , res = 1.780357532344608e-07
it# 442 , res = 1.7290775892418325e-07
it# 443 , res = 1.6792746716270367e-07
it# 444 , res = 1.6309062350842178e-07
it# 445 , res = 1.5839309630431889e-07
it# 446 , res = 1.538308728051071e-07
it# 447 , res = 1.494000558782948e-07
it# 448 , res = 1.4509686047338516e-07
it# 449 , res = 1.4091761070014672e-07
it# 450 , res = 1.3685873665495154e-07
it# 451 , res = 1.3291677098868567e-07
it# 452 , res = 1.2908834637893084e-07
it# 453 , res = 1.2537019269001405e-07
it# 454 , res = 1.2175913345699953e-07
it# 455 , res = 1.1825208434115628e-07
it# 456 , res = 1.1484604931719646e-07
it# 457 , res = 1.1153811885454406e-07
it# 458 , res = 1.083254673225213e-07
it# 459 , res = 1.0520535030204047e-07
it# 460 , res = 1.0217510258931045e-07
it# 461 , res = 9.923213553203039e-08
it# 462 , res = 9.637393543321519e-08
it# 463 , res = 9.359806036537822e-08
it# 464 , res = 9.090213942258633e-08
it# 465 , res = 8.828386948285739e-08
it# 466 , res = 8.574101417944709e-08
it# 467 , res = 8.327140098270293e-08
it# 468 , res = 8.087292052105102e-08
it# 469 , res = 7.85435239546566e-08
it# 470 , res = 7.628122130819316e-08
it# 471 , res = 7.408408031921845e-08
it# 472 , res = 7.195022385854986e-08
it# 473 , res = 6.987782923755689e-08
it# 474 , res = 6.786512615906708e-08
it# 475 , res = 6.591039529181838e-08
it# 476 , res = 6.401196685792243e-08
it# 477 , res = 6.21682192379275e-08
it# 478 , res = 6.037757727924875e-08
it# 479 , res = 5.863851162842097e-08
it# 480 , res = 5.694953656287484e-08
it# 481 , res = 5.530920927349143e-08
it# 482 , res = 5.371612861423911e-08
it# 483 , res = 5.216893374509203e-08
it# 484 , res = 5.0666303039814737e-08
it# 485 , res = 4.920695283710884e-08
it# 486 , res = 4.778963648778357e-08
it# 487 , res = 4.6413143468722836e-08
it# 488 , res = 4.5076297661127835e-08
it# 489 , res = 4.37779574391854e-08
it# 490 , res = 4.251701340416732e-08
it# 491 , res = 4.129238853596638e-08
it# 492 , res = 4.010303683907488e-08
it# 493 , res = 3.894794217813425e-08
it# 494 , res = 3.7826117920355016e-08
it# 495 , res = 3.6736605752768957e-08
it# 496 , res = 3.567847501321533e-08
it# 497 , res = 3.465082178379079e-08
it# 498 , res = 3.365276820593264e-08
it# 499 , res = 3.2683461807860314e-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.09414079549704508
CG iteration 2, residual = 0.12203689059859807
CG iteration 3, residual = 0.0863034962644804
CG iteration 4, residual = 0.054765061373381614
CG iteration 5, residual = 0.02703987161164137
CG iteration 6, residual = 0.014941028668260055
CG iteration 7, residual = 0.006694630914310554
CG iteration 8, residual = 0.00284839806274339
CG iteration 9, residual = 0.000987701508151471
CG iteration 10, residual = 0.0003968124927947859
CG iteration 11, residual = 0.00012407544574560332
CG iteration 12, residual = 6.330209216165865e-05
CG iteration 13, residual = 2.9785257662869824e-05
CG iteration 14, residual = 1.339030629221604e-05
CG iteration 15, residual = 5.363365741544918e-06
CG iteration 16, residual = 2.0462661924717143e-06
CG iteration 17, residual = 8.231874195477862e-07
CG iteration 18, residual = 3.4544567760902373e-07
CG iteration 19, residual = 1.2951230160193418e-07
CG iteration 20, residual = 5.355359552432075e-08
CG iteration 21, residual = 1.7611013168353534e-08
CG iteration 22, residual = 5.970141704993225e-09
CG iteration 23, residual = 2.467425255803206e-09
CG iteration 24, residual = 7.865934453395217e-10
CG iteration 25, residual = 3.1722300643176497e-10
CG iteration 26, residual = 1.1006534206433618e-10
CG iteration 27, residual = 4.0954503446546655e-11
CG iteration 28, residual = 1.599634501034386e-11
CG iteration 29, residual = 5.326653784450982e-12
CG iteration 30, residual = 1.988704684131753e-12
CG iteration 31, residual = 5.843515407395414e-13
CG iteration 32, residual = 2.2678569699661713e-13
CG iteration 33, residual = 6.065405904122058e-14

BaseWebGuiScene

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

20.311617625163528

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)

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

35.052324628378834

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.8770481008864492

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: 11111111111111111111111111111111111100000000000000
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: 0000000000000000000000000000000000000000000000

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.993158
 0.997949
 0.999962
 1.35426
 1.82282
 1.90685
 1.96916
 1.98517
  1.9997

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 0x7f03dc135370>

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.826237
 0.849126
 0.897069
 0.933508
 1.00501
 1.11965
 1.36946
 1.43186
 1.54355
 1.65858
 1.76409
 1.85156
 1.90736
 1.97186
 1.99987

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