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

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.0144743
 0.0565951
 0.0940742
 0.111584
 0.148079
 0.204136
 0.279083
 0.379019
 0.45319
 0.53846
 0.662266
 0.764709
 0.854406
 0.966403
 1.07186
 1.18706
 1.28597
 1.40124
 1.50818
 1.62481
 1.74747
 1.88402
  1.9919
 2.07749
 2.18043
 2.29565
 2.38514
 2.45388
 2.50384
 2.56324
 2.58875
  2.6348
 2.81085

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

max(lams)/min(lams)

194.19561167462462

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)

1678.5881226918377

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.05505568711828939
CG iteration 2, residual = 0.09357794070949911
CG iteration 3, residual = 0.11255359311920529
CG iteration 4, residual = 0.0883578866157642
CG iteration 5, residual = 0.09465752852240354
CG iteration 6, residual = 0.08675791614352733
CG iteration 7, residual = 0.09174564400414802
CG iteration 8, residual = 0.07950174974020313
CG iteration 9, residual = 0.05558110481441834
CG iteration 10, residual = 0.044799001641920515
CG iteration 11, residual = 0.03885061649304077
CG iteration 12, residual = 0.030408675499253898
CG iteration 13, residual = 0.02288741507787512
CG iteration 14, residual = 0.02051092486797348
CG iteration 15, residual = 0.015405453021112554
CG iteration 16, residual = 0.011893431227990001
CG iteration 17, residual = 0.010292148243983716
CG iteration 18, residual = 0.007888479867611086
CG iteration 19, residual = 0.005481703824765564
CG iteration 20, residual = 0.003612015168020746
CG iteration 21, residual = 0.0023840872403187567
CG iteration 22, residual = 0.0015057839927216509
CG iteration 23, residual = 0.0011613533060643182
CG iteration 24, residual = 0.000869420245521821
CG iteration 25, residual = 0.0005329860738693037
CG iteration 26, residual = 0.00037108669875231834
CG iteration 27, residual = 0.0002859501917016282
CG iteration 28, residual = 0.00022783224021485763
CG iteration 29, residual = 0.000166594323733906
CG iteration 30, residual = 0.00011944094829848384
CG iteration 31, residual = 8.147974841738306e-05
CG iteration 32, residual = 5.5020653070099096e-05
CG iteration 33, residual = 4.190089073735408e-05
CG iteration 34, residual = 3.369706727003576e-05
CG iteration 35, residual = 2.7766665129423223e-05
CG iteration 36, residual = 1.9736038657175716e-05
CG iteration 37, residual = 1.363994533705487e-05
CG iteration 38, residual = 1.0506189808377937e-05
CG iteration 39, residual = 8.983446326952995e-06
CG iteration 40, residual = 7.660521208331822e-06
CG iteration 41, residual = 5.801953313066017e-06
CG iteration 42, residual = 4.051197313069691e-06
CG iteration 43, residual = 2.9961283895352596e-06
CG iteration 44, residual = 2.3762446752543294e-06
CG iteration 45, residual = 1.8104198799132735e-06
CG iteration 46, residual = 1.2473066060491126e-06
CG iteration 47, residual = 8.123842440232735e-07
CG iteration 48, residual = 5.602213594389701e-07
CG iteration 49, residual = 4.297300852190805e-07
CG iteration 50, residual = 3.046923636128493e-07
CG iteration 51, residual = 1.9532592484907017e-07
CG iteration 52, residual = 1.3346842701507282e-07
CG iteration 53, residual = 9.387214810153153e-08
CG iteration 54, residual = 7.042267578742155e-08
CG iteration 55, residual = 5.0876599295554645e-08
CG iteration 56, residual = 3.6287644199056415e-08
CG iteration 57, residual = 2.4168652132655575e-08
CG iteration 58, residual = 1.6475086046502934e-08
CG iteration 59, residual = 1.0716935189887449e-08
CG iteration 60, residual = 6.890521019254731e-09
CG iteration 61, residual = 4.550452935048494e-09
CG iteration 62, residual = 3.135967176609285e-09
CG iteration 63, residual = 2.033380757491356e-09
CG iteration 64, residual = 1.331589593370284e-09
CG iteration 65, residual = 9.48052142950352e-10
CG iteration 66, residual = 7.191889610955045e-10
CG iteration 67, residual = 5.584212662146425e-10
CG iteration 68, residual = 3.963364255887182e-10
CG iteration 69, residual = 2.642865235172716e-10
CG iteration 70, residual = 1.7946264132246952e-10
CG iteration 71, residual = 1.1724755382242814e-10
CG iteration 72, residual = 7.742596210186727e-11
CG iteration 73, residual = 5.2131525221998164e-11
CG iteration 74, residual = 3.62410036123939e-11
CG iteration 75, residual = 2.3341859703488752e-11
CG iteration 76, residual = 1.5421560226282028e-11
CG iteration 77, residual = 1.0252664873785537e-11
CG iteration 78, residual = 6.927119148775129e-12
CG iteration 79, residual = 4.713860806170039e-12
CG iteration 80, residual = 3.1545465764504343e-12
CG iteration 81, residual = 1.9657493943078776e-12
CG iteration 82, residual = 1.2237497611515653e-12
CG iteration 83, residual = 8.299364158574324e-13
CG iteration 84, residual = 5.756065313208486e-13
CG iteration 85, residual = 4.107851399532092e-13
CG iteration 86, residual = 2.6545116923689193e-13
CG iteration 87, residual = 1.6948927421017946e-13
CG iteration 88, residual = 1.1040106281723313e-13
CG iteration 89, residual = 7.528580384049475e-14
CG iteration 90, residual = 5.3907824314899504e-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.08119734356602216
it# 1 , res = 0.07696557260572953
it# 2 , res = 0.07344471412858082
it# 3 , res = 0.07025007673116736
it# 4 , res = 0.06735911057304736
it# 5 , res = 0.06472122579225087
it# 6 , res = 0.062294223622381836
it# 7 , res = 0.06004424891557193
it# 8 , res = 0.05794435368653123
it# 9 , res = 0.05597310993397294
it# 10 , res = 0.05411339643696332
it# 11 , res = 0.05235143129498116
it# 12 , res = 0.05067602643031652
it# 13 , res = 0.049078019166521124
it# 14 , res = 0.047549839165887446
it# 15 , res = 0.04608517709824829
it# 16 , res = 0.04467872944058139
it# 17 , res = 0.043326000417520126
it# 18 , res = 0.04202314712044809
it# 19 , res = 0.04076685752716331
it# 20 , res = 0.03955425380652286
it# 21 , res = 0.038382815214659936
it# 22 , res = 0.037250316285503696
it# 23 , res = 0.03615477704179933
it# 24 , res = 0.035094422710862015
it# 25 , res = 0.03406765099651578
it# 26 , res = 0.03307300538722954
it# 27 , res = 0.03210915330716986
it# 28 , res = 0.03117486816799722
it# 29 , res = 0.030269014573674553
it# 30 , res = 0.029390536082149768
it# 31 , res = 0.02853844504667718
it# 32 , res = 0.02771181415333076
it# 33 , res = 0.026909769345597963
it# 34 , res = 0.02613148388613044
it# 35 , res = 0.025376173353034902
it# 36 , res = 0.02464309140603591
it# 37 , res = 0.02393152618837694
it# 38 , res = 0.023240797254965688
it# 39 , res = 0.02257025293720382
it# 40 , res = 0.02191926807110181
it# 41 , res = 0.021287242028414446
it# 42 , res = 0.02067359700122097
it# 43 , res = 0.02007777649909595
it# 44 , res = 0.019499244025139206
it# 45 , res = 0.018937481902962742
it# 46 , res = 0.01839199023151119
it# 47 , res = 0.01786228594851358
it# 48 , res = 0.01734790198658913
it# 49 , res = 0.016848386508685472
it# 50 , res = 0.016363302211716505
it# 51 , res = 0.015892225689077555
it# 52 , res = 0.015434746844209748
it# 53 , res = 0.01499046834862759
it# 54 , res = 0.014559005138852558
it# 55 , res = 0.014139983947550387
it# 56 , res = 0.01373304286488472
it# 57 , res = 0.013337830926694156
it# 58 , res = 0.012954007726595163
it# 59 , res = 0.012581243049534415
it# 60 , res = 0.012219216524660239
it# 61 , res = 0.011867617295679266
it# 62 , res = 0.01152614370711084
it# 63 , res = 0.011194503005064697
it# 64 , res = 0.010872411051340005
it# 65 , res = 0.010559592049799267
it# 66 , res = 0.010255778284096096
it# 67 , res = 0.009960709865945666
it# 68 , res = 0.009674134493220828
it# 69 , res = 0.009395807217235921
it# 70 , res = 0.009125490218650528
it# 71 , res = 0.00886295259148325
it# 72 , res = 0.008607970134778353
it# 73 , res = 0.008360325151510997
it# 74 , res = 0.008119806254357436
it# 75 , res = 0.007886208177987818
it# 76 , res = 0.007659331597570535
it# 77 , res = 0.007438982953202166
it# 78 , res = 0.007224974279999439
it# 79 , res = 0.007017123043610654
it# 80 , res = 0.006815251980920437
it# 81 , res = 0.006619188945737578
it# 82 , res = 0.006428766759272515
it# 83 , res = 0.006243823065218345
it# 84 , res = 0.006064200189265575
it# 85 , res = 0.005889745002886657
it# 86 , res = 0.00572030879124145
it# 87 , res = 0.005555747125053994
it# 88 , res = 0.0053959197363275336
it# 89 , res = 0.005240690397767275
it# 90 , res = 0.005089926805785061
it# 91 , res = 0.004943500466970432
it# 92 , res = 0.004801286587914051
it# 93 , res = 0.0046631639682756165
it# 94 , res = 0.004529014896994936
it# 95 , res = 0.004398725051543881
it# 96 , res = 0.004272183400128201
it# 97 , res = 0.004149282106743879
it# 98 , res = 0.004029916439004409
it# 99 , res = 0.003913984678649816
it# 100 , res = 0.0038013880346601174
it# 101 , res = 0.003692030558892565
it# 102 , res = 0.003585819064165283
it# 103 , res = 0.0034826630447189954
it# 104 , res = 0.0033824745989797455
it# 105 , res = 0.00328516835455907
it# 106 , res = 0.0031906613954232324
it# 107 , res = 0.003098873191168193
it# 108 , res = 0.003009725528337407
it# 109 , res = 0.002923142443723854
it# 110 , res = 0.0028390501595974916
it# 111 , res = 0.002757377020801111
it# 112 , res = 0.0026780534336611917
it# 113 , res = 0.002601011806660501
it# 114 , res = 0.0025261864928203598
it# 115 , res = 0.0024535137337449242
it# 116 , res = 0.0023829316052758204
it# 117 , res = 0.0023143799647147513
it# 118 , res = 0.0022478003995643862
it# 119 , res = 0.0021831361777465536
it# 120 , res = 0.002120332199253914
it# 121 , res = 0.002059334949193202
it# 122 , res = 0.0020000924521809734
it# 123 , res = 0.0019425542280515152
it# 124 , res = 0.00188667124884025
it# 125 , res = 0.0018323958970042887
it# 126 , res = 0.0017796819248471253
it# 127 , res = 0.001728484415108924
it# 128 , res = 0.0016787597426922676
it# 129 , res = 0.001630465537489084
it# 130 , res = 0.0015835606482767284
it# 131 , res = 0.0015380051076530215
it# 132 , res = 0.0014937600979809229
it# 133 , res = 0.0014507879183114146
it# 134 , res = 0.001409051952260031
it# 135 , res = 0.0013685166368059222
it# 136 , res = 0.0013291474319896562
it# 137 , res = 0.001290910791481944
it# 138 , res = 0.0012537741339998386
it# 139 , res = 0.0012177058155444591
it# 140 , res = 0.0011826751024382805
it# 141 , res = 0.0011486521451370662
it# 142 , res = 0.0011156079527963382
it# 143 , res = 0.0010835143685682465
it# 144 , res = 0.0010523440456102749
it# 145 , res = 0.001022070423782998
it# 146 , res = 0.0009926677070188294
it# 147 , res = 0.0009641108413425826
it# 148 , res = 0.0009363754935217705
it# 149 , res = 0.0009094380303344205
it# 150 , res = 0.0008832754984305099
it# 151 , res = 0.0008578656047740584
it# 152 , res = 0.0008331866976479217
it# 153 , res = 0.0008092177482045541
it# 154 , res = 0.000785938332547231
it# 155 , res = 0.0007633286143278866
it# 156 , res = 0.0007413693278441148
it# 157 , res = 0.0007200417616235552
it# 158 , res = 0.0006993277424800745
it# 159 , res = 0.0006792096200284055
it# 160 , res = 0.0006596702516446203
it# 161 , res = 0.0006406929878589201
it# 162 , res = 0.0006222616581690505
it# 163 , res = 0.00060436055726136
it# 164 , res = 0.0005869744316283719
it# 165 , res = 0.0005700884665716965
it# 166 , res = 0.0005536882735783817
it# 167 , res = 0.000537759878060613
it# 168 , res = 0.0005222897074479053
it# 169 , res = 0.0005072645796222377
it# 170 , res = 0.0004926716916858104
it# 171 , res = 0.0004784986090509398
it# 172 , res = 0.00046473325484572476
it# 173 , res = 0.00045136389962282847
it# 174 , res = 0.0004383791513645787
it# 175 , res = 0.00042576794577677067
it# 176 , res = 0.0004135195368599968
it# 177 , res = 0.00040162348775344846
it# 178 , res = 0.00039006966184131563
it# 179 , res = 0.0003788482141164032
it# 180 , res = 0.0003679495827897905
it# 181 , res = 0.0003573644811442805
it# 182 , res = 0.000347083889621145
it# 183 , res = 0.00033709904813420744
it# 184 , res = 0.0003274014486056376
it# 185 , res = 0.0003179828277161132
it# 186 , res = 0.0003088351598637657
it# 187 , res = 0.00029995065032554954
it# 188 , res = 0.0002913217286154457
it# 189 , res = 0.00028294104203361865
it# 190 , res = 0.00027480144940109265
it# 191 , res = 0.0002668960149748972
it# 192 , res = 0.00025921800253807676
it# 193 , res = 0.000251760869659655
it# 194 , res = 0.0002445182621203887
it# 195 , res = 0.00023748400849752467
it# 196 , res = 0.00023065211490659756
it# 197 , res = 0.0002240167598942669
it# 198 , res = 0.00021757228947727304
it# 199 , res = 0.00021131321232509877
it# 200 , res = 0.00020523419508077932
it# 201 , res = 0.000199330057816347
it# 202 , res = 0.00019359576961870548
it# 203 , res = 0.00018802644430348806
it# 204 , res = 0.00018261733625069224
it# 205 , res = 0.00017736383636165128
it# 206 , res = 0.00017226146813114382
it# 207 , res = 0.00016730588383317808
it# 208 , res = 0.00016249286081648773
it# 209 , res = 0.00015781829790597314
it# 210 , res = 0.0001532782119084258
it# 211 , res = 0.00014886873421832
it# 212 , res = 0.0001445861075214198
it# 213 , res = 0.00014042668259322085
it# 214 , res = 0.00013638691518935573
it# 215 , res = 0.00013246336302542817
it# 216 , res = 0.00012865268284443935
it# 217 , res = 0.00012495162756690073
it# 218 , res = 0.00012135704352520126
it# 219 , res = 0.0001178658677755353
it# 220 , res = 0.00011447512548853386
it# 221 , res = 0.00011118192741384818
it# 222 , res = 0.00010798346741852885
it# 223 , res = 0.00010487702009640282
it# 224 , res = 0.0001018599384444072
it# 225 , res = 9.892965160873351e-05
it# 226 , res = 9.608366269287981e-05
it# 227 , res = 9.331954663063205e-05
it# 228 , res = 9.063494812001635e-05
it# 229 , res = 8.802757961551883e-05
it# 230 , res = 8.5495219379298e-05
it# 231 , res = 8.303570958838247e-05
it# 232 , res = 8.064695449507518e-05
it# 233 , res = 7.832691864221051e-05
it# 234 , res = 7.607362512815404e-05
it# 235 , res = 7.38851539222216e-05
it# 236 , res = 7.175964022889903e-05
it# 237 , res = 6.96952728988565e-05
it# 238 , res = 6.769029288535838e-05
it# 239 , res = 6.574299174563177e-05
it# 240 , res = 6.385171018519757e-05
it# 241 , res = 6.201483664378394e-05
it# 242 , res = 6.0230805921899757e-05
it# 243 , res = 5.849809784785605e-05
it# 244 , res = 5.681523598151747e-05
it# 245 , res = 5.5180786357110256e-05
it# 246 , res = 5.3593356260585004e-05
it# 247 , res = 5.205159304330611e-05
it# 248 , res = 5.0554182969542835e-05
it# 249 , res = 4.909985009658866e-05
it# 250 , res = 4.768735518807469e-05
it# 251 , res = 4.631549465744327e-05
it# 252 , res = 4.498309954286539e-05
it# 253 , res = 4.36890345111116e-05
it# 254 , res = 4.243219688964488e-05
it# 255 , res = 4.121151572765292e-05
it# 256 , res = 4.002595088329978e-05
it# 257 , res = 3.8874492136984886e-05
it# 258 , res = 3.775615833118976e-05
it# 259 , res = 3.6669996534226535e-05
it# 260 , res = 3.561508122802995e-05
it# 261 , res = 3.45905135193489e-05
it# 262 , res = 3.359542037477192e-05
it# 263 , res = 3.262895387554767e-05
it# 264 , res = 3.1690290496133045e-05
it# 265 , res = 3.077863040160895e-05
it# 266 , res = 2.9893196766707928e-05
it# 267 , res = 2.9033235114061888e-05
it# 268 , res = 2.8198012670473635e-05
it# 269 , res = 2.7386817743212115e-05
it# 270 , res = 2.6598959113254317e-05
it# 271 , res = 2.5833765446684974e-05
it# 272 , res = 2.5090584722258722e-05
it# 273 , res = 2.436878367582749e-05
it# 274 , res = 2.3667747261157187e-05
it# 275 , res = 2.2986878125337466e-05
it# 276 , res = 2.2325596099946585e-05
it# 277 , res = 2.1683337706852214e-05
it# 278 , res = 2.1059555677893493e-05
it# 279 , res = 2.0453718488728933e-05
it# 280 , res = 1.9865309905590275e-05
it# 281 , res = 1.9293828545663746e-05
it# 282 , res = 1.873878745013851e-05
it# 283 , res = 1.8199713668585545e-05
it# 284 , res = 1.7676147856358916e-05
it# 285 , res = 1.716764388322107e-05
it# 286 , res = 1.6673768453234718e-05
it# 287 , res = 1.6194100735217993e-05
it# 288 , res = 1.572823200465937e-05
it# 289 , res = 1.5275765294734605e-05
it# 290 , res = 1.483631505894191e-05
it# 291 , res = 1.440950684171057e-05
it# 292 , res = 1.399497696002339e-05
it# 293 , res = 1.3592372193153448e-05
it# 294 , res = 1.3201349481631834e-05
it# 295 , res = 1.2821575635203542e-05
it# 296 , res = 1.2452727048636231e-05
it# 297 , res = 1.2094489426362285e-05
it# 298 , res = 1.1746557514090899e-05
it# 299 , res = 1.1408634839154064e-05
it# 300 , res = 1.1080433457791323e-05
it# 301 , res = 1.0761673709773286e-05
it# 302 , res = 1.0452083980019454e-05
it# 303 , res = 1.0151400467229204e-05
it# 304 , res = 9.859366958950267e-06
it# 305 , res = 9.575734613805014e-06
it# 306 , res = 9.300261748557444e-06
it# 307 , res = 9.032713632887942e-06
it# 308 , res = 8.772862289351362e-06
it# 309 , res = 8.52048629847047e-06
it# 310 , res = 8.275370611002951e-06
it# 311 , res = 8.037306363685937e-06
it# 312 , res = 7.806090702165665e-06
it# 313 , res = 7.581526607646251e-06
it# 314 , res = 7.3634227293118e-06
it# 315 , res = 7.1515932208259704e-06
it# 316 , res = 6.945857582274746e-06
it# 317 , res = 6.74604050629952e-06
it# 318 , res = 6.551971729182755e-06
it# 319 , res = 6.363485884694585e-06
it# 320 , res = 6.180422364198965e-06
it# 321 , res = 6.002625179413754e-06
it# 322 , res = 5.8299428292598825e-06
it# 323 , res = 5.66222817129389e-06
it# 324 , res = 5.499338296012727e-06
it# 325 , res = 5.341134404672039e-06
it# 326 , res = 5.18748169216689e-06
it# 327 , res = 5.038249230868083e-06
it# 328 , res = 4.893309859858835e-06
it# 329 , res = 4.752540076507952e-06
it# 330 , res = 4.615819930754405e-06
it# 331 , res = 4.483032923469462e-06
it# 332 , res = 4.354065906944336e-06
it# 333 , res = 4.228808988352778e-06
it# 334 , res = 4.1071554363812995e-06
it# 335 , res = 3.989001590112266e-06
it# 336 , res = 3.8742467705605776e-06
it# 337 , res = 3.762793195206251e-06
it# 338 , res = 3.6545458947066984e-06
it# 339 , res = 3.549412631384941e-06
it# 340 , res = 3.44730382126366e-06
it# 341 , res = 3.348132457516503e-06
it# 342 , res = 3.2518140361437767e-06
it# 343 , res = 3.1582664843194334e-06
it# 344 , res = 3.067410090307078e-06
it# 345 , res = 2.9791674351533595e-06
it# 346 , res = 2.893463327586205e-06
it# 347 , res = 2.8102247390729675e-06
it# 348 , res = 2.7293807420582952e-06
it# 349 , res = 2.6508624493284514e-06
it# 350 , res = 2.5746029555447406e-06
it# 351 , res = 2.5005372800047275e-06
it# 352 , res = 2.428602311397095e-06
it# 353 , res = 2.3587367537913716e-06
it# 354 , res = 2.2908810751191085e-06
it# 355 , res = 2.224977455242853e-06
it# 356 , res = 2.1609697379225345e-06
it# 357 , res = 2.0988033820302652e-06
it# 358 , res = 2.038425415820471e-06
it# 359 , res = 1.979784391192782e-06
it# 360 , res = 1.9228303400422064e-06
it# 361 , res = 1.86751473182168e-06
it# 362 , res = 1.8137904323959279e-06
it# 363 , res = 1.7616116631852195e-06
it# 364 , res = 1.710933962500022e-06
it# 365 , res = 1.6617141480405906e-06
it# 366 , res = 1.613910279576796e-06
it# 367 , res = 1.567481623597086e-06
it# 368 , res = 1.522388617989059e-06
it# 369 , res = 1.4785928391996155e-06
it# 370 , res = 1.4360569688617811e-06
it# 371 , res = 1.3947447621082154e-06
it# 372 , res = 1.3546210169659847e-06
it# 373 , res = 1.3156515438445063e-06
it# 374 , res = 1.27780313688417e-06
it# 375 , res = 1.2410435455193557e-06
it# 376 , res = 1.20534144690881e-06
it# 377 , res = 1.1706664192937224e-06
it# 378 , res = 1.1369889161115048e-06
it# 379 , res = 1.1042802408675386e-06
it# 380 , res = 1.072512522340481e-06
it# 381 , res = 1.0416586912721046e-06
it# 382 , res = 1.0116924572609038e-06
it# 383 , res = 9.825882860519578e-07
it# 384 , res = 9.543213779212043e-07
it# 385 , res = 9.268676467172414e-07
it# 386 , res = 9.002036990374693e-07
it# 387 , res = 8.743068146792112e-07
it# 388 , res = 8.49154926748328e-07
it# 389 , res = 8.247266035600695e-07
it# 390 , res = 8.010010294992781e-07
it# 391 , res = 7.779579880918734e-07
it# 392 , res = 7.555778445297398e-07
it# 393 , res = 7.338415284689586e-07
it# 394 , res = 7.127305185586999e-07
it# 395 , res = 6.922268259952688e-07
it# 396 , res = 6.723129797348643e-07
it# 397 , res = 6.52972011145823e-07
it# 398 , res = 6.341874397022707e-07
it# 399 , res = 6.159432591872167e-07
it# 400 , res = 5.982239236226157e-07
it# 401 , res = 5.810143344804237e-07
it# 402 , res = 5.642998274814528e-07
it# 403 , res = 5.48066160098755e-07
it# 404 , res = 5.322994998415674e-07
it# 405 , res = 5.169864116029918e-07
it# 406 , res = 5.021138473580336e-07
it# 407 , res = 4.876691342245933e-07
it# 408 , res = 4.736399636878286e-07
it# 409 , res = 4.6001438155557564e-07
it# 410 , res = 4.467807776072234e-07
it# 411 , res = 4.3392787531824345e-07
it# 412 , res = 4.214447228242701e-07
it# 413 , res = 4.093206832812781e-07
it# 414 , res = 3.975454257726523e-07
it# 415 , res = 3.8610891653779513e-07
it# 416 , res = 3.7500141055927573e-07
it# 417 , res = 3.642134431479294e-07
it# 418 , res = 3.5373582178470617e-07
it# 419 , res = 3.4355961871592467e-07
it# 420 , res = 3.336761625887141e-07
it# 421 , res = 3.240770319307513e-07
it# 422 , res = 3.1475404704815144e-07
it# 423 , res = 3.056992640780088e-07
it# 424 , res = 2.9690496729704544e-07
it# 425 , res = 2.8836366314048594e-07
it# 426 , res = 2.8006807353371656e-07
it# 427 , res = 2.720111297267254e-07
it# 428 , res = 2.641859666903443e-07
it# 429 , res = 2.565859162847454e-07
it# 430 , res = 2.4920450273121613e-07
it# 431 , res = 2.4203543627645496e-07
it# 432 , res = 2.3507260804371125e-07
it# 433 , res = 2.283100851126359e-07
it# 434 , res = 2.2174210510631445e-07
it# 435 , res = 2.1536307145382084e-07
it# 436 , res = 2.0916754861423567e-07
it# 437 , res = 2.031502573761954e-07
it# 438 , res = 1.9730607042775195e-07
it# 439 , res = 1.9163000789331896e-07
it# 440 , res = 1.861172332360944e-07
it# 441 , res = 1.807630489240887e-07
it# 442 , res = 1.7556289285271858e-07
it# 443 , res = 1.7051233376704244e-07
it# 444 , res = 1.6560706824401367e-07
it# 445 , res = 1.608429163935476e-07
it# 446 , res = 1.5621581875407193e-07
it# 447 , res = 1.5172183250744885e-07
it# 448 , res = 1.4735712843300102e-07
it# 449 , res = 1.4311798736849158e-07
it# 450 , res = 1.3900079692954408e-07
it# 451 , res = 1.3500204922810624e-07
it# 452 , res = 1.3111833662057514e-07
it# 453 , res = 1.2734634987053666e-07
it# 454 , res = 1.2368287493932198e-07
it# 455 , res = 1.20124790127106e-07
it# 456 , res = 1.1666906358766115e-07
it# 457 , res = 1.1331275074173177e-07
it# 458 , res = 1.100529916738983e-07
it# 459 , res = 1.0688700855493636e-07
it# 460 , res = 1.0381210399306976e-07
it# 461 , res = 1.0082565762606738e-07
it# 462 , res = 9.79251248214122e-08
it# 463 , res = 9.510803392718951e-08
it# 464 , res = 9.23719845903885e-08
it# 465 , res = 8.971464531215289e-08
it# 466 , res = 8.71337519116926e-08
it# 467 , res = 8.462710516719772e-08
it# 468 , res = 8.219256903327963e-08
it# 469 , res = 7.982806920203638e-08
it# 470 , res = 7.753159096448283e-08
it# 471 , res = 7.530117723681155e-08
it# 472 , res = 7.313492769685997e-08
it# 473 , res = 7.103099633872838e-08
it# 474 , res = 6.898759056141798e-08
it# 475 , res = 6.700296896013302e-08
it# 476 , res = 6.507544066154603e-08
it# 477 , res = 6.3203363048093e-08
it# 478 , res = 6.138514099161166e-08
it# 479 , res = 5.961922529860007e-08
it# 480 , res = 5.790411099885679e-08
it# 481 , res = 5.623833684638744e-08
it# 482 , res = 5.462048343618774e-08
it# 483 , res = 5.304917195076123e-08
it# 484 , res = 5.1523063722207655e-08
it# 485 , res = 5.0040858239050835e-08
it# 486 , res = 4.86012926604021e-08
it# 487 , res = 4.7203140117398695e-08
it# 488 , res = 4.584520943575234e-08
it# 489 , res = 4.4526343338363224e-08
it# 490 , res = 4.3245418146443594e-08
it# 491 , res = 4.2001342484177105e-08
it# 492 , res = 4.079305600114823e-08
it# 493 , res = 3.961952933918804e-08
it# 494 , res = 3.84797624760545e-08
it# 495 , res = 3.737278413472251e-08
it# 496 , res = 3.629765116928812e-08
it# 497 , res = 3.5253447403659214e-08
it# 498 , res = 3.4239283263301304e-08
it# 499 , res = 3.32542941094594e-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.09386193808775962
CG iteration 2, residual = 0.12289904446040574
CG iteration 3, residual = 0.08558846385984858
CG iteration 4, residual = 0.053786375091255334
CG iteration 5, residual = 0.02799513368494736
CG iteration 6, residual = 0.014969155114208049
CG iteration 7, residual = 0.006689317102568868
CG iteration 8, residual = 0.0028248604314505604
CG iteration 9, residual = 0.0010332508846612996
CG iteration 10, residual = 0.00040318420091610234
CG iteration 11, residual = 0.000126153484125304
CG iteration 12, residual = 5.134068138544936e-05
CG iteration 13, residual = 1.86267529855379e-05
CG iteration 14, residual = 7.86571512990075e-06
CG iteration 15, residual = 3.310916766942548e-06
CG iteration 16, residual = 1.5461433004860725e-06
CG iteration 17, residual = 6.174156961107636e-07
CG iteration 18, residual = 3.2008170496573533e-07
CG iteration 19, residual = 1.243696230710687e-07
CG iteration 20, residual = 5.906688350434197e-08
CG iteration 21, residual = 2.207041609152557e-08
CG iteration 22, residual = 7.324210399036425e-09
CG iteration 23, residual = 3.0122986220535976e-09
CG iteration 24, residual = 8.584333680651104e-10
CG iteration 25, residual = 3.6570806010889224e-10
CG iteration 26, residual = 1.2692446936158052e-10
CG iteration 27, residual = 4.461371781815046e-11
CG iteration 28, residual = 1.975092094382798e-11
CG iteration 29, residual = 5.5343243734924135e-12
CG iteration 30, residual = 2.307643414710973e-12
CG iteration 31, residual = 6.724334725524573e-13
CG iteration 32, residual = 1.9582752398605173e-13
CG iteration 33, residual = 7.628337311332365e-14

BaseWebGuiScene

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

20.294597459170067

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)

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

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

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

35.59893732672581

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

Add a coarse grid correction

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

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

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

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

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

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

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

The matrix coarsepre below represents

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

coarsepre = a.mat.Inverse(vertexdofs)

This matrix can be used for coarse grid correction.

Pitfall!

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

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

       1

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

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

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

Additive two-grid preconditioner

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

twogrid = coarsepre + blockgs

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

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

 0.993322
 0.998163
 0.999953
 1.34169
 1.80171
 1.82729
 1.95319
 1.97402
 1.99976

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

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.81817
 0.846572
 0.896024
 0.932504
 1.00024
  1.1144
 1.36948
 1.43102
 1.54233
 1.65822
 1.76373
  1.8514
 1.90737
 1.97204
 1.99987

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