From e1e59d2fba8c07a281599c69f8869c8d305eb826 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 1 Aug 2017 22:33:27 +0200 Subject: [PATCH] add Utilities::LinearAlgebra::Lanczos_largest_eigenvalue() and Utilities::LinearAlgebra::Chebyshev_filter() --- doc/doxygen/images/chebyshev8.png | Bin 0 -> 29790 bytes doc/news/changes/minor/20170801DenisDavydov-c | 4 + doc/news/changes/minor/20170802DenisDavydov | 4 + include/deal.II/lac/diagonal_matrix.h | 7 +- include/deal.II/lac/utilities.h | 318 ++++++++++++++++++ tests/lac/utilities_01.cc | 97 ++++++ tests/lac/utilities_01.output | 21 ++ tests/lac/utilities_02.cc | 299 ++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=3.output | 7 + tests/lac/utilities_03.cc | 194 +++++++++++ tests/lac/utilities_03.output | 123 +++++++ 11 files changed, 1072 insertions(+), 2 deletions(-) create mode 100644 doc/doxygen/images/chebyshev8.png create mode 100644 doc/news/changes/minor/20170801DenisDavydov-c create mode 100644 doc/news/changes/minor/20170802DenisDavydov create mode 100644 include/deal.II/lac/utilities.h create mode 100644 tests/lac/utilities_01.cc create mode 100644 tests/lac/utilities_01.output create mode 100644 tests/lac/utilities_02.cc create mode 100644 tests/lac/utilities_02.with_mpi=true.with_p4est=true.mpirun=3.output create mode 100644 tests/lac/utilities_03.cc create mode 100644 tests/lac/utilities_03.output diff --git a/doc/doxygen/images/chebyshev8.png b/doc/doxygen/images/chebyshev8.png new file mode 100644 index 0000000000000000000000000000000000000000..9eb23a1a0b4ce5cdd31e07e7dace9b1ffbda5a04 GIT binary patch literal 29790 zcmZ@=2UJtp*4CNlIEsKl0coOCktQHTN*G6~G!-z?rGxY?5MXTdB3=53ARVZxeqx)?aJn(@PPFg+&Ok=$M^Em5qWR=c^Myks;P3OlodfqT>cm}+GeP+LN8`L< zZ>x6U=vw=BdAm!lC7V_C>Z(JHSF(Gb6lQSoF>+r26!F%7uT|1sVhKNI#&^HGBpC48uG)I4Hxou&@xA?Z^C}cIzY6lwXIBg_YHkNA<$O zsle>UM#{^VFH`v4lA))kXieD-YKj~aFnyEW?2BvJSf-P-?stc zzh~<1{U=%Xz8Ri%dU}ShUtcWDecS9qt0g3gtHN`>eyyjc z_j)zE;mlSI4gP|F$P`%^1lcW)ZdYCKV0!&}#7yr_H2GNU3ttb(iy@-$(qYG-5v5C3 z0!G%GYd9PhijnE6wOM~1Z8o0PqSfSgRVXsD*~Q0S^+ru)Wn}uR|5HH2 z+HhY(!(@%>vo@e#TYh-(Do4AMY{)cTFA>WGjcu89ET7Xqq2&)THbwwSBy(78ZK6MX z;5*0v&X$aqCtGRu7Q^e;^tk&?w80m292^|n+}4K6o#uZ$T*tFb0k(FFLE9~YcnyPm z1^z)Tz9_VRnLMR}b>zdfu--9LcSWQ!;$dC?D}G5iBFRI8Y}bq2f}LG(c*_w*4ufk( z){K~1<5NT?=nHrCY<_C_b-glJh;!xR+p!uA{`mq1kAp8H_d*hO9?-bxSX8(Xim+^v z{D708&Zzw7C`%@XJP#fvH#ax`K5h!X5v7Ho-XxR76RZsl;e3)Y*T6QhvaEWLQ{*ex zg$Jtt^f@4+mG&a5Z-z@Z*A$+sxjfGn8Ft8S41%aGjYCZW2%DdK;Om{XfG;VNJct3!^Ho-z)xgi%U4bQfKf{?8vRjVZ+V^+ zDgJd0fp1T`2-~nm;=sNUO8dmGjfUrj{uej$B#JJJyAck)Cz-dXkt@MDb)RmBG_o*Y zykq|5sP2HDDN$kFAsAzZV1Kt>{6v(s7mzhkJ@#aM9f%EBdJU?(6q$iry0|A%=BhKq z#k&1@2k{Qz4i#{R1h|7$?=i6Ycsc8Oh~C{AWMeoVSI~i>tFU6znGj`SdJ=H)lJKUp zKlG12UzQ2lBMR^A?CLTCtWLEAtZwwr-zT}--xKqiJk(iH3+U1wi!A9K_?a-@D^cA= zz{-v5`uCOAvaZ@Wtwp^ejpqn0zyn!YP8r-^3Ae5Bc3jPVEa{qu6g!}Je~)Fc-hOID zN#7rBRLL(g0Zll`k0d&W>u7!PuIMEO#Pk0E3BG&_^n~~;$)>@Ms+piFH4Pd21e3hP zTgcifMT)DlHEZN$K(QLh&486Q9?*2UA`dS9u{)dsteZULU~M=`zzDE#*i*hD9uU51`S9@Swp%9%uRQE$9?KxIvAD)pBAg@bIlqbWFa8tl&RvZrK zVxqu|8N?^BU;czmzwp|aGwzJxp*9%TRaYQ6Jycy>=GgWmK{Brk4tda`%O19bUnxZE zOLb98iNgN}iQWuQCW$7KKY{4?#3L>@JS1On>eQaTp6S|Wz40T$&(D$4D$o)(thB`leNf)KiC+-0jByyXibJY1g3{XI~hYex* zJ6k04A;U4$0;2F~i&9i~R;02wz$^RAgVl+#>Zy$&(A3-_)_IG%`0gn$ofNHVk6T(Q zo9RA!mE4@b7RSho1N5dXwj5c5O?qP+T2!V z5X_YM#<2HITwZ#5`rEf=sM;5lQp;aLnIu?P^>}%^Iy*ZHbzr;C{GDA~1`ubwKVWZl zeE7I(ihm}{Y}OX9CGiSNOQU$G@rHd?%E=UVkuuYGfwxyw2?&HPjGMy7pZ*@i zzho*~-EOm76#gbPbq}*2MBr)`U-x2dbnvt|r|5MhXCvJ-|A`pS7Fijy^KRa_zBbzc zrQpiBduVi8`MRTz0v&zQ8J682c&wV`4-9@~m>X*}-v42ytj;3E+aQnPP!CBlv6@BB zog#&9c@ST&%NIg%GMuREq>yvPC=WI`m@So_)t1#GQ~E=LNq2^@fJETc1qU&nsu!vZ z&}sxS1I8J*_3;n~yz$S0=3i#SXv@C=G12c}S0OAV8)WcQn4el6d5Rn`3d%!o62da& zi~Efc6LIg*;4IHOoE3qKhd7U;({kFdvION{76U{O9$VSHR`=tF>eJJMnALF|n5@sy z6^vQsD^kue^Q));H)Ku^OkEn7Iw=pi9$1#yv`CR+#yZk_AA>FNWI}Yli&WkLafy_R zya5QT@-GE|?;bhXvR=Qg@y1MW{RzaFn~zVw&@5ym>k=6;p*WyZmhJLQdIn;Eb0Pe- zPwPICm4lo#w0ooOr`B#qDEjv5Sg>&rn3k!!N_Nkd=KLzKM%LVA3A4?oI5oQU7M z5Psd00HqB)T#f-6BHlE?eBcP31Nz57cFe}nFj)BkX-<9cBPn=S$S^)%O`a)3Wxf00 z2|3*5z6VpP`!ukJ4fVi^qCZxA#e@j@#-(@fi2J!*v5w)JA5PBfM ztg)ccMefD%2fysO`U;uWaICCF5vn^KIL#mMPpNJguO|b7YN?W2>GLaN@n3)aT}8!l zr6*=gcHOJ{1tSqUgmI+&x-<@Cs(k7rU%uSYi{Il#PFdOGws>LJyr(>MAOb7ZF=F%+ z`52%F8#~}eE5%qnX3>P^;QPH;;vDd*163)eCl1#x3gX&{j1Rx(Gf{zKxOL6ip@#|z#YAFjC=$-!V%7RDze9UKs%2-FmUkL+BqSTq%bnqM! zlA!^*>&B%G4^azJ6+H3u5L1Eg(>%Gfe>{FaC5H$|z++SEXPt|GCJWu2Vfyi(hr*MjsCaINwz!O)+qJ*^?tBUXKJMKlf$yZHqq42@B-O0Qlt`Z+n3F@|Fd$`xI*96 z?`{!Xe;vAxuVG-{E#DMkHGtD2?Rf*(u^S!O+j%=y^>06T=uK30Um6DJibT`-ONV&+ zyZmtf6-&#T6|pg_Bd;!i$k?7-?8y|H<_!dPn7nR4mH=lmyd%}NkvI45CKUoQnF0EX zjg4Xp35kd>S`$aaZzP$VK6_}8Dl>BSBS(VD*y@KqpBDr4{q%&FmxwBhj`+HNKS^L*7ad@+z;<_cZdN2UuXAgnholW*Dp2Y*$c!D-r9JvcVe)*K*sCSoPPa4$W`10a-5Lv@;sVf= zAK4QcA}V|WSLX|S8B;xA9Fed+D?JZ54^06{u6lWh@>?LkeRqilDmQ&@h0^BAsN&m? zhfrHx5`D>zAl6kJhz1N|1Pt2qOM&vO;nJH3<_v{OO0a(e3qJsbO zJU1`zCr_i-cLRa+?GL5DM}q+uv`4cF=n(K~q;RnSbZ99W%>X?hnB;B(Pb=f|6B-xQvYQkwUn6Szdv=8a3y4?Ix@s50}I~I^)x~Wt#qeOg@6A1Z6urln)S;x!|^~nf(9s&aRK!pi8qrGXxATX z>|E7=%aMCrChc(2KOy0nhYYL;i)ZI(4mbxiBP(-Qkwg85o!vmIXI$^cyFa?RR;5$1 zYwu-c-Q_>q;AQOWk-9#+c+)5~5z4EdC&pMZyes(YV;yci3#O<3DnZYH5s9*iNB_yM zyWsdn#@xh2CPhkjtB_5;Juyjmo+bxJRf}X7*N%^W>SLbrZl8jKt&)>RTo`tZdre{B zX{uRgH&}(~Oxa=c`QGZD^?h}_I931N8~%j^8pYcY8ygRO<~!xDg7W0~L+NUSK3E_( zSbme z>tB|(TWZS~YlyrKlleM||An<80F}`*^gz=S#5}&Y%-$9;;nRw;>br>h?dq=*dnVZ{g9E{$F_X=zo^ zT+$rHg^%HU1HOk|zW_4ZglC{mvq?-gk|@CsQymogYZ~; z+rxd|jz3EXK}C1%uMvurW>8}!Di%qs)^&r{gsUYx3z(Gz?zF8*Y3TNHug}KvSbE0y zX)br+TXI#>(uMZ7dNQPQht|J+PS{OeU3oRsm}dQ6K$7LYJ2JzXGp-dBF7_kgfWONC zURE`})Oqq;P-9D)sS+$VF&5XNQR?2&P)#_gK23NTpg ze9O0!v(V);WYz|MNIzIM?9J9%8}q3lehs=GgzeW^S==lvwX3G^^;_zZ+*}y$U1deZ zrc}^%;M1LP7T3|%u46;k5C*gEw z3NH9pIy6!keow{t4LK9kuWz<`j`^yN)Y8a84We*z&27HirDdMlN*f_%mzW=YUrJ~) zccXrCNiL|;`ZhMyWpW8EQ>QL0CAqX&oV4C#SjpU~l~S4`wp-{{6v8ZxsI06^OPlC| zT@=2(vjL+MZvNu<^!1xJx6bc0Tp_xB&$;|YPynz71%PI`tv0*nD)+@UJN>D!^)kp} zxEd1z552d9TBvhaqDj_UPf^1K1Yd)fh_Y|`2F_Do3}4-rC{eI}^*Ej(UdoV5)Kn?5 zg_R|V#VhYm(JJ_hng&`3;GB2w{41qM56Z^65em$9#Pe*2R_QyQ%Z7Oi zA?zHpQNa9XYF)Tqpl;j{FfOqha#bs;T#iZ8xjpp8=flo=H0v7 zN|aN%?jukx<2gWnTnzcm9jKs|6e%?s87eaGyG6K%l?okq$=wF9+ptiEd1QG=9bJcu zs;x!LOkJ4*q=|PeCdl7Jg%`J_daICAK%hNuC94vpzGQMXjtJ{2_v-(|Wf>41yS}pv zb487?>h+oBI+wAQMCQVFVcVjx8alVmnwk$Yf9;e+kKi*0fVP%k0{crF(7dr#3Jx6`=*f2lO!=v2Gs^87f8ziFeuL-FTpv zQk;${-x3Rd`gZtaY?KJC$@==XgwT~ClzL&LGRKc1lv!PdTFXDjMTYwx_~li5z){?mj1e-D5{I|1y>|^zFs>7~x<7t0TE1ehyLMRm z5@zi+M#gGhvl*IQFYil#P(&~NzaC9aE@<(Rd`BSL)_^;REdYo3Ut}z*i}}33;p`_11Cie8T0xKrEX~kTmR+-XgH{0pMNpc z%JLp)l+0mR`WpNoY#L8q&1bjs>S<3xb8VM-Rk?{L99y&!Lc`noPa2!0#zQ)Tf6X}K zln?P^Jb+P{FDRwtqM{P#%`pnM3fTr;kK*go1F*!*yF~I774N~!*H{oTpF{-XKnGwy zfp{dR3*%H&1k#v!7MqZD^=fFbDP#8V4xwnMkK=sx&Sqvwhx1ADh)7O=wPW86DDJ7+ zKGwd5UUdZ!#e3VjtT;d1cExw#oT!pnir&CLX-1x;vxVVSokZV&=WUnn{a!RN~+T-tnvtEEb|MY1hRP5<6p9QAbN`<8m3 z68p&XLcRzC;&>daJ!SP`m*)5@-j!jmG|2k~)|LxCQpd1#0n5Kj7#!qYx4Jv~%(eglSG$pzI6sv44(t zb+KFwd%nKzW_;>blHiMJ!Aw>iIUxw8($}AOS*V09=5W1UNGv}tY1|%crCE6WJg^Di z#Yv?b5HNcr;KD;eL1Q*d=eaq@l4`QiqkW~h@h`<~r_Kb$n$>;#) z%n*f6>bE_Gq==`l`tq}{iTBMLhgxS7!fHRidpl&SaM$pq=ab|<%vIT6@966rT8bxM zZ~z@VQIGUK7xVGm&4Id;8)ujUP@@-qf-?wg2TXPW2y#$R$3tRQ}SD!983*5^i2Sy`+>jM z-^LmjfKsnup^8~Z?7Ol2Bq0(ll)}7Q;kI;!8u2CLn;YS0mhu{nh`_t>t;p_C5pD0oufqngGH1)99TD+{28N%9bZOLvnZDjMRbWAn|oWDK^ElVpZ zztpl?n#EOyQqKXf`bRlLyKdE$5Nt~rv>NQ5YHJZ+OF=uv=C%rZD(>8|NA~>>5E=9! zBiz<=CAT)}6^7FG9?AqeP0n~FLk(oeWSPB+9Nj-uju#|*L&+69auPrJ%2KEo8;`QK z4!cM%Ha#&z8@P=hjl5!iwWNsIywZ7WgE(C4M0ZKG_rI9%1)3a=v+q@0G$j~ih}4rD z;NG=Qez&Qt>{PE%RAk=Gv0qd{QOK*dKfl2((o@E(6~p`oMU;BRoPOQ(yI*KCN}{mQ zu-jOFwSuwB^;mCueg$K_IUTdBSUEmkk|dPMtDOXS??TW3T@DP;@%{XSCSR2a(S4GF zPfOBM_r8ZC_&k!^nUklT8-(8p;YfJ|nv|efA%1|lqXNsal;RUiKOP&Rd6%jy~AN*``0w9I5va-rbPwa5HEC2KCT2N5Z zr~|ph%S)HNt8le!4I=P>@hM@Z3!7cNHnU3f>2|gKIGjyOhw#pvnbH}?E6nhRYb;Rt zVqjxG65D%_q0JLe0k{QHZ|U9wgJ`S%R(GdYAG9PJpHI8m=bo+g$`EWg3KTM_4~_b& zktcxXa-#5PE3%Z}8q=3clU19kS05!Lwg@i1l}i77>y_pZKGRjqpD4rqu1`8kNj1{k z`l9FC;)S8o5S24y((dbF1rg(H>+S3zWJwpB?mf(wuy~m$X-OnAQ&mt4)He5B0U&b} zB~+qb>Cv26H78I2uZ>w>lS*8Mr%1W+w+|Kk^<*UZT^~w)9&{*`Vg z#r9FejKwL+QYUins;|5t%c|4!wrU=~y;VxM@$3(-WfPioP=Y6Ja=-SbTDP(SKuMbw ze!T2OM(+1J#>nI0Z+$d1S7V>I2AjxI*qotfv#YULV`*2&!cl)&Wpn3hE;3=i{6miwe)ntDe83pCwM5xW0xQt9XAe|_-y|U2P zfCn!qQq&j_{g6|*7)_-Xzh=kmcsS&Cc*=7f@W zh1d43(~8upv=im9M%iOC-=3C5oep1|l|i8=(^b$kVEF;V0v^1%#Xt=FCMr}Ah+7+? zxnejw`9)=mE3U?I*g{qn6L^~m>H>sG+~Zz9MV$?ewY}leywg~4!3Xd73ZVBWOqQGI zby&`;v%*uN)Fn`{jUbP9n|FUD0W18`gJyxGNOA=PP?@c<28E1C6G8jZO?mCE-y^Cg zTyUwUXtJAqo}E;(%p3_i3yB6EkI0Yxs*?=B4sCXw^QgwG9{S|CAU{e3zLL z%@n5#Py|ruCMp7{Kz%-dFrm0pJ9P;7odcCQBDY0v66|G@$&+|_^%i|SC=Bk>|yS&P^TU0TzLAe zjpz1X-K^uh+A7~zSlGn1=^5tr{RDwkQXE>PEK6?zk`vn=Ut=wK@uI%AWFyL6jH7h` zHl1%~xT9|Jo&jWn3PN!_pb|h$P2H9x(YHUbbw1!h3%SAMU@hmk9o+={1`hEqaFAoV zAh^m>tyrh%x(_>NOQyOj|KIx?AQ`qS3lzw;&)UC}d@QcF!EQ`MM4-xRago`Sy)k`i zv^*`pHS%h8Pn8(YSx7qj9yo!y=hp#)XzZ1s{rQ6i6 z2hOi)+<2F@|434a`ZR>wQ3eKZ27R+@Op0PlibF@vEo~ZpovK9)j6&p!$QRAG$c&$< zeBs|SfTy$SVDdI?j~QB>Y*sSHy3V`wuJZa^-OSAxbZrkVv2gLeYl6)9n=_8*2hJkh zIE#O=MmJhxmsX~#roYlX{Bsn)F}G%qey!zY5g;9#H$evRUTYg{DYQf5%a^3zVR; zJXW7@o(oKDuWsZtr-6B}Wf#KjsvR4hW$)dan7NzK_wTO^N>Zg0|2oQ^Xl$J}q&?>9 z_K7Jo(@>33n-;p*$y{curg&$L0nWlA==_dR!fc8DXW9oY=J#;}M{nX^oh8=zz|U^R z#eJ#qeDCrjeTg`DI>?O$indWUBg4!kn~_4($}Y?94zy`lMa~P|C=209=qjt z+cLAw>O*yi@_m5!0)v7;9cQ{d$!dS-wL~61_v$L}6&0v+;S-&aX2CRWJ-%Grtom1y zZwv;?p0J*|Gf*`!CEP!%$dPpjNeiWnKq&b`4U^pGK1GNr_KUW8LCIjY`D2NZcgStD9r`rwG$_gx2y3#wO9MK`Svtp^sqRuCzR9kxjn?9U61&r-| zzxdU3As@L`j582*Z+-p>)wPOEgblXjY2j0QT2AXa7JThdG<$19_Mtj*kE^6Ethr~M zi%UeJHHmP3apUjpUTz*NtWxB{#ceqm9j?i{+}MJV%{)d??N9tm1EZ+IU+-@W%L309 zu)i$|;02qYfUuWL-S_}kHwpU|4h<9^7+wb;Pm%J*>a(>&tDDUb}_oO zg~1o6WOod96**bn{4}@twtCUkve2$(I&%nCWmNAwubIrh&F;3+u2N}rCT}9b_e^*s zXIpie03*%EgA4>vua>MaLk2g3Bm4BblSS&Z$0$cPd9OV-E!CD*cxAogzLU0`& zJi~WfPzG6A*BsGnYp<##lvGRBO&Lgw?F+)LPCo3rBq&~uYKl2JHa`W!9=z9_1^4U1 zJOGIx%G<(=TvomYtCa=}G4l-HfOj@58m z)zFZ&T^8fY%hR5kl`_4eau=kW9Zfpk2EO|f5g)V=&CE+?ZgK2hkb|w%m^CY121D;wj~0+-M&r3t6gM|`a7e00Ms5hfbVL# z`l@7p<#&~^U{Q|*gd(zVb9otGJK5cql(CqbKWq*lH%Cwp&jRl;u47@b0royrHjUi5SAF_aUS&u8^PxK9g(;VfFp zr%Sl+4+}fw3|G$A8>lsx0R%kQa=&iN{N8jsU2Xxjz>9pF$+gGygKx$qJ@rSrdAc-; z9VVrH`FtL?r)Afv`Iv4I7Dh%!*xA_|cHg$g3Cemeo&z&H+lXV}OhUMHL_}J_y+0_H zcTYdDdk*dBcLQ~S)?~tl3xl19smau!rrkwHxvGem^-T_3RSTAoDo>lMqHq3O5^_^|CEE#DUVyHYNdDGbHq)It)E)*&&y?^jTJrIQp>HNEDG9?!3&) z8xeJWOTvfUm;uq z;CZAyMd6ZLC&GY%kczzn7XUD3;|xL?_HFam_WSNRc%ck1S5=vA)aOHE4u^@f+f`+@ z^}%eO*R#hDBP)vxkrTyzPu-DZ_$^vV;73VV7p$b0bcX3sUmJgW|C(kbC$9sx+H={B z+47P2^^%=vY1oLyLW7iVatD<%(zhv`wPd7v`p;K3rJS?CXYWgOpyV_&6ZOy;M#&D=+*tkuzN`dj%gX%K_Yi?%v zR}CHPR&;QPb5^p~x2EST8BPOTMXl(-_P97=xwK~M{dL*2GzG6}XV=?hO_VOHr~FP3 zk3f*pzJUw*-CTVJcEnwP5^*tN9{W_{p1MR6C!ZkYrA`Ik{=)$O!uR&=KpXxC z;~$P){2{08hkhOVO#EL_a!NrZI-cIsw$1Q3k#VGGA|f~Ed|Nx7*KOx^gIq;ZWK3_2 z?d}jtwe)kDY+lTiO`sn>_Bh1GGYD*XL)8FOx8NZ&V~+~8IZwo5$VW&4u!%*>x^l`guz zC>!GsXXscwSNol;R8JG1v|rV)WP!RWm;2!ZP!wZv#!H-|dIC8XP4BM6D#_?}?dphz zw%l303%}*pti(UyWc3S?&&dtYGnxE(+ygJmP44Cvcvg{}VY=EFdUz#C^su@=!Py(d zl>G0s7hRa^7f&F42UeMW+VQx;+R_=Rm~=1XGbfQ~;TK-su``x=!xiXH=FQ6bv5L86 zf7c?@%LjaQzzHOK{eM6h;J$hlZA-}j>Unw)W@)}#zW>=Ri8Zu9(r!O9XwBwMroq}P zzY_GCYK`mWoDi!uzUXGw5%(5rIdb{d1i^s;+am7)4mgm2V@Ixk&m5kP0pnYruDmlD zqZL+WmDs4765fBjX3Oi5V_Jc^v#bo8ux2Q8Y3E?lKn$qP;J9;beQP zk#%9Oesz}p#6;8yh7?`xm$sdOts*A=Qv5^T()N^zG`?U>oVJSOG{t}WHnuJpY(q%I zY0guHH~72sW;dL6O&a>P&QE;69!eab0(J~?_0uI`A=OHvSq}$G`Y5kU367@rJouDL z&OFxj=Mkf)TIhUpwcV+nV5}8!54zJbhvDxnsjo{rlX7j{@5PDt(%~Js@(T_RjuG!6 zm2G(;paTVfT-Lxw&|IVT?l?-`8$KDs`@5h~GfiFV@I=0aB?rhlT9yp@4^q-Xb|c}i znLv=WG%YSpW=&0*g(NS`_TiY|Bf;rDI^On5FNRW!i4CsdWwPs?eb04OGI}67@nA;$BcS=;_W=$ zL1__~)WrvJOH*Y0rK<3y7j!V-1K&G0Ry%;!ge0rn=hwGBA-G6RU@wy{WJiC0L?^)V zd2Kl$pFG7sV`;pwF*M3VxIV$?+|KRwmA{&wZU0qBi8^ID0`%f3|oA|jwie$*1FsF{oDdm#z)p^q)q2Lp^)B(MB)0!M z8cy2cp1GX3fexjxw2u1p?3t4&@(iVMw|Vx)3;A`Ya)YjJ6dA`);&7}J76$3j^9x>Y zWIq6|3*+EZYj&lk%5=;Nhpj61b6vtStLf=YmAF4vLXToS+6jNFM{L5vuSBmm^ve}_ zAm;)S*b6&y=jMjGVHd@Vu6d)AW(Rq-<1>9JV~rcZ%m%dU++j?)8A_PBD?64tgnZzK zb9R|u^h@FX{K#OZ_nYOpuE=zY^47sIzVAM7_=Pnom`;CZ`bif*Lzto2Sco1-ofuxa zPi|RKdT-a?SV~Ud#QJ*q=-sGc$x>V&N%m)BZ;BJT>2;p@Xeb0cVwE@PN_#l^bu!1DNuQ zt;GX7_ikVRGt7MQVSnIU;c1-2NV->^drIe^5(0jpkUwiXnK8QfmK>s}6j?fYdUXDm zQqLANd=7WSKYWR$ zrl)UqKQ6M@O*k@}N?QbmKH08opuc)Ay#f{UgZ?p}PJE1(Zs!f5LSpftF*RN5+xx7E zmA+@r{6BF)yN=oW3QUhWQD?p5{{utL9x^SOuqHp|KQQ`7(j|Ti$nUlUU;S87ALUQq zRZI;rLZi9`Qc$7{C~@TqPr6GP&c zOD(wiqYVGVh}oIo{8S#j8w?d3h>Y@@EwVm!ZlDxx&Th=C7nNesRKXJ&$HnY*KA4=D z8FLE6wY}0>)xWY^o$?XkhAS$+nOdjmly8*ai5~$9KUTY)jdfAZ8yx4YoNvQqy1Jh^ z-}xz|KRxH*R>}y?-VoBNv+(1$>J92(E?>d27CsW~u-26_T`I0IxUO1AF$RP*bLX0E z-A@Du0JEkOyR)%;sy630Bm51xVSqX;hPfh7xL&}tUR`HH#$thLge)yX$;AcSyJT>9 zOqWcx?z^;Oxu3K|$+F^-_-so|TgXL73mC;kli&mCAP2at!kQ3_&~8&NaaoBvbpxen zETZ`pQN2=3U&`K6p!!7w4pNy!+K0rVOJHEx5VC>Fmp0L#k|`pbH3t!J3b;~G7=yQV zV|NR`Pl9vKHaaPc%;oI!6iNCy94re9RPRAR>bk6Jz_a2jVt;Te>Sx&@C{%d(dz zz~^i5)<76vfARXizAA_?Xu%erSkn6?-5xeG&h_Rxer z=ZAAjI;A=MmQjdP#x*siywb4BXACfXx?uT*h4ADb{QMh)TP5AZIPax>{W4I39SNbGk5;y;tL% z-ky(+m~4qn*(Mcv*Wb7d;R5`1y^80??B)+VA*mq$^z6b(=Fh@?d-p?&7UpWTBFS(P z3E?#sm{hvbCUji9)Z_9-;l9r&nICwPk$ejY?z`04*nqMC(pSsvq!8T@UQ=rGW8^7W z)#=gkNTDEZDXLyAM5%sh9zGxuht~?6t|6u)q{-Qi@BY$+zj6^^hvxz+-FmC9m$+sed1_n(_&h9@_R(IDTb-l1axYb zE2`FwPv5t@N~V;uFcmSsP~P1us$F!=%GBh+=ug<|zPzqxdK)nK8DD<)uCBqZO{27K zUP0T=taQo_ZO4u)Ow9L2YigD;tQo3WkN0;>8dxb5M^j3For>c@H@*(3kzB*R^6B8> zPae0I331Qlb}t{{G*pj)4LII6PlxUhoi{Bq|g)n)v@(2bjme zL&=3i0Mo*69Pz>(UH+976Q3b;H@Z?=67ob5mT%>sMV8m0xR<1ovU@Fhz%lWp{gIC)7oN74Z z;vB#I^P+G1NIo1R2mkhuiJ7*ujnajo98IW@l&4egYuc$vE025QPF>yUCHIGI1B#|_ zFZLiI!`R+Y1wb?<<4Xp+BJAQ6NH}J#BQ0U2!*r96?U7KGjz{pT*EQrqi}3&Qb@+S& zu4{kZDIH?fhZ5R&5VM;7w0<(F{VtBr-$CrHbeYjC_K{SFwZFfvQ;3G9mW#ttHMFeL z)O2hUrqLUs`a=S4KGg1wj`;AEj@!*u#Ur0gs^ z2sm(XNE)1oMGxK=SY<-o?AsRO1nsEj_35ftUP-r4<;{lD$12!C=$5$Rp}e6cyXTTl&h!^e6-<#b|)AH<=RjiGaZaP zwS%$i=_&R~hE*GFks21OsFF%4NcqMi_gycg+}0x$^91(K$HdbojRN@E@5}!`hMwlc ztnSk=_UW>MIfka_qu{k8)j#XK$%MkG?s+L=xju#%VN0_Ze}9>`wAQ26J5Qg^#am6J z$MTd5{$$wafk(0cRr`xzaDrU3QfZd$Z}E=Zf8Vf5&i6XdRoFO5JFWD9GkfeXtR<_w ztgM!!%e`X<^&1-5VebF8kq=*|Ny)GsQ=e5zEnM=shFM~eSwpCF#4yI>O^@BU%-Y#+ z)OyYSUopme2K2xYlr!n0*Wgm)v&QOArTULnYH2G2MA%fCQI^9pE-oAp-oR(&Yn2<`d^=98bE3uZXm&0OfoXiM*-oIQEs6Q$ zh6zmTc(4*A;M*}Nf*UAv`;?%RLBlhWZ?iq3(;>Z{TR!3L?$wgIwKYK9(pY2Q_UUf? ze-(i(TR{={)edUVN-sKk%AmZsI43i6AXRwj*=c57KA4qN0&%c72>RhU3cOhbMSRyv zr7oXS9(*eJx~N-P+JiDeL_`FbE#}oN`+{*UDJt?;1WiFnQ=PsYdgzJSji4r7rOO8j^7H#&9ZEn1l=SVZKk6RdT1aI<)jf%Rk z$smYj-hFw)$z%+gg}Vh_bmQTY8=vxaXDqeeYhY&z2ww6OTWO zC7zMhf-FejKED!}(icIQ~x4V_m#DGW^AyN>0Fv-%*y7TRKMHTWc+YaaBNkW!+y7TtgwumiejS-hVMXk7| zzKNp&tecgOH+nWtbiv~$^%rnZOo+_NO3J}3rye>ma8HX%PtOP~ePKUSB+cKLd(`2o zGiC$G7iMq3voe&$b{P@!2D8B@`*BGUr3)xOGyO;s<{-fpkGkEBSopjZhd%$i{}&dA zn$O+Pu_{Sq@oC#U%2u+=IDc zZ+$zOzI6=yk9n%?^z8~RczmMy;zaoPYG`8oNQp9zVZjsERNNQuBEVpeT#FMhZiR0N z{uIOmAfng(hj=qH-|DtrykPKAsBxxqaoFjtUdT~iC?{wI>C1sZ(Az_Y-qFgm&^ zN$R8>6TP0$4}*o7LaW7G*eTHkkt(l z&l!TwsL$6Z)RuyaKAkuA@1YH1ypN__G2gDzNuH{)Po6_YgX zv%KpKqFY-oM`cRJ24-t5End!9pNlmXawQIiZM>*epVgE zd!;AXMpS#&k1s6s9mQ1qFa1AFmDh%2*_Uzr7|cEGtL5c8TjEDWy0=ueo)ubqu90Fa zeb~ElYbDQnhwq4_mL$;QjR3pnR<JGw!zp3TriZe^=p(DsLtJToV;i5Q&;oX^^{8;N+*K#32; z&7}#*%^{VVn`mmW4zoE0#(vT@pN~ppM|P`SLJ(fY#-W|hA_1ODS31I>6uYq$bkqW- z4eBJ^nr%ws5Hy?l;6yD3Lq)}+sp;!;v`D;8yI_@fyxr*q1JnS%l&}L2`Ke02U1lV;4ug=6t z>1%lIqkA;JjxHUISat^EKdF=b!QenE5N6fJduJCq0yT;MBs;&3uGx~smLV&^RbX22 z%tamyX0cc25AUD!G+1t{nR{{GPN$SH!69un3kR3i_tHdwQ;fTkN4HBC2-*{aLzA-x zrM5RjjQu)b>+6ft(hP(p5^xZ`WqlDgmxFceX$DC?Cy2(cW6tDYWZs63-v9qstx8rOG_8=6w z&XxVwYNH=+zEvf!b+m2>&oxxfXb$3UE=9P+NjDr-ZG=pfTOK%X>FRnmhUf^+RL=6b zubY-zNAB-o&)ck1Eimu)a}RpflBW200qbVE>%ERgDs9spQ6q3=Uv~H= zj&cy}f;>gK7CC7y9a(I}zhZt=eQ{k7AHr)!aW?h0|8bIaMJx2Ftqrpp>{jsTUG>_~ z-J-NVlGlu;v+j)G(mver=txdsf=Z{SE81ey^>bQ|Kit83z8yY= zW=-FF4Ax;sJs*w=BE4G`%rvAEYPI}Lq0h%{__EW~=OMe65K;Y!yrap~C_9Vy^4;}$ zyJzJ@=P)hpcWbIG(qQ{L+9&#mwR&k1RT#MRa~%KiAwHP@lJrjg)~vzer<1hnN=KQv zrU{W=^lwdJ-U1fJI`WT<%x9l1zieGD9|AIVm&)%1@gzTK{y^_K+HEZQ9Op#4#IqVk z5SnS9C_028y;u{A7;2d8&n~rf=C$A5c0@l10cd)M^@uHRUI*Q~ zIBA3XIp3@2rBJ>dGM6teCX6ZE4R?0-w}0@3?C_*GxN>(B4AsipA_lN01_D|2guG^T zA`-ZRXj}4pz!(Pfa128V0@G}Iueq4#Gf+)m368NBG+vC%NRC(=10^lygA#!hWXv;f zYt2W=v0x-CpFtGPp`;v@t^`#5oF9nD_K7DANey#-3yT?F^4CmEj$>Y~{IbUxG4Z%$ z&xv|k)AdDve#yA7kM zm9gin`C^NT!PE@kM^Uk&QczJzU?VH1S|IXr3Y-N30WX20@O%Cofaf=IiWcu=c|Kf z4}I=kG2wSHSI>tUlomYQ?x z);!cz8jmmq$Xct_q#+V@-CmN0=I4`=Gn0^z((?d#iFV(6`&DIwW2qSM@JMNw_k|=Y zDqP^%m8>}l-#OjEL@aliL!?k>jfedXB!J1=o8K-Sn0&>*yDywa?W{oHA8M|l2_A?a zTJM>hVZOZ{`I!$q)7z`;R@Xc5Qa<$ec64NxD42v~;{?Ft%+Ck-hnka;FxcC>AWH$4 zgqBvrfCI>nu`LJW8OURaBuPqQ?1=!#BxL|-f&3;x%YkP=?X+ovMn>k>VnrL9qQA}| za07BZfJt-&B5yzn1Mq(q2cQK3w`yLXG;wFgRhod9GD0q#5I}q+Nkll1BpS9Fo%*yo zw%1o%MgU)#J55mdrE;X3=ZN(EGX{pqc{T@B^ZgUQKi})>l;uJI=5^D7X#t$r5YSKp z1N?3i7@1#FQ4zAgv7uZgn-sC|gI0=j01`9xObMd;I!+`TsJmlgQI!t`HbN5sNT-=g=?}-SnnqApU-S?sA-0*P!ZJ`%VYg$pl=I+}jvIuGc^^ zad9To^YbO?vjmILq)bPmZ2mrkJ6q5C;bksp@3C}%9}p2yQHff{KUv+IIi1$5v4jA* z&nS%_PTR3JXSo`0oe}q1+QWG=E{FB}E=5JfIoD}}$Lkgyh1F#$MD(eJN?Y@M&SK-Y z_s^dL-2(gTv9ugzWl1^Q<^!sy_p2+}N=kY;_GmB`1_pGNSCM!IS+6K~AJ)8hK|=Z| zG}VUE8uB36Cr?u7P^RpVHMg~X(pNwh9+S{_ZIUo?(7OG0uboyh-b@hmUo84b%W2P8JKLFXha29-S7+@^n zO=@UJB|!y+M#jp&dn7B%>%#Pkg@B$Y1RaE_81JUpGr#`ASNyqI~X#wv>k*VT0=cBW# zcU#Gu6-*-HqGUcWRkl<3wD0eW#+}wLPE35tq5FDgYxO2WOhEznfyT-~oFXt@AwoVr zD$&XM)m(1ZaSFuF3qx8WT3K(=Z{=^_g_3|Dl7IEnXkQ2Q_A@fdPBHssw}`#d0ZY?I z5Er)5$jB(G4`hk*QYWTvAHrw+QiNAJhf~- zeP#uVQBO>$?vGiyty-4MlCc5IQLl%DggXdDKZ&HIWRJDj z$Gp@F;~s=u4b^9M6SvoI-voVQfA_A_PSK9D5=IAnQwDU+5+C<*Ac@G^d%ltQ{9MEh z*V?-2(&9%(*zXVORd$`T6+N2jrbkA-eU~@ob)~^o4llKrhO1F@-G?#bG>w%Iop$DH zpDxd#l3cT&Pm4;-Q4U9vma}GET_J-*U9Cq2 zGBOg!gc#AeGWuNbto#rEtid)uT&Xc)sYypt2)QzQf>1E?J9 zyJ0>!nS@(grKK-k;@j8)TsGP=9igxv3{2kJN2eGzClIIuGo)ot~tsQr?aon0}dcYU3hakX-D zIU7$z%yZBTuQ@yY4v`3knoAE4RyQ<)nSRgg+3X{haq&llzI+$UZoGNYa-3H+F6J+$ zmh&TEm)X{cF2`#l44&IKmJ4ZI@j;c9d@oc2b2}_few~Hq^wfOZAJ3uvnR{?_1oJr5 z=efl98ItjAtpD}^c6-RNurO;V<50sL(9}0c2gbBNqo-0Yq>i=pfX4#J-Q&t}j z*kRc34ij;D>U(TbrF%n;1MCh2TyCXFcEs{XM#eV)DI7cdBLU}H`;~?loAkA)qegDM z;M7!M)B#lLhF3B6;=?49tp&^W*fJ6ls7p-WpL>5(kjA;cpEP0Qb+*lhGuBGIxc36qUp&Ja? zer#*ShJXbcxYfCzm@Rq3_e<1$1HuF>J6avG>q~LJeoYH3uR1D#!C+}OPl#@w^6)sW zw0T3M@4V#X-}e^wd3PpQg^DWeqZ|7Mq#$bJKd)K}ygiaWMxi`&gx7AO)gt<0 z490T@eh)R>{>eueLzWS!yX(QxKG?k^pA=o=QgAE6Cg4RBGqGs`t;VBx2c3nwPUK@?LXS~n7KNgnBG@H4nu!SQ8HHwBtw*2>ufbsA? z(u=QwY9=OYUl+x^hE_*LaVU8ofwLh)`LM&1X?o`1+l!dhy!MzSh&zqciUBG_t=@geblupr-yY76^yuO5fGy zE>Hgqs<)M1?x0OcDXoZGt-%H>cCH4#a^jtwnwS$Fd_r;56;k!9+e4h4-D$1!IiOw|bIWKm_J|6rb@8opb9X8)V zTh&&euGW8X@^c>?aHZMX*X%1v$!*ew;pT=#D~jOumwLi@6A@%UzBv^ zjn;w&Wmky8eJ)`NR7|k%305Pqs{o$5@3lo#c z52~WCT~@vGP()MPvhRBUX_i(AQc(8@xHVB}Y4LvHb!D*_SCqxr*bH0Je<|MH4yg69 z?D-}V85!L+sR(BQBHxE(E5!m?0X|S~zlos;bkaL$t%8dd>lGq}M*EuDs zp#GusIZ> zeD@NOk?}?^>ePisA~J*eowgaka1qkbla>*TIWqn4Vr6B%u+X)zOO@9`5XXQB(x{go@S+GIZ%UDs4F5-j-JMZI1%> z?dl zEG#mvFH004C}3%Zv`V{S-r0Xo;%rE}x0jqz*8KSV+{(T$3_mhj_D#ZRtO#6y9O)x3 z|CFVl_FU7FmD@;OlMHZq;B%s@u<;zN75tv=o{5?gQ)=!!dDqZ5{z43V0QFKD~W}$O4OU za=?YEM)&;#$mEhMF{0x-^9i(rjt^gI^|(T8Xn7gXv8#9C-LhC9Mix+68I{;OZzEr0 z`pRaITzFmuAt2zcAZBI}e`;svbTM*ZwnfRTt&`}ZcCtaF(6i=Fq0Cdz!q1=vAvgE@ zc3LhjDhXZhR9j_S0#In^Q^7|dOb}tPT8Tu|{CdzR=*6XJldPxkST2s#{)`@HwSEfg z_s0758*c*I2Z4d#Mjd?&)vkm#6%ahX4mF&)1qrh(IbzA6{2!-cXulX#h3V3 zS56yPc0zTHU`C8jpKu8kHGH|_DN5%qo*D6N()>&OG>p}j}IS~Ga! z<7+A;S0WcE{p(jUy92jk|Ve2u_d5n_rJR|GoA7Qks33 z!EBX&y;Nh(ls$(t#2?7sT5KGjDt}7IXSSgIUH4=ygxJT2yokZU6cB@qR^aum8tBv$ zuz9?cpy_V$@JnOv_wFOedfJnd_Ag|NLNg0XbySUO<6pdDb^jJGL`LGOR|O09{$pBP zr$qg>gHc2hNyaT~w9@nWF&8GI1s8A7c;W3KS?%f{cHaIum<_DcBYE;W=e((_FNRXG z9cz*ToEmo-JrLPv-fN>z>F1F-N<{FDXI_KiMh@*~#7c%+**E;#)yem0*;|=ObidAx z?5nR}2YY*BF25s(K46&eQW1bcde(uwKQzCMJSbG9q$fwxx9x+Zgwv&roLq=n$zrQ# zRoeFmoJvU%A0%5R7w2V=f$oirxp`sF`usm_mO(fcwU@wAUue`%o&W$CYmC1r{8CQ| z8B}0h75hZv5VbvAG&C%m&bJyrq}!X)b7pA4Gj@a z-qh|I8BTsa=mi|OnD^H|jSvIKZ|E0G>MX<{FB3uP=CF_SPWUm>yi~I25o&**dA#;1 zH#b+U1k;573jn3M`7P-c+)3p9xurkh%ad;jIAHE=ZYV|713j^}R;11!wHT@QuyaIE z+!qCdv{%kO)PE1dt5kDnv67SixVvwj24Z53C<1JeQ z0gTR48Z6aP0jFRTlt*WZ`Jth;Hct`wdzEJ43;d~L29ShHKQoyZ06 z|3pH+YMLCnhh09mK=YS-UABAjLmwAdcj0l~S4FwYCxu>8DayFMi>ZVLzmT61@00U? zm2P)0$qHah%vq{tKcOG1Q`qGkv!`>PO6Jq?{q8pz-XOOBKxa8Kqxm+pgi^NhZR3{b zl#WpUdsjj0-kRYN9r}8Q^t)^IVIunBaUScPKWnFhDIFV#{NE#X2B~p|3sj^oNEC^u zV~ZO_V!Qo#mB1XWn4gHWG%P=p*dn1GjW%|eAsDx$g_~cA#g&#cdk7G|(ESh>% z==%cKsl(@@%liVzLDA$Jqq57WC$i|Bvqp(^UWccvd--!`vDnZ4tG6E!+8w7zhk>Lw zH`V%{bBXmaljmqA_EiMFSAp7`>M{>Kt2rBHvls7L3D+`t=hvsgmn~`+>^z>bf%l_; z;$?Vjsb%E*xMX=Q)`qp%Gxr`|v8&0WU-}S74&R+h-yudbxJ_~cTp=M389W3t#LMWa&kHjFZ^Q_ij9GW@ecdN&v^p_@fO$BmA##L zXOk164#PoSCPp#qrekO&#k*J zCYtZ~CBJ%A%KMLYld(0Bk%OSGW^f?__%06y))$SWghdk^PI3cJQ!6yi>}(2gkNb1= z#V4F-K4X8peB$Kw*JIMRF+dTa-adrC8Lz}+j3DU`+bI499ipykw1&#OmYt7?)gj8N z0k&!jm4-aqji;}89wRkVJEVjdpI%}ohaFno5=+BW9wuL1?SnsKBz>wWaiuQC(woIn znvWjiZg*Fe*BCT#P4YNY9X?V%N3sv~fg(Q55P2C;6gZAJ2xM~M62l5AYcs`1Lu7l- z{j^BJJc-E%Ef95neD7NO>>Boav)KOGzHdY}aQStQm+>JM2?{wo&qsTuS^UZam%0E^2m<6kSx~+Szu6 zMh&9=cGjYCa&=jj0&{NTMa5mBV|&h2!nkIq0J!EWqyJ&C+${>=G4wjHjizBWi5uYw QoDhnpqM#;UDQgn^AHCA{WB>pF literal 0 HcmV?d00001 diff --git a/doc/news/changes/minor/20170801DenisDavydov-c b/doc/news/changes/minor/20170801DenisDavydov-c new file mode 100644 index 0000000000..8d0d164fda --- /dev/null +++ b/doc/news/changes/minor/20170801DenisDavydov-c @@ -0,0 +1,4 @@ +New: add Utilities::LinearAlgebra::Lanczos_largest_eigenvalue() to estimate the largest +eigenvalue of a symmetric linear operator using k-steps of Lanczos algorithm. +
+(Denis Davydov, 2017/08/01) diff --git a/doc/news/changes/minor/20170802DenisDavydov b/doc/news/changes/minor/20170802DenisDavydov new file mode 100644 index 0000000000..3aee291080 --- /dev/null +++ b/doc/news/changes/minor/20170802DenisDavydov @@ -0,0 +1,4 @@ +New: add Utilities::LinearAlgebra::Chebyshev_filter() to apply Chebyshev filter of a +given degree. +
+(Denis Davydov, 2017/08/02) diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h index 621ece0709..83480f965b 100644 --- a/include/deal.II/lac/diagonal_matrix.h +++ b/include/deal.II/lac/diagonal_matrix.h @@ -179,8 +179,11 @@ public: const VectorType &src) const; /** - * Initialize vector @p dst. This is a part of the interface required - * by linear operator. + * Initialize vector @p dst to have the same size and partition as + * @p diagonal member of this class. + * + * This is a part of the interface required + * by linear_operator(). */ void initialize_dof_vector(VectorType &dst) const; diff --git a/include/deal.II/lac/utilities.h b/include/deal.II/lac/utilities.h new file mode 100644 index 0000000000..8383159838 --- /dev/null +++ b/include/deal.II/lac/utilities.h @@ -0,0 +1,318 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__lac_utilities_h +#define dealii__lac_utilities_h + +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Utilities +{ + /** + * A collection of linear-algebra utilities. + */ + namespace LinearAlgebra + { + + /** + * Estimate an upper bound for the largest eigenvalue of @p H by a @p k -step + * Lanczos process starting from the initial vector @p v0. Typical + * values of @p k are below 10. This estimator computes a k-step Lanczos + * decomposition $H V_k=V_k T_k+f_k e_k^T$ where $V_k$ contains k Lanczos basis, + * $V_k^TV_k=I_k$, $T_k$ is the tridiagonal Lanczos matrix, $f_k$ is a residual + * vector $f_k^TV_k=0$, and $e_k$ is the k-th canonical basis of $R^k$. + * The returned value is $ ||T_k||_2 + ||f_k||_2$. + * If @p eigenvalues is not nullptr, the eigenvalues of $T_k$ will be written there. + * + * @p vector_memory is used to allocate memory for temporary vectors. + * OperatorType has to provide vmult operation with + * VectorType. + * + * This function implements the algorithm from + * @code{.bib} + * @Article{Zhou2006, + * Title = {Self-consistent-field Calculations Using Chebyshev-filtered Subspace Iteration}, + * Author = {Zhou, Yunkai and Saad, Yousef and Tiago, Murilo L. and Chelikowsky, James R.}, + * Journal = {Journal of Computational Physics}, + * Year = {2006}, + * Volume = {219}, + * Pages = {172--184}, + * } + * @endcode + * + * @note This function uses Lapack routines to compute the largest + * eigenvalue of $T_k$. + * + * @note This function provides an alternate estimate to that obtained from + * several steps of SolverCG with SolverCG::connect_eigenvalues_slot(). + * + * @author Denis Davydov, 2017 + */ + template + double Lanczos_largest_eigenvalue(const OperatorType &H, + const VectorType &v0, + const unsigned int k, + VectorMemory &vector_memory, + std::vector *eigenvalues = nullptr); + + /** + * Apply Chebyshev polynomial of the operator @p H to @p x. For a + * non-defective operator $H$ with a complete set of eigenpairs + * $H \psi_i = \lambda_i \psi_i$, the action of a polynomial filter $p$ is given by + * $p(H)x =\sum_i a_i p(\lambda_i) \psi_i$, where $x=: \sum_i a_i \psi_i$. Thus + * by appropriately choosing the polynomial filter, one can alter + * the eigenmodes contained in $x$. + * + * This function uses Chebyshev polynomials of first kind. Below is an + * example of polynomial $T_n(x)$ of degree $n=8$ normalized to unity at $-1.2$. + * + * + * + * + *
+ * @image html chebyshev8.png + *
+ * By introducing a linear mapping $L$ from @p unwanted_spectrum to + * $[-1,1]$, we can dump the corresponding modes in @p x. The higher + * the polynomial degree $n$, the more rapid it grows outside of the + * $[-1,1]$. In order to avoid numerical overflow, we normalize + * polynomial filter to unity at @p tau. Thus, the filtered operator + * is $p(H) = T_n(L(H))/T_n(L(\tau))$. + * + * The action of the Chebyshev filter only requires + * evaluation of vmult() of @p H and is based on the + * recursion equation for Chebyshev polynomial of degree $n$: + * $T_{n}(x) = 2x T_{n-1}(x) - T_{n-2}(x)$ with $T_0(x)=1$ and $T_1(x)=x$. + * + * @p vector_memory is used to allocate memory for temporary objects. + * + * This function implements the algorithm (with a minor fix of sign of $\sigma_1$) from + * @code{.bib} + * @Article{Zhou2014, + * Title = {Chebyshev-filtered subspace iteration method free of sparse diagonalization for solving the Kohn--Sham equation}, + * Author = {Zhou, Yunkai and Chelikowsky, James R and Saad, Yousef}, + * Journal = {Journal of Computational Physics}, + * Year = {2014}, + * Volume = {274}, + * Pages = {770--782}, + * } + * @endcode + * + * @note If @p tau is equal to + * std::numeric_limits::infinity(), no normalization + * will be performed. + * + * @author Denis Davydov, 2017 + */ + template + void Chebyshev_filter(VectorType &x, + const OperatorType &H, + const unsigned int n, + const std::pair unwanted_spectrum, + const double tau, + VectorMemory &vector_memory); + + } + +} + + +/*------------------------- Implementation ----------------------------*/ + +#ifndef DOXYGEN + +namespace Utilities +{ + namespace LinearAlgebra + { + template + double Lanczos_largest_eigenvalue(const OperatorType &H, + const VectorType &v0_, + const unsigned int k, + VectorMemory &vector_memory, + std::vector *eigenvalues) + { + // Do k-step Lanczos: + + typename VectorMemory::Pointer v(vector_memory); + typename VectorMemory::Pointer v0(vector_memory); + typename VectorMemory::Pointer f(vector_memory); + + v->reinit(v0_); + v0->reinit(v0_); + f->reinit(v0_); + + // two vectors to store diagonal and subdiagonal of the Lanczos + // matrix + std::vector diagonal; + std::vector subdiagonal; + + // scalars to store norms and inner products + double a = 0, b = 0; + + // 1. Normalize input vector + (*v) = v0_; + a = v->l2_norm(); + Assert (a!=0, ExcDivideByZero()); + (*v) *= 1./a; + + // 2. Compute f = Hv; a = f*v; f <- f - av; T(0,0)=a; + H.vmult(*f,*v); + a = (*f)*(*v); + f->add(-a,*v); + diagonal.push_back(a); + + // 3. Loop over steps + for (unsigned int i = 1; i < k; ++i) + { + // 4. L2 norm of f + b = f->l2_norm(); + Assert (b!=0, ExcDivideByZero()); + // 5. v0 <- v; v <- f/b + *v0 = *v; + *v = *f; + (*v) *= 1./b; + // 6. f = Hv; f <- f - b v0; + H.vmult(*f,*v); + f->add(-b, *v0); + // 7. a = f*v; f <- f - a v; + a = (*f) * (*v); + f->add(-a, *v); + // 8. T(i,i-1) = T(i-1,i) = b; T(i,i) = a; + diagonal.push_back(a); + subdiagonal.push_back(b); + } + + Assert (diagonal.size() == k, + ExcInternalError()); + Assert (subdiagonal.size() == k-1, + ExcInternalError()); + + // Use Lapack dstev to get ||T||_2 norm, i.e. the largest eigenvalue + // of T + const int n = k; + std::vector Z; // unused for eigenvalues-only ("N") job + const int ldz = 1; // ^^ (>=1) + std::vector work; // ^^ + int info; + // call lapack_templates.h wrapper: + stev ("N", &n, + &diagonal[0], &subdiagonal[0], + &Z[0], &ldz, &work[0], + &info); + + Assert (info == 0, + LAPACKSupport::ExcErrorCode("dstev", info)); + + if (eigenvalues != nullptr) + { + eigenvalues->resize(diagonal.size()); + std::copy(diagonal.begin(), diagonal.end(), + eigenvalues->begin()); + } + + // note that the largest eigenvalue of T is below the largest + // eigenvalue of the operator. + // return ||T||_2 + ||f||_2, although it is not guaranteed to be an upper bound. + return diagonal[k-1] + f->l2_norm(); + } + + + template + void Chebyshev_filter(VectorType &x, + const OperatorType &op, + const unsigned int degree, + const std::pair unwanted_spectrum, + const double a_L, + VectorMemory &vector_memory) + { + const double a = unwanted_spectrum.first; + const double b = unwanted_spectrum.second; + Assert (degree > 0, + ExcMessage ("Only positive degrees make sense.")); + + const bool scale = (a_L < std::numeric_limits::infinity()); + Assert (a < b, + ExcMessage("Lower bound of the unwanted spectrum should be smaller than the upper bound.")); + + Assert (a_L <= a || a_L >= b || !scale, + ExcMessage("Scaling point should be outside of the unwanted spectrum.")); + + // Setup auxiliary vectors: + typename VectorMemory::Pointer p_y(vector_memory); + typename VectorMemory::Pointer p_yn(vector_memory); + + p_y->reinit(x); + p_yn->reinit(x); + + // convenience to avoid pointers + VectorType &y = *p_y; + VectorType &yn = *p_yn; + + // Below is an implementation of + // Algorithm 3.2 in Zhou et al, Journal of Computational Physics 274 (2014) 770-782 + // with **a bugfix for sigma1**. Here is the original algorithm verbatim: + // + // [Y]=Chebyshev_filter_scaled(X, m, a, b, aL). + // e=(b−a)/2; c=(a+b)/2; σ=e/(c−aL); τ=2/σ; + // Y=(H∗X−c∗X)∗(σ/e); + // for i=2 to m do + // σnew =1/(τ −σ); + // Yt =(H∗Y−c∗Y)∗(2∗σnew/e)−(σ∗σnew)∗X; + // X =Y; Y =Yt; σ =σnew; + + const double e = (b-a)/2.; + const double c = (a+b)/2.; + const double alpha = 1./e; + const double beta = - c/e; + + const double sigma1 = e/(a_L - c); // BUGFIX which is relevant for odd degrees + double sigma = scale ? sigma1 : 1.; + const double tau = 2./sigma; + op.vmult(y,x); + y.sadd(alpha*sigma, beta*sigma, x); + + for (unsigned int i = 2; i <= degree; ++i) + { + const double sigma_new = scale ? 1./(tau-sigma) : 1.; + op.vmult(yn,y); + yn.sadd(2.*alpha*sigma_new, 2.*beta*sigma_new, y); + yn.add(-sigma*sigma_new, x); + x.swap(y); + y.swap(yn); + sigma = sigma_new; + } + + x.swap(y); + } + + } +} + +#endif + + + +DEAL_II_NAMESPACE_CLOSE + + +#endif diff --git a/tests/lac/utilities_01.cc b/tests/lac/utilities_01.cc new file mode 100644 index 0000000000..f747722c37 --- /dev/null +++ b/tests/lac/utilities_01.cc @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test estimate of largest eigenvalue of a matrix. +// The matrix is the same as in slepc/solve_04 test which has eingenvalues: +// 3.98974 > 3.95906 > 3.90828 > 3.83792 + +#include "../tests.h" +#include "../testmatrix.h" +#include "../slepc/testmatrix.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char **argv) +{ + initlog(); + deallog << std::setprecision(6); + deallog.threshold_double(1.e-10); + + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + + const unsigned int size = 31; + unsigned int dim = (size-1); + + deallog << "Size " << size << " Unknowns " << dim << std::endl << std::endl; + + // Make matrix + FD1DLaplaceMatrix testproblem(size); + PETScWrappers::SparseMatrix A(dim, dim, 3); + testproblem.three_point(A); + A.compress (VectorOperation::insert); + + PETScWrappers::MPI::Vector v0(MPI_COMM_WORLD, dim, dim); + PETScWrappers::MPI::Vector y(MPI_COMM_WORLD, dim, dim); + PETScWrappers::MPI::Vector x(MPI_COMM_WORLD, dim, dim); + for (unsigned int j=0; j(Testing::rand())/static_cast(RAND_MAX); + + v0.compress(VectorOperation::insert); + GrowingVectorMemory vector_memory; + + for (unsigned int k = 4; k < 10; ++k) + { + const double est = Utilities::LinearAlgebra::Lanczos_largest_eigenvalue(A,v0,k,vector_memory); + Assert (est > 3.98974, ExcInternalError()); + deallog << k << std::endl + << "Lanczos " << est << std::endl; + + // estimate from CG + { + ReductionControl control (k, + std::sqrt(std::numeric_limits::epsilon()), + 1e-10, false, false); + std::vector estimated_eigenvalues; + SolverCG solver (control); + solver.connect_eigenvalues_slot([&estimated_eigenvalues] (const std::vector &ev) -> void {estimated_eigenvalues = ev;}); + y = v0; + PreconditionIdentity preconditioner; + try + { + solver.solve(A, x, y, preconditioner); + } + catch (SolverControl::NoConvergence &) + { + } + + deallog << "CG " << estimated_eigenvalues.back() << std::endl; + } + } + } + +} diff --git a/tests/lac/utilities_01.output b/tests/lac/utilities_01.output new file mode 100644 index 0000000000..fd5aca82d7 --- /dev/null +++ b/tests/lac/utilities_01.output @@ -0,0 +1,21 @@ + +DEAL::Size 31 Unknowns 30 +DEAL:: +DEAL::4 +DEAL::Lanczos 4.83338 +DEAL::CG 3.72123 +DEAL::5 +DEAL::Lanczos 4.95881 +DEAL::CG 3.90336 +DEAL::6 +DEAL::Lanczos 4.81862 +DEAL::CG 3.82808 +DEAL::7 +DEAL::Lanczos 4.98166 +DEAL::CG 3.94033 +DEAL::8 +DEAL::Lanczos 4.70266 +DEAL::CG 3.89841 +DEAL::9 +DEAL::Lanczos 4.93451 +DEAL::CG 3.97020 diff --git a/tests/lac/utilities_02.cc b/tests/lac/utilities_02.cc new file mode 100644 index 0000000000..b754902120 --- /dev/null +++ b/tests/lac/utilities_02.cc @@ -0,0 +1,299 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2017 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * Test estimated largest eigenvalue of M^{-1/2} L M^{-1/2} using + * k-steps of Lanczos algorithm. Here M is diagoal mass matrix obtained + * from Gauss-Legendre-Lobatto quadrature and L is Laplace operator. + * + * Largest eigenvalues from pArpack are: + * + * DEAL::1014.26 + * DEAL::1018.29 + * DEAL::1018.29 + * DEAL::1020.72 + * DEAL::1020.72 + * + */ + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +// debug only: turn on calculation of eigenvalues of the operator by Arpack +// #define PARPACK + +#ifdef PARPACK +#include +#endif + +#include +#include + + +const unsigned int dim = 2; + +using namespace dealii; + +const double eps = 1e-10; + +const unsigned int fe_degree = 1; + +void test () +{ + const unsigned int global_mesh_refinement_steps = 5; + + MPI_Comm mpi_communicator = MPI_COMM_WORLD; + const unsigned int n_mpi_processes = Utilities::MPI::n_mpi_processes(mpi_communicator); + const unsigned int this_mpi_process = Utilities::MPI::this_mpi_process(mpi_communicator); + + parallel::distributed::Triangulation triangulation (mpi_communicator); + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (global_mesh_refinement_steps); + + + DoFHandler dof_handler(triangulation); + FE_Q fe(fe_degree); + dof_handler.distribute_dofs (fe); + + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + ConstraintMatrix constraints; + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction (), + constraints); + constraints.close (); + + std::shared_ptr > mf_data(new MatrixFree ()); + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.mapping_update_flags = update_values | update_gradients | update_JxW_values; + mf_data->reinit (dof_handler, constraints, quad, data); + } + + MatrixFreeOperators::MassOperator > mass; + MatrixFreeOperators::LaplaceOperator > laplace; + mass.initialize(mf_data); + laplace.initialize(mf_data); + + // Gauss-Legendre-Lobatto mass-matrix: + DiagonalMatrix> diagonal_mass_inv; + { + LinearAlgebra::distributed::Vector inv_mass_matrix; + VectorizedArray one = make_vectorized_array (1.); + mf_data->initialize_dof_vector (inv_mass_matrix); + FEEvaluation fe_eval(*mf_data); + const unsigned int n_q_points = fe_eval.n_q_points; + for (unsigned int cell=0; celln_macro_cells(); ++cell) + { + fe_eval.reinit(cell); + for (unsigned int q=0; q1e-15) + { + inv_mass_matrix.local_element(k) = std::sqrt(1./inv_mass_matrix.local_element(k)); + } + else + inv_mass_matrix.local_element(k) = 0; + + diagonal_mass_inv.reinit(inv_mass_matrix); + } + + const auto invM = linear_operator>(diagonal_mass_inv); + const auto OP = invM * + linear_operator>(laplace) * + invM; + + // Do actuall work: + LinearAlgebra::distributed::Vector init_vector; + mf_data->initialize_dof_vector(init_vector); + for (auto it = init_vector.begin(); it != init_vector.end(); ++it) + *it = static_cast(Testing::rand())/static_cast(RAND_MAX); + + constraints.set_zero(init_vector); + + GrowingVectorMemory> vector_memory; + for (unsigned int k = 4; k < 10; ++k) + { + const double est = Utilities::LinearAlgebra::Lanczos_largest_eigenvalue(OP,init_vector,k,vector_memory); + deallog << k << " " << est << std::endl; + } + + // exact eigenvectors via PArpack +#ifdef PARPACK + { + const unsigned int number_of_eigenvalues = 5; + + std::vector > eigenfunctions; + std::vector eigenvalues; + eigenfunctions.resize (number_of_eigenvalues); + eigenvalues.resize (number_of_eigenvalues); + for (unsigned int i=0; iinitialize_dof_vector (eigenfunctions[i]); + + std::vector > lambda(number_of_eigenvalues); + + const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 10; + PArpackSolver >::AdditionalData + additional_data(num_arnoldi_vectors, + PArpackSolver >::largest_magnitude, + true, + 1); + + SolverControl solver_control( + dof_handler.n_dofs(), 1e-10, /*log_history*/ false, /*log_results*/ false); + + PArpackSolver > eigensolver( + solver_control, mpi_communicator, additional_data); + + eigensolver.reinit(eigenfunctions[0]); + // make sure initial vector is orthogonal to the space due to constraints + { + LinearAlgebra::distributed::Vector init_vector; + mf_data->initialize_dof_vector(init_vector); + for (auto it = init_vector.begin(); it != init_vector.end(); ++it) + *it = static_cast(Testing::rand())/static_cast(RAND_MAX); + + constraints.set_zero(init_vector); + eigensolver.set_initial_vector(init_vector); + } + // avoid output of iterative solver: + const unsigned int previous_depth = deallog.depth_file(0); + eigensolver.solve (OP, + mass, + OP, + lambda, + eigenfunctions, + eigenvalues.size()); + deallog.depth_file(previous_depth); + + for (unsigned int i = 0; i < lambda.size(); i++) + eigenvalues[i] = lambda[i].real(); + + for (unsigned int i=0; i < eigenvalues.size(); i++) + deallog << eigenvalues[i] << std::endl; + + // make sure that we have eigenvectors and they are mass-orthonormal: + // a) (A*x_i-\lambda*x_i).L2() == 0 + // b) x_j*x_i=\delta_{ij} + { + const double precision = 1e-7; + LinearAlgebra::distributed::Vector Ax(eigenfunctions[0]); + for (unsigned int i=0; i < eigenfunctions.size(); ++i) + { + for (unsigned int j=0; j < eigenfunctions.size(); j++) + { + const double err = std::abs( eigenfunctions[j] * eigenfunctions[i] - (i==j)); + Assert( err< precision, + ExcMessage("Eigenvectors " + + Utilities::int_to_string(i) + + " and " + + Utilities::int_to_string(j) + + " are not orthonormal: " + + std::to_string(err))); + } + + OP.vmult(Ax,eigenfunctions[i]); + Ax.add(-1.0*eigenvalues[i],eigenfunctions[i]); + const double err = Ax.l2_norm(); + Assert (err < precision, + ExcMessage("Returned vector " + + Utilities::int_to_string(i) + + " is not an eigenvector: " + + std::to_string(err))); + } + } + } +#endif + + dof_handler.clear (); + deallog << "Ok"< +#include +#include +#include + +double cheb2(const unsigned int d, const double x) +{ + if (d == 0) + { + return 1.; + } + else if (d == 1) + { + return x; + } + + return 2.*x*cheb2(d-1,x) - cheb2(d-2,x); +} + + + +void +check(const int degree, const bool scale = false, const double a_L = -0.1, const double a = -0.01, const double b = 0.01, const unsigned int size = 1000) +{ + deallog << "Degree " << degree << std::endl; + LinearAlgebra::distributed::Vector ev(size), x(size), y(size), exact(size), diff(size); + GrowingVectorMemory> vector_memory; + + for (unsigned int i=0; i> mat; + mat.reinit(ev); + + x = 0.; + // prevent overflow by not perturbing modes far away from the region + // to be filtered + unsigned int n_in = 0; + unsigned int n_out = 0; + for (unsigned int i=0; i= a && ev(i) <= b) + n_in++; + else + n_out++; + x(i) = static_cast(Testing::rand())/static_cast(RAND_MAX); + } + + deallog << " Modes inside/outside: " << n_in << " " << n_out << std::endl; + + // for x = x_i v_i , where v_i are eigenvectors + // p[H]x = \sum_i x_i p(\lambda_i) v_i + const double c = (a+b)/2.; + const double e = (b-a)/2.; + auto L = [&](const double &x) + { + return (x-c)/e; + }; + + const double scaling = scale ? cheb2(degree,L(a_L)) : 1.; // p(L(a_L)) + deallog << " Scaling: " << scaling << " @ " << a_L << std::endl; + exact = 0.; + for (unsigned int i=0; i::infinity() ); + y = x; + Utilities::LinearAlgebra::Chebyshev_filter(y, mat, degree, std::make_pair(a, b), g_, vector_memory); + diff = y; + diff -=exact; + + deallog << " Filter ["< std::abs(diff(max_i))) + max_i = i; + + deallog << " i =" << max_i << std::endl + << " d =" << diff(max_i) << std::endl + << " ev=" << ev(max_i) << std::endl + << " x =" << x(max_i) << std::endl + << " y =" << y(max_i) << std::endl + << " ex=" << exact(max_i) << std::endl; +#endif +} + + +int main() +{ + std::ofstream logfile("output"); + deallog << std::setprecision(6); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + deallog << "No scaling:" << std::endl; + // no scaling: + check(1); + check(2); + check(3); + check(4); + check(10); + + deallog << "Lower scaling:" << std::endl; + // scaling at the lower end + check(1,true); + check(2,true); + check(3,true); + check(4,true); + check(10,true); + check(30,true); + + deallog << "Upper scaling:" << std::endl; + // scaling at the upper end + check(1,true, 0.1); + check(2,true, 0.1); + check(3,true, 0.1); + check(4,true, 0.1); + check(10,true, 0.1); + check(30,true, 0.1); + + return 0; +} diff --git a/tests/lac/utilities_03.output b/tests/lac/utilities_03.output new file mode 100644 index 0000000000..5c0d203ca7 --- /dev/null +++ b/tests/lac/utilities_03.output @@ -0,0 +1,123 @@ + +DEAL::No scaling: +DEAL::Degree 1 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 1.00000 @ -0.100000 +DEAL:: Input norm: 6.16127 +DEAL:: Exact norm: 37.4711 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 2 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 1.00000 @ -0.100000 +DEAL:: Input norm: 5.92779 +DEAL:: Exact norm: 586.339 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 3 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 1.00000 @ -0.100000 +DEAL:: Input norm: 5.89974 +DEAL:: Exact norm: 10657.5 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 4 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 1.00000 @ -0.100000 +DEAL:: Input norm: 5.60617 +DEAL:: Exact norm: 157052. +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 10 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 1.00000 @ -0.100000 +DEAL:: Input norm: 5.36039 +DEAL:: Exact norm: 5.05476e+12 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0.00146484 +DEAL::Lower scaling: +DEAL::Degree 1 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: -10.0000 @ -0.100000 +DEAL:: Input norm: 6.43748 +DEAL:: Exact norm: 3.99850 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 2 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 199.000 @ -0.100000 +DEAL:: Input norm: 5.74697 +DEAL:: Exact norm: 2.70598 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 3 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: -3970.00 @ -0.100000 +DEAL:: Input norm: 5.64201 +DEAL:: Exact norm: 2.34171 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 4 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 79201.0 @ -0.100000 +DEAL:: Input norm: 5.70416 +DEAL:: Exact norm: 1.92446 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 10 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 4.99312e+12 @ -0.100000 +DEAL:: Input norm: 5.64307 +DEAL:: Exact norm: 0.989370 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 30 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 4.97938e+38 @ -0.100000 +DEAL:: Input norm: 5.82031 +DEAL:: Exact norm: 0.745548 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Upper scaling: +DEAL::Degree 1 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 10.0000 @ 0.100000 +DEAL:: Input norm: 5.51882 +DEAL:: Exact norm: 3.18201 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 2 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 199.000 @ 0.100000 +DEAL:: Input norm: 5.83880 +DEAL:: Exact norm: 2.78386 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 3 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 3970.00 @ 0.100000 +DEAL:: Input norm: 5.49421 +DEAL:: Exact norm: 2.11010 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 4 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 79201.0 @ 0.100000 +DEAL:: Input norm: 5.55782 +DEAL:: Exact norm: 1.56672 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 10 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 4.99312e+12 @ 0.100000 +DEAL:: Input norm: 5.85002 +DEAL:: Exact norm: 1.24742 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 +DEAL::Degree 30 +DEAL:: Modes inside/outside: 10 90 +DEAL:: Scaling: 4.97938e+38 @ 0.100000 +DEAL:: Input norm: 6.25978 +DEAL:: Exact norm: 0.896899 +DEAL:: Filter [-0.0100000,0.0100000] +DEAL:: Error: 0 -- 2.39.5