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

2.1.2 Building blocks for programming preconditioners

import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
%gui tk
from ngsolve.la import EigenValues_Preconditioner
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
Draw (mesh)

fes = H1(mesh, order=3, dirichlet="left|bottom")
u, v = fes.TnT()
a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v) + u*v)
a.Assemble()
f = LinearForm(fes)
f += SymbolicLFI(1*v)
f.Assemble()
gfu = GridFunction(fes)

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 & I \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.0146975
 0.0567568
 0.0954973
 0.112599
 0.149013
 0.198016
 0.266185
 0.354506
 0.452415
 0.535142
 0.63906
 0.751009
 0.836039
 0.932691
 1.04077
 1.16473
 1.24579
 1.36685
 1.47433
 1.58784
 1.70789
 1.83511
 1.92562
 2.02092
 2.12651
 2.20727
 2.27824
 2.38882
 2.44463
 2.47644
 2.51031
 2.54019
  2.8097

An estimate of the condition number

\[\kappa = \frac{\lambda_{\text{max}} }{ \lambda_{\text{min}}}\]

is therefore given as follows:

max(lams)/min(lams)

191.168278526202

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)

1585.7752877045514

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)

it =  0  err =  0.05522230048185231
it =  1  err =  0.09382744174195908
it =  2  err =  0.10587866002039004
it =  3  err =  0.09033871153608707
it =  4  err =  0.10086392272140011
it =  5  err =  0.08517119355086304
it =  6  err =  0.08906575044224518
it =  7  err =  0.07684062913929467
it =  8  err =  0.056047175061895105
it =  9  err =  0.04472071103380234
it =  10  err =  0.03419741252004591
it =  11  err =  0.02682484049645099
it =  12  err =  0.022322983522404753
it =  13  err =  0.020056199326844725
it =  14  err =  0.014372559087884462
it =  15  err =  0.010772925966275729
it =  16  err =  0.009633149024762148
it =  17  err =  0.0073015321836821315
it =  18  err =  0.004811737390491996
it =  19  err =  0.0032534005521665317
it =  20  err =  0.001902580694234345
it =  21  err =  0.0012857479602116576
it =  22  err =  0.0010097069801732902
it =  23  err =  0.0006971063184209084
it =  24  err =  0.0004415360534270248
it =  25  err =  0.00031313626798158353
it =  26  err =  0.0002459345482987088
it =  27  err =  0.00018102122685671873
it =  28  err =  0.00012566996433136174
it =  29  err =  9.081609378020966e-05
it =  30  err =  5.920492010745938e-05
it =  31  err =  4.0326988735642384e-05
it =  32  err =  3.115865171601253e-05
it =  33  err =  2.2328835017506844e-05
it =  34  err =  1.63735443518241e-05
it =  35  err =  1.1514375356105962e-05
it =  36  err =  6.8946264410375216e-06
it =  37  err =  4.683506893773661e-06
it =  38  err =  3.5825089011129006e-06
it =  39  err =  2.859816977428613e-06
it =  40  err =  1.8589434122159337e-06
it =  41  err =  1.2888774288482613e-06
it =  42  err =  9.608284726114791e-07
it =  43  err =  8.138518532407976e-07
it =  44  err =  6.559651432751399e-07
it =  45  err =  5.088391213652928e-07
it =  46  err =  3.7554232878879976e-07
it =  47  err =  2.866653641968371e-07
it =  48  err =  2.3813379098727406e-07
it =  49  err =  1.784776125351117e-07
it =  50  err =  1.2713130281505307e-07
it =  51  err =  8.879135461994969e-08
it =  52  err =  7.189655596299818e-08
it =  53  err =  5.3603595340687565e-08
it =  54  err =  3.682678853570715e-08
it =  55  err =  2.5543483270413777e-08
it =  56  err =  1.7986625269364435e-08
it =  57  err =  1.2658159613085641e-08
it =  58  err =  8.74794105831229e-09
it =  59  err =  6.004571771303868e-09
it =  60  err =  4.240187119793316e-09
it =  61  err =  2.9124607839179425e-09
it =  62  err =  1.95454042937834e-09
it =  63  err =  1.28442596904015e-09
it =  64  err =  8.403738519702617e-10
it =  65  err =  5.901727212281454e-10
it =  66  err =  4.443364083343974e-10
it =  67  err =  3.1599407725917446e-10
it =  68  err =  2.0604768404034608e-10
it =  69  err =  1.2843022850328627e-10
it =  70  err =  9.063493172245921e-11
it =  71  err =  6.766709792706894e-11
it =  72  err =  4.746461848492383e-11
it =  73  err =  2.7809057936360353e-11
it =  74  err =  1.709795688455029e-11
it =  75  err =  1.176897616349203e-11
it =  76  err =  8.51804968461161e-12
it =  77  err =  6.277337591794344e-12
it =  78  err =  3.9132090932737885e-12
it =  79  err =  2.2135436745856382e-12
it =  80  err =  1.3679858558673982e-12
it =  81  err =  9.755276306296867e-13
it =  82  err =  7.137682424651503e-13
it =  83  err =  4.4876451218901374e-13
it =  84  err =  2.8581994514363945e-13
it =  85  err =  1.8635801952125264e-13
it =  86  err =  1.3655897098879978e-13
it =  87  err =  9.36583145099722e-14
it =  88  err =  6.003735062164544e-14

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.08192679624402945
it# 1 , res = 0.07726762555864525
it# 2 , res = 0.07377040495808554
it# 3 , res = 0.07062220895606912
it# 4 , res = 0.06777582421551728
it# 5 , res = 0.06516770753241384
it# 6 , res = 0.06275317789329397
it# 7 , res = 0.06050089690216452
it# 8 , res = 0.058387526692409986
it# 9 , res = 0.056394947347621995
it# 10 , res = 0.0545086993995133
it# 11 , res = 0.05271702122673418
it# 12 , res = 0.05101020777328185
it# 13 , res = 0.04938015983450581
it# 14 , res = 0.04782005439140267
it# 15 , res = 0.04632409596387132
it# 16 , res = 0.04488732463071546
it# 17 , res = 0.04350546518198168
it# 18 , res = 0.042174807039847335
it# 19 , res = 0.04089210774233109
it# 20 , res = 0.0396545147959703
it# 21 , res = 0.03845950204060288
it# 22 , res = 0.03730481759400025
it# 23 , res = 0.0361884411063836
it# 24 , res = 0.03510854854333058
it# 25 , res = 0.03406348308437501
it# 26 , res = 0.033051731008142335
it# 27 , res = 0.032071901655931344
it# 28 , res = 0.03112271073987291
it# 29 , res = 0.030202966400259202
it# 30 , res = 0.029311557527407415
it# 31 , res = 0.028447443952510115
it# 32 , res = 0.02760964818388455
it# 33 , res = 0.02679724842336503
it# 34 , res = 0.026009372645006548
it# 35 , res = 0.025245193556920734
it# 36 , res = 0.024503924298637167
it# 37 , res = 0.02378481475221903
it# 38 , res = 0.023087148366535692
it# 39 , res = 0.0224102394114777
it# 40 , res = 0.0217534305931838
it# 41 , res = 0.021116090973105302
it# 42 , res = 0.020497614143416644
it# 43 , res = 0.0198974166192704
it# 44 , res = 0.01931493641498571
it# 45 , res = 0.018749631776714237
it# 46 , res = 0.018200980048633027
it# 47 , res = 0.017668476653454283
it# 48 , res = 0.01715163417113754
it# 49 , res = 0.016649981502262958
it# 50 , res = 0.016163063104665278
it# 51 , res = 0.01569043829370675
it# 52 , res = 0.015231680598052692
it# 53 , res = 0.014786377164053767
it# 54 , res = 0.014354128202872912
it# 55 , res = 0.013934546475367345
it# 56 , res = 0.013527256810457532
it# 57 , res = 0.013131895653335023
it# 58 , res = 0.012748110640368406
it# 59 , res = 0.012375560198007516
it# 60 , res = 0.012013913163348851
it# 61 , res = 0.011662848424338797
it# 62 , res = 0.011322054577854098
it# 63 , res = 0.010991229604123473
it# 64 , res = 0.01067008055614577
it# 65 , res = 0.010358323262924123
it# 66 , res = 0.010055682045473519
it# 67 , res = 0.009761889444681583
it# 68 , res = 0.009476685960204578
it# 69 , res = 0.009199819799668862
it# 70 , res = 0.008931046637527158
it# 71 , res = 0.008670129382983446
it# 72 , res = 0.00841683795645999
it# 73 , res = 0.008170949074129685
it# 74 , res = 0.007932246040080845
it# 75 , res = 0.007700518545721331
it# 76 , res = 0.007475562476061317
it# 77 , res = 0.0072571797225452084
it# 78 , res = 0.0070451780021285534
it# 79 , res = 0.006839370682320229
it# 80 , res = 0.006639576611929876
it# 81 , res = 0.006445619957279532
it# 82 , res = 0.006257330043655511
it# 83 , res = 0.006074541201789725
it# 84 , res = 0.005897092619174911
it# 85 , res = 0.005724828196028781
it# 86 , res = 0.005557596405733547
it# 87 , res = 0.005395250159586561
it# 88 , res = 0.005237646675707954
it# 89 , res = 0.005084647351957455
it# 90 , res = 0.004936117642721644
it# 91 , res = 0.004791926939438961
it# 92 , res = 0.004651948454736169
it# 93 , res = 0.004516059110055665
it# 94 , res = 0.004384139426658663
it# 95 , res = 0.004256073419894378
it# 96 , res = 0.004131748496628999
it# 97 , res = 0.004011055355734409
it# 98 , res = 0.003893887891538068
it# 99 , res = 0.0037801431001420736
it# 100 , res = 0.003669720988521344
it# 101 , res = 0.003562524486313286
it# 102 , res = 0.003458459360217854
it# 103 , res = 0.0033574341309254496
it# 104 , res = 0.003259359992497481
it# 105 , res = 0.0031641507341227644
it# 106 , res = 0.0030717226641792865
it# 107 , res = 0.002981994536530998
it# 108 , res = 0.0028948874789913454
it# 109 , res = 0.002810324923890657
it# 110 , res = 0.0027282325406816797
it# 111 , res = 0.002648538170524317
it# 112 , res = 0.0025711717627896324
it# 113 , res = 0.002496065313426207
it# 114 , res = 0.0024231528051331455
it# 115 , res = 0.0023523701492869746
it# 116 , res = 0.0022836551295689427
it# 117 , res = 0.0022169473472448743
it# 118 , res = 0.0021521881680458556
it# 119 , res = 0.002089320670604589
it# 120 , res = 0.002028289596400842
it# 121 , res = 0.0019690413011714345
it# 122 , res = 0.0019115237077420857
it# 123 , res = 0.0018556862602383012
it# 124 , res = 0.0018014798796365335
it# 125 , res = 0.001748856920614219
it# 126 , res = 0.0016977711296619743
it# 127 , res = 0.0016481776044206479
it# 128 , res = 0.001600032754206714
it# 129 , res = 0.0015532942616918012
it# 130 , res = 0.001507921045702434
it# 131 , res = 0.0014638732251068592
it# 132 , res = 0.0014211120837572047
it# 133 , res = 0.0013796000364567705
it# 134 , res = 0.001339300595921035
it# 135 , res = 0.0013001783407047141
it# 136 , res = 0.0012621988840657815
it# 137 , res = 0.001225328843739243
it# 138 , res = 0.0011895358125941567
it# 139 , res = 0.0011547883301479946
it# 140 , res = 0.0011210558549129947
it# 141 , res = 0.0010883087375511042
it# 142 , res = 0.0010565181948121436
it# 143 , res = 0.0010256562842346165
it# 144 , res = 0.0009956958795847124
it# 145 , res = 0.0009666106470129384
it# 146 , res = 0.0009383750219080032
it# 147 , res = 0.0009109641864258245
it# 148 , res = 0.0008843540476760172
it# 149 , res = 0.0008585212165446786
it# 150 , res = 0.0008334429871366016
it# 151 , res = 0.000809097316817568
it# 152 , res = 0.0007854628068396603
it# 153 , res = 0.0007625186835328571
it# 154 , res = 0.000740244780045596
it# 155 , res = 0.0007186215186194956
it# 156 , res = 0.0006976298933805504
it# 157 , res = 0.0006772514536345653
it# 158 , res = 0.0006574682876495266
it# 159 , res = 0.0006382630069117587
it# 160 , res = 0.0006196187308428548
it# 161 , res = 0.0006015190719613711
it# 162 , res = 0.0005839481214803695
it# 163 , res = 0.0005668904353229085
it# 164 , res = 0.0005503310205485652
it# 165 , res = 0.0005342553221744136
it# 166 , res = 0.0005186492103826856
it# 167 , res = 0.0005034989681007053
it# 168 , res = 0.000488791278944353
it# 169 , res = 0.0004745132155138423
it# 170 , res = 0.0004606522280309772
it# 171 , res = 0.00044719613330827107
it# 172 , res = 0.0004341331040406897
it# 173 , res = 0.0004214516584100991
it# 174 , res = 0.0004091406499930594
it# 175 , res = 0.0003971892579639203
it# 176 , res = 0.0003855869775833213
it# 177 , res = 0.0003743236109654287
it# 178 , res = 0.000363389258114413
it# 179 , res = 0.00035277430822268883
it# 180 , res = 0.0003424694312235002
it# 181 , res = 0.00033246556959025305
it# 182 , res = 0.0003227539303752261
it# 183 , res = 0.00031332597748125803
it# 184 , res = 0.0003041734241585562
it# 185 , res = 0.00029528822572135686
it# 186 , res = 0.000286662572476714
it# 187 , res = 0.0002782888828602794
it# 188 , res = 0.00027015979677244453
it# 189 , res = 0.00026226816910930074
it# 190 , res = 0.00025460706348211725
it# 191 , res = 0.0002471697461206412
it# 192 , res = 0.00023994967995462916
it# 193 , res = 0.000232940518867916
it# 194 , res = 0.00022613610212036573
it# 195 , res = 0.00021953044893315316
it# 196 , res = 0.00021311775323138693
it# 197 , res = 0.0002068923785416417
it# 198 , res = 0.00020084885303698165
it# 199 , res = 0.0001949818647277607
it# 200 , res = 0.00018928625679283826
it# 201 , res = 0.00018375702304660783
it# 202 , res = 0.00017838930353898983
it# 203 , res = 0.0001731783802837913
it# 204 , res = 0.0001681196731116986
it# 205 , res = 0.00016320873564443518
it# 206 , res = 0.00015844125138680747
it# 207 , res = 0.00015381302993245984
it# 208 , res = 0.00014932000328108505
it# 209 , res = 0.00014495822226236616
it# 210 , res = 0.00014072385306511352
it# 211 , res = 0.0001366131738674067
it# 212 , res = 0.00013262257156536697
it# 213 , res = 0.00012874853859732883
it# 214 , res = 0.0001249876698608756
it# 215 , res = 0.00012133665971992045
it# 216 , res = 0.0001177922990992615
it# 217 , res = 0.00011435147266356294
it# 218 , res = 0.00011101115607990808
it# 219 , res = 0.00010776841335868004
it# 220 , res = 0.00010462039427341043
it# 221 , res = 0.00010156433185559654
it# 222 , res = 9.859753996242925e-05
it# 223 , res = 9.571741091609476e-05
it# 224 , res = 9.292141321128581e-05
it# 225 , res = 9.020708929067876e-05
it# 226 , res = 8.757205338429242e-05
it# 227 , res = 8.501398941303418e-05
it# 228 , res = 8.253064895263055e-05
it# 229 , res = 8.011984925731201e-05
it# 230 , res = 7.777947134147551e-05
it# 231 , res = 7.550745811731941e-05
it# 232 , res = 7.3301812586311e-05
it# 233 , res = 7.11605960842311e-05
it# 234 , res = 6.908192657717135e-05
it# 235 , res = 6.706397700723422e-05
it# 236 , res = 6.510497368653672e-05
it# 237 , res = 6.320319473849986e-05
it# 238 , res = 6.135696858402421e-05
it# 239 , res = 5.9564672472634836e-05
it# 240 , res = 5.782473105574829e-05
it# 241 , res = 5.613561500249084e-05
it# 242 , res = 5.449583965495543e-05
it# 243 , res = 5.290396372371141e-05
it# 244 , res = 5.135858802061855e-05
it# 245 , res = 4.9858354229300375e-05
it# 246 , res = 4.840194371106144e-05
it# 247 , res = 4.698807634598913e-05
it# 248 , res = 4.561550940747264e-05
it# 249 , res = 4.428303647042051e-05
it# 250 , res = 4.2989486350428265e-05
it# 251 , res = 4.173372207456912e-05
it# 252 , res = 4.0514639881982265e-05
it# 253 , res = 3.933116825356327e-05
it# 254 , res = 3.818226697063176e-05
it# 255 , res = 3.706692620003855e-05
it# 256 , res = 3.5984165606928974e-05
it# 257 , res = 3.493303349294157e-05
it# 258 , res = 3.391260595975471e-05
it# 259 , res = 3.2921986097031015e-05
it# 260 , res = 3.196030319400556e-05
it# 261 , res = 3.102671197426804e-05
it# 262 , res = 3.012039185265958e-05
it# 263 , res = 2.924054621418891e-05
it# 264 , res = 2.83864017136738e-05
it# 265 , res = 2.7557207596285808e-05
it# 266 , res = 2.6752235037222646e-05
it# 267 , res = 2.597077650136112e-05
it# 268 , res = 2.5212145121445153e-05
it# 269 , res = 2.4475674094500226e-05
it# 270 , res = 2.3760716095100676e-05
it# 271 , res = 2.3066642707015354e-05
it# 272 , res = 2.2392843870726235e-05
it# 273 , res = 2.173872734695419e-05
it# 274 , res = 2.1103718196545164e-05
it# 275 , res = 2.048725827450787e-05
it# 276 , res = 1.9888805740331945e-05
it# 277 , res = 1.9307834580586206e-05
it# 278 , res = 1.8743834147660917e-05
it# 279 , res = 1.8196308710329865e-05
it# 280 , res = 1.7664777018027522e-05
it# 281 , res = 1.7148771878047257e-05
it# 282 , res = 1.6647839744876715e-05
it# 283 , res = 1.6161540321541195e-05
it# 284 , res = 1.568944617237857e-05
it# 285 , res = 1.5231142347775923e-05
it# 286 , res = 1.4786226019085057e-05
it# 287 , res = 1.4354306124536763e-05
it# 288 , res = 1.3935003025872397e-05
it# 289 , res = 1.3527948174279261e-05
it# 290 , res = 1.3132783786710763e-05
it# 291 , res = 1.2749162531343342e-05
it# 292 , res = 1.2376747222128253e-05
it# 293 , res = 1.2015210522471506e-05
it# 294 , res = 1.1664234657852027e-05
it# 295 , res = 1.1323511136057308e-05
it# 296 , res = 1.0992740476262984e-05
it# 297 , res = 1.0671631945903052e-05
it# 298 , res = 1.035990330474014e-05
it# 299 , res = 1.005728055716698e-05
it# 300 , res = 9.763497711510873e-06
it# 301 , res = 9.478296545492734e-06
it# 302 , res = 9.20142638007839e-06
it# 303 , res = 8.93264385864601e-06
it# 304 , res = 8.671712733319411e-06
it# 305 , res = 8.418403657416411e-06
it# 306 , res = 8.172493983483656e-06
it# 307 , res = 7.933767567766574e-06
it# 308 , res = 7.702014580399288e-06
it# 309 , res = 7.477031320963506e-06
it# 310 , res = 7.25862003899917e-06
it# 311 , res = 7.046588760873389e-06
it# 312 , res = 6.84075112049103e-06
it# 313 , res = 6.640926195685272e-06
it# 314 , res = 6.446938349203836e-06
it# 315 , res = 6.25861707457476e-06
it# 316 , res = 6.075796845556455e-06
it# 317 , res = 5.898316971400937e-06
it# 318 , res = 5.72602145524725e-06
it# 319 , res = 5.5587588569724775e-06
it# 320 , res = 5.396382160124041e-06
it# 321 , res = 5.2387486429536435e-06
it# 322 , res = 5.085719752475799e-06
it# 323 , res = 4.9371609832884464e-06
it# 324 , res = 4.792941758633401e-06
it# 325 , res = 4.652935316508108e-06
it# 326 , res = 4.5170185971830715e-06
it# 327 , res = 4.385072136057004e-06
it# 328 , res = 4.256979958117818e-06
it# 329 , res = 4.132629475959835e-06
it# 330 , res = 4.011911391118441e-06
it# 331 , res = 3.894719597675873e-06
it# 332 , res = 3.7809510893390857e-06
it# 333 , res = 3.670505868684667e-06
it# 334 , res = 3.56328685916848e-06
it# 335 , res = 3.459199820193475e-06
it# 336 , res = 3.3581532637970477e-06
it# 337 , res = 3.2600583745142916e-06
it# 338 , res = 3.1648289314628957e-06
it# 339 , res = 3.072381232041618e-06
it# 340 , res = 2.9826340189667907e-06
it# 341 , res = 2.8955084082790716e-06
it# 342 , res = 2.8109278206005745e-06
it# 343 , res = 2.7288179132925763e-06
it# 344 , res = 2.6491065154456158e-06
it# 345 , res = 2.571723564207027e-06
it# 346 , res = 2.4966010434604014e-06
it# 347 , res = 2.423672923871059e-06
it# 348 , res = 2.3528751048535465e-06
it# 349 , res = 2.2841453582609813e-06
it# 350 , res = 2.2174232737614964e-06
it# 351 , res = 2.1526502055997703e-06
it# 352 , res = 2.0897692210266675e-06
it# 353 , res = 2.0287250505702893e-06
it# 354 , res = 1.969464039117375e-06
it# 355 , res = 1.911934098781118e-06
it# 356 , res = 1.8560846634527472e-06
it# 357 , res = 1.8018666438437067e-06
it# 358 , res = 1.7492323847821698e-06
it# 359 , res = 1.6981356230441588e-06
it# 360 , res = 1.6485314469923352e-06
it# 361 , res = 1.600376256649062e-06
it# 362 , res = 1.5536277256803486e-06
it# 363 , res = 1.508244764407703e-06
it# 364 , res = 1.4641874830384717e-06
it# 365 , res = 1.421417157300543e-06
it# 366 , res = 1.379896194015483e-06
it# 367 , res = 1.3395880979599996e-06
it# 368 , res = 1.3004574401699742e-06
it# 369 , res = 1.2624698266073e-06
it# 370 , res = 1.2255918677765647e-06
it# 371 , res = 1.1897911495989067e-06
it# 372 , res = 1.1550362048340911e-06
it# 373 , res = 1.1212964855609552e-06
it# 374 , res = 1.0885423358983803e-06
it# 375 , res = 1.0567449665072173e-06
it# 376 , res = 1.025876428871224e-06
it# 377 , res = 9.959095908646414e-07
it# 378 , res = 9.668181131464157e-07
it# 379 , res = 9.385764254222216e-07
it# 380 , res = 9.111597045982461e-07
it# 381 , res = 8.845438526144883e-07
it# 382 , res = 8.587054752484367e-07
it# 383 , res = 8.336218617996018e-07
it# 384 , res = 8.092709647933546e-07
it# 385 , res = 7.856313809619995e-07
it# 386 , res = 7.62682332102173e-07
it# 387 , res = 7.404036470599436e-07
it# 388 , res = 7.187757438579928e-07
it# 389 , res = 6.977796125536263e-07
it# 390 , res = 6.773967983881534e-07
it# 391 , res = 6.576093858627596e-07
it# 392 , res = 6.38399982732111e-07
it# 393 , res = 6.197517047638082e-07
it# 394 , res = 6.01648160936633e-07
it# 395 , res = 5.840734389872142e-07
it# 396 , res = 5.670120916469611e-07
it# 397 , res = 5.50449122686387e-07
it# 398 , res = 5.343699739852417e-07
it# 399 , res = 5.187605126956174e-07
it# 400 , res = 5.036070187690837e-07
it# 401 , res = 4.888961729790659e-07
it# 402 , res = 4.7461504513396437e-07
it# 403 , res = 4.6075108286727444e-07
it# 404 , res = 4.472921002565832e-07
it# 405 , res = 4.3422626753263947e-07
it# 406 , res = 4.21542100204526e-07
it# 407 , res = 4.092284497445619e-07
it# 408 , res = 3.9727449282075296e-07
it# 409 , res = 3.8566972244087746e-07
it# 410 , res = 3.7440393860630063e-07
it# 411 , res = 3.6346723919561035e-07
it# 412 , res = 3.5285001131485737e-07
it# 413 , res = 3.4254292287366425e-07
it# 414 , res = 3.325369143416729e-07
it# 415 , res = 3.2282319100270845e-07
it# 416 , res = 3.13393214867448e-07
it# 417 , res = 3.0423869744928535e-07
it# 418 , res = 2.953515921734483e-07
it# 419 , res = 2.8672408788377657e-07
it# 420 , res = 2.7834860134073113e-07
it# 421 , res = 2.702177706588409e-07
it# 422 , res = 2.6232444952714855e-07
it# 423 , res = 2.5466169985204204e-07
it# 424 , res = 2.4722278645321504e-07
it# 425 , res = 2.400011709010443e-07
it# 426 , res = 2.329905057260869e-07
it# 427 , res = 2.2618462877411483e-07
it# 428 , res = 2.1957755802810812e-07
it# 429 , res = 2.131634862099809e-07
it# 430 , res = 2.069367755911583e-07
it# 431 , res = 2.0089195322001486e-07
it# 432 , res = 1.9502370583299227e-07
it# 433 , res = 1.8932687573653563e-07
it# 434 , res = 1.8379645543772082e-07
it# 435 , res = 1.7842758409160176e-07
it# 436 , res = 1.7321554260859983e-07
it# 437 , res = 1.6815574985223712e-07
it# 438 , res = 1.6324375852972194e-07
it# 439 , res = 1.5847525120331187e-07
it# 440 , res = 1.538460365826503e-07
it# 441 , res = 1.4935204577373734e-07
it# 442 , res = 1.449893287671978e-07
it# 443 , res = 1.4075405091746439e-07
it# 444 , res = 1.366424896713164e-07
it# 445 , res = 1.3265103096566273e-07
it# 446 , res = 1.2877616666977125e-07
it# 447 , res = 1.250144908844406e-07
it# 448 , res = 1.2136269731172844e-07
it# 449 , res = 1.1781757605355208e-07
it# 450 , res = 1.1437601128620873e-07
it# 451 , res = 1.1103497789963403e-07
it# 452 , res = 1.0779153926650003e-07
it# 453 , res = 1.0464284469397545e-07
it# 454 , res = 1.0158612644480159e-07
it# 455 , res = 9.861869787720146e-08
it# 456 , res = 9.573795077326162e-08
it# 457 , res = 9.294135313771748e-08
it# 458 , res = 9.022644680141832e-08
it# 459 , res = 8.759084540422157e-08
it# 460 , res = 8.503223247641311e-08
it# 461 , res = 8.254835915869222e-08
it# 462 , res = 8.01370421416458e-08
it# 463 , res = 7.77961620073764e-08
it# 464 , res = 7.552366120809265e-08
it# 465 , res = 7.331754236129858e-08
it# 466 , res = 7.117586639524648e-08
it# 467 , res = 6.909675081918946e-08
it# 468 , res = 6.707836818986148e-08
it# 469 , res = 6.511894450425933e-08
it# 470 , res = 6.321675746324933e-08
it# 471 , res = 6.137013513206768e-08
it# 472 , res = 5.957745439495801e-08
it# 473 , res = 5.7837139561921364e-08
it# 474 , res = 5.614766108888818e-08
it# 475 , res = 5.4507533852471614e-08
it# 476 , res = 5.291531629186205e-08
it# 477 , res = 5.1369608995322195e-08
it# 478 , res = 4.986905328823826e-08
it# 479 , res = 4.841233020610855e-08
it# 480 , res = 4.699815949500824e-08
it# 481 , res = 4.56252979541951e-08
it# 482 , res = 4.4292539149818404e-08
it# 483 , res = 4.299871142157872e-08
it# 484 , res = 4.174267767994465e-08
it# 485 , res = 4.052333390190662e-08
it# 486 , res = 3.9339608322299644e-08
it# 487 , res = 3.819046044416984e-08
it# 488 , res = 3.7074880325813256e-08
it# 489 , res = 3.5991887391963725e-08
it# 490 , res = 3.4940529730331197e-08
it# 491 , res = 3.391988323546586e-08
it# 492 , res = 3.292905081126291e-08
it# 493 , res = 3.196716148699039e-08
it# 494 , res = 3.1033369958910765e-08
it# 495 , res = 3.012685531488535e-08
it# 496 , res = 2.9246820885573157e-08
it# 497 , res = 2.8392493149548758e-08
it# 498 , res = 2.7563121083519314e-08
it# 499 , res = 2.6757975813604662e-08

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

it =  0  err =  0.09421429683611413
it =  1  err =  0.12242092674801901
it =  2  err =  0.08627911974342342
it =  3  err =  0.054294636877542546
it =  4  err =  0.026834953641369612
it =  5  err =  0.014442946534311274
it =  6  err =  0.006828162088958496
it =  7  err =  0.002947538932807923
it =  8  err =  0.0009884874866274685
it =  9  err =  0.00037122112200935675
it =  10  err =  0.00011770383684169319
it =  11  err =  4.7294009691925314e-05
it =  12  err =  1.795956404345207e-05
it =  13  err =  7.73981616121209e-06
it =  14  err =  3.688537792536532e-06
it =  15  err =  1.7853513785615997e-06
it =  16  err =  6.793337817618327e-07
it =  17  err =  3.267271476258251e-07
it =  18  err =  1.0509766947753396e-07
it =  19  err =  5.2632976611271436e-08
it =  20  err =  1.7518834053967413e-08
it =  21  err =  7.1749421044469985e-09
it =  22  err =  2.5718606236788454e-09
it =  23  err =  7.674623733376304e-10
it =  24  err =  3.2317029440816587e-10
it =  25  err =  7.901568725512893e-11
it =  26  err =  2.830256119137346e-11
it =  27  err =  1.0594770382019746e-11
it =  28  err =  2.983199288638775e-12
it =  29  err =  1.2174564854809225e-12
it =  30  err =  4.3306940793653186e-13
it =  31  err =  1.6903295709170732e-13

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

20.165532778155335

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.

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

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

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

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. It is also possible to experiment with coarse grid corrections using NGSolve’s python interface. We now show 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.

vertexdofs = BitArray(fes.ndof)
vertexdofs[:] = False

for v in mesh.vertices:
    for d in fes.GetDofNrs(v):
        vertexdofs[d] = True

vertexdofs &= fes.FreeDofs()

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 imply that the condition number of this preconditioned system is nice.

One well-known and 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.

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

0.993081
 0.997768
 0.999961
 1.34206
 1.80314
 1.83857
 1.96023
  1.9828
 1.99987