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

2.1.2 Building blocks for programming preconditioners

In this tutorial we will discuss

• 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 netgen.geom2d import unit_square
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(fes)
    a += grad(u)*grad(v)*dx + u*v*dx
    f = LinearForm(fes)
    f += 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()

<ngsolve.comp.LinearForm at 0x7f64c528cbf0>

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.0144926
 0.0564694
 0.0914842
 0.112088
 0.141509
 0.196298
 0.24898
 0.350942
 0.446016
 0.524543
 0.631951
 0.755926
 0.836108
 0.928214
 1.03829
  1.1669
 1.24882
 1.36421
 1.47177
 1.56997
  1.6888
 1.80786
 1.93984
 2.01178
 2.14276
 2.22343
 2.32306
 2.37764
 2.44145
 2.49227
 2.50466
 2.53857
 2.81153

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

max(lams)/min(lams)

193.9981173490965

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)

1570.206119819098

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)

WARNING: maxsteps is deprecated, use maxiter instead!
CG iteration 1, residual = 0.05512594965944649
CG iteration 2, residual = 0.09391567400128714
CG iteration 3, residual = 0.10246000044758205
CG iteration 4, residual = 0.08873850258940057
CG iteration 5, residual = 0.09874906344132314
CG iteration 6, residual = 0.08747939898167575
CG iteration 7, residual = 0.09128016616311646
CG iteration 8, residual = 0.0791956254476097
CG iteration 9, residual = 0.05330659844275126
CG iteration 10, residual = 0.04307660394188752
CG iteration 11, residual = 0.035617827784361906
CG iteration 12, residual = 0.02789752363830375
CG iteration 13, residual = 0.022718179737348052
CG iteration 14, residual = 0.019727678781107868
CG iteration 15, residual = 0.014290616919723395
CG iteration 16, residual = 0.011091139866485028
CG iteration 17, residual = 0.009639213829158201
CG iteration 18, residual = 0.007125510641526397
CG iteration 19, residual = 0.0047624262344669556
CG iteration 20, residual = 0.003226667693700545
CG iteration 21, residual = 0.0018108143713303159
CG iteration 22, residual = 0.0012193549503749512
CG iteration 23, residual = 0.0010106346909277552
CG iteration 24, residual = 0.0007755358613273704
CG iteration 25, residual = 0.0004810987056609594
CG iteration 26, residual = 0.000340838076605625
CG iteration 27, residual = 0.000281424337955284
CG iteration 28, residual = 0.00022326399718238923
CG iteration 29, residual = 0.00015673369168887007
CG iteration 30, residual = 0.00011473436199699396
CG iteration 31, residual = 8.176258331617258e-05
CG iteration 32, residual = 6.909747711413725e-05
CG iteration 33, residual = 5.818927681384233e-05
CG iteration 34, residual = 4.2700302819676935e-05
CG iteration 35, residual = 3.08401438150349e-05
CG iteration 36, residual = 2.339272507941274e-05
CG iteration 37, residual = 1.753259412068235e-05
CG iteration 38, residual = 1.402732973402308e-05
CG iteration 39, residual = 1.0983127415964254e-05
CG iteration 40, residual = 8.006229771250064e-06
CG iteration 41, residual = 5.17309362019899e-06
CG iteration 42, residual = 3.845892074500316e-06
CG iteration 43, residual = 2.844359902943899e-06
CG iteration 44, residual = 2.020854613531387e-06
CG iteration 45, residual = 1.4192478693100715e-06
CG iteration 46, residual = 9.72381969099622e-07
CG iteration 47, residual = 6.408805372619777e-07
CG iteration 48, residual = 4.6008294229236534e-07
CG iteration 49, residual = 3.2875521350168574e-07
CG iteration 50, residual = 2.0937320695325325e-07
CG iteration 51, residual = 1.3531555501210955e-07
CG iteration 52, residual = 9.877177665547215e-08
CG iteration 53, residual = 7.07097797324112e-08
CG iteration 54, residual = 5.372229001089104e-08
CG iteration 55, residual = 3.796030574369298e-08
CG iteration 56, residual = 2.7372591402441928e-08
CG iteration 57, residual = 1.8704231204220455e-08
CG iteration 58, residual = 1.2090336245670144e-08
CG iteration 59, residual = 7.3006475756812826e-09
CG iteration 60, residual = 4.495414383915852e-09
CG iteration 61, residual = 3.3246736071415413e-09
CG iteration 62, residual = 2.392299791391423e-09
CG iteration 63, residual = 1.5947684052573668e-09
CG iteration 64, residual = 1.0435355977677207e-09
CG iteration 65, residual = 6.959500480557338e-10
CG iteration 66, residual = 5.312559273824841e-10
CG iteration 67, residual = 4.1118439010046324e-10
CG iteration 68, residual = 2.8962302224610927e-10
CG iteration 69, residual = 1.7271886422841232e-10
CG iteration 70, residual = 1.07963546966214e-10
CG iteration 71, residual = 7.488799066846355e-11
CG iteration 72, residual = 5.532064035922531e-11
CG iteration 73, residual = 3.9698286235580665e-11
CG iteration 74, residual = 2.507107753052256e-11
CG iteration 75, residual = 1.526175973459143e-11
CG iteration 76, residual = 1.0325343132723367e-11
CG iteration 77, residual = 7.578565822870134e-12
CG iteration 78, residual = 5.196599298169043e-12
CG iteration 79, residual = 3.5518900795842666e-12
CG iteration 80, residual = 2.1685214799947498e-12
CG iteration 81, residual = 1.3136253620509309e-12
CG iteration 82, residual = 8.595049185746301e-13
CG iteration 83, residual = 5.984329049331948e-13
CG iteration 84, residual = 4.3198485382294164e-13
CG iteration 85, residual = 3.0347543373935257e-13
CG iteration 86, residual = 1.9156538617009383e-13
CG iteration 87, residual = 1.1336925960575171e-13
CG iteration 88, residual = 7.716438928027766e-14
CG iteration 89, residual = 5.620731319496194e-14
CG iteration 90, residual = 3.7348468950444507e-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(). Its 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.08160965515786975
it# 1 , res = 0.07705906465506132
it# 2 , res = 0.07350156794600012
it# 3 , res = 0.07028688539903526
it# 4 , res = 0.06739530161684565
it# 5 , res = 0.06476364944760173
it# 6 , res = 0.062342780807217495
it# 7 , res = 0.06009661735499908
it# 8 , res = 0.057997960786708075
it# 9 , res = 0.056025810660797135
it# 10 , res = 0.05416362287952465
it# 11 , res = 0.05239814445301902
it# 12 , res = 0.05071860990359418
it# 13 , res = 0.04911617036914302
it# 14 , res = 0.04758347710615852
it# 15 , res = 0.04611437055532574
it# 16 , res = 0.04470364366504399
it# 17 , res = 0.04334685884586829
it# 18 , res = 0.0420402045743434
it# 19 , res = 0.04078038189946108
it# 20 , res = 0.0395645138747252
it# 21 , res = 0.0383900728015662
it# 22 , res = 0.037254821458334995
it# 23 , res = 0.03615676540439092
it# 24 , res = 0.03509411411463134
it# 25 , res = 0.0340652491941432
it# 26 , res = 0.033068698295959674
it# 27 , res = 0.032103113650757346
it# 28 , res = 0.031167254338726247
it# 29 , res = 0.030259971606919605
it# 30 , res = 0.029380196671723764
it# 31 , res = 0.028526930554168488
it# 32 , res = 0.027699235581922253
it# 33 , res = 0.026896228260764148
it# 34 , res = 0.02611707327371169
it# 35 , res = 0.025360978410640445
it# 36 , res = 0.024627190267318525
it# 37 , res = 0.02391499058202732
it# 38 , res = 0.023223693101684383
it# 39 , res = 0.022552640888706964
it# 40 , res = 0.021901203995600544
it# 41 , res = 0.021268777447111694
it# 42 , res = 0.020654779480294923
it# 43 , res = 0.020058650001448762
it# 44 , res = 0.019479849225931947
it# 45 , res = 0.018917856472666587
it# 46 , res = 0.018372169089894697
it# 47 , res = 0.017842301492683413
it# 48 , res = 0.01732778429590299
it# 49 , res = 0.016828163529074326
it# 50 , res = 0.01634299992169302
it# 51 , res = 0.015871868249463992
it# 52 , res = 0.015414356733397727
it# 53 , res = 0.014970066484981082
it# 54 , res = 0.014538610991684586
it# 55 , res = 0.014119615637940796
it# 56 , res = 0.01371271725746021
it# 57 , res = 0.013317563713361317
it# 58 , res = 0.012933813503103173
it# 59 , res = 0.012561135385639008
it# 60 , res = 0.012199208028571919
it# 61 , res = 0.01184771967339634
it# 62 , res = 0.011506367817168376
it# 63 , res = 0.0111748589091667
it# 64 , res = 0.010852908061288415
it# 65 , res = 0.010540238771083841
it# 66 , res = 0.01023658265646672
it# 67 , res = 0.009941679201251958
it# 68 , res = 0.009655275510770996
it# 69 , res = 0.009377126076896773
it# 70 , res = 0.00910699255188732
it# 71 , res = 0.008844643530514068
it# 72 , res = 0.008589854339999876
it# 73 , res = 0.008342406837334343
it# 74 , res = 0.008102089213577695
it# 75 , res = 0.007868695804799018
it# 76 , res = 0.0076420269093261375
it# 77 , res = 0.007421888611008521
it# 78 , res = 0.007208092608225099
it# 79 , res = 0.007000456048382657
it# 80 , res = 0.006798801367673498
it# 81 , res = 0.006602956135877871
it# 82 , res = 0.00641275290600751
it# 83 , res = 0.0062280290686041895
it# 84 , res = 0.006048626710517259
it# 85 , res = 0.005874392477993166
it# 86 , res = 0.005705177443922282
it# 87 , res = 0.005540836979094875
it# 88 , res = 0.005381230627326925
it# 89 , res = 0.0052262219843237905
it# 90 , res = 0.005075678580155175
it# 91 , res = 0.004929471765221525
it# 92 , res = 0.00478747659959895
it# 93 , res = 0.004649571745650129
it# 94 , res = 0.0045156393638008715
it# 95 , res = 0.004385565011378318
it# 96 , res = 0.0042592375444169646
it# 97 , res = 0.004136549022338745
it# 98 , res = 0.00401739461541948
it# 99 , res = 0.0039016725149549547
it# 100 , res = 0.003789283846045805
it# 101 , res = 0.00368013258292044
it# 102 , res = 0.0035741254667195957
it# 103 , res = 0.003471171925669908
it# 104 , res = 0.0033711839975732504
it# 105 , res = 0.003274076254544013
it# 106 , res = 0.0031797657299277625
it# 107 , res = 0.003088171847335115
it# 108 , res = 0.0029992163517303665
it# 109 , res = 0.002912823242513506
it# 110 , res = 0.002828918708536509
it# 111 , res = 0.002747431064998835
it# 112 , res = 0.0026682906921642346
it# 113 , res = 0.0025914299758499387
it# 114 , res = 0.00251678324963187
it# 115 , res = 0.002444286738719753
it# 116 , res = 0.00237387850545045
it# 117 , res = 0.0023054983963557985
it# 118 , res = 0.002239087990755443
it# 119 , res = 0.0021745905508333135
it# 120 , res = 0.0021119509731542707
it# 121 , res = 0.002051115741577662
it# 122 , res = 0.001992032881528742
it# 123 , res = 0.0019346519155891016
it# 124 , res = 0.0018789238203657833
it# 125 , res = 0.0018248009846049175
it# 126 , res = 0.0017722371685114887
it# 127 , res = 0.001721187464242257
it# 128 , res = 0.0016716082575360116
it# 129 , res = 0.0016234571904511813
it# 130 , res = 0.0015766931251749085
it# 131 , res = 0.001531276108877098
it# 132 , res = 0.0014871673395754816
it# 133 , res = 0.001444329132984478
it# 134 , res = 0.0014027248903200585
it# 135 , res = 0.001362319067030615
it# 136 , res = 0.0013230771424298157
it# 137 , res = 0.0012849655902035762
it# 138 , res = 0.0012479518497675256
it# 139 , res = 0.0012120042984483132
it# 140 , res = 0.0011770922244673326
it# 141 , res = 0.001143185800702178
it# 142 , res = 0.001110256059204391
it# 143 , res = 0.0010782748664499757
it# 144 , res = 0.00104721489930529
it# 145 , res = 0.0010170496216823767
it# 146 , res = 0.0009877532618685547
it# 147 , res = 0.000959300790508967
it# 148 , res = 0.000931667899222562
it# 149 , res = 0.000904830979834465
it# 150 , res = 0.0008787671042067093
it# 151 , res = 0.000853454004649727
it# 152 , res = 0.0008288700548978281
it# 153 , res = 0.0008049942516337092
it# 154 , res = 0.0007818061965434139
it# 155 , res = 0.0007592860788903442
it# 156 , res = 0.0007374146585890037
it# 157 , res = 0.0007161732497681308
it# 158 , res = 0.0006955437048060461
it# 159 , res = 0.0006755083988268753
it# 160 , res = 0.0006560502146421262
it# 161 , res = 0.0006371525281273097
it# 162 , res = 0.0006187991940190667
it# 163 , res = 0.0006009745321212831
it# 164 , res = 0.0005836633139090658
it# 165 , res = 0.0005668507495181527
it# 166 , res = 0.0005505224751093039
it# 167 , res = 0.0005346645405966025
it# 168 , res = 0.00051926339772934
it# 169 , res = 0.0005043058885168522
it# 170 , res = 0.0004897792339871674
it# 171 , res = 0.0004756710232694172
it# 172 , res = 0.00046196920299054354
it# 173 , res = 0.00044866206697756196
it# 174 , res = 0.0004357382462561546
it# 175 , res = 0.00042318669933805604
it# 176 , res = 0.0004109967027872326
it# 177 , res = 0.00039915784205896105
it# 178 , res = 0.00038766000260181835
it# 179 , res = 0.0003764933612158773
it# 180 , res = 0.0003656483776612427
it# 181 , res = 0.00035511578650648774
it# 182 , res = 0.00034488658921321346
it# 183 , res = 0.00033495204644805836
it# 184 , res = 0.00032530367061594227
it# 185 , res = 0.0003159332186092286
it# 186 , res = 0.0003068326847646606
it# 187 , res = 0.0002979942940243076
it# 188 , res = 0.0002894104952922365
it# 189 , res = 0.0002810739549839236
it# 190 , res = 0.00027297755076024626
it# 191 , res = 0.00026511436544309993
it# 192 , res = 0.00025747768110486054
it# 193 , res = 0.0002500609733296287
it# 194 , res = 0.0002428579056391199
it# 195 , res = 0.00023586232407863246
it# 196 , res = 0.00022906825195920068
it# 197 , res = 0.00022246988475259878
it# 198 , res = 0.0002160615851308381
it# 199 , res = 0.00020983787815067326
it# 200 , res = 0.00020379344657610972
it# 201 , res = 0.00019792312633504285
it# 202 , res = 0.00019222190210814327
it# 203 , res = 0.00018668490304326287
it# 204 , res = 0.00018130739859450478
it# 205 , res = 0.00017608479448029662
it# 206 , res = 0.00017101262875869426
it# 207 , res = 0.00016608656801500392
it# 208 , res = 0.00016130240365937125
it# 209 , res = 0.00015665604833185309
it# 210 , res = 0.00015214353240933828
it# 211 , res = 0.0001477610006151907
it# 212 , res = 0.00014350470872478939
it# 213 , res = 0.0001393710203666236
it# 214 , res = 0.00013535640391589525
it# 215 , res = 0.00013145742947715664
it# 216 , res = 0.00012767076595383356
it# 217 , res = 0.00012399317820263158
it# 218 , res = 0.0001204215242691085
it# 219 , res = 0.00011695275270380696
it# 220 , res = 0.00011358389995453007
it# 221 , res = 0.00011031208783554954
it# 222 , res = 0.0001071345210676783
it# 223 , res = 0.00010404848489032917
it# 224 , res = 0.00010105134274261876
it# 225 , res = 9.8140534010042e-05
it# 226 , res = 9.531357183737666e-05
it# 227 , res = 9.256804100396043e-05
it# 228 , res = 8.990159586011692e-05
it# 229 , res = 8.731195832316028e-05
it# 230 , res = 8.479691593091948e-05
it# 231 , res = 8.235431995221594e-05
it# 232 , res = 7.998208355025056e-05
it# 233 , res = 7.767818000014814e-05
it# 234 , res = 7.544064095728857e-05
it# 235 , res = 7.326755477571974e-05
it# 236 , res = 7.115706487515235e-05
it# 237 , res = 6.910736815415546e-05
it# 238 , res = 6.711671345031177e-05
it# 239 , res = 6.518340004392133e-05
it# 240 , res = 6.33057762046314e-05
it# 241 , res = 6.148223778103521e-05
it# 242 , res = 5.9711226829432764e-05
it# 243 , res = 5.799123028281277e-05
it# 244 , res = 5.632077865903524e-05
it# 245 , res = 5.469844480425775e-05
it# 246 , res = 5.312284267444853e-05
it# 247 , res = 5.159262615083581e-05
it# 248 , res = 5.010648788976716e-05
it# 249 , res = 4.86631582060384e-05
it# 250 , res = 4.726140398821711e-05
it# 251 , res = 4.5900027644755054e-05
it# 252 , res = 4.457786608101311e-05
it# 253 , res = 4.329378970559823e-05
it# 254 , res = 4.2046701465521205e-05
it# 255 , res = 4.083553590819315e-05
it# 256 , res = 3.965925827214525e-05
it# 257 , res = 3.8516863602228475e-05
it# 258 , res = 3.7407375891037503e-05
it# 259 , res = 3.63298472458565e-05
it# 260 , res = 3.528335707763255e-05
it# 261 , res = 3.4267011315743694e-05
it# 262 , res = 3.3279941643038925e-05
it# 263 , res = 3.232130475445511e-05
it# 264 , res = 3.139028163683184e-05
it# 265 , res = 3.0486076868672698e-05
it# 266 , res = 2.960791794071321e-05
it# 267 , res = 2.875505459638052e-05
it# 268 , res = 2.792675818997187e-05
it# 269 , res = 2.712232106501221e-05
it# 270 , res = 2.6341055948650285e-05
it# 271 , res = 2.5582295365879216e-05
it# 272 , res = 2.4845391067949227e-05
it# 273 , res = 2.4129713479176295e-05
it# 274 , res = 2.3434651159147662e-05
it# 275 , res = 2.2759610279838648e-05
it# 276 , res = 2.2104014118698835e-05
it# 277 , res = 2.14673025659638e-05
it# 278 , res = 2.0848931645813766e-05
it# 279 , res = 2.024837305180328e-05
it# 280 , res = 1.9665113695454848e-05
it# 281 , res = 1.909865526805474e-05
it# 282 , res = 1.8548513814734272e-05
it# 283 , res = 1.8014219321085125e-05
it# 284 , res = 1.7495315311398608e-05
it# 285 , res = 1.6991358459179657e-05
it# 286 , res = 1.6501918207681143e-05
it# 287 , res = 1.6026576402815935e-05
it# 288 , res = 1.5564926935398194e-05
it# 289 , res = 1.511657539412831e-05
it# 290 , res = 1.468113872912425e-05
it# 291 , res = 1.4258244924004207e-05
it# 292 , res = 1.3847532678949667e-05
it# 293 , res = 1.3448651100829162e-05
it# 294 , res = 1.3061259404692407e-05
it# 295 , res = 1.2685026621154093e-05
it# 296 , res = 1.2319631315426521e-05
it# 297 , res = 1.1964761311252656e-05
it# 298 , res = 1.162011342456552e-05
it# 299 , res = 1.1285393204730265e-05
it# 300 , res = 1.09603146831795e-05
it# 301 , res = 1.0644600128266956e-05
it# 302 , res = 1.0337979808542306e-05
it# 303 , res = 1.0040191762326845e-05
it# 304 , res = 9.750981574042902e-06
it# 305 , res = 9.47010215609978e-06
it# 306 , res = 9.19731353884433e-06
it# 307 , res = 8.932382664546687e-06
it# 308 , res = 8.675083189227387e-06
it# 309 , res = 8.425195288377078e-06
it# 310 , res = 8.182505469929897e-06
it# 311 , res = 7.946806391231912e-06
it# 312 , res = 7.717896682316495e-06
it# 313 , res = 7.4955807738521536e-06
it# 314 , res = 7.2796687295951765e-06
it# 315 , res = 7.0699760848554805e-06
it# 316 , res = 6.866323688273475e-06
it# 317 , res = 6.668537548903354e-06
it# 318 , res = 6.476448687787794e-06
it# 319 , res = 6.289892993426421e-06
it# 320 , res = 6.108711081771182e-06
it# 321 , res = 5.9327481595787554e-06
it# 322 , res = 5.761853892354251e-06
it# 323 , res = 5.595882276567648e-06
it# 324 , res = 5.434691513979811e-06
it# 325 , res = 5.27814389075226e-06
it# 326 , res = 5.126105660315429e-06
it# 327 , res = 4.978446928441887e-06
it# 328 , res = 4.835041542587461e-06
it# 329 , res = 4.695766984065512e-06
it# 330 , res = 4.560504263355307e-06
it# 331 , res = 4.429137818540112e-06
it# 332 , res = 4.301555416472654e-06
it# 333 , res = 4.177648056660295e-06
it# 334 , res = 4.057309878945815e-06
it# 335 , res = 3.940438071839337e-06
it# 336 , res = 3.826932785904063e-06
it# 337 , res = 3.716697047560082e-06
it# 338 , res = 3.6096366767627906e-06
it# 339 , res = 3.5056602062022766e-06
it# 340 , res = 3.404678803410105e-06
it# 341 , res = 3.306606194696833e-06
it# 342 , res = 3.2113585915398934e-06
it# 343 , res = 3.118854619068975e-06
it# 344 , res = 3.0290152460879916e-06
it# 345 , res = 2.941763718282942e-06
it# 346 , res = 2.8570254922225025e-06
it# 347 , res = 2.7747281716497917e-06
it# 348 , res = 2.6948014455257386e-06
it# 349 , res = 2.617177028338906e-06
it# 350 , res = 2.5417886016222495e-06
it# 351 , res = 2.4685717569773773e-06
it# 352 , res = 2.3974639415254237e-06
it# 353 , res = 2.3284044041299157e-06
it# 354 , res = 2.2613341437348796e-06
it# 355 , res = 2.196195858606616e-06
it# 356 , res = 2.1329338977876953e-06
it# 357 , res = 2.0714942132542484e-06
it# 358 , res = 2.0118243139889334e-06
it# 359 , res = 1.9538732207587686e-06
it# 360 , res = 1.8975914229620042e-06
it# 361 , res = 1.8429308362751796e-06
it# 362 , res = 1.7898447610616418e-06
it# 363 , res = 1.7382878434375732e-06
it# 364 , res = 1.6882160353151628e-06
it# 365 , res = 1.6395865580301005e-06
it# 366 , res = 1.5923578647415308e-06
it# 367 , res = 1.5464896054634081e-06
it# 368 , res = 1.5019425929123035e-06
it# 369 , res = 1.4586787676901176e-06
it# 370 , res = 1.4166611677204324e-06
it# 371 , res = 1.3758538950531329e-06
it# 372 , res = 1.3362220858499662e-06
it# 373 , res = 1.2977318806862105e-06
it# 374 , res = 1.2603503953027437e-06
it# 375 , res = 1.2240456926902698e-06
it# 376 , res = 1.1887867560698385e-06
it# 377 , res = 1.1545434616570472e-06
it# 378 , res = 1.1212865539070226e-06
it# 379 , res = 1.0889876194694675e-06
it# 380 , res = 1.0576190638783379e-06
it# 381 , res = 1.0271540870225026e-06
it# 382 , res = 9.975666614174736e-07
it# 383 , res = 9.688315088719502e-07
it# 384 , res = 9.409240794155068e-07
it# 385 , res = 9.138205302531375e-07
it# 386 , res = 8.874977055977215e-07
it# 387 , res = 8.619331163644192e-07
it# 388 , res = 8.371049213740427e-07
it# 389 , res = 8.129919085967485e-07
it# 390 , res = 7.895734770958201e-07
it# 391 , res = 7.668296193251278e-07
it# 392 , res = 7.447409038024537e-07
it# 393 , res = 7.232884593098371e-07
it# 394 , res = 7.024539576502217e-07
it# 395 , res = 6.822195989393709e-07
it# 396 , res = 6.625680959753308e-07
it# 397 , res = 6.43482659328939e-07
it# 398 , res = 6.249469833530954e-07
it# 399 , res = 6.069452320813209e-07
it# 400 , res = 5.894620256640086e-07
it# 401 , res = 5.724824272442413e-07
it# 402 , res = 5.559919303736298e-07
it# 403 , res = 5.399764465129149e-07
it# 404 , res = 5.244222923183113e-07
it# 405 , res = 5.093161794958599e-07
it# 406 , res = 4.946452018597014e-07
it# 407 , res = 4.803968254557409e-07
it# 408 , res = 4.6655887695056554e-07
it# 409 , res = 4.53119534078273e-07
it# 410 , res = 4.4006731458663145e-07
it# 411 , res = 4.2739106751887404e-07
it# 412 , res = 4.150799628862171e-07
it# 413 , res = 4.0312348268801676e-07
it# 414 , res = 3.915114117158501e-07
it# 415 , res = 3.802338294894691e-07
it# 416 , res = 3.692811007379039e-07
it# 417 , res = 3.586438679275675e-07
it# 418 , res = 3.4831304333306013e-07
it# 419 , res = 3.382798006040081e-07
it# 420 , res = 3.285355678568031e-07
it# 421 , res = 3.190720201673138e-07
it# 422 , res = 3.098810723489597e-07
it# 423 , res = 3.00954871913938e-07
it# 424 , res = 2.92285793057158e-07
it# 425 , res = 2.8386642902687236e-07
it# 426 , res = 2.7568958681308237e-07
it# 427 , res = 2.6774828065288935e-07
it# 428 , res = 2.600357256698156e-07
it# 429 , res = 2.5254533266927417e-07
it# 430 , res = 2.452707024310038e-07
it# 431 , res = 2.3820561968893028e-07
it# 432 , res = 2.3134404829197146e-07
it# 433 , res = 2.24680126216696e-07
it# 434 , res = 2.1820816005825902e-07
it# 435 , res = 2.119226203866321e-07
it# 436 , res = 2.0581813724652912e-07
it# 437 , res = 1.998894952923866e-07
it# 438 , res = 1.9413162931360107e-07
it# 439 , res = 1.8853962015966607e-07
it# 440 , res = 1.8310869005620007e-07
it# 441 , res = 1.778341992532099e-07
it# 442 , res = 1.727116414562915e-07
it# 443 , res = 1.677366403377984e-07
it# 444 , res = 1.6290494511532158e-07
it# 445 , res = 1.582124282727928e-07
it# 446 , res = 1.536550803714497e-07
it# 447 , res = 1.4922900808186128e-07
it# 448 , res = 1.4493042991789724e-07
it# 449 , res = 1.4075567333873202e-07
it# 450 , res = 1.3670117154515981e-07
it# 451 , res = 1.327634607095345e-07
it# 452 , res = 1.2893917672178983e-07
it# 453 , res = 1.252250520510957e-07
it# 454 , res = 1.2161791367541325e-07
it# 455 , res = 1.1811467987344487e-07
it# 456 , res = 1.147123574534058e-07
it# 457 , res = 1.1140803992375061e-07
it# 458 , res = 1.0819890390463315e-07
it# 459 , res = 1.0508220791689381e-07
it# 460 , res = 1.020552891726667e-07
it# 461 , res = 9.911556146170513e-08
it# 462 , res = 9.626051336249763e-08
it# 463 , res = 9.34877056737119e-08
it# 464 , res = 9.079476940135878e-08
it# 465 , res = 8.817940376585874e-08
it# 466 , res = 8.563937434023701e-08
it# 467 , res = 8.317251118871242e-08
it# 468 , res = 8.077670655436459e-08
it# 469 , res = 7.844991374192955e-08
it# 470 , res = 7.619014467713193e-08
it# 471 , res = 7.399546887305388e-08
it# 472 , res = 7.186401120385897e-08
it# 473 , res = 6.979395069784048e-08
it# 474 , res = 6.778351889053394e-08
it# 475 , res = 6.583099798747893e-08
it# 476 , res = 6.393471993273108e-08
it# 477 , res = 6.209306466042533e-08
it# 478 , res = 6.030445863425405e-08
it# 479 , res = 5.8567373947445594e-08
it# 480 , res = 5.688032634053713e-08
it# 481 , res = 5.524187448146306e-08
it# 482 , res = 5.3650618752398984e-08
it# 483 , res = 5.2105199573544157e-08
it# 484 , res = 5.060429649387206e-08
it# 485 , res = 4.914662733929086e-08
it# 486 , res = 4.773094666047237e-08
it# 487 , res = 4.6356044948394464e-08
it# 488 , res = 4.5020747732457e-08
it# 489 , res = 4.372391400955877e-08
it# 490 , res = 4.246443598175326e-08
it# 491 , res = 4.12412375047599e-08
it# 492 , res = 4.005327352550195e-08
it# 493 , res = 3.889952910586058e-08
it# 494 , res = 3.777901872432093e-08
it# 495 , res = 3.669078471195969e-08
it# 496 , res = 3.563389774371637e-08
it# 497 , res = 3.4607454444519045e-08
it# 498 , res = 3.361057831176517e-08
it# 499 , res = 3.264241734233499e-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)

WARNING: maxsteps is deprecated, use maxiter instead!
CG iteration 1, residual = 0.09414953738541333
CG iteration 2, residual = 0.12217425692967859
CG iteration 3, residual = 0.08601476799518139
CG iteration 4, residual = 0.05488212824534889
CG iteration 5, residual = 0.027459914932220164
CG iteration 6, residual = 0.01523316222565969
CG iteration 7, residual = 0.0067962047416292165
CG iteration 8, residual = 0.0028909296245792616
CG iteration 9, residual = 0.0010063644603491206
CG iteration 10, residual = 0.0003940448891274729
CG iteration 11, residual = 0.0001221088115779579
CG iteration 12, residual = 6.060285452707124e-05
CG iteration 13, residual = 2.827209240249563e-05
CG iteration 14, residual = 1.2739595841174722e-05
CG iteration 15, residual = 5.2941742436082805e-06
CG iteration 16, residual = 2.0472849339416365e-06
CG iteration 17, residual = 8.04952183815169e-07
CG iteration 18, residual = 3.5475789106359884e-07
CG iteration 19, residual = 1.29144633434917e-07
CG iteration 20, residual = 5.438554003499729e-08
CG iteration 21, residual = 1.7535915194415497e-08
CG iteration 22, residual = 5.84564970760941e-09
CG iteration 23, residual = 2.4957367825415158e-09
CG iteration 24, residual = 7.701980016761407e-10
CG iteration 25, residual = 3.1356081094504934e-10
CG iteration 26, residual = 1.0993771377265834e-10
CG iteration 27, residual = 3.912508493816494e-11
CG iteration 28, residual = 1.5718876754805824e-11
CG iteration 29, residual = 5.040743053318683e-12
CG iteration 30, residual = 1.9766256334138664e-12
CG iteration 31, residual = 5.671566563797397e-13
CG iteration 32, residual = 2.1155291673720475e-13
CG iteration 33, residual = 6.067815798067473e-14

BaseWebGuiScene

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

20.327559946070362

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 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)

34.93029206976627

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

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.99316
 0.997965
 0.999962
 1.35407
 1.82288
 1.85073
 1.96899
 1.98535
 1.99969

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

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.