From 88087ea6a62e23b7707397320c4488f90da72ff4 Mon Sep 17 00:00:00 2001 From: Bausager Date: Sat, 13 Sep 2025 21:44:20 +0200 Subject: [PATCH] Finishing up and starting lu decomp --- bin/abc_lab | Bin 80104 -> 16512 bytes include/core/omp_config.h | 65 +- include/decomp/.gitkeep | 0 include/decomp/decomp.h | 4 + include/decomp/lu.h | 61 ++ include/numerics/inverse.h | 95 +-- .../numerics/inverse/inverse_gauss_jordan.h | 94 +++ include/numerics/inverse/inverse_lu.h | 0 include/numerics/matmul.h | 87 ++- include/numerics/matvec.h | 100 +++- include/numerics/transpose.h | 22 - include/utils/numerics.txt | 11 - include/utils/vector.h | 14 + makefile | 67 ++- obj/main.o | Bin 130544 -> 2816 bytes src/main.cpp | 559 +----------------- test/test_all.cpp | 2 + test/test_common.h | 52 ++ test/test_inverse.cpp | 126 ++++ test/test_matmul.cpp | 164 +++++ test/test_matrix.cpp | 142 +++++ test/test_matvec.cpp | 237 ++++++++ test/test_transpose.cpp | 88 +++ test/test_vector.cpp | 211 +++++++ 24 files changed, 1502 insertions(+), 699 deletions(-) create mode 100644 include/decomp/.gitkeep create mode 100644 include/decomp/decomp.h create mode 100644 include/decomp/lu.h create mode 100644 include/numerics/inverse/inverse_gauss_jordan.h create mode 100644 include/numerics/inverse/inverse_lu.h delete mode 100644 include/utils/numerics.txt create mode 100644 test/test_all.cpp create mode 100644 test/test_common.h create mode 100644 test/test_inverse.cpp create mode 100644 test/test_matmul.cpp create mode 100644 test/test_matrix.cpp create mode 100644 test/test_matvec.cpp create mode 100644 test/test_transpose.cpp create mode 100644 test/test_vector.cpp diff --git a/bin/abc_lab b/bin/abc_lab index 05cd2c8ebf1e5e0f547290cd7bb49080cddfba43..b4fd476313e9fd37b7424b33d343f1635e2a0d2c 100755 GIT binary patch literal 16512 zcmeHOdu$ZP8K1Lxm;~?_A%cY9z* zh?qdyT7)WTtEOraiK?xN)JnA~R3%l_8c|5omqMyaNbMgaDvFWP2AZ^jR^|HpX69S3 zH+Qs(s{B#ySbM+ue)BzMzn$Hi_02rqwRLk4x!jx&nxabz%-a zmx;T@49L?Yr`lr*(5lMCs;O9``D!59Ev8BtT&ifn$UP(^yGmv8Ox40DqDm&a=~UHn z5Fge3WaM@)oMBioTkc4OLZOxrOE*4N1S zq@2*PPUka-jcB977_=MEc0ue^W`dDoP}z0aO}w9YJ=$)XA>D`|avVmE_ZiqxUjEy{ zmwBg-SLr%Dtn<b&8rjsykE6 z(HFM8KaDwk_=L-feawnNRLd72IULciX0PE}QIi`qD|&R=#C#BHeQ!BkaAq z+~%fKu3+~$g`|yZo_F%Yb}E~44Na1CQAKPAUxo3UM7eaS5s%_cY-LdKZ`|TF9)of`NsJ3@>+f8PRxn9UNyV-2@YAuu1 zW^*kip&Ix0bX@O<{G%}t!AN)kNcnEW(Q*A*&@;r3G&g+V2e;D`;T&Wt5%~B%ko2i6 zNTue`=zR@!;!Z#FH4LX}@kcEe>r%65w7m(L3O~R2xxs_zlzu)jKd!$wYGTq%j#a9$ zd~u`VY3#BMF_(NijTJ5vKA!G1E?0fLx(`&#bsw+43S`4CG*2AmkISn*p6^|53o*<< zn1L_@VFtnsgc%4k@c)>BcWUqchjsL=8teG<*VhPPJvr`Hm0q=u{=DX_Jl)dTFM}?v z_$oeYm&F9>4^YdMH%g^a|I<>Y3EP#|Jh`;uFtqz?mvsuMUs|y*An%zX?>=|q?Qv`L zvUT)Nm#iaiorJf9b$r^2m8e*w$LhZhY#KjzAjslWCKf>k%I(W};(zihP5T(7nw=dBAj-34&4O*^P5T~hsG zd&>XFhAuRV;_BVj(G68pMKayl*6|JZLYPS1ER`k_=ZJPvxs{zYK^{ZtvY{YEqus)oSi&XFF$4Ur+fP7jr2)l%o|l8!xCm7%s`ldFau!* z!VH8N2s037Ak09RfiMFfXa?waq{x!$PsR|Hej8bOwN!F}TY=959|WET9s<5g?bk}B z2Z6KyS}Oe-cpvZ*kbc{w5_xQ=hzvy|OJ>ikc^Yrn#4pA1BJA%ah0dAsC#ZT=Pv$q?d=C^M#RRt^&f4imp}iJJ;gf(J!yzK72;hSWN!eukj_&$)}f#G zXp*>-IC_1Y0r@nJ5!inND7xlGKVP+F&Wx%r-iEHX7XmD7+T$qHe!U(c&GA{9d(K%L_F!@EW`^y)KW} z>dALRDA9W#6<*u7NQq)X2WEbymbu^KS`UjFW#RaDYnk=_SyJbRyh;kcW6JY=QJ*M& zDU~w!3$L|O{?BY5p^OOQw-lv`y~G%!S{h`72-)$C@!ymW&RndztCGJ z%lsHB<@#UgynZJAED%fYVGjx``JvWw{(mFQ=8I{Dtau-nP>YB~A~=q&NPe;4`;HCh z0!GE1UVfOLfke%#nlEY$W)#+eUk;0ye%|nDEW!kPc)o&p`luSeM`Y|Fpz%QSw zDy4_O-vR$4dYsGu-Jqq9{8ZiPtutoRoaEoz&nGJI&yxSDd4ivp?B@mIS@7YH!PEIh z1LwO0=S%%G^f;9N`-0u?D)^bGz*m!>Dm>No`Eux4;OTs$F=iAl1W(sX1bO(6Vw^7y z@XMv2_vWV!e7W=5A^lI@FC6^=@a6LJ1@Hzd`S4rXe=Lxcu(x6RF1PJKGLAFr&g^o!GhJPs&0W~n(7gkD7!s+hT`VM_ zH)Rh)jJG4u&I7&2AyLTL15P$UTPnOj2OQ)}#_=eMjVmCxUQmvp9oR78dB+xqDg8^J ze+l%jiTbBDXT+xLk_c|b@I1j0ZP!37?EUcaNLf)EHKul%yXVgh$CkOV`w9g%yv61a z6@GF@$1S&^_=E9gV4`a7&q1V0AsH8@n;e2D2ZNc<$+2oC2lOp6kbpOqln-^UC=0G5l#@xpnKH)( z5+-ia4DLPO1Lb{IY5t2xBfisB3GPSa`LAx4x?CsxFb&=b1_pq74ws%tf(<)|5@;V$DEb>&+~dlj?eS@VE-N1)102|dEW1# zj!S){JhMH|zf<67D#-Ty9>K`(4N#GVeE7?w5NJ-$_B;<~|~@ zN_8nS}-(wgXefDG%wEreJ$_e}T<_UVf z!0&h3jqEGe37^N|w>PvMBgbbO-Ry<=Zos}rD=@NtCHt{}{lnUT(O^T3?2qxKfc=PW zV0=9=VAzh`s`>V0{&?QU$nVdrAH09g>G(X(grPMx2K8BHqvO#^88Lz>%)@RS}OQ#w@qOhKlCS~UT zg$5eC9G~Ayc>gmehVDNqY{%E}d*JB&J^PKMh5t{iQP`uv`(%4YdXMqji0SQJ0i5QfKB0+*iK@&*iwsa5$Bq$QLASjzj2T&1^q{XJa zijKJJsN;g;zA%alCLx^!+yF-fQQ4GK8<0hX1d#miTXk>u?M|XIely?qJpVJ{A@|m) zs#B-VR;QL*S(M|xEXHoPnLn*;m)Mv}zEwdojf|Cj!(}pU7uwp}GHs{Zl5MS#@4$bv zT=a8u3^MU!KI2r<)(TJ7OX5p2&#nq=GoNPbe6n6kKYbHbj`_4Ds5I-v@->;i@RKZm zhn{UeP-*jNttZ>tZj_IDKGRdd&8M{ zXFmHH{4}2?KZ%b2FEscW{p^c8{+Q3`dee-0(a&sE&1OE?27gBW;FY_ge8v%M=(k9_%8e++iot}@zdNe=ssew$CDdE0q2rcWBs|GXJf&YdxRR`L9E z=U+JB+yVXj6wU6_Pm(~olJTG78#;V69;BmJgm#166bG!xqe9=Pr4lTgTTcGEK7VwX0p?pRQ_-UWX_#geVs)h2`x4{4N zTPQ!c1^5qWuw!(iNYf7_Hv_-A1w7}pP<~wt@;U@OJ(|H^48EOiOR`;Pj892M83)RD z#HgK_Ov>QT)Wh~7_`veXQRQv<`Pa{!JuAP+GqKQo;j)(e%k^uyFdU>s^yadeYJ+Pn%en?62#^Of1YnPQQHD>>~d6O!R>8Yeu<7dNK;;PcE3_ znLc|K08?g+m}?+SJHN1amS_6R0!8tVw4AGk{Rx6|{#-vVC@d%{D4bi6V^H(&@(JW# zF|t1d_h*Aot87#fH&b$NU)JrK#t6TuqP^YfQIpX z5;SU7{^W_C$O=R^Dj`G^b(Cl=nEKYiA856H{E ze*XOYIR%A9vu915G2L@>{@lLJ(rk1F?GDJ#hhh{KOq>ZEV7{kWWx9$KDt$9byGG@Y z>YJZ4azOrwVI$AapMhpT*oCIrPntL-e_~2o2)g1H4V z!X>S>hhC98B7e@rLa0K)j1Yc>o6XqrN6(r!9YkJLP*gk<=1uQ4`NsUo({9Y4I&u08 z+suNQGiT2&utCW1MaYNB(Bsg>&73%$5US+V>C~883RoW0@-cbN%{CM+@=TfByEmiF z0Zh&ey#e!ajZyV)N_ z+1)hbX|38nR5~sEdcN_RmxpHh_c_0xDr6)xY(w($vIe_=@BBVSYerVTh@_=l1cH8P z95rI(zgC#hv_gdKK&lm@-8Qy3TU(^t;r|J?c+5d#@YccRMARIEbc`&|8m%w}*%8mh z*iHmI7Vq)+&-_^VkK?A@mVwyL9M4CbdGiU30Ivj$KX%~gVdQ5bU3|*)wl-iJYruT> zK7jDMqpiS5FO_MhZH}e>IhOhx&NvezV+Y%Ncw4ENK9lLTw!1?09k!uH-Uzhc zdHVGBHuU;_19rb`FV6Oak$y|2TiaeR(u!V(q!-oEpEkA=jXb6Y*FT66Ji+#s0rM?N zMY_Fht&zTY|2N_R6JBPlJH#ccz`=x@=>!AL+&&V6`t@3ORi4e6b0OxuV}YYXp`S7f z+!|NaSm0ekr7@?nz|XM2H(TJT7Wj4xytf6OWTC&81)gkycela~{(!BaxjVDXxQlQr zAI!LmaLxtIpNCFU{AmM#dQ&-1THsKC(9cQ>9HCq2=QRr)3LE-)#{y^D=1-XgZm5FH zsPY9#fu?F(&)6n z<1BFc5c4O&0=Ld5lPvHP&BAylTi{$DFn@Yk;Kq`c%<5}_pJW!s^Fj;UjK!IqX@Q?& zDX&@Jr&{16EO2w}joEn?csEP=2^RQi7Wgy^JlO)DV}YM;fzP+VdsyI$E$|cz{7wtJ zrv<*;0zcCNf5-wq%L0GW0zcaVUul7J&h)>Z|8d}d9QYpx{>Op;ao~R(_#X%U$ASNG z;QtRC_}O{d5v^oToaT>R*J!h8-ZD>2;|8syGH$KxUSs+sY|S+GoPyuZ?wK~0pT?Zv z&c?>Zsmo=W6UyMGaGDdzU`06131sk%aGDdx;M3tWCyc@S!f8$vgG<6`PQHSL;WQ_P z!AapXCxpSP!f8$fgM-6qPW*%Y!f8$jgFV7&P6UIU!f8$bgE8SWCw{?$$3y+)gfF-= zoaO>Ta8o$ViC?fHoaTft_(nL*iC*yOaGDdm;C=#aRVi)WYPIE#R>=aIOA{UGar@0^zJosCv|H(}645vAf3vLRh zIe`mSgwl;YFZ<1Ao9gVIZIk-b*z>ewk?Af`=?+op)=}x>MeWKE5Mx{@SN_UA$cL=4WgK_x+f%uH~`}T>cu-y-7rQ#cCp{A$Ls$%(Sm#-$AaAz>hJoCGzB zQPG4tyHG`oHX-P{YeEEl zj#m=&kpwl0wM<8_;Z_oJPKcPN*|eU12;+d{TtG}{Ld+l#b4vpX=;?L8Dq=nY9ip{C zjOOzSyn;x`Pb^1$xBq2Ag4dG4+W;WVbzg~p;#D!>Y__7OPYAawjyA`H7}ai((e6mL zd$)lwLMmIvhM4dvUdeFtEUi;Oq}_!L18vw1yz$?6az#zYh3^u}Ne z12bmUHyfRT{u>%N6je!MV;-?> zbNlMVFLYiOPm3B*hm3!FR9_Ah<#zI4sTs0(LF3=MAW8swHh(XuuxkfE<`bf524R z%h-~hzSwL#AhafJUxT(MskT2qsEB_v0x5l2))>j9T;Cos0I%d}9|>v{XEJSNNpXAf z^B&+3e@AQLS^R3+lk)&IdJ-i=@z1Lo@0Wr7<0{Yc)3c)7_6c=75@d72qF8s3eK=ZA54b&N)Vk$C? zj%a=_GqlDUr+2DdVW~X~1cTXGu<^)f$wqMyPxA0bJOf7(VrRjNu3_ZKLb2$++jmHo z2X93IDhG~?i#>UjEZno}JIA8gC-O@I^2 zNn00a*jEN>laLBDEb?aJ>GtU(@Gun3l^SsSgjQZ30|0yDEWHruwRd2N{8{AA5Y8oY zjo!V0-VN6sy)<7riYshhVsn-J=5UryW=;Pp3;|#i2ec@<*|@ib{9~G)8sp6os?Xd#I%%xXnLD@u!$U(ZT()TazWe zMlpn3UWmG3d9j&~uT-q7#FL$QpHFKao!;*;{iF;}2L@ISM}D87;kwm|tP&ETr?)rB z3W%&`$+}jN_3<8s^d(RvYEe>Q8Kkh>9>#L?U~`>IR9NmXuuLVEr%;#l!2m2lR&PLC z7y~wpP<^k#lYK9f=|=Gm(^ei7!>ZALHZm|o>U|_kPtOgLI-*%pZ-Awu+}1kP;rD<* z970LO%6%BDX%3};1fMd-V+h8GY-SRiXIQd=1zj|nw~!I2?*)&jv$NtEpp~A@FVlhg zKION>cuq?zGhCf#EG#$v_Q1jMH89_J`n$UC70=flkR5phu3sutl~_th1Oup0Pr*RY zkxsHBhs8lqAkM`vw{JHG<2Rr=q8rJ;IO1I+{c5+@Us2=83RpA7kPxJPcJtk z)hh9FQzUy;Y_%44g0lT3lEb2N2+7|Sk{i(z$)13i))dJj3dur+WM5EG@)uzJlMe0!yq=SeEU(Al3kgLSgCtNyjEAse$^cm=b&jz7bC-#SBXGPaE!2}?mSg^ zg|bjn%`Za>#sE{-o+Y5_YV?(H&OB&5EwJ$?boLswNUA!q4y?8vo{Vqqh9Me?)hK!b z==L3F{@}g>nwd+&VY`y&GE#A zkO8=B#pwpF;ENVqFQFDPiR<(*uBj%jMCnUw-TKRmIerdTEsM)!twwPPo+Pj0-4M!s zZr{rrXyI3}tUtCl{f(ac5`!#bVuR6N^Q~v7^_mP-yk&S)=d8qOhfKpiXug~$H68zb zRdA0hiLY=d@vY#@$nE=ybEPslcE@f+~uA}L^Gt5eV2}S>RkSsmD``-;ht3`G+p~aj&jEAkE+-Ee3|geWD|gyMb4c__V_EbQsHoCRj+~EecEB_eu~Sk-FCWWj9EKco5!l z(7T5SV3k8%iXWHaX-?M9SMO5w&yw%f!0Q?+%5noC| z95sMf5U5=HVX<4*#`L&WyoX>=8NgFQno+;!*@E;&`Gi=_$m#7yzWNi3VM8$V-8)Fd zK2a4^Dy8&37zZ^~x=1q!wnb7kafxc;^w*6hs>HF{s3wZ90{Sir{Zfhku(&>iUabpb zU22awm_sC(|Kmwc>}1+XZ1HIJZy4%|sc4r3fs;ny`3u^uTcu#8xo#)yC7175aHz3l zrutL&P`ER7_luK(Om$~AiPzIhUo*N>EzXVZj???B^4~+?FHrWa?TXs_4T=;wlyDE( zk6Lkx(XUQT`qP=LzBQsEW@F3Hq#|egZb{BbCOHA|1R)VTxfa}gXO|=;Rg$t#v<>y= zuL{B5a8P<(GyKNsWbRHUG6UK~mqqqZYDjt0EgXF<)uY^SM2q6(Xpcj8@r{qCT zNj}#9Nq#2Dzt|ujq!_BQ4hVwRAmm_^Dljxg(8y@;DcR5mj4MoMD`UM#KMU)!Aw04% zJaVFjxB(5ieFx<%W1skVisENyxHQ$Hc0{YEfAlw_M*&gueRPkmM~`;GSWs@{cM8Gz zMn6o(!oLnife1$Qrv029DrgxFO8H%&<~vHo^r2qe8qtHqcEc~lc{vA-WWbY>uHR=s z@!!SCC?-5WWVpIi_$KXfMi6`Brr|+aN4F_;v!a z;QfXtiVS%U1FI48m`Z7GI6anyxxi4=68%tLb=ATbuV4hKzno0*jQ|?)B8r8@_gjG? zXEe~LJ`H0J5Vt)|eoC*}Ho8y6!_c2IRDbUOMj`M3t#}nBTjX{>iN?z9;58T<#Vx;A z>9?Kar?rxrk3{NdEB3AsMa3sTW{9aUB&PX9qNne8D%7O`60OBfo!K0< z>ek~S;;c;NUKW^|$!3r4z^|)iv!+cvbiJx<%QTsK2!(+4o-hS(5fY)?Wa_6d29&$h zXuUhoifce?i%k7iqOtNvF}3IKRa)7kQgcjwg@JlE6cyy3!PNGl1%We&L{C5VWQeId zwnud-#MHx3dX$~^wc?8wBF@UxOVA#eI)O}m4yc4(Hfu8V{OPLli3pyN$2-dSpvG~;@Q}2{$nlg0*;4M}f*|gGRt}$=jhF=5F99`2-gGfWI z*p26hNzgbk4%HX~egBdX0~tp5i!8ZMypQJ^mgHD^sZscC7LJj?TCpA_h~eol@UU0~ zAe?$8gU^fQLc6C8e2ON6(XY+A)pFKovn{?3){Ie?XX z!1q>8?~iB*srg7boZfGk^5-n(!W5RxY}J{aZMG*foy*#y$4+l6{02Dl9OtrnlyrIz zl1#OzK_2-BN|Ji+!rX67YssU1imWg>h|8aqbR3^eWT6GHMsXWzgKyiRgqpwo73h}R zcVrD(ASX0mTye*cYZveay{QanvaM3Y?P zd>>YeCeM{k){6FMQfMJ8Pbe(I!dS*Ni)F6DQu_r7uNBomAY$-K5m{>xxm?zjb}zCo ziso4^dG>Os=`&Q*d!Cg%0|rsPHL3$nFQPh*{f(PtCq|%oA|d^v1ohdx3fUP*Zp>g| z)qFs{TLp|@x>octXg1>fjbe}_S(mEy4m62B-L;IG!hujNLf!f~UJxPWwI2jtu6OQL zxuap$=E_|=dQ{l8j0(@tY8HSxmpSV8iMpn}c$yZe54!RxSo}a_EA85*h!CfM)Zjq0 z8>R8VlhHhj1|KW8qB>re5#nXV`S%sjBf#49@7HRmU6xy61eu%nez-Y{3#Elm~ zTslM}>_&O0q#??CP(FB*A-1}=#Wtjsu&o32sbGX~AWP)4%TVZ(fRh^R3p#PZFerkbkdOPqRo<<`qRL~ z#iO;N8zGT$4JwgbOXJkjU)spB+r%^M-!_z%8vKHx!A%ycx`9Y~Ic03YZ%~f&@xO%1 zOST(Ey+55jvKsX$8ueH-1q2+BI2NIGwJJ>W6?A+%BhO~Om zWiXEzxBW+s7k6rn<#oG7FN5r`Yncly!D7@G@kU-)=#P8_CHXqiAU!!sr!PSsm@*xf zZy3h0FPkv1S~0NGKP3Z8Hb{<$$!G|q(VAmVe1?XLU$D<7@H=GDon(ezX3$MBLt1ou zgyXl&E%pEq>q4QBC4bYO=nYH02a8VR;DZ1vf_0LoryqYHMDVbt1V5xKc_NGlWjB4M z2z~>76?;)yVH%<^-Tq_lsii#RAh%i0q@g|f1rQ!Gs=tEjG_>dV3>(^mNWRkw{8NDr{d?J(wn#%Ytx?a6 zGIK>zIIc%&L|KS(rqew9tL8&J+E$yx_D zKbh({WL+SBp`YxSX>%J5n|mBjCAvHRDV5;#6`+dlyi%M= zNQC7vwwOj+*V9*HETI;?PEPMZNl+-&b_3MjxvV;dYRh|j09OOidotq{nX%HykZRjc zW~`7I_o$3+?Er%joZ?*OgovEpw)j0gAx3f8=L{!i$MykYZM7W*Syin$s{ zabE%&K_Iib941fr_x68&ZX3rqRs@YDj4$eTkmwZZ&rT@HMY4_(KMX2; zin!~l6h&Wcklh*_rf3Br5twMOufkH%rDTI58}VLj;0;a(FBQ5w6}qf@LmfP`S#;AB zy4C9p`v)em2PLKZIVFTHL!tZ7j1DWs;m?|*uA@S?#6UNR=pHxGeFT9j>V9NpJ$+di z-O6UsReUTvd%l70M4}sRqFVrTVLod52h??ue{K~2o!$vB9u7ReuoFrRu@5g+3{k3e zAF`!RYr1|MDvFgsx#=QRThqOA>RK&h6DTM!>35|Yz62ROLQHRE1$Z#V^^U_3AbmGF#qOP)iAJNR;)%6sMp1rw@D7m;mpuXv9h33D>jp91==)Z5NQq?3r&1)#?_m$6StSx; z)jy~A1Z7seKVqlSAH(RzL1O`AC*DTs2x-teh>g85)=%WygGN3o7*oxjfT+k^0U$wX z@{glV-5Nzz`6rU9^>>G;T0}@BRUc9+{=^k92~gGAplUbJ2aYU!mSTYq;x??|dMnp% z=@GD#dvnOWy@OeO1;>6MUW^qWFfMR3A$9@O)2+|g5La1#YhD37`K=ubY=NU4pM4Uz z#F?m8TIRW25t^h3y**6mmmfDF^l3%t9;h=2y{<~JZx9F-^HHB|Ca_Jl_u12s(c>W* zy-AdM`nYf-H#Kc!Ahq3}xKuS#A5e{K27Zfp{Bjdo^vhz|59HC*Gw~jrPf_gwQkD&i z>J~sEL}iR&2Qwk6H$eB}WQoDKU3sxrJt^Mg^u=> zwmnaXRj$Riq8H^gl2!eTUUV|J6a$O~ml{DfR_eS6bPP?_74GN4xQ9jKc6!@j*aZHG z3jd}`#icjUPq7!Jg}^WJAb1!YkuO8g!!QSf0dlIsC>7?^2ta}i!(^Rz@nR^QsJaww z_+OR_5|}iA&2?f`by#7qR^e zx`Np~g^d2O5_!;^b%s3(H2+!YH(HG&4bnrO!f(sWquW{Rk20u^c#oQe)HE{-IWPkU(v5ZCh1W;| zWlKCiA(v3OEV_!`$?c!H80)4qFD>Js)BAU}jRk1|%VX@BLr%p01KRIg#tEX+%ke|z z$T&7ZcI+9M^8|hazr>|?=3LX}(z6rE9vLNXhuEx!)TN9thha+iRW0g996~pM*XlCxEzpX{bdUL<@1W_-Bpt6RpAKj z8A7a7%Mf%83+Vvx9#OcwVO)95;wn_Qc73RNzn0V=M@iZHK~Oh^ zZkR%MRTy2{X3?FY(A{sK^AO#uCc2;Cw_N%H-Ii8{DYYCQGHN`XU{sE!fg^dbhZw-) z#P3cl@Vo8ng597`s`-xtk~4$uL;)qwWrRd6ak(1loT7C zU(!Qz=<4Ns4Em)U1H^~j$PYQ>SC%)-`qu5+5!f+ftr_%c5yOOvX%~#4dG_l07SOnT z4Ki$HY|AufsCN;)ifqj#H02IM@g}fGA(rFM8ANqC-F2Qe$tHQPVWFJRy@yiABT}DrW?>Nc9&`ck+48|%( zca`kvcAX5VLoSUD&U9-mV7RrDt|qda#qtmsxEni{)v*7D2bq6?Od3(D{B{s?@dOlL zIfK#Ym(RT+AArG`&eIoYh!BY z_U5{&*5GI9kyWUkatwjUk zI7%vm;?9gI+oS=cOGU=<^4&tC(+dWoC7eHOLgE#&VYL?Uc{<$v%ZkLqWs=0p!z8XE zB%FP3M>&X1^Aw^0lR5sI9L*tnkF z1?%VrIafC&r$Tyz#2Jd5Z{Jbmd;lb(9;L%t!1d)Wpdpywh`eu%yzmHd10Z!@tMV$I z;{MSW^ugP(4vlVbybGNuQCZlTfv-LnG>G+UfVL;vWkWBk8j=Z5sdt`76D#9R?aYhBdHGI){CIX()e~j~6fP4#%qwk6s+=@cQ?gbhsQXz%;hQCLrW1|DigZW^^s&B&AmPRWowXStQ?^m_|}W zMVQ7}AVQII8ynZtPY;Jv7dItmv~t>4Dsn1TDRQ0$67elcN1DdbzNTrE8+l>VxB?Kv zG!$cu*^6NwWvKdl=J7s6CtiFH&9+9{G>^Mg4atZl>RmDdS9p_iuQHm_2PJ`Uv8%I0 z&UvyC+Z<$UBh3SGPLz32oyC*GDd=9Y8n`1BTdry&WOxU;hR~_zEYlxl@QzW|zhWiY z(DaH7$YZ2CcEH!EQ< zX!>NM=>|!oGQ{ipkT?~hU8&w>h?b+?my@Ur^-fI@=ki@o|8xN#gFXXd@kl44uo#)s zyN&JMUO54Fl2&$Y)-_qF4rft+xL(q zZe6Y@P0{@7iWg{oOefwZ!1w*KrQn)0}YXBIg z?mbM3z$~^`hQ9FN_%gRqJZb{bjmQoP#H+n-OlaQ4BSJkKS(KT|WbOrb#l?3{aSgjl~%qZR)b_W!_$3;6h$(I~B)-iLtmQh~(6HAuts6SNE##!1LC^*sc zZO~F{upNz6=OZz5&%i3?-IY6UPKy0{5%cq~+uspafuP-;$29xKm90FQraN?r|42f& zw@?aIS_gKbk+{k+fhm2o4M~w-^oK7;!EsCHlGEcPkt{0i ze1qua2fK9o^^kF(d`2~5YPsp=)*T6MACBIX9Km%dtz^Yj4C^rP<5m=j;#Z(>s~bL0 z+ok6~?d9UQR#M;EdB;njL#o2wS>zdRS*g*=e}pRRget6`I<=3pJNzkrmPhGBMKQiC z7Fx;n`m93%C`2p}K_Q3&M=|w6D8%uYxhb07|8}VmI`BgwT0LHEzhXJgbX#iQ{ zDFEN)&}o=&`5_V?|#=lwF;gZ+MulFi$%nG?3(ug6ld zieZi73-FTC=?wj$?IyNY(rJ$;c^xa$jp8><1t^_CGc1BgLMDOnpu-SZt!TxXex0!g zOfw0OmQ7O{r1a-B8)cgD2vksC+EM7IflhwpEuS(OD4asYa=7ZK;Sg1>y<~gWD5yxJABU$oe*e)x0KWX6DScI0+KWp&M z(e~(te*T_4+JO`TZpI$rp_M(7vtQo;W=Ojn%R^xuOnZcf4%ws1p&img(+-t{)qo6l zBYV3C0xZQQhD6uRZ*= zpTBdO~cyLvX53_4&@EtH+g1fumK)< z=UHsOs{u?1Gp<&&kr4pT?1*DzED&~;2zww_c6z(wHxZQ)5fo5LyH3?YnMqaCHOUH<_9|da;xj{FS#G{B`i**MScoBFs6gp1qE$=Mq#SydGvAIF_ou zpOB$FpkgP8EG-Kh>!nP=8^=;8Jq*~=)&*A{l8=c;c>for7`zQ_rIiJfQ4WuX>*;x{ zv)Jpzd4HE#9gZSIBzm^C@d4E*X*ev9huiZ#(hT3-D7tfLwIL zFV1lh$dJ~FFQWS;_@!W#gdK`Z?lY3XkuOCQn-G~C8%pxJq~HZcGT7r~Rm)hYZwr4b z=h(yX61*0!)6(+aLrIBapi$PS`lW%%ApI>P8EhS@CA;(hdKjz@!6c?vB8l4x2}vR{ z-!fhX!wMB45{KV}e;Q(!WX(+>JTe($v=W?SAXG=UgLj80mBp4K8T{?-h!)a}7L3k3 zYjjmzHWi!}YF3gH2LrNDcOX<(qW>NO3I}*Bm$ujC+c;kH?ZleJ{t=_p<#YO=w0#(u z5{sU(*?tL^$12Z8Kr}seJs%nMo8EN4_6C@tS|o znW!dieh!K{1+-ShUWan_>rmHbtsk6+OoN`$BQ@V|@1V=Qd z`Sx9V9i}e9A^n5yu(-NcHQdFx@KVk9qxieb zsRk3Yi*Jqdtk8V7#?{>exvrZ8=FfPxv9b77xE(n2j&16mE&c%SIJyScQ};&k`cMX( zjVxjdo&{lq`A7L$SI0(j-3;8Ep!tU;yZ!T$-2U0UaBk3@GD7oZr{v{&58(zUIour; z_mCfxw~Ja3k+uozqbxPK#_ikhI{XWV%O%q>{%X!wYupHC@abA#SH_pllKHIezZD-R zw957Uhz|_Pe~fW)>f>90`uqtzVP&#Uaa{^jfep`a0RprjT zD;=!^zjbu9X7KCI*igKioW2m8z6J!bSx8|s`f2yvME3OsBl>B2*R;UVPUTnVUAI)0 z-<)u2WyOw%9GmglvDHeuv$Q*_mTW(cDXC_!Kt#wsI?7yn=Z)k=42X^i9^n+t-~#Si zLrV5E*aej_7fw4=cR!iXx9)wM$VuBIzm03!QL{hPce)igOZdHcP_=3;wJP}R)AUf; zHrv`TYeV_>~cLmnbV)mT;uT_=I03IcCCD0j9b6Ljun6; z=N+u%_I(N>X8izv4+-9u^s!dHI~J=KtxI+%d%9qSDGuw}$*znt?wrECWC8_Xn* zb<~m}=S*Z3|CH;iadEw&WJjCJb;lCVWH8sATCM5x&{Sj0+|$^`7XeYV-d6e4{EqQU zG7nX?^7NET4Z7nstWr6$eZAlfJRR{Km+hN|z;&Ysz-fru;W~Z85-hC1=vVh7nyA|a z7x%WGR^EU)#9^m*C|!>im-Ke$OGzUp5UG6geU^4e%^&t* zDX2D9%{LQL`USoj8{bzesY(v)ui7#izcwPC&MU|NwSoPKTPk@~n)F+(thRatU z*f2K!HOT*%+g=6FEAbi)G)3;o0^9n(=GHH>n_rw&?sTnoyy42IDe8&51s7rZRuxmZ zGi!DBg;0Y}i<}Y)*GEX)+T~W>%Kjq~x57Z`TE!)KZ|7*aO^k3PYm_uH;eG<%LQmg2 zk^BVX*NS$8_!G~66c9cAnFpD%O_csMTB)u=IV#mME4Nm)$-L z;9MZAgm@Q2ybLk8I%v$!2SB#{OCqM>ftsYKw$Fkh2jy(GW>IetC;lw5=P@sGNUd2oD5puVaGs-BP$C>sxYu&xlX65$Sa zF#beP-`2R>)mVPmQPdmD52}~eTzAIyA~?3}6XDDf;*2^tcRY^)ZD_Am?X*(*m`d_Q z2p)MD!062~9>;v+N5GI!ryiC*DwQ}dKcvJtg!mI@JtUQQ+x<$3KYT1&iFK5t5(khI z8owAL%f&eV_&BZiX$ZhwKCC;0ha=B)YnZ&*xda>N(n=PxSzkHFnsu7}qria~<(lt0 zM!Y#Tm%o3C3j;44+Eg}c zT#X_WqerafDF0lCjQM52KNQw?mgWn%Q!&W>7URmO^jsX+F=lJp*1)c_ zKb8I?I0k_!!dH9kL+&kA8qsvhoYfUKa??sBrnlCxC!OtzG_N zu=6@_@@GU)j)k{iiwzSY@ha=++1+fp{&P}XzRThyzul<eJsHM_qy|4v0mB_!C#I00ieB z?qkM2(SRCq*z~QfYiquofE*gf`>2A&D8Al=6&kwkAj*g6Rpaf@IOv4tKM&%mjJzwbOe82Xo9DhhS`vci!P*QqOX5|0T5wcgrF`n|Fa>51yi?>zIl76z3CkEq0wOJ(Fj>XE3t=-^SNhg^< zDu#X(`FYuPymgF;WB|^u1lK87xi+macsc4c2?1OWNWHs7o^z8SKpR5*iGBY!n8z2{Qo81nfq$$Q*?Q{FLwU&ev?|3~u1IN7E0scGXJZ81_}>+p}| z%_jdq-s65p-pKot@*aD?ly^lQ#VX}}5+VM?%zLH0U)Pnq?_Cis?|hV_ydOl)e^=gP z|C{oLeUJT5*|)hh^~dra`^WMg`#bVR-k+5B^!uc|JBxkvy$$e#*uPJ|CW%~0&^ykMT{P|ay+tHslHu2}<|8;pE zviS3Aw1D}7^5+~XFmvamJ<^}UzQ;${_kQ|wY(?#cKd-j>a|}!Su>}Nwj)^+`xm-)Q zTp!pD_I(oUd%U#oHnbZqvDuGX{JAWPc^sW}o4^q^i6ur&UkHB=f89^|>jm`ZRb~+- z?>5rkcZ&4q**VPO5U&E8!hieC`Wm}jU0=- z<@$`8FKQQG;4I-cB5r?-<{ykP{fanuD)`Q!yxi%%Mh@lcJlA5??|@r(q1iE5t#C-* z$cs_BnyFkp&p~dvd`EIqzjl=zNt)LWt6FJgmhfZha>RQHj>TZ49P4?&XD=2%_G2B0 z%jTJw05#xz$Uiaax_s->%3L@gX;k~jPL-=ASL2ItlbYAB2HMp=5i={|R?46~>kQXw&C=9L+ zY+ZRfD5Z%y-h~R`fY4{ok}FY$@a#AT#(oe0!80ExXiCyUdAP!)itsl=JwxlkMCt>U zILd-Wcmt&FUNzf!_--lRcSn+P>A%Mi;!j*mL)Oy|EM>+vvGE==vbESln5J_aema$4>xpytKn|kd{go$6~T|o zat330f4}Cx0B%BQ2^_?0QcDhbMo~M!bGy@lwlsq}*@gD&5JX7oeFdsg+>%digWp+Y?zH3u^fN9f1ffPpD+37u~r zuE)z|yw1zuN?LiAD7R2=DuzoRh-@eZH|9;nIP5GLL*WCS-PYl@!=OLenGH&tWmm7yA2RSk!#QEb#0jT+5s(y72e zO&W<3=J{j46W@(q+`$ur^>zOk%u%n7!xEr7V^i@@x!!f2Iaya@r#ovD)^=Bb^|sMz zdq?70Anztm7A6!4n*ZhmtkOQ@wr|5)Xcvq70DieXS>Z{_@%LYxozZ`(=O>)#=7f&Z zIp~Mh7|6Bna@#+5Wdw^~l0!?~Hu+HbDxPh+{*7L@xc=?_gZ@m&!k*>U(df<>?9Rxv zy{b3oTY9ru^(M}wA5(BgHCB%MILh#Mb%LrBYq8&ArT<+;|MyY!|4Zu`_0S6ZEW(5} zyq@vOGAhr%w4M=XUC%g^N`aG}|LS_iQLJZh7YHj4_PU7qr#R}fdK|f)@dx~Mmg1Zh z3}a#9s~7@X;Pud^ypEUiMXqO8c(8?+5D$JAjf*b}92kW^ZgO7D++8Nu!yYh_|IaD@ ze;v(#<0QSg^NK-YHX^7#ZvRy1xxnUFW#z7$oaKOHjR`Kjd>?i*zGEk_hT_5vA=N{Y zkwTvqVZRnvs+K;Fl=Z61Sa1oaku z;yo_#43}xTGDI<6xGlX|nwb1%Vl9YQdpJfX($?i*CVC@u4I38eIoPn+yp+bI_ETtY z;FDY@&X03vbT3|%i`A1N{8_)t^{$U|(F@b-$H|SXcbIZDmb-ihHT}l=+>9+nxK#Rz zgv<3&vACaYJ$$9zZQp329_z<)vW$z|Krs5y4*zkI$bOKK#V+EqZ+6)a0#)HExgdGn zTIDKJKqEZj3U|oy?YsjJf8q(C8yO20GXsBL7W|ezgU9vc(0Mj%V%cUAa@61F=2sP? zxv=%!NNz=Kf&JQAn^C*4xd4DPTKcsf6`LtLNpgY1o&R{n&Z$4yn*YD!$RU1 zNeGr#!}UBJLT!&h4tA&C9GFv-=&mFn^Z5x((OL+5{iV)UCGbfuWa}~ z*{_7W=>5UAQ| zfz5hMG5qJbe4pi};*U`s(aMitu5qxiCwxy)2bZ3`%7%E_gXO}Z=yb$xJ;#6lA{<4* zLZb1(!C~i;AlMCd#mK1zb~*;8X&DX9(zggM`3ZOqVgm(R!QnBn)X0P9 zYL1DJqXRdFahfzM854}QSct@x7uNk?KoG%1><8cOEbRk)WV7!}cWRZUU(IG>=Hgf( zwp-+m@a=~xTH%BQKSt1TCpI>O8Qp0aUwKaCjxb#o_5#21;G|;5OgIIc42kTwMIzd2 zIjr-RPuZ{K442q3~qwU6ta!${RH2&B$)1y;XXIfFPD*%AvyhMdD}0ZMZ#$CVurqm%0_XWZgR;u~bx{AVvvfBa zB|lTK({Mo3=i24mp&NDw%hA7qtI1QxrE-_Ox4Fyi#_9qQ!+WYeRGOjH?{}Iv5G!u-eoY(NsG<# zXtdnyf&t2H+L0I~>=%U2Dx)DRE{#E}*sySw3#=?`kF_YLoPS}z0P9O!yvO=6w7rZ6 zM7X{r&n2t%rFfu_8j3{*-}TG?14;hxk*>%Gj;wXh(E)I%S4{!ITm-p|~F+OlRp`s6%1mBd?&# zT3_^DC>Jy2rNBVWZ<17mRdE}sOaD2GTEWlP#9~7Hi7(tjKb3yUJX+#DaRdlBrgM=( zZh!IR1mu9j;saJuaaZA1HPrlw^5GuI`0_uf?{~=t|6P4=4VwO=`p$DV|JU>#>*wRt z`uU&Gcctm`6Ot;QnlR4M?qAmTasP_GkNZ7+N9O+neZTnD7W947EmGevzBHuo%L(x( zu7;vR-_Q0aeQ&qmcl7-eR#E!CWV6!uBWUoyr|*x+2LD}s?*f|sqxwGf|E9jfzjO2E zKj7cT{_FmI?7yP#V}DQIk@^2X--pd_LEm@cR6Ep~f3L>y?@tnfzkd!zN5)0?I+_0c zl$(D?-@j!QrSJEBq4d2S%K!KD{UzDpzpL+Of(dGDM1YK7Mh=m;A&}wG2kw)(GLETu z#x2hzM8q@eT&d-ZXBsT=%=ka6@44QxEV&%zcIFtTAFf6;;+CsS^3-HRK!aFYo_`o` z@$WdJQ7qip;**Eu%CKK%#50JyT;^t3Z!@03W*N4?6S~8_k9T0(tYe$3LuH6&c=SM> zi!5sOTI>TM z@@N{*9P~(i->3)QMWe5ZNhR#oK5VVIE$MP<{(T z`5iMpGD7(icmONcUl5Mvd&>GavH6*SB86-*QP2kD86(&M!b7Q#IFM6`Icg9ZnO~2jPqE$KWjNS5!QEtB* z-}K`kx;dG^R3a=4tT$jj!K=?h@iJ7S?lVgfPY=e;bO>CDT}7@0gQ$>pC^!HrFXFM; z2fpD8Rvv&Ff&-FyWuVVqw+tw=ZSp)-OY=2I+?cO@8986$6eddGynBfYF&`IF!K8na zLaxUjA^SZ+yNOCo+Oj|@K5vmnbh!|q{G9X71_-6H{K(l=?3|P}x*8$w3nh}~#JxPQ zC*phID0#2>LwS$=LwTb_3-VSuR(TuWfVVtBS`K-aS>-J+JKC>p{HbKSU29uaxdVKF z-joIRp%mBOpbxKF^Z^WXgrz3-t;^z00Z@*BWn7v}f4p$x|pwT(Z@ z#wyrY1skh-(^R1>%9(!!nyQ3PgM++vJ*=?CdZt`7{tx<(d7v$1{9Fq7_x;EBRuun7 z{^R%b>5u)#e@UPIYyRWoS<(JuY*YVXtha`S3;FfTbZCm&Z#rfU71E3F12&#H$Hz%( zeGlIYVx!|PXh_qGj&V667cK!?FaxOs{UMu1b@K7QAM z^*=14qGT%!0*=>Fe7PqX!-NB-fOWs{dMelb)Osq`{XjSuR72kyh0TTQ$<=Mp93JY! zBz!9Od&M)Qy}xX+cbM-;_<1bCU34}U;l7<=HAmo=>l01ogWUoSu`S7`iX`b>nNX}ff@3t=O_>H7OUeo!o6YS!clR2Es#;ob=9l{BwXNjHP%(}mhig_`!oEFMhM=1*9*M89ea#& zzZfS4oTWuz-)J5@j@W4UMZ@22QX37KZh1y|DAsZ7<=4I$8=Omc34E^aICg{Wh^Snt zxF{8$1LITRj0$H-Us^@^E_{@E+?DzrPNK&-?>;9MFsMXddFvDJ;6~E~C&& z>ah0zwfxwFVmz8KHxFNEz`{~aDB6Bmz;N8362=riEcqMNWs*(;|DTT~}5kLV9Tj2Aax=5h>|XSj4nwX}cgCwijnpZW>*OZ&$g^m9Je{*OZ# zi_sA-`9owBp&V^@gBxP#4Ns?X1HazWKExg#e-`2Fy43XxsR1QB zndiqhu5t&2d`os?FQhLmL=;RT$C%=#hznn!z!Z|PlMhcj1pKLhg15zP&}kt))hEmy ztNFg7I%CHOMD2q**Wz=`t`~5}=p;?QEE&vc?aEl|=^A#?0EPV%@(g`WCvkE^t^b$O zM_oo8=-|Rp5BR9dr~~+`Vo1?P4d?f&L|iS~{XexnmYM`g^q9CFxE2#)?1fne%#e7$*;`%w6fJcTNF zMSXAp03rP=#tQJQ2_2Z_gdzsTkq88+|1@?dRegxY8T<>9Rua9>gF5CrQjp&&_)t_A zyyJ}g)@&ok7IY)cmVN;l7%Fb;jiKUd?5cqzbADwVBf{Iz_!!JLFeT!gMiS3y;1q0c z=q3jReqKcR`}?So`WK#skf&d)6Ef`26`#ek3Q3`J5l(LjHPnwY8IhuOVYaKg5lQig zhRfF|dN+ViQ)M>5c|mu5w_^Iyt^4q+ngV>jJ^VuLD2lZRKaP z8|{Hzd0PSp64&90zsXr%{3+FVIC~R^9pE;w<8!cM+&KWmw#p?NSI3vD9jeVmsV+PA z#L8kSZ@)F6e|hE8>sGrv!S8M^N~nD5)`SI}5;S{tWkvbw_y$TU{%5VR%9Zh})9Zr9 zQ4QBaaLZlQj-QDdnYeHb0~LUgu;OomPXefV6ojc7`%9B1Q{L&dsi@3jE^9?HArWH0 zk&&$csv{$7>$+pLruPJt$KSt_UuA6*3vh9Vx%-r&`3I(eGV#?{(TIxQW(+nPh4EyO zlqPKz-_QBC8%4qfXPgg``>^s96b?6@5}|Q%#>t@Gt(C@Yc6zUpnsWpGjvfz-9b^1g z?MNY$tV)N|U&SC8VzMaR5LU%i92lz16I*M!?~D7Bsc&M+hTb1>6{P zmFzl>1E>zEWBd@x<ep1WEX%ZyPc8l-F6&q0Mx`_J$p)34afchg#|EToG0)zg10ifEH zJa;0H>FJ-G6;hhI2};&sHbQBhLout;9J)>_&Gp9B;X@~x*@3z`fgDURZI zs52L%tf4den>1Ctf@hMFP{j1dvQ@^ne}gD-kS26(FHX_$vsLN77~^s+kl)~$@h#5K zU_hUV5Qp(?!z_C!xJ3?vTe-hpzAGja++xRmpA2px4H?%qEaYj3i_vvd$1At$Eg>xe zLtglJ_amregtW2-L)zBiO8P*Avp1E3S(nEAx~u%6te!sXjkxV$bZT*m1kc>+t*Z!s1RHj%|R@kSL0 z_tzm~=072S`Q;C8j6?U|g8Ugur(*enttsIj2@Ahg3I8o6{L@(05Wb7TZ?6DK_!4py zb|~2aoeF+~M5O$q;_dD#-i|TjZT|ULOT7Il{qu*2c=!{Fx052`?IgzA@Zgx%L4U;2 zaSS$*kVC((%!skU)C!nsTiTWzqE^SUhNC&R46}*BGdd|@7-x+*GMDS%uy|^y(rBpAT$7+CTy#(YkCdct+c%k?d@vhNfe1v9(;%n8} ztr4&uCb!&ST-6~RTy-x`18fSnnC$1r?O3Z8?3C#t*|Z=^@EyB?>8#AQ7IK8}lM&hsq8Hyv+z?%8pK!v5SQV;~pc$=r9@o&2$8m}1} zrk~?m)UYPEIdB&}|3lcb1-rgc`%A108P*s47k2j1ECDZV&KP zy7iY*B6?keA=#9}rCM(ztEB#9=N;>CP42s=(0S|W z--CGWV84^sY4-4w)ubrmqEk`?U{jHP@Y}i)?Z4}?4gjDoeQGsp{ zbc1~2JUO14U-jA&(w)y1{E42?lB6H-p$@~$018WaX#UpPyP>~bk&)o^eg&c2K0l6s z*}~~%jEyfRu4kZt>CRDkT|{&!lWeH}UBLrY@!OogXex4RKgc48#ofnUuw# z-02CyA@+z-%4_6QIjr7?o$#-Ci&(O7$vz(8w+rWB@c>7^s>uH3 zsN?p%l=38=A~1%X9fE_vXRyWj#LP@0t}adKj0b--y!SQUK~QN*G7>zH=EbLvcsb;m z)C|jPl3sbNy_B{; zPJ`N!I!(Vzl8;MgxJ`qJ2^@jwiNM>fzn_wc66`BWCE;C`LNtZundFb`C(Dth5e86v z|0`9iG-Uz+CW5YbH=99OCPN5VjtJ%$K$4Llo&Xlqz1OpiH-%oiG5NNVI5Q!u%c0r{DDf4k(0`hT{EUutKjKJIO7(Vj01=Po) z9S$M#sytk=rupBJdMj_l8E#y^a>8N#1dVHcc6`cq#B)_v$$}L&yJs&xs`;1~#w|u2 zTaG{5Ys>ML6(@sD4qVJUYHi@>iVZRN;Hcc~9C8rS7QZs2 z5qERsX8i2*u4+woD4vroc_f5qv?pZy`Y+A)b-jaw*9chbTDN^Wnu1cO>xAw+ffDMN z%-ZZlVC(oY7Sl=;&oMHUro=%+&O5$t%SPUkLPYFT$+bDejh$Ot#(=JJypNup;vW-kMfzi!Hubm zeKWq8;+?>D9Lgw7S?eS3-SO2rSnGgvii=sq9ctEP*7INVlrA;7#EF=&CTNK|eo81lBcC#CH zHzAOOkEr68smy9+YOWMSD$%F#LMXbe)lrhFM$o~Q<5Qc|{Txrn7WFtL1@ z3`aIRMWX-;ka*WYt+t+|FF8V%&lOl~TUo9`*DKA3#7lG_gx5me8j#!i{V$z=nYmyc zby!4_lOzLE^%>>xv}xLR%aKF6>)v)eL(k{Rnm;oIZm>B!J#BM#eQM<15cWH#ya-b4 zbpF~D%F%at7xxBXv3cVSL%QDV^Xdc=pBPwAfECLIc`l zKXoqr@CUGP(k4APe?&zG7m&S#6T=o-MXCgbF3de!m*Z=C5r9Z_(HJ-fcD6#!z+WOL z!~IO$AN9e3>HJX@W-=E}AkARx?W8Wa$x*4N?#_EP>aTRZX4i{$9{$IJU1@Eys};V3 zFYe!<^EFkR$(+9ktwGJdSm?D{^b3kMne#8Chi8;QZ9x>&Hj`RHokf?P)H37@w%A&S z;T6=*9edVwdYwo%lNvaIGlQbhK#E4lCjE#P>n~jmt^&QOL%Ou`X0lRb8LO7 zU;RPb)Deg@%6=4TA_HP&cu&t+WGh$rwZ3%)A{M0hhO=9 z%4_r%4zJfIfFY(1F*f)TQ+%;$B_?|nxf8hwp5AEGsmS~jCVzj$2OqnDginXz)A|lW z%{PR_y3;p?fOxHgQ*&2{{U}(2T0!3jO~8&MjBz+}6AzuQYk#K5B>b6PLE{;7NBu`? zJ{&@Tr*9m#!XGK-XD3$$>*U`@6;hBM&sUtG1Y%61SmGsJa-sig?dh|>{-dK;k_J&H zc~yi>Iowi(Oo#LKU_FoynfS(v6G^`5auqGzrgX&qf7szv{Sji`f^sqT^;Hv5*#+>$ z#94U%-9n-qgP{`(LV`lIpFpaWrZ-{BEFd3$=f5J2etk~NdF;i^p{8E=KqtLf@~+o zV5~*L_*nIZ^I^Xk+8>^ZRu=Y@YUe}cyr|lzFnw;LHbUd`F-;Iv`+1^<=TlHP+9$rJ zicjZgPcX*Q6L`Mgc{k4Y&E9$QJZs-7$P>F}! z@+Ajscp05JKOeYqI6D=B-eW@NI}T`;J!y^=%Tsf#ZV_oNr*Kos=)~X75P{_#gq-bV zChVujbgR?mX-8T3IVgsM&~E6Gf=a>eeC1)^GUxlUDI`ct9JqpXfRX1;qJ# z&xxeimh^;3e%#s)?^DtW5iBf!-Tt;Ug;Ok_( zz5G|qOB^gh<0y~;)%Xg7_8gA0;zKI~?N3Wtnet@rvMn_8#q}Gp@ldX}T6Uin+Ky$$SQ9M z(YsX_UjAWbjunk>GBX#6m`bOZ2WU!-!_gPeGAa?^efJkHq+}!gLF?!f)L8Tp+{lkZ zIPQVoli<%pwf<@x;Zf0u_hl|PO#IC|@c$&c!oEykcTeUx^pa*N>xCP_2?wXl;hu$U6~8axI@QH z-I_+0#u*&qGf|O`_CEfXm;pOp^ZJ3U(;*n8lOOD-GFomBKC$!Z#9D^Bwf+VK(w}x6 ze^V!R(2`auj3Up9#K{k^MB4MvUG1i3m&Y08T^Z7%DGvIcF+(t0uDlBYpJ460};dKONx&ub2#^Ez+9gX&&XRbRN zA3ey+Q5`~NZGFfQ^7wrN4Tlk|@z;A@hRai3jhyO$zt*T?%whC5s*G&6*Ihf?6Lj~} zWL0is`mllchZSzR9LC^w=?VtopQ3X?SyfV~Cg7iK_zhrgLKJ1ecr8PrZWT29i!U7t_|ry1f4Jgpqn_ zIp0;D%hTX-wQZy}rqtyPZl^pcr@hh`3%Wyga+gGJNGyoxoPO_j(p5M@0Z(ILVQ{vi zjsypdK|eL90FTe*Zo~t9Y(4~vd=ea&Z;TpcR69Igw~KrLhVX?vP6rW%(C1Ks(;bu! zE#(=#W9emA15f+(?Z~h8`W^8lDEXjQ79zDsDR(1zpVMDk=Loos*`8310hbQd2dYHR zSKUq3*YLR=0q-2UdsaPap=u_0LVhngg3swz463|1`q+|Z*SLK~gb69Xf&T5qAB~XL zIr!@w!Jt&s>-PuUQ7Oc!_^KR1kJAXbgQ1{7^k{&?>lWbeq3M2rPUj2O5qq^K5DbZ? zogMJ|W~%br(god4zt5FO*p@C-;|a7OOHyv4r~3T9K|$CRtxzkROK5^5f>$A$sSaTq z)riC{m!~!+D1r#2-%V%pHMj$)^%Z(;kJj_9NMn;4b7&3AwIHMbscxv6g zAclpY;qtqK2)6i|#)I9txQISjs2QSB8fT z1tZZ(wyU>Lb9|5E0!sv#`)nSLgP*n6JLpjS!FSkZ3ndhgA2Mk+dp`pw1+bP9O(h z(xfP9Mx|y-z>kV%)CK&mdS`6zplD_s*PwH9^-q?o=5Uz4CCiQKV+@VeQ*q@(|42>a zG=cv4njx9~k<2)*{zidjBD^Z408`MoQX)E|MiGWcDUsa4Q5E7!K}{pYN4Qi?gY~s) z*kcq7BY|I4cklt+`dV)D#06Xvmh+^9GDz zVvfOr{ga}dUX0p<18xyGyq;>ea}H3Nso3i^WPDFD3IuB$brH&XUt|=b`94idYW(E+AMC0qv_ksSCD~lCWjBHuENBbPe#00Dus!oD(;UVAiKC9O;~N1+)%oIvAmla znOe1-f2k7NW~f5iW{A;^F?x@uk^j zRt1@cGMj09SB3APprN3~XB8!uD+OqN&8OvOfxn4#Avu&*KPx@MBVDBK*B3pr^{ux* zxbfXlcf5IdYkE$XpML$osaGqiDu32;;_>zWd|~rXZ|M5)xI=U2t9tbXA6Qo7jG-E!qR6Jza0*9x z_n=(BX&a#r;F|kDPk6lFe;u#}uodtuV2@PjzX`7<0UI9(hnoO5{ZBZ&4{-4x!{M(1 z%N`1cuflfTH4lfwb%0ra4u>}a(%++f5Af`h;c!>H5^cr1%yPg3Tf^b`fV-a!hc^Re z?FffI1Welr`{DbNO}oP3BEYjRgu{)1*)N8}y9xg?>e9( zm6jDLEqr-ua|g>urjAf=bh`Ot0EhtDcKoJ7$5KF6OgzOIS5= zW$6<$vXFvG~*XerJxnp4JQ}+~OD<*^4gI+{$uc%?@k#3}mdkwkg?->^C1gv%s6O5|z;P9Nl86T%OX* zr&ARtq+6J%D@vDTZ36v19EOiq1>!RBe;s}!R^spZlB~&)WvMff14t}NH(La1E4=jWsXOgmfPd(M>KTy<3_{{K{GTt4aKs$yA+L) zoTdkKyK1eXNpDL-zD)fo?f2pE1%N4EAW&%E5UpumPM%SkZhnIxRbm8qLZCf^cr}W6 z-izUxKptRTM6_em%`Xugo$l4Ff|qu4pGAzLvqsrWvUh1l)^dw2!&sOyCL_CjXn!F6%YkbIZkdv|++0Ha3i6s$8Y1Q&qvS|GX#&sd z;Q2m&MP+OTeF*XRk5opG+BjN9OdD6|Ldrn)A$U$w3`aYWU#%p+T0wr*k}^6YdpY^l z!qno75zQUOWfYk|=$J7AE?S(C4T^}DsebSwZ1(zwaQJse>ugQg(+jn&f_yynarO^w~p9*bUM{{D)>_{-b^6= zLR}u>VnVzyze-IzQSi{1WMC}13A#%C#&%A%MdPTsjA%>J&1VRfr&C;_anFl!>hI)7 zhj1h6H4c$vg$3J{IQrcYq8^`aIbTGkiRoe-I0MqY4=KLqzdzqV==k2B6yLqypU(xp z)fg{@o{P@OXBB!@rn79+uN7UR9OLUT@&R2pTZ@H4uEi@TCs7{u|LvN>}8Ql}fA}<5hGB)lQ`DqJ}YEShfEryuNGtUc#P50QmvUZMq{9^qk+1H2I?7;XrQiS zWamckAIuJiUn~XxaM-z6auI=$A{ZpA5JKWCXG&;@uRtYfETJnLHlCOnj$_J!_JRby;Pi zL!J+q`Drs-ZbAPPReRI)h?%{hs&>(&`w%_Wa{dHok8u7hXNR~{{y&@3Wxa_nNMXk< zr0*Y6ussI4bZIg1r!8!~g}-fKuUn=Hyt;?up@W($9&Cp$~V zZb{+4@5$~-;VXJ#0Ok+&WG|)g1wGlFWU)VT{>7ziJLex?%3k99n87~b{L5bKu~cO4 zN#$D%_I@gV0O(YHvL{>Gfxp(19qE9~vmJQ2Cu>RLJA1PAX?#U5wl|Icz88Bm9htk+ z`K<=~TRK0~lRe*wziY5vnLxjj$-gw%u}r?WH*2x-L%rCM^GVeC{Jvi7{R{YOm$GNF z9zY*Da|z$mi-j-AIRJs(`Fg|hNe_~;^HQS0Bz>Ow27s*Ic+SBuymozO)5^|6A?s?s zz{DE3i5)iK?px7P8_fJ6GZ~4oWfp$OLe>GIIfcKH#x|w!Kc}#>^r*O+^G7q;R+uD} zz0LVssqB_iO23-wXY6nX{#+(omd2m%$nH<$i&ELXG=4UXZ8(pw?Z|eYNBHF(A*Hn= z-HKf$mY34`M;VqEIw9|sPW)jjyQ4Fbr#ti8y0DF878vB2E^MWh^4_!p zamGq{Pjo@@^Dg}F=d*R^^WR>;ww%vj?ZXb8PYU050WzPxfFJ6{S})*dyRqB5@^DwS zr7Pdjm3`F}ng4Skf22Em`9h8-_>&j%KU~D_xQK7Nh<$YtGS_zFd%Lmi-S~+M*oRsC z{T^(|#r&Hd?5>OXl1tf(7xN>RTHd&rKX<9+%ZvHyOD*?YLYQrr@TdE*k1rvb=I;E# zKJ3BnB>q@;ey|T)*@M4$9^2f5FS?Aq-Ghp{uP2Inq$j^~5PPlX;H87uAA0e-2C^@E z0lmlQwr2oa)0;m%fNk&1Kfc`ZYHz;va?7{9`Krq;_xB;p{g?53`?7tPQJq>ZCXo6cLuW614-Zg zG~d}ah_B0I9}gn2>vMp3F^B(U2s@oa0#^*?uMJ^O4(6wau#LIMd?lAJ$YWpR5_LF_ z?=P@y$mbghEW7h5eJUR$rwjOsp_Uax_^rb%PY$8m%Z^C%?m9L4V}W`7;Un?|#5M(skFdg*8UbP4<7XMAf3TV7m% z=v~b38qNMyOeA|o^FNGbpN{66#oL2ud4VR z)og>4(mP%JfSVn4QF^_b@2q8e+>~B66Yt{wF_XVo%huNL`98L@hOhClt-nI4pZ|(4 zo5fCek^G03p9=6>Yx#SP?CDy5vXT9x7U;kF_!mC*iH}%r_wyZow#85B!+w5i9a~gK z=|gpV$1JvR7Rh>a7MM=V;_CwJwg4fX3h>qddtam<3}QOACB%=#~gNN1Keiy zY<}WKwsSUr>_+yl+5GQKCbqGW@1DnAZG^1b=kV|xwq*|gYA*X~4iHbw&_2fX_h6AuHw7m)Isn`wO5jQ{T;ET0bC3C3gAk>TNjylGvNDx&jK!547p3e4cG#B%MyA2$`aE<=u&c+3Spmb zBYr<#KU&~N3;bw-A1&~s1^%B~;6_=(c6vzqH!d-FfsdXs=+gIn67FeU9GgEf;eMRD z?Gi0K(y8rq;nYdWpTYJcox~qG<8WEQsKC@JQkA(bxd>h6S8isScKA}q$%bJt` zdi6wCtIAi)XY2x%pO;Ro&zDlDN7l*ot>lbZLAI z-VxEI>xBu8xUkG8E?s_4b+75L(Mr@0wQw)uU7fkMClzG9nYgS98Vcqr zSfpU3f-@AXQ?Nz?-nrUfEmPU0c6E;JYy$uOL=j^oQ3b@ztai7 zCU!|AK57nHf9-?%p1f0H)~bwC3ctNcB6=%)M>AvRt8^?V0N%>Fvup#5F;^k*B)@W# zMAEtvUB6WL?Q}*3u4NLi8x($3lSJrpR|DV0+_t`(fKR4R5Cez-eOQ7%Uy~f-I4I@g z9ksa9(H}_9x_jik*0aCD?|(<)X3icWpg@4{6>FHYxx|Rap8~PZkef2neUpsPM06v-jyi4@v zE=Ow70^~M*lh=^oEMV zmn;0$S`G~axOPdvdJfog=vt_+<9EyB+)!?cm>4^6QkJG^%pXDtyyciJ-Mr zx~_!7P<>}SChx^sF_PB@<`?b-v#4l zk&BaW}NNxb174oO4aarc##U7QLB+fss zeUGk3{)kx=sV1OrP`MBBW1Zzw4v$t8SD0^lcJdl?q??XL+yf)~N6^mP))=0=5qL zWaU2C4*ozp_$%<>mrTxR;7R}e^jBGMiEn)X>)WA!wH^H5+rgXg9F(lw-oR72+ncCV z__hSBP~nUIC=oh-(iw%x`9BAq^km9E-V*m_5gi{h0P@UrsKVr-Zyg+mujVp0maguwxs6Cls_fK@41h+<4^2 z&Lg*@*5-teDMRg*6P!7oV89{RCr$|s!`JoLqin05;;>bhmgW^;R~0rzdS?3Uj`~Jh zmA9Fo#tK z9JOw{tG>2&4p?IDuqVtDiX-(`)5ceq*2``mJGr=`)LuHV#BK-gxQW-;OUo3ptYk8? zms~rsxWYDu*~gVn8eLp&pEP#tl+vm8sl}tqOTls7#0l|j-aXk8+`VYu2z*j{%)Z-6m>sJf-94GIEy z>w)clc6T7)546(=u!O^21^@NgopioH5n8Btw9THEQ$S)=`>U35m5#|Ttwj|D7hF~B z^pf_l=jY_-3}xtvUbim}MV&^mV286G#XzM-r@M}hRyb#Q)unlPdIwmn1*x0b?P4pi z9i`0-)yO&~cp0m8*J9&noHFQjLBAb`KDfN@wiQHyLv(?1cdBl+h%Yrb(3h1f#ajo= z!!)zrTC3TqQ~c0&M2k&=2efgwXt#`Xm;&s^Q$>~hA}8Rk&avBPHa6O^l|AUERtn9r zH&Do-)~*h?+iQSGEvd&K@KuqKC?tW^>#1_W&mDn~omwrUFNYE&C0H#m~sq=9dCT$}>N<*jU}jp=SV!lX2a z#>ELLJ=XN3VjH{14g-;BMq_(}+Zh|oR71(JeImvUJ=#wx(8GPa2~=e$%%C39{&w;@ z=dbO~nqS*#r)Q*jqpgeS1%(4Q20ziRdDN~nJ|>r*r|oE;S}{hCkC!Mvoxntgt_f z8HmhSY?UsUx>UGuy!R%wwY1I;jZ1mC(8lEv$4A8biWt_DxJzU_x7EV$q1#lNKuie@ zablW~=xl8)g;7!-N)tIwK@9FRdyw&ns?o+pX~?xzM~jPuHUw}T*?#cER2nRPc9X2b zL|ae=y&v3y9^HZny^ O_i%lhpTg3sEH>|u@8#`L^by7*r5DJPR#68!C=&2(Qb=A zi=NwVx5W@dg(UP{Hjd&a4o$(3?uNh99hOsKl#S>_ojEor#@7PdkdxP zwuuy%m#_VWHj=|t5OJ#`F(ZFCRK(bm=kW*anEqin@KgmHfjJ^v#6`ZGf_&NDKeU4k zM<{j&8r-FEA(%P`jYhPQTgECo!d5g~#W>^0$YyDJUeS{k<%I7V-)tP? zWS{KDaSt%zm@(0L(huzl^ynlHUn`wl+on4ti`zrR1W6Cjv9lmOHxP3(H8-Fx9@Dwu zqq>VFof0Gy7c;5Sj;MN?t28e^_5>%X6g}a^&jjNwYoA%`_o<1pJu?4v&G9*EJoYt$>0 z8W*JLI~3z|s%KI)HxN?gLy453(bkD*%jM zL8{l-EDP(KbZOW_-bVgcuPY<>mm!5En+X31#jjy*giG_f8rPQT$kj6IY`C4q&l_5UIcZ9!f(xM}aP? zBiT;#>-WkY;Hb*FeEofuhWdLdkccbM{1T0S#yLzf}`rxmWaE75<}Avc!aFeJH#MY>{23SO6sOT+td7t24RkIdIl(@89? zP035c&A5%_uT%UQ#-0g5oSIL=rxN(}_c|I*-$Go`OY>>Wvq;32uh-Qz#PUzvm1sX& zCnf%fnGPy)1*atxm~hKp1bHmK{ys*-b(9%%C6@m#@?!bVo|42GVoPV-rTs@EKSCmw zKSh1#qv3k2!^U2T<^L0zSb7qDnCSHPM%GSJ@x=1ufK1~`kwrm;HQ#Q zH;JrU8Hzubl#02Gm^3A6nT7g}=GWgLRce0Ku+;zQ(t7K5q;;g&^7Z$eO&Ns6Suhfa zH2tUfHT((kWBDWN(P{xZJHjawnopw#fSGVwf4%OTt@yJN$}iH43g)Wv6YZq=2PDvI zaj`Js`!TfBRQV%R2f{{N@hz;@U&nuPFRgnFJ6m -// Configure OpenMP behavior at runtime. -inline void omp_configure(int max_active_levels, - bool dynamic_threads, - const std::vector& threads_per_level = {}, - bool bind_close = true) -{ - // 1) Allow nested parallel regions (levels of teams) - // Example: outer #pragma omp parallel ... { inner #pragma omp parallel ... } - omp_set_max_active_levels(max_active_levels); // 1 = only top-level; 2+ enables nesting +namespace omp_config{ - // 2) Let the runtime shrink/grow thread counts if it thinks it should - // (helps avoid oversubscription when you accidentally ask for too many threads) - omp_set_dynamic(dynamic_threads ? 1 : 0); + // Configure OpenMP behavior at runtime. + inline void omp_configure(int max_active_levels, + bool dynamic_threads, + const std::vector& threads_per_level = {}, + bool bind_close = true) + { + // 1) Allow nested parallel regions (levels of teams) + // Example: outer #pragma omp parallel ... { inner #pragma omp parallel ... } + omp_set_max_active_levels(max_active_levels); // 1 = only top-level; 2+ enables nesting - // 3) Thread binding (keep threads near their cores) is controlled via env vars, - // so here we just *recommend* a good default (see below). You *can* setenv() - // in code, but it’s cleaner to do it outside the program. - (void)bind_close; // documented below in env var section + // 2) Let the runtime shrink/grow thread counts if it thinks it should + // (helps avoid oversubscription when you accidentally ask for too many threads) + omp_set_dynamic(dynamic_threads ? 1 : 0); - // 4) Top-level default thread count (inner levels are usually set per region) - if (!threads_per_level.empty()) { - omp_set_num_threads(threads_per_level[0]); // e.g. 16 for the outermost team - // Inner levels: - // - Use num_threads(threads_per_level[L]) on the inner #pragma omp parallel - // - or set OMP_NUM_THREADS="outer,inner,inner2" as an environment variable + // 3) Thread binding (keep threads near their cores) is controlled via env vars, + // so here we just *recommend* a good default (see below). You *can* setenv() + // in code, but it’s cleaner to do it outside the program. + (void)bind_close; // documented below in env var section + + // 4) Top-level default thread count (inner levels are usually set per region) + if (!threads_per_level.empty()) { + omp_set_num_threads(threads_per_level[0]); // e.g. 16 for the outermost team + // Inner levels: + // - Use num_threads(threads_per_level[L]) on the inner #pragma omp parallel + // - or set OMP_NUM_THREADS="outer,inner,inner2" as an environment variable + } } -} + // ---------- Helper: may we create another team? ---------- + inline bool omp_parallel_allowed() { + #ifdef _OPENMP + // If we’re not in parallel, we can spawn. + if (!omp_in_parallel()) return true; + + // Already inside parallel: allow only if nesting is enabled and not at limit. + int level = omp_get_active_level(); // 0 outside parallel, 1 inside, ... + int maxlv = omp_get_max_active_levels(); // user/runtime cap + return static_cast(level < maxlv); + #else + return false; // no OpenMP → no extra teams + #endif + } + + +} // namespace omp_config diff --git a/include/decomp/.gitkeep b/include/decomp/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/include/decomp/decomp.h b/include/decomp/decomp.h new file mode 100644 index 0000000..eb167a7 --- /dev/null +++ b/include/decomp/decomp.h @@ -0,0 +1,4 @@ +#pragma once + +#include "./decomp/lu.h" + diff --git a/include/decomp/lu.h b/include/decomp/lu.h new file mode 100644 index 0000000..af13042 --- /dev/null +++ b/include/decomp/lu.h @@ -0,0 +1,61 @@ +#pragma once + +#include "./utils/vector.h" +#include "./utils/matrix.h" + +namespace decomp{ + + // Stores PA = LU with partial pivoting (row permutations). + template + struct LU{ + utils::Matrix LUmat; // combined L (unit diagonal implied) and U + std::vector piv; // pivot indices (row permutations) + bool singular= false; // set true if zero (or near-zero) pivots encountered + + LU() = default; + + explicit LU(const utils::Matrix& A) { + factor(A); + } + + void factor(const utils::Matrix&A){ + + const uint64_t rows = A.rows(); + if (rows != A.cols()){ + throw std::runtime_error("LU: factor non-square"); + } + + if (rows == 0){ + LUmat = A; + piv.clear(); + singular = false; + return; + } + + T big, tmp; + + LUmat = A; + piv.resize(rows); // piv stores the implicit scaling of each row. + //double d = 1.0; // No row interchanges yet. + + for (uint64_t i = 0; i < rows; ++i){ // Loop over rows to get the implicit scaling information. + big = T{0}; + for (uint64_t j = 0; j < rows; ++j){ + tmp=std::abs(LUmat(i,j)); + if (tmp > big){ + big = tmp; + } + } + if (big == T{0}){ + throw std::runtime_error("Singular matrix in LU.factor()"); + } + } + } + + }; // struct LU + + typedef LU LUf; + typedef LU LUd; + + +} // namespace decomp \ No newline at end of file diff --git a/include/numerics/inverse.h b/include/numerics/inverse.h index 031c904..e53acc1 100644 --- a/include/numerics/inverse.h +++ b/include/numerics/inverse.h @@ -5,100 +5,39 @@ #include "./utils/vector.h" #include "./utils/matrix.h" +#include "./numerics/inverse/inverse_gauss_jordan.h" +#include "./numerics/inverse/inverse_lu.h" + +#include + namespace numerics{ template void inplace_inverse(utils::Matrix& A, std::string method = "Gauss-Jordan"){ + + if (A.rows() != A.cols()) { + throw std::runtime_error("inplace_inverse: non-square matrix"); + } + if (method == "Gauss-Jordan"){ - - utils::Matrix B(A.rows(),A.cols(), T{0}); - - - uint64_t icol{0}, irow{0}, rows{A.rows()}, cols{A.cols()}; - double big, dum, pivinv, temp; - utils::Vi indxc(rows,0), indxr(rows,0), ipiv(rows,0); - - //for (uint64_t j = 0; j < N; ++j){ ipiv[j] = 0;} - for (uint64_t i = 0; i < rows; i++){ - big = 0.0; - for (uint64_t j = 0; j < rows; j++){ - if (ipiv[j] != 1){ - for (uint64_t k = 0; k < rows; k++){ - if (ipiv[k] == 0){ - if (abs(A(j,k)) >= big){ - big = abs(A(j,k)); - irow = j; - icol = k; - } - } - } - } - } - ipiv[icol]++; - if (irow != icol){ - for (uint64_t l = 0; l < rows; l++){ // SWAP - temp = A(irow,l); - A(irow,l) = A(icol,l); - A(icol,l) = temp; - } - for (uint64_t l = 0; l < cols; l++){ // SWAP temp matrix - temp = B(irow,l); - B(irow,l) = B(icol,l); - B(icol,l) = temp; - } - } - - indxr[i] = irow; - indxc[i] = icol; - if (A(icol,icol) == 0.0){ - throw std::runtime_error("utill:inplace_inverse('Gauss-Jordan' - Singular Matrix"); - } - pivinv= 1.0/A(icol,icol); - A(icol,icol)=1.0; - - for (uint64_t l = 0; l < rows; l++){ - A(icol,l) *= pivinv; - } - for (uint64_t l = 0; l < cols; l++){ - B(icol,l) *= pivinv; - } - for (uint64_t ll = 0; ll < rows; ll++){ - if (ll != icol){ - dum = A(ll,icol); - A(ll,icol) = 0; - for (uint64_t l = 0; l < rows; l++){ - A(ll,l) -= A(icol,l)*dum; - } - for (uint64_t l = 0; l < rows; l++){ - B(ll,l) -= B(icol,l)*dum; - } - } - } - - } - //m = temp_m; - for (int64_t l = rows-1; l >= 0; l--){ - if (indxr[l] != indxc[l]){ - for (uint64_t k = 0; k < rows; k++){ - temp = A(k,indxr[l]); - A(k,indxr[l]) = A(k,indxc[l]); - A(k,indxc[l]) = temp; - } - } - } - } + inverse_gj(A); + } else{ throw std::runtime_error("numerics::inplace_inverse(" + method + ") - Not implemented yet \r \nImplemented: 'Gauss-Jordan',"); } - } + } template utils::Matrix inverse(utils::Matrix& A, std::string method = "Gauss-Jordan"){ + + utils::Matrix B = A; + inplace_inverse(B, method); + return B; } diff --git a/include/numerics/inverse/inverse_gauss_jordan.h b/include/numerics/inverse/inverse_gauss_jordan.h new file mode 100644 index 0000000..4fd5db5 --- /dev/null +++ b/include/numerics/inverse/inverse_gauss_jordan.h @@ -0,0 +1,94 @@ +#ifndef _inverse_gj_n_ +#define _inverse_gj_n_ + + +#include "./utils/vector.h" +#include "./utils/matrix.h" + +#include + + +namespace numerics{ + + template + void inverse_gj(utils::Matrix& A){ + utils::Matrix B(A.rows(),A.cols(), T{0}); + + + uint64_t icol{0}, irow{0}, rows{A.rows()}, cols{A.cols()}; + double big, dum, pivinv, temp; + utils::Vi indxc(rows,0), indxr(rows,0), ipiv(rows,0); + + //for (uint64_t j = 0; j < N; ++j){ ipiv[j] = 0;} + + for (uint64_t i = 0; i < rows; i++){ + big = 0.0; + for (uint64_t j = 0; j < rows; j++){ + if (ipiv[j] != 1){ + for (uint64_t k = 0; k < rows; k++){ + if (ipiv[k] == 0){ + if (abs(A(j,k)) >= big){ + big = abs(A(j,k)); + irow = j; + icol = k; + } + } + } + } + } + ipiv[icol]++; + if (irow != icol){ + for (uint64_t l = 0; l < rows; l++){ // SWAP + temp = A(irow,l); + A(irow,l) = A(icol,l); + A(icol,l) = temp; + } + for (uint64_t l = 0; l < cols; l++){ // SWAP temp matrix + temp = B(irow,l); + B(irow,l) = B(icol,l); + B(icol,l) = temp; + } + } + indxr[i] = irow; + indxc[i] = icol; + if (A(icol,icol) == 0.0){ + throw std::runtime_error("utill:inplace_inverse('Gauss-Jordan' - Singular Matrix"); + } + pivinv= 1.0/A(icol,icol); + A(icol,icol)=1.0; + + for (uint64_t l = 0; l < rows; l++){ + A(icol,l) *= pivinv; + } + for (uint64_t l = 0; l < cols; l++){ + B(icol,l) *= pivinv; + } + for (uint64_t ll = 0; ll < rows; ll++){ + if (ll != icol){ + dum = A(ll,icol); + A(ll,icol) = 0; + for (uint64_t l = 0; l < rows; l++){ + A(ll,l) -= A(icol,l)*dum; + } + for (uint64_t l = 0; l < rows; l++){ + B(ll,l) -= B(icol,l)*dum; + } + } + } + + } + //m = temp_m; + for (int64_t l = rows-1; l >= 0; l--){ + if (indxr[l] != indxc[l]){ + for (uint64_t k = 0; k < rows; k++){ + temp = A(k,indxr[l]); + A(k,indxr[l]) = A(k,indxc[l]); + A(k,indxc[l]) = temp; + } + } + } + } + +} // namespace numerics + +#endif // _inverse_gj_n_ \ No newline at end of file diff --git a/include/numerics/inverse/inverse_lu.h b/include/numerics/inverse/inverse_lu.h new file mode 100644 index 0000000..e69de29 diff --git a/include/numerics/matmul.h b/include/numerics/matmul.h index f7026e3..cb5916e 100644 --- a/include/numerics/matmul.h +++ b/include/numerics/matmul.h @@ -3,10 +3,12 @@ #include "./utils/matrix.h" +#include "./core/omp_config.h" namespace numerics{ - + +// ---------------- Serial baseline ---------------- template utils::Matrix matmul(const utils::Matrix& A, const utils::Matrix& B){ @@ -19,10 +21,8 @@ namespace numerics{ const uint64_t p = B.cols(); T tmp; - utils::Matrix C(m, n, T{0}); + utils::Matrix C(m, p, T{0}); - //#pragma omp parallel for collapse(2) schedule(static) - #pragma omp parallel for for (uint64_t i = 0; i < m; ++i){ for (uint64_t j = 0; j < n; ++j){ tmp = A(i,j); @@ -34,6 +34,85 @@ namespace numerics{ return C; } +// ---------------- Rows-only OpenMP ---------------- +template +utils::Matrix matmul_rows_omp(const utils::Matrix& A, + const utils::Matrix& B) { + if (A.cols() != B.rows()) throw std::runtime_error("matmul_rows_omp: dim mismatch"); + const uint64_t m=A.rows(), n=A.cols(), p=B.cols(); + + utils::Matrix C(m, p, T{0}); + + #pragma omp parallel for schedule(static) + for (uint64_t i=0;i +utils::Matrix matmul_collapse_omp(const utils::Matrix& A, + const utils::Matrix& B) { + if (A.cols() != B.rows()) throw std::runtime_error("matmul_collapse_omp: dim mismatch"); + const uint64_t m=A.rows(), n=A.cols(), p=B.cols(); + utils::Matrix C(m, p, T{0}); + + #pragma omp parallel for collapse(2) schedule(static) + for (uint64_t i=0;i +utils::Matrix matmul_auto(const utils::Matrix& A, + const utils::Matrix& B) { + const uint64_t m=A.rows(), p=B.cols(); + const uint64_t work = m * p; + + bool can_parallel = omp_config::omp_parallel_allowed(); + + #ifdef _OPENMP + int threads = omp_get_max_threads(); + #else + int threads = 1; + #endif + + + // Tiny problems: serial is cheapest. + if (!can_parallel || work < static_cast(threads)*4ull) { + return matmul(A,B); + } + // Plenty of (i,j) work → collapse(2) is a great default. + else if (work >= 8ull * static_cast(threads)) { + return matmul_collapse_omp(A,B); + } + // Many rows and very few columns → rows-only cheaper overhead. + else if (m >= static_cast(threads) && p <= 4) { + return matmul_rows_omp(A,B); + } + else{ + // Safe fallback + return matmul(A,B); + } +} + + diff --git a/include/numerics/matvec.h b/include/numerics/matvec.h index 0f9d50c..0c137e6 100644 --- a/include/numerics/matvec.h +++ b/include/numerics/matvec.h @@ -3,11 +3,13 @@ #include "./utils/matrix.h" - +#include "./core/omp_config.h" namespace numerics{ - // y = A * x, where A is (m×n) and x is length n and y is length m +// ================================================= +// y = A * x (Matrix–Vector product) +// ================================================= template utils::Vector matvec(const utils::Matrix& A, const utils::Vector& x) { if (A.cols() != x.size()) { @@ -18,6 +20,27 @@ namespace numerics{ const uint64_t n = A.cols(); utils::Vector y(m, T{0}); + + for (uint64_t i = 0; i < m; ++i) { + for (uint64_t j = 0; j < n; ++j) { + y[i] += A(i, j) * x[j]; + } + } + return y; + } + // -------------- Collapse(2) OpenMP ---------------- + template + utils::Vector matvec_omp(const utils::Matrix& A, const utils::Vector& x) { + if (A.cols() != x.size()) { + throw std::runtime_error("matvec: dimension mismatch"); + } + + const uint64_t m = A.rows(); + const uint64_t n = A.cols(); + + utils::Vector y(m, T{0}); + + #pragma omp parallel for schedule(static) for (uint64_t i = 0; i < m; ++i) { T acc = T{0}; for (uint64_t j = 0; j < n; ++j) { @@ -28,7 +51,34 @@ namespace numerics{ return y; } - // y = x * A, where x is length m and A is (m×n) -> y is length n + // -------------- Auto OpenMP ---------------- + template + utils::Vector matvec_auto(const utils::Matrix& A, + const utils::Vector& x) { + + + uint64_t work = A.rows() * A.cols(); + + bool can_parallel = omp_config::omp_parallel_allowed(); + #ifdef _OPENMP + int threads = omp_get_max_threads(); + #else + int threads = 1; + #endif + + if (can_parallel || work > static_cast(threads) * 4ull) { + return matvec_omp(A,x); + } + else{ + // Safe fallback + return matvec(A,x); + } + + } + +// ================================================= +// y = x * A (Vector–Matrix product) +// ================================================= template utils::Vector vecmat(const utils::Vector& x, const utils::Matrix& A) { if (x.size() != A.rows()) { @@ -38,6 +88,26 @@ namespace numerics{ const uint64_t n = A.cols(); utils::Vector y(n, T{0}); + + for (uint64_t j = 0; j < n; ++j) { + for (uint64_t i = 0; i < m; ++i) { + y[j] += x[i] * A(i, j); + } + } + return y; + } + + // -------------- Collapse(2) OpenMP ---------------- + template + utils::Vector vecmat_omp(const utils::Vector& x, const utils::Matrix& A) { + if (x.size() != A.rows()) { + throw std::runtime_error("vecmat: dimension mismatch"); + } + const uint64_t m = A.rows(); + const uint64_t n = A.cols(); + + utils::Vector y(n, T{0}); + #pragma omp parallel for schedule(static) for (uint64_t j = 0; j < n; ++j) { T acc = T{0}; for (uint64_t i = 0; i < m; ++i) { @@ -48,6 +118,30 @@ namespace numerics{ return y; } + // -------------- Auto OpenMP ---------------- + template + utils::Vector vecmat_auto(const utils::Vector& x, + const utils::Matrix& A) { + + uint64_t work = A.rows() * A.cols(); + + bool can_parallel = omp_config::omp_parallel_allowed(); + #ifdef _OPENMP + int threads = omp_get_max_threads(); + #else + int threads = 1; + #endif + + if (can_parallel || work > static_cast(threads) * 4ull) { + return vecmat_omp(x,A); + } + else{ + // Safe fallback + return vecmat(x,A); + } + + } + } // namespace numerics diff --git a/include/numerics/transpose.h b/include/numerics/transpose.h index efdca2e..7e6cb5c 100644 --- a/include/numerics/transpose.h +++ b/include/numerics/transpose.h @@ -43,28 +43,6 @@ namespace numerics{ } - - - - - - - - - - - } // namespace numerics - - - - - - - - - - - #endif // _transpose_n_ \ No newline at end of file diff --git a/include/utils/numerics.txt b/include/utils/numerics.txt deleted file mode 100644 index 69cd323..0000000 --- a/include/utils/numerics.txt +++ /dev/null @@ -1,11 +0,0 @@ -#ifndef _numerics_n_ -#define _numerics_n_ - - -namespace utils{ - -double random(const double& min, const double& max); - -} - -#endif // _numerics_n_ \ No newline at end of file diff --git a/include/utils/vector.h b/include/utils/vector.h index 753c5fd..3786271 100644 --- a/include/utils/vector.h +++ b/include/utils/vector.h @@ -66,6 +66,20 @@ public: bool operator!=(const Vector& a) const { return !(*this == a); } + //################################################## + //# VECTOR: nearly_equal_vec # + //################################################## + + bool nearly_equal_vec(const Vector& a, double tol=1e-12) const { + if (a.size() != v.size()) return false; + for (uint64_t i=0;itol) { + return false; + } + } + return true; + } + //################################################## //# VECTOR: Scalar Addition # //################################################## diff --git a/makefile b/makefile index 052f8f7..314f2fd 100644 --- a/makefile +++ b/makefile @@ -10,6 +10,11 @@ SRC_DIR := src INC_DIR := include OBJ_DIR := obj BIN_DIR := bin +TEST_BIN := $(BIN_DIR)/tests + +# All test sources + + # @@ -23,13 +28,30 @@ SRCS := $(shell find $(SRC_DIR) -name '*.cpp') OBJS := $(patsubst $(SRC_DIR)/%.cpp,$(OBJ_DIR)/%.o,$(SRCS)) +# === Test sources === +TEST_SRCS := $(shell find test -name 'test_*.cpp') +TEST_OBJS := $(patsubst test/%.cpp, $(OBJ_DIR)/test/%.o, $(TEST_SRCS)) +# The single file that defines TEST_MAIN / main() +TEST_MAIN := $(OBJ_DIR)/test/test_all.o + + # === OpenMP runtime configuration (override-able) === OMP_PROC_BIND ?= close # close|spread|master OMP_PLACES ?= cores # cores|threads|sockets -OMP_MAX_LEVELS ?= 1 # 1 = no nested teams; set 2+ to allow nesting +OMP_MAX_LEVELS ?= 2 # 1 = no nested teams; set 2+ to allow nesting OMP_THREADS ?= 16 # e.g. "16" or "8,4" for nested (outer,inner) -OMP_DYNAMIC ?= true # true/false: let runtime adjust threads -OMP_DISPLAY_ENV ?= FALSE # TRUE to print runtime config at startup +OMP_DYNAMIC ?= TRUE # TRUE/FALSE: let runtime adjust threads +OMP_SCHEDULE ?= STATIC # STATIC recommended for matvec/matmul +OMP_DISPLAY_ENV ?= TRUE # TRUE to print runtime config at startup + +# Export OMP defaults so child makes or tools see them (not strictly required) +export OMP_PROC_BIND +export OMP_PLACES +export OMP_MAX_LEVELS +export OMP_THREADS +export OMP_DYNAMIC +export OMP_SCHEDULE +export OMP_DISPLAY_ENV @@ -49,11 +71,19 @@ $(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp # === Run with OpenMP env set only for the run === .PHONY: run run: $(TARGET) - OMP_PROC_BIND=$(OMP_PROC_BIND) \ + @echo ">>> OMP_PROC_BIND=$(OMP_PROC_BIND)" + @echo ">>> OMP_PLACES=$(OMP_PLACES)" + @echo ">>> OMP_MAX_ACTIVE_LEVELS=$(OMP_MAX_LEVELS)" + @echo ">>> OMP_NUM_THREADS=$(OMP_THREADS)" + @echo ">>> OMP_DYNAMIC=$(OMP_DYNAMIC)" + @echo ">>> OMP_SCHEDULE=$(OMP_SCHEDULE)" + @echo ">>> OMP_DISPLAY_ENV=$(OMP_DISPLAY_ENV)" + @OMP_PROC_BIND=$(OMP_PROC_BIND) \ OMP_PLACES=$(OMP_PLACES) \ OMP_MAX_ACTIVE_LEVELS=$(OMP_MAX_LEVELS) \ OMP_NUM_THREADS="$(OMP_THREADS)" \ OMP_DYNAMIC=$(OMP_DYNAMIC) \ + OMP_SCHEDULE=$(OMP_SCHEDULE) \ OMP_DISPLAY_ENV=$(OMP_DISPLAY_ENV) \ ./$(TARGET) @@ -77,4 +107,31 @@ info: @echo "Source files: $(SRCS)" @echo "Object files: $(OBJS)" @echo "CXXFLAGS: $(CXXFLAGS)" - @echo "LDFLAGS: $(LDFLAGS)" \ No newline at end of file + @echo "LDFLAGS: $(LDFLAGS)" + + +.PHONY: test +test: $(TEST_BIN) + @echo ">>> OMP_PROC_BIND=$(OMP_PROC_BIND)" + @echo ">>> OMP_PLACES=$(OMP_PLACES)" + @echo ">>> OMP_MAX_ACTIVE_LEVELS=$(OMP_MAX_LEVELS)" + @echo ">>> OMP_NUM_THREADS=$(OMP_THREADS)" + @echo ">>> OMP_DYNAMIC=$(OMP_DYNAMIC)" + @echo ">>> OMP_SCHEDULE=$(OMP_SCHEDULE)" + @echo ">>> OMP_DISPLAY_ENV=$(OMP_DISPLAY_ENV)" + @OMP_PROC_BIND=$(OMP_PROC_BIND) \ + OMP_PLACES=$(OMP_PLACES) \ + OMP_MAX_ACTIVE_LEVELS=$(OMP_MAX_LEVELS) \ + OMP_NUM_THREADS="$(OMP_THREADS)" \ + OMP_DYNAMIC=$(OMP_DYNAMIC) \ + OMP_SCHEDULE=$(OMP_SCHEDULE) \ + OMP_DISPLAY_ENV=$(OMP_DISPLAY_ENV) \ + $(TEST_BIN) + +$(TEST_BIN): $(TEST_OBJS) $(TEST_MAIN) + @mkdir -p $(BIN_DIR) + $(CXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) + +$(OBJ_DIR)/test/%.o: test/%.cpp + @mkdir -p $(dir $@) + $(CXX) $(CXXFLAGS) -c $< -o $@ \ No newline at end of file diff --git a/obj/main.o b/obj/main.o index f843db1765dadf41dee09e9e1adec6fdfc157fff..589799b2d52074b912079c939ef312c0a1cd7ae4 100644 GIT binary patch literal 2816 zcmb_e&1(}u6rVJ;ty55@}<%P)e{Aw3l_0Y)!D6gw3QMB7UIO zmLe1c5l=mO@aDy1P*TxsiGjHbi-tWDcH?x!6ObniH z2m}~OfSqNlGeSW(t&Zbj7-t=9r|?BnGc@G2K-&zB6aW_WgZ~y&^}cH{w51UH)3p=RWFpKD^g&e_4c38GXL-@KG4(b9W<;LB}f9<(qwx zkm$prkLwSvN3vpA6cJwcMWW!!mhMJwp|1W0MzZ>2@)OR`J6Rir-O;=oJ)<9U)0G=z zl{c@1jrI5To(>-!NtN@w9FBEs-CFd-sj`@M-RSDkqCK4whUksZLeMrr>mM3n4UAO~ z3qnKba1{$&yTAfdp}?Nj=9UF0ICAz%+X1N&6B_|dSOmsUZXr*wPqc+*8v3_2HOy{- zQ0KO~rn)_zZ6o`(;^z=lFNfj~6t_VE0< zX#Wo2e*V|$_&*eQ3-qo~oM$jSsO*dID&RYmeTfu)$#048y!DmMfHulv%Qj@m$I=ooOSN&v6fn!gYyonepU= zy*4a4#Omd9naT=jA;`czFvBQuGkw)aXRjJ#W^Nod46@5r;$Gu>@N31Ddia>)->Uc* z6&&vf#ldq!!EZ#x9|_dO+h9!oVTt2SQg2g_g44H#kgnjEqb6Eu+U918X(jHsI!hyN zO>w4WOl~qQRVp#fvBpgeQO);|v9|SmsLUQ8D7C}OMW&?-wr%AhUB19A4PJ~B;uPfk zO1`WWokG!a__R!6WsNb%v@NDxaS9OOKgmmov)w0FTw-{ZDCmQ48`Y9@J^mn)KvX{d z&Jd6E9>DO`P=20Z#@I#CukX-T75sh?Q_mjl|H?lCzW7RPlph1GSigEH0lN$h1=pvF z_boXNMyy5L!@npCv|Ey~{#YY=(Fh*_%rF0flJo zqlA0+&%q$_lY`W=GXG^c$Nl3w;mv;^Fx-C>nku6Rf?$1Bk+9$(j&&kT^_%eneS8-I z!~ag6<72au9|oH$w0^=~zyMvn{HKziB}b{f{O5rC)$fz6pDuY=zjkQEr~BiV8z~{V tpsZZS4#04qUOxUdc)uNh*pvj{;Y1M`u7=>c{PU)jvM_G9f_BY!nm~VZ~^Ppc18O5Njfd%)$myp$0(&AqbW?lU+cO7uhVj zE{pg=A6xY)R;~C{TVE?$#Ux}&po*Y^sIBs%!V*+K#3Ug3eb2dfW@k6Myja_R|Ig={ zd^R(4?)$mtoO|xM=iZrD;=Rglvnj@9Q+}bu-wddt^m!)^Eo5+^vP!AtxcPX$5Wm~; zy92-3_)+Y>`2AcyUyS=D^0@%_e)wG`-}~c!xqQ9?_d@)NasMTL*UR@2xZi-^Ncnyv?xXR$Nxt8V`#Ai@ z%lBJwpNQWi`F^Xs{|fiX@_CB9PsM$je4dW`ZSuJS_e%Nf$9)EVGv)i8xX;0Fu6&<| z`+WI)7w!x23&{5>+=KWnl<#-reh+>j`F=0%Vf+@!_XzHb<@0^GFOkppax z_rJ*JS8-n#m2s*F)d#qzHPTE>jPEs)xQAd7~bx zjWI0TP7lx8soTB@&vciKyfLuO#_%N*&eysV=0In|D93}y(WQs4(8B{7ie_qAl|rsJ(0CjYXqY>EXf1>H|HrR?lnE3wBjr ztVepEtOr`0bIW9&h>*A26FB6Yd#xUN*BjcVhrSj;p&GAks~#Gh+4yHY@Pl*i6?$lw z45BtZK@u+6Y0aV(S>*B>Ym7)R(dG(qQ3$=4$*TT zJO_3-=<(9{V^8##9$K%5S|%*f19#?@DOyz%nx}^=P{trveyd0O^8s(9|K7Tq*{-bC z(%@HG)iYcPpJD1u0#k{r#Pb8~h z!V)0Ot_Nx}^+0uI-L}2Y6gC0bg`nNF{?E1Ug^J!AT$|G@U*V>7VOc2dk(j@R;|&hzD59Y9o#3{G*ge1Zys4Rs_4d|(M31ixQg*r zXxA3Lh&WeAMs;bb+cxLa4f)juf7R8@aXB4z+rC=v%!vMwqm7px4TZStKIsrb2A=u9LsH5Rd-4n z*m796?a^(ux~+lQ=&S~gHV`iCsI)!Sf>sz<(1QK�GVjKtp1&in(74-pdjG`^v? zR_DJLi~P{q+P0zN#!j6&&F3ChQ9jmJ?ymHW@7Hfe`8eOS3isr(6}OaExP24HPIXTi zJL{&(Nq3a@%mX_38au&Fvr_7vGS)OM^ySdyw^5Wk3qd^JRSoZ~K<6YDne=^>9R92Scn^-YzrhD2{ zH`UBNWm4sovA%H=MRv;6X%$n(PNsn5*msmyOmkB($DBTGW_d+V(PXo7N62~7)ajGQ zjw`=u?D+9*a!yr`Y17Lq#sZza5@YwIBK)^ZbR(a&p@ebeq)g!w`MJ-u#uu!HH0k4~ z$$!G+X=4*hkohCN#1K`BDpfv=|4CzBe=*AdcBA6r=| zDmr=Ew94{$E|gHFyLoKoq;YOvd8Mz?P4UqI+a3*ei1PuMZffOpN6Y|uQMN_9uJ+BgYMK4TY;^KPX5yGn< zx@kJZHkx6^A3te|MJOT?q8>e-XX=dd3bguX1IPL+E6*D=tz!JxsXt4FDK6gpNBWWu zwKQ1m>!gQgLRG2P^H#G~s&y|`FsDkqjv8NxSx0- zYR-)o4U0B^Q_V^m)XsK#{XiE*g);5$jG=ue5A8Ruvgik)J)|Dj?7UnH9z;xUB$w26 z#T(keYMUqG%RMwGlI5=3+B|hL6mO{Hb=$Z0sIloa-`Ks8t8IGVa2IVMCGv*Wp%PQS z&}+ViMlh$#+j`BHsyE`y2z=@Gog6)slc`4v-JXJK-{}(!{ehtP>2Bx^g$VM0T^efe z1hzC41-7;u6mc%_jl*E=&8ycVGf`BleMUDbV`FsZ>a|MU?eLwg2fjQ6 zCF^-TGSpcd>Iv=AmxJfb;?P7Vy#P2dlU4E1V3!3ICHd8j%TPq)cB{%vXobeFhdx6P zBf>yPYpbo}$V0Xd?aFSuswkrWG(CK>9=-}n;JnQI1A5-Odfs}ypvgZ#53SZ~4%>+w zP!M$N69~DDIcEEzljsF|DtjWxc~`V$lygC}X;@9P=>p`9S>u-c1ED>uEc)C^(&7N} z`#=o+e#t7z`J@a$>=sX&Qa-g3YDlGf{Iv2)4DLQ~ZJbspM}yu_QLgihOdYrlHS2ks z^qOY7UhusZ97w8fuoe0jJp@|2Y76t|DfDjsfoRk9)#!PLUb90b{WE&Vm3v|I?<2DM z=z-Pl=$_Rd4adI?t2rL2~Mhs$8Rblw2K8j4>_K)v9w z|BrfT9lEpWKutTX>W`>2scRwN^m}h)V5=w4WS?AY^6=QZbxo`QzTvk`V@KPQ@DuiIAiTJfbIn{|~M zL!o*TBy-ggYhQ5tfe$MTRb3^p6!Sm<9sU;aAI+b{%!x@!kNljvf@dYG1SH&HJp=%KQoK@f4R4`h*`=uw{D$go^YC_dj()LPNRlRKW;6;9($IIFynmF2c6(|L*?C6P2Xi${Wzxj z>;_4^{Zmjwt^YZ12op{SgadUThxuqL6`fbBM+RYX&}^S^F=(N0Ko4!u16vPak`UPD zh&GH^cN9|Sr$;&kk(r`LbkL?N>!iAEx4l;5E!bMw9k1S#Knc~L1cq&??8q?C0p6o~ zu3r;vI=?}?`vkB>wB?#=38}dwL#xUF%v12GR@H#a1-?ELl_r^pS?FfH#zYw5t0&{@ zgVAeA4sOv4YW=&8*5q70k`;(Ht!@BoxC)}xqfwun_6aEPdE*VtItjIJAKf>)0rUR$ z)L?36FFhh#X$|!jY^j9cEItXs(v>OGWD*C$05h`a1~E6h(G$73X%+Q;esydz=3Kay zf!+;i8QTRIUYPuULZV-&s*kJ6qpBc!pq1#EY&EY3tdarIUAbCRjt;3&hlGV(! zSzAEXG{|;Qu}}?-ZEJMf-=cdb*XW@e$zoHYc*1>iJ(xH`KU*J;Go$5Zw4P zJrwokVNm_S?kTAAT^8MX{U`aKM7N*2j^$76dKgwx%CF;qnsaGS~8;tZ>X-T zxg*EwU90*N zT9t(QMgp0rj}q!P>bAycY;vt9G}?s<&&%+H{XqYX!2Gwc7&zzNp@%+(RZh}T@fk{K zxOc86bX6wzk2h~8hH4U1dci6!NLJ7MPfY!y#Z~;6f?QqR$AmQueEZB$}{o?_4$84Zq0Upa*oD?OaeI#%9?829zO22V{Bdc0l>-b%^|?63o%t7K^l zJ@gtWa_XT?P`&f3NkN@?iR4eT2P#%e87W&i&IM4vyam-2G4#E#7hA`2q6g%=0{JN( zt)!ZiIo5{u0f5mPSxZeQu1`arXW63r%ba#db0DvLSt*cY}=+3eXNcIda8Tzg5`{fR;iUciGU5=|q z`q~ABmN_wuAu9A_Z^3FUa3RO?*)fjQ8&QgFPKx6N)&eJT_%J)kd4!i0+whcQPegzM zxxPk-@>Th*u|Y`lFH~t>cBt<G={){*} z{fsyQ{UGc=o_?f}_DDmQkT+tW3PbE*E`^aQRE_zy7_U|v^RD(v9Wz786(Siuc@YeIOk$(@FN9~JFyF?2-4F$v-w(H?5p#xo=>CHnwlIl~V1#e*0XN_+pRvq_xBOa`@T@UM)lX2~# z*>P9ntu#{RB z$H+57TlfKL>ItpMulCf$9G;r3s;Mv6)!gaIf>kIds~zb+-hu}IZavgM6L>b~kPiae za_V6_ge|8TYm2n*33aHh1}e{d5LE`*^*Y__2{n5#qsLMx3_DttAG!wu<($2LM%`Vg zDTxf{>mW&+Uhu9KxCU_uUYzNKu0#`f5nj-s1%7UZW$L!o1b@+N8A^|JMmU+cs$-RS z0;KW&)O!1oOg6`#f&>P%@Wu>!X#6g<#qs1P9qL|(@;WK>n;dG3L)q>4St3!c|mJrCGr-m(*lnWp@lmT4SJxqhA)9w6I2L+_B>PJ0pBp9 zKJ=+g3+58dMuxP)gv&i9Xwi7{Mn)vl;yFxUc<{S6V;>8>G>EH*25x)-9);aecYSd&3`o=ATL zgYYmdGzsr}EpR;%7vRP1on9McpGkOow7_L%Sf_r!&vhI&*xiQ zFmuLX6&bE4GwPW4Hr z@@pq3)x0PZZ&=H(_Fk6lv%|k2d_8O&zs!UQ+FS6R|LfA=8sGGyUzX&5QZx)}5lboI z@cgf?%l||VzUM2#1kk02Z+Br;;bE`sGuZA=wumpSJLdJUr`P(rl!W`vFD~eNx9{ux zYMQ0;^d9RkdTX@Qw%u!c-%}9tKg&a9<7a#mD>GET<6XZwT6@M64MTf|qd6O?IoIWX zC0cWlr8S#HYch@c!-V6k$4bZf&gwhtzL;plD)^gN>;Jl_|K~~dKPG!v6Q~9k+@5}% zJuLH|v4^wN^#0b>2vC7Nu}(Vk~*# z8+Bkvq<>F6L|$$6-q36}2k7t<3(fAL2fy--DGh!O?}{c zp&Ha@EnI3s@GXN^QV6LnECqQ?{zoK_Y`U#kFWBVY?E%gL)lN^rPCp*VEpIck82b#Z zLgmGpqijU}ni5zJ$Sujbz-E|&;_JrKro!7|27Fb1el<6(# z?w?l*>sw_>{+dCNu8)zC73P!22zx~<=*`nw;|cB4Bfn}YE%>N%Fcc`)pvW*)(QRuX zdu?9Z21{$Pe!<2N7Y4p+YI6tt3ZcjDBTZ;KYd z<0dUGK4VksY|2|U7(uX7jVWtC_!7y5JQBp)m)W&oM;^(uz2RHYaR({c*1GK<^{BQ%g{%=2S)UMWR@ZXthGy`sl`S7PBmoKa>4sdxx;{;_!o2PNy zY0ljIHTmxW-bKd(vs+SzH#9j{M=8a*mvXLd^0@1x*Tcn*N)p5klYH5Dy_aa5e=>Su zqDg=lPmC3BS?7}c&xQY&CC#N)LOgQ!S2KN@Y&3*S-B#_~0-=QkPa z`B;n+>-jOUp6?((y1-%Aj8T~Wkr(SfFk!?R*ql@UW?jlZ9gcqHQoh^5J{W%~2tUhu zKI>6B>rs(TEl^?f$d#Gi(7E6!TA{xZ)B#GwixUsD^aA(IgcW*A0pL&Q<)2vlj}Uf* zmHJJ#7wpPvTB#S8#P1&6;?m0gxNs3>$5uB_vBDR4LhqL5ed7(a=rt{vYwfEz8?vu5 z+Y>2%MbQfyd@wH$1QHUwB}xAX&x0#44A1ZCHLXbZy|#dyKY%_vPpzpQZe=!U3mGBJt*!vQ0st)YG4~VP%I$_3erCHf)9Ns zkp~5dK&*v-=u_zlLIa+;iS4%Zkp3F06|Jq^^}G#w!Cq~_A(A%*`?Ljj(L=;}qGx5? z<$8U^sru_?ih}jJHn0}Hz;c{}^S3Xzi)4E|VQ7Mr@TG*vtX$yZ73xz@1Ol4eEWQQ~ zhoEsya7P)0vb>>p35Q_2^@g3^8|idPV29gxD(H@IXosh&ubZX%SdY!2*=wO213)49 zv3`wgK)mdag4(~XRegy<37>gb`P!>TX4rUca|TwsYtX*_uMtk2SFl$~FT*Rv3-{H? zbiyl0eIMa(a-4Swcbu9kD>jR+n zYI3L|dtH;5NX7k}91A{8THt%4D#QS# za}LG*h~~i~Nw|rt*Ate20-r%L6F+`ROzW>UXlrLupFyE8dI^-^zYrR(R<#Ies52pR zV;*aI?rPv7b}OLBd?#+HhI~&b)pD$%YRf@JEXR_yrXxa-s}M~ooy1PQ`VUY)A5qo6 zsLBhf`j)C}RWSg9c-$wsf-xwvQeXkm{HCS$4vjDj)ne5 z9f|%C_s4>|Zn5BLCq7$=p=xH8fv9gz;|6s8aE)lNo5dWi{Shmb*S9G8n z+7G$52tVTL6xz?|Wa>XXh4wQ#f%==FX;`!$1Rjn0NxwT*{f9ftTd>2T|CENdP}lwc zLH{A667+v_KWPZs|97xjVMpnNOw0hF{%eLM2mVeVh$QvDifLpAbdKWOu0qFyqgrlP zITiSbrUF|h|4mK>I+81Lle@CDRMrcXG4)xlf>9lf2l z@GxoYp-s}E*PNp3wsq0nWn?>~$)&ypS}@kH;D$!d?6n|F+C@!@IcB7^kUkXpc8h&R zcUt@r-n@Ith$K6pAfArfHZ6Df3XK(9M8lFfV!$p*^LTH;w#u{g$S|8+G8ZAaASTTCL`GU5m&ciufTldZ&rO7wQGFW62;ep^ls9*$&&mEVI>wLRr)~e zVz7XL`fQyY?gqeNxU3pw+8SwNI<6zhJM_4m?`=w)@6i+}sl(q89WE7ap6_udH(g(S z#E$k-R=kUvlVW=^yR4HT0h2y7OS`8T)vT*&Id?TRj{6!-P3Z1RIg?b~y(GGO#1Z;? zC9*zpe~&n7ed=CJMVjQ<9SG$IfEoMV&^|6~7Z+Ab zh1F7FjW1!6$0Srl{j6XM4SEL6ZXcsP%7B=cpr`-M^5aPKbaeUgFQTV^qx=}rR({Y5 zzIC|Z-D;CSDPq0skTF!S6rlqMXb%Q9mkawK&U-NY%_EhE$ed9Mk$@$k{$x=j+p8z! zggg12*xVnvlB}aJAUcr)bRS_rG%c`Myp`>A2zO*(NZu{jVU-|HDF{Eqe%L_vL+CYM z!hW~~zD~{;w5s2Lm`Nz)?Gwde1UTIrqJ2@YAHoQWm<|jA&KD4UiO-E;!UAHBXE@s5EZ9CdJ{=!MMh&smxT_Pq{rKGBG0a-)>*ZTknS z4>ZK4>z~zBrASv^nSiy=(mzeK#XqTk9A^LMtr!=vJ~3DfzNe!1zP4MJCR|v1kX=*l z>nzqIXidWC57@`a*nPV&c4uv%euxxds|2>i6x5Slw#bDS53pSiyDXkuM?@G9d))X< z&M*PDB;VlapjA}@eZ$GoJz=3@zaGLa4qm7z^s-O^dYXO*^Was4tw%G#mI}3o{|Ub=zid-ZoDN%K(8EB!34sE^i2?w=YF9%pf?WOcd$c6~mKOb;yBI0*zQZ z_?S08Q8ylR%_zf;C2-iO4Yc2>CWfBeszN_MGq4r=MX}P7>&EdCIa#^zp!aowNw63c zWnEt)&wA=lt?CsnVB4Y6Nbe8Or=Ca&_QTAkO&f59ubP8UfUvm)AapH2w*thkuGmd@ zH0_Ov>QSMh=dD8fv?^Zl1n%Gf)ZJ460NF>|Hc4A}4KO096U{=;2m*y)(JyHAM)O1l za;r2Cmb^|{lF%N2NMBr!nam;VV)I0t_00c7pXPcD(I?Q){EsKl=Paxp97bo%M?+}I z09{61$fNC-z((7?oCQy$?J* zbZK7b(L~=^*}~B=3qE1`-&&ftm)U`K%j$2Lt}cv3Am<&} z{@nEvT5&o>k6h&jaxy#xt9)JKA{wB$d_tI{=buyMctNi3CYDiG5e;N};MxQkbrsRT zkW6|b88w7_?$H0$NaPKzWj4mjrO}{xW4IAs+uus_c6e=#Sj0i1UvVOI{+E!#Og(XO z@!Ag3D)7d-?Xe!Ku*`G7ceW=y9SgSKV|kvl2eEnc4sZ&r?GuPZhL_Nafp(r!YWt=% zG|PqI$yf`aQ)MjU3pXxZAg!m`*2g;=^C>J4VySNn)`&sjvHs8wEtv^@HTE^2#(o}s zg&>qx4V{j}m~}aM*bx~V#A=r%g9%Hp`}rWPW8oa2Tq41%n_?FOVA7u-7US8jYzopq zh_Rh52^k>$@3Fcm(xYp%*n@brl%NHlM(rrvnTv2o?BS#)1V_p63@W2UiFxr(CGAHw ztV4xy;=vJ|Lm;A9U$%}R@nvpomJIdJP4E-xLVhCfGVTT1r4fa;(-D1y≤4e)7d& z?;h)f4eNjYpOLde7t=383kHaS!|(}7?0OGo!x1kuEpid@gv7i+4`CnW=F)I?DE>4I z(01;~>{c2MsGD#W*e6=mduUB5xY6YapX>>JS`_%wf&TV@)ab}1ES|N(m3R|Hia49G z*M9HyYk2&3GSYd7lKJSec=}u8(tkf;o@LE=lSOMEyndh-#%?*2m@0#uETT}P@HLrP3gk0WS z>8g8tw(HIlUAnEluD0g2tY+%TtZ((!)t-Xyv>@#ZC%Lv39lOL+@BEgsMxeLgYYb!n zCc^xGk1YdGv?w+aiBqj#nY&UGx*@@3WmRMM)GjMTBebZMtB3pNqGrs_iT&w+)zs1A zt6fP^9x{AP$iHQoO$|DgMhAI-{?0iZM;XY2H;zoM!k>xPh2cHNhqj%TK%tq^lc`sR zpwP`)@LHy{TQjAyN%~~HGnbHRB~p09`O!dhtA$7hP*NPthM~{!U|U7&(4npD!o-O1 zgg)|wJ}C`%$_2r0hwO5x}DnSsabWPGoR;Ydi~7Vx_ToKK>PHXv`92snv%nVpCstD9aq`M-Jr! z$1rLZ(+}&O}0kbQ)=kX@YR9-JOwHg>CeqmZKm>!5Gc>Q z{$I1*ZYt04DENkzXG3ppo@;n(AUPh%mB6AhL-E8Zf7;^TB9T;?Ifbb*i$szBP?@2W zy+C?opb_ptq%y-lyaF@RgP58AjSJliMWamW%`lU1Rk5zgdUH3alg*t67a4d6CC^U} z&&!M+o+tdX*k4G86H8!gQ{c;X&^H&s5RQUrRL)x~%47p&XkWB>dIPH@kh?UrXDy=X zsaDmV>Pb`7=DDPz;$%Ib5CE{_Ty`6R!wewwmKpmyo-EV|`E1n%*dJHw+iww#0^O#* z&0ERi_z+WP!pGaFfS9QhK1ku_El!5Q4Ua8SxaZE(XtsVL+HP%s;`|&4c-5z4V(x(s z9=R5dJveHl3%8t`{AycADbZ@_C?#7ArHoqZ3C$tJ99@35KmceWQ71R(WR6hB=LmIt zO_Dke?hdTcR0T7+fkb~euMh6Xff5%H=|0n75-Gx4HH}v^tLi#);I;M@FAF0gZg;j@8wbbE()LNjwJ@J+oT`kdCm^r(0ju2v&}(>^_`{LK9zWdMzP5xFS7e=V)7z?v1X9gv&fuv zez<>QvY(Y}HW-@@4OM@Y@GTq4aZOd-Toj`?|z>V<1Y3Z|L)gZ0PfB<`W6p z=<`H%?;#gsKA;8f2l*61F%7kgElkqQvHqq7?|?>6;YFc=hhkF_)%yuc{24{zjG|CG z(#X3bPz#bdo%3_VNO?=BZXwuo4_xOI3SfHCg12yd=cM(f zBKc}F|KNo_7_jczt$^7qx?NRn6X+A_U0`lzmR5BVUgA0kod&YHDwk%sm;~T-kWigB z@@y_|Jc|gp9;)JZcxE9JywFmP9;rNyIl`O=UZMArL4l@IwS{XmoI!jVa+Pk}hR0Wm z__R@`7SRn$nQyH`zlezOzb1gV3o2oD@*G(EVbm*$whp9zmm3$sq;H)Q1zu_tH zA)SgV9Q)y{N_$ay!C}#?!jO(j%JK;u` zj^{q|85OF^b>oE`{DO3@HL1WiQ{{u}C*`?q3IB~_+mghX4zi=#>5)ZTeH=fP9?YGWQIMLI-)3!jTeLEvI)}n4+9bcRln_ZZ+=M%jOOL z3e<&ThjnD;rlP|mF&v2MQomBk6k&`fa%V35hu_Fu1_&sb9Zc~`=Zur}Z5GZeK{yV( z!I8?)Y3TH3Qg`VzKct5d6Hc9U^fQhO*26C|-SP>~!*UBLTRKY|ad7aCH}ow6JUBWI z;Cm?hnlC}8U`IS-)c*=M$PH|}c{z781qW864RpXW+CcTC9Vt9vvTi9P8>kIi6I*=O z76tBHs@QyAVUyI`bO7>vq)|%3#X+ScSnYQMP0rQP=3%R%-`1|TV*^=@R~xVozs=a} zHMqeI*I*bo{p-ve{eK$JinDb}3%=EYuVfH91kMR9WmFM5CwWY9sPEmyp{@&QupJ7H zUFEfHK~W%z^ZB^3dpc4N<+v#|wIce-=xU0l2L#SBl2+wrqKmYJAGfDMUgl1Ou2)C4 zO6a+rGAv_O^eegkh&>$eUx5ERWkgt3eAiJK)N!;cA`iCV7CP7(*rV)mj&lP}C4&-~ z^>%dci1yLuoc0H5ZN6?$Y;b_{dMcRWVeU$0M-Pr5izeihDT`2L&;!V(Dp#TUR?(UE z1k+*IIh*T_)4uhHbCpr?OJQ$LtZyyVm(HE`T}uc9$qWODS?fw!8+Ao}apqU*x>7tk zcN)c^u7qO79f7?duqLWo=p%T{8wl`%mZsAVrir9Ee!t!@8j)u<-oQM2C=G0v9;`an zSCds5N~jksNK(}cJn}tny%$q2+h*)&E5Q~i8=O?Ca_3VFERrCXC3uY6k4s^rtD=MZ zshg-$HZndrHxH_JZpNG9>j2#kXNwQRV5^58Vs6JqLQqAm&|k6dv4sv8uiMuA+#ot! zd{C38V51+5D)ZBV`xDc`J7ukD4+SY^?Mszj@&}c<9`(SUL9b+A-6m0C+d!dIP3V@T zwtWjB-@!6jZW$8;KlB-ePugn5e98%mTOUc|B# z2Y`FSLL7N8>qh}3AYrSaURy@Fhqe>u`!E#8@-LF2YeQfwFDd}QREzu`Ms!c!-{?Gd zE%+!gEF{Su$^%o2FD{E57~U@8*i{ZZX_r+!)$c(eWJS*Unq#URRG~#A2|gNr4eTqNCk6sO#k_O1{MwI zll8=0_&KqdRE>R-)&4K2=*l4Bi-{PxFsew6z|aLNDZEXLk9NMZQJme^L?@+V!2^35?I<` z!9g>oB)kd861;@SiHHWXO&WDdwLe4otA|%m6hP2)D`?ty!%*`BL2+i;4{pB0427>i zEav?^Sd&{!{!)Xn!wX{^j{L-}HEE~6U6iR3*Nn_{@72HCddufh5# z9a{2(3n%9Mkkg;GN;*5>ya(W6uYWrwUV(Bk_T@@NRImU%NSZ+Xk8_S<3`T{pYzh*O zbi$YmL^Zo?`1QlMTY#S(`_GW#AD9n7mpV))BwH;8UMxzq;8thBZ3l^3d+g5~TMW2) zT$T?5F64{tllE%>w!G0s!%0{~r1ul!{(v ze2Qvch$`pDtG)L?ifThNuCjDNyxOZMHYlHjLaaw3r?M*{?b*h7`W5UEb&unF6+P?j zn&ZmtiRT>rziUp8v!^Ezg*(sJ{?DinNHU-qLoq4$hBkX(Oh=ml^Mcr^5O9wT%-4v! zU2xLvE-csG?kZBAT*&co|_1^jA2p1Eedyfy{KTLFXkOrOPl9AHJ0Y>T!K2;NrT9$)mv*&F&soVJI7h;wXMaG zvHq=%f5i}hBIaV<6#58l7$iqQ?rc1F<7e3uixG#!17R$7m_qQ`o<1BSX$t2rvMy|5 zdSo4mc_zY2UPX59?Gj3>ia#EA2QUY9?G9cCAs^|^ALUKnRr~(L%r)P<+UE_ zy=RH_81d~On#Iz|BszO~>vSdtyuD$qJ!mlkVNoAP!}OW}uNK+^I?x5jH8&v!FPivL zjlP03XI%4mBdG9S5a@zCU(!ZUH`3v74A?UDa7h<&k}sSd*gsI7<}3I8(P_RhPk0Je zeaGjzOT%-Zc?|C2jrhoxJRZhQX!cma`9ag}zK1%BzUos!DmgwD6r^QmiU9lmfz{5Q zc>YXwHaQg*?C|0E53;jFrrRVO4YgP>wsDSim;hT((%}67Yhrn5tEq1nX@6;V53)kZ zjSMBQ^*39ZP!rVuLL)({`kR3Xsm^s@ zWF#WEEaf76hbW;Gv(9LdnZPW?480&;g@jVjQnPq7m27FHe~SF}i9%l}5T8NeFA`aw z#n*{sp!f_?<+!o1LFXB<46Aj@z{C$1$!scbt1!zmQp*Enj2OLU(8pZFhnr^L>n8Zj zOZj-E^4yf*ag*^qmkTP&`RKIClO~jpyB#1SQPJedZt*3C6yFD_oH%y68Ph-2{5}bN z35q^hGi4e+Uo#0s8z%_flIjFO1IQH+-(N96I5NkNk39&WHIp_c!chtWCS^jPgiPY! z^6^cZZuyW)9DegOpAX-)q7x29xg3OyD%*S}CLS_*Y{K_r;wS$I6p4)DP_B}@XP@5Z z=E*|ZgykP0?80~$Jf-*;o#Zdr09k-@pCx?!#)#kMOEw6A8;ccx5|V&E;pt+HB_Rp1 ziq6IIOwXb``jn9epS_{ybBn~=F!Q53ZNJPTOKck;3uzm$gd{BtpBU?iy`LU0GvSzqPVSC-)ewHS1~O>Ub;UpMCt5xjY$jVr9szj2;T^ zWM$htZ@?RVjv%p-z@Hhd@k4(64?uzzJP}jJ<3b&9rLnv6-IWFE2Qu{8$X;uwfMfJw#iHNbyEa@p}-H2IUWL@(kiHKO)M zk{13GjQ-}h{O^b@N$NvR5Nkl|QX?>xSiC;D5yy*#HgeS>9gT^in$=fu+qWn!Qtn+&_twp9?yuqA?%MOqXDy+Coc+HK72m0a>1?9 zsxC!Sth+?<`Q90p*1;m)llFIoWNtyq4aR$?R>i)2;a7xQrLXyVL&41k*?PnLjHH0l znK9;7qiJ&-#hBV7J=7vWCX3lGydodJMCCo)U0o89x|Z@pfB>BRC0K zP-^yS>~}Eg#NURaepE>8FPQ$Ora;5MV9gxs<*5t{`PDGpKkyg|!25 zInUSx^hmteV@8A)ygak5Kos&dbtWk^2QBymc6%jbG%lDDt7>#>V%s_yO^zeFaU(Y< zbv07RWi3db;!3WDpdo?U)kK5~K{2WGOeLz?C`A@5$NS4@lOalt#w4^+3x1B(vBW;K zYJQEcrLI;AjQiVE%OU`Onp{Lb9NRd~{gmAHxXxf@JaF~Dq^igk4mKt`JPsj2EkMUu zG#maFSG?EwQ*>PSXj4x}#0W_~@lKL0%$xnb~(b;B= zpQpxo3LF#)+_T2R%3L;|`}s*~Fu%EsYn5$VK#KZK*>t!5J?}kwsr45c&`|2M8z%3G|m#;p6eI@(k#(QRYgh4Cag| zWnPGxF0alkOLur9&&@Z+Por8!wUe^M-oNM88yv2OStHXUrO&b-S(qm_DAM(~45f%NXUyrLPmO%4s zh8L+|S>~hAE6hNGw@j`|WT*K_>?HWXXs#tUDEtW-frW@Kr!$L|&S&175B`_D$+JC? z8S&;wQaa5>z}NLNZWjf;kCrC$qY-YhGvZ(c{W&2$L@yg*7HMJu895}$h?`LlGnSDb z#OK#s9W#Ukl%Xb=4%k%LaV(sTAkjXfC_@%mNy^AksvshP2^6B;$jvM$HIt`cBQ1Ca z;Wed7B8ZRZK9&hGwY;{+Yw99WC@g;XLMBcEMRg@@BYV)wvDr|D8&t)(rVM~&i{jb% zihw~!HQBZCSx@LtY#@YcyaA05Tb>`{%u+X}p2L!+7f+>--*H*o#yljlkmQ<>mA{jr zBZ%sdjr@qSu@oKXADN6@LEEX@lOh8xxIHZeaMk%=yn>d_GE7v2YmhLU89A;1i&=9j zicZ3TLHtyNNYuE2^GV%=6kXC*il-WAP!$WQ#BPAfK}8D1Fs(KSRwO0Q5VoxACCb^S zNw7qoqfIt(e2JZQM)w-`)GbQUDgIkh_i7wR1}7vUvsWTAu^qq|4yqo4sGuzwQYtmB zj8|6V%F1k8b@*rq@c&OdBa6sB=B8@FFM+`H%?0i(+#GvIN}{wXk3p2y1*h0Rq~dCd zrdse4q$8QpeeR#M`|zQUM3HAvD~i}Tko^A@-;?#E0h{{J+j{>f0(xrcrn?H>I9b^ihLLd6r znEciAtP$T5yCS)dNtN3IMZl1!raYD|cNT{}98H4_T&Q@YWIi@+U=q19i(W(Z^;|xO ziFF3(!}uPvdfi7;?lcUCfE_y{X_}G{Nhb zGC%)46y-u9y?`e`TGbZ_!Jc1N@lrDC<8vXcwT;_Iq8cP%G1eyzjqpeNjBCaJgb)_% z41(GcW6v;RAG>&hUuP+`usDi-qIiCCQ*FsR1{B9W4{GAi8MqDvEMLrE$xCNV)Iy*7NjpnLWnWfg1g!l znmHzFKc@o`9K`<&f;K4l-z$o%g(ya$Bz$@I*xj@zDOH!G;ZgpAR6@TRO)}fhUrixY zVPrO^5GjmRvA=_~X+Hj-q=zIS~BjFh^Gbll-)#*4s$>d{}VX+jv0O*06#q^Kb)4689m3`^dqkS z5V<6BqCr;twHyj={7ldkv8HfxNMy8gj)&1!2<4IVFOu)oE3v1Sq`p<_nOsdW+SGi| zkDx1|yoFyO2J(hkKa*Mx^c^EO(*G4?2SY%o@90h&jD|#t@g>I)zFI7Du(i||Uz({X z)rGVfG}NJp=)zT$zlOiyRj%Vh$5p=v2wm#Dhm@R-p1>+6$To}Kat0hC;{C`Z(AZlh z4;`83Hb|a(4ENMscf7g-Nmm)%PwvbLfN2b?GaEUjqeV(E+^8p@poQt`cl8vTvL17MVUF?~h1H z3N58wKU4%)nwY2>fP{L}l))&Z$t)!N3R&dQ?Dcdr>MwP#&{73Fr@u~|NjWBe! z*Sz2!BL_^xg2ib3F;g~9^8jX&<*`s5k2#aJlFg!e3(0nDN^{$3Nn;(Q7aDit>YuJg zu}Q$?Png|lKW%pO9f+T>c6TDXI3D%Ln?&o(Hc5l?7#izL>pYX8v15ry>j4BaP3w`! zVUQ>cyS$#xq^(DSJaw(-Ap-=8oFexi7V#<%RqV`1^h&c>$+bP1^`p%_S@I#@R*^3> z7|IQ*8YvoF3bLEe+2y$F9=|=O>rPRXczH%Oh9X70P*RMj(86?uit~0+AfYH%K}#G&a_FYRfg&g#`4V$?}aWkE2$6lwX4O zYN(Vrbml4Q7OZR1GOQH6b*nz=5yf8S1T|nOjs2tzz{&?=!EQB{yxz#<^>osA8ST8C zol{P`*PdF2yd(;JgN~!3zbEoZSy9f-ER$;>3($g(9bX^yCu&bwh^2ePA~Y8Pq?oR$ zv02C1Rbt76(>K0|dOt=ek(k?F#+Vb^{f#!s*7%B!8D0L_sas*>C`vQ2wDkC})Q;9l zQeug7aSX_55{MMgfx&_Q4c{w)GtKJ5`(}1 zB|Z|+hvQI-tV?N)WckA;zkP-QFfLI*`A?a0LrMq|7lyI#dnqp3sKE8bOEO8h)k9fJ z7(<5rGTmG$h1RJBNtsNIdT9W)`(GeNOjO;IIw@~5U@OXjgu{R2 zq-Vtg8aSX<24FX-k%H7SE!ZRjIzg2T-U*z?Bhc3ma3C7~Cr$Culhm<{l%!cmLyB3* z-pO!Cx7dMEK`x2>C)k8aX3s(`ydo#CvfB>-MHZGi?&B{f)f`prRMR0&&P$0=5BnGiD{NRez(EE@3pOoodJ3}_TmSm?#sA< zwu;>!(|d0ndtHuXqo)crKKaw-6&ns75Us=-*+x3-RSP}?c7pvRPK^nTZ#zIBn1+1% z^uj;9p{@A78NNXxlwLhym}slhdPB{uwi>#vIYSkitzkkF;?UGLk$Y0c=J*2@jY$C+ zo=sAOeX%eCCF!j=jxHw`hn$H|YxD3RH8h-4r)c-?ht{J7$v%zu!c6Vn?F{$_D?_Z? z8F~Z4ab{H{Ya^AW-8-5EQW;?|o?&q&a-b=v7W@VNCyA;=!4TJ;j6JNDWdcL0SwtEM zI76&aZ7pJrn#E3FK!)HCF(jL~2}>=N9s4jB$1w~%*hxgAT_Nv!oThtyHPb>o8|j4$ zrvT_9oCH0!g$5O5g4ELfZ>;()C6@`YFyYdwo(2++nS!<8^T1eY?&M0D!ZpnzPuuBo zl3cY;o>v@my^bnN&3c&wxnvVjYwG}jk5(|He#Sw>jfqvXvR~tKKV21tQ=6!0Y_ur| zXho?pEP6an!=UZ0MBzNgv(Pp%M*|1c;=d6f=c1xT#sRd#;9*Rrj6v_1+XIUsp(&-7 zogk=*5G9t=tvsgM)pR5uvXq)t=@@QptXO8Ir#u7O8e zY%;_VVLNs^9?ZUq{R(de0gGXv1G3SA^ld#0oG5-TTCD{ex$Q#H#F8=DaX{-lLy|lT z6q3OJ@sdR|x0nVp!;)m&ExO%tG`I#O94Bj71d~-8{YhA3_>qNnoLvoyJqFg0ewAYO za*XWp2D=8#qO7@_6vea$%k#5860X-*)J)Spc$slniY%$*sFI|C93z-{nK=K%6 z!%LSF=|oDNuP<3a2V;{{uPP5XQQik(Uv?+%a$-xG)<3wUH za3RhTqf&UbSX6j8uzGkTOd7P{XT+~})r|q4c8wR@C1hrSsVuYY719Yyx}mbYke zsi?K0$@-z2qHC}7M%XxkGxiA|-UyqMET^oMG~^2nls>Ho7QnHJg-)W__v;uY_3K4GClGb31WUTru#5Auv&bit#d+WNr-7A7bbDe5w`Ct)C1$cb(=|$ zLE-19U=g{fz62;(R%j#=F|Zkr;q$NXM&6`FFIh`bw!zX#br6vffN5eH*?i2>z;M=q zTdjCB>TcP!g320l*$6da8*FBjWyDLUNJ^lQA}Qwo)s75iaXS2u#ho7ES)XJU@~DL& z(bO`JzYV;HO7y;Em@??i$WM-$2KkBUyum)qbC}!9!Ief}Iza;;)5KYcHLds%7QUSN zF$sv^niVWMw32$9>WSt4uG+%?!U8c4MeZWzub~b60i76soJzOtmFIy+_m`0ypPRms zyoJ8Ehwq(kpzlu|sBMSuHX{pY>DV-(=ixN@p*TAgU(3Na?mI26+je^w?$u~h8J7BU z=mc~ALZ)rQatA)z6mdpKJfNL&QLj@*bQ@x}AZFwmcoDW#BA;Oy_$FJ8-3%#ct~^1H z6lY`tg0ytTyPUj#Y9i$eRLa0tI(C)GIE{g7xnJzK1xq+WG7H^6vwYqJSYI-q+a{O9 z2sa^itrA2-Gj&f5W-nEbv_s#B?fmhRQQ&K-EjTliLafs_-mObQzX^RDt(ja8TL-VA z;^=W4c?rK8`dnF+gA=|&Lm{S0*eOh?z2W9C?4z`5v(IwFX&hMvVUIJU1iVj3U5GpbPK~NSr#6MC2f;eoSfsab|y~|1v z86b!Ze5IB@8)4gnZ|C~TKr}oYAZ`x21YRh{Sv81?gE(>8Z{#Z)tjP??@|@U4ux1t%71IqG`$*`T$E^7& zhgg+OcC{zg#HS)_qH8AOvm&WjGsp?~7$d=&roIZ+6bwse7}0E^6da1L*?~2CBN25m zdWul9E6$hrAO!mNko^9~A^!eHGCu}EpGW7%5h?f)5?b)1!3=^gN03KYw!x2(&?Y~M z5aL6U(8PzLhgq5W9_=z)TZI)2e$%QRK;6 z2=b@dKBJe+&lEwhSKk$F8exMNf-4q{XNPoK1Bs!lG2pV4Y}5LzlhLB~XnlJXI`L;G5qP~xL6p$2yUjf?UnSSk|oprx$vbG#Jpq<_2V-yZt6pZ+z` zze99rKP=!lz<)~O1wJw6^H(RPXc?E zlV7vw6$8O6m0w--+8L4@avo_r;Zt&z;t}U!8!NCx@>)<3zCcF;j?y5`(5b#O!&78v zW&{GX^oTiqiBYNa<)qM!^JA~!6bSmZ3r!0^$;0Oxqe<`#$ZL)YqKi$(Iqr1tiFHr9 zFF-Q7H$kXG59Nu?6PyR2GK=V;o`MyEnG_4LsbG9zZaJL|za};or@!MsGi+Rn%|Tf> zv)zq&xc#Py1e?`xbt7kqplpaM$UNdVG!v~8Yv9$OU0ZA^r(;bZArL2waJb5et@(D~ zeB=LUf(u2)s4G`k!0HTGiq~U|CeBJSBbyImifz8~gZU=h&MmKW-jey6l+;pno~jeR z(>Ez0r*xGwi%OsOKZ%-jB_tJ?RKYf(W~U`}kfO9c^03`jW%o=}lp34*wp-a`Q+IV& zKBhxHzp<R-DlO%4P`oa$rUl;1nm#od(Wo$9mQlyy!u z*iHGhioidr>VY$qYE|8FhO$moce$0_s@mLLc{u}tYckXoZspSq^%;OO)IHsld)uk+ zcT={vLts-owY8fP$yDpQDa$g|#od*)nd)QRm6tjqu%V-Rzgzi7M|EpAWmT5?iCc+k z0Dr8h&2D9vrrvX=5^SeYE3K!UzZn^x zt}b&sc6X&*>dv4z_$JdF`)nc$d9^;W>;&6=C;S|FU84qV$_&+}Y_s9*Fn`ea3A_57 zosgs`VTZcaLD&Hx;8ZtcD$AYfi%zABJ(R5{QiUHl34w2&q;5S$*?E%Mbc*s|C$+VcvZ9k(+evBZguwsl ztiEu%@@{AKQJb=-v--P}m8B=ED^FHhPDbD(r>JXBQL0Z-zdA|zyo>s2S0!|+`a@Uc zH>awhGn98uRkxqv*m$a1bB3e&RP~`V9KSz}AS+K(U+g28uuw{ zdMOW|slMJzsXkNPagO7?Gu1cGaU4EVU2=}&=^g}m`YiRST;=Vvs7*W1QXkD#BDwS& z%~ikdt?bHGk!AR7^{um&=g+2m!sj649Y1qDb-wa)Pj&Zs%CGYPcqUJMrjPPQp1S`6 zC6=fD;{xSSp1S%1W&ODb{Nh~ooePyc=c@1KI~Mg)f1mGor58Pa(@XuhxAIVLs_**V zh`91R_0fD~$9a_PG8+HhIbVIFud?ra%5d=o>id0_S1(Za^;MqigTM`a)L_1Hpby1v z%~#i70jA{QTghYD0nJ)l1b!3LJmu=dTOY75yCd_EZ0Qnd7B?^t`j5 zdf-CGqRZ4hmpNXzjGo`>j{-mJuWm0??z$Y0zrI}kb&*nkxjO$!<%i3oP)62Wq3$bI z4qTzWS*$E7x=2wr6{)|uQaM;eQPy0k{_ZN}>nqjeS1I=mP*)C6ULHWtUk^ZD_ZF*f zd6bumDdkti0Q|AUv%f@nYM{E+<9KnPx@)jw6P^bvtpn93Jjx@wdY?zB)9HDOu5R)u zK@UAQ3_{Y62C2XGDt8T5A1+n?G+2GPRQZsfYrP2k%uC>JO4SuZm4}C@D~BrI3{lq( zbsQR^etNBA+12XCYn3%u6MWaz>TiZB_YI}zSFb_Rr>;|IZkQ0SzS6od14$r*Ns;=0s8&r?>L#iWPgUN$6{UWAs~Wyd`D!vAzniS?tx)fuqW*oB^7<5Y&n)G; zDFD|`RS!&6c2A`w4^C5SrztC@(et)x>iyG|yQkCh*6C{PZOX#iD6f}pL!z&4Qy;BR z9;hI|YZdCw3guINezx*b_>ue6#yga+d}{6O%CBdDY#y4aes!l(H&cE2PUYZC^~?D- z<;huU!yM(kS;*_b+tt?Fl@+(EEq5p_w*&B(JJc`kP(HbX0B_G$e>+$Gc(&R!M_qWQ z`uQB?&v#PFop&PQ(mCpW99?C3-?EG z-+=plxZi)bO%34wDeiCKe)m1d_g^oqz=uPpn1_4B`&-RVyQx86}D+ntMaboD5t z|Cq8C@@@WTaaS&p_afZ4&@+{n`WiT3KtFfS;WzuI`uy(vi_X93{66Ph?B~}D=U&+N z{62m2B$(4vKhY8os7V<(WVm8Sn+ay%4F?(UWxL~A#Wu@jJ3Xsi<`SEt^Z~T14BLrB zN=deBp?$y!8Q0wBSnQ0bi!uW39ADTvAo`UZ?K^N2X@<%)C3G*MG)3860*-;%?qMCW zy9~^B6=i4Y9fxFhK^PvnoGO{-Vw9^~neAFcX%;$)vfTma!0bNu`5m&|fFQ+{*_lNh zxjidnnlt5{+W2-7%&5QZj-g&|UAspp#_2{I`{E^6rnwDe6JE@AUKw!=UE~a?BdLjl zJ316@D@qsmq{FD`nk8dl5jWtp7n*fsoVpR^d0A%6DAVy^ z0{$q@2geYn=@{a4m32#37SVs9jB}FFR|XGoZ`v0TWt4QZZ=^dl0oNOfJw2LHDGL_y!_%6xrvdH1db}w`e$nF_X(XYe^0?k)-9G;!I zYW3Uem>(Ks+A~oHx~Tqtm#{Dakl!MEG4U1h3pi((=pP{S5p?2WqC7UbVjRi%f21-h z6U(slaj8*=j%YrarYH3i=;}VAtHnfD5$BcJJr@yOEzBs&zBJHoV0NMXvkuvpf<%k5 zdm@HOWs(m*mu>tMS#h@9M=WzdXZE0ICL|5w_eY#gvh0^pHb0{L2$jExlv4nf!PmS?o> z7}B}8ozAWR7_?$q~Z&lQJ%O0Gn}5JV3}2B1${j zr_&w%O!;p<#{5Y(IW(t{wX#Zz z7x_CUCP4WW@Ckq#W=Y)A{1erMggL$fkmfp!b+i!vhWV&{Dh$1n$FlT#YQ_wW-gcYdL<{R9o9UjU<71wJ6tpf>4xI>xlG&3(k81s z0Y4|olleABXIn~%!s8anh)ZSr{@{}Jp)vXz9-{|gj4sP?X*5P-pq@?x^+Xz|$6QSV z^-zWIJWHm32TBn4OOZ%L5jYo4Ad+$>l$~RLF%9vhXe5m#=`JEUr?R&-Rk%YrFBKTN zUVa;R2SDr=IO}+T4P-Q(TTT>a#-`kRmqK$7p3@utT{j5g^|jik=2OI1Ris4SYx%_%&(Z*QJ3EPXoU` z4SYlz_{cQy8`HpVN&~++4SakW_$_JRlhVL%O#}Z`8u*kn@TqCw)6>8!(!l*`;4{;} zZ%+fCod$kq8u;8a@Of$Accp;`(!i_Iz)2dky&THjZGnI|lu#P@y=mZ!(!e8W;I#jv z?d4GJZwmy(p*)xdzBCP-yj2Lt`ObpYi_}dF!c@35TkZa9u9p^An(oQq1rdm5E7kn$IZ8i^`ZJ__Y%La1wlkgfB^g)8QVe z^0`IAlPPhkgfC0V$IozBY?H1#7*1Ma;x$j^lL%FmDu%z53|2x6|3w@4VuoMY2F`N> zWl$UV!yJET8~9@kAJzu`dxn>_fiGuxGSw=|pBO&2O?;ZaDs6g8QC??wd7Jpv44=>j zPPQuQuVngpo8i7T@!w)OcMLERqKPtcFdpy{Vc>>mel;oCWN_cV)WJ`E5rF3HWMan(#toi9Q8UOht zxLd-%OoE>y;p>v%y(L`8t$5&t622g*U6(TaqBfG4_RA{DOG=k`7W~9>Y2aa9Qa*!Z z{)T5b@Ks50VQ0HY*_@Qm z9WwsMN%7}1{Gv7-vrxwGniPMLgpW^xKfv&pOxw6o{KGPTB`N-YOZcTmM0tNw=8q|F z!bJy45Pdc#F6by8)F9!*lH#wE@Lwgt zKag-iU-5iCW;pPiRR1j!F63oAKG}^FU!#WwPxRL1X>6n~h6>q+pD{Qi=# zYs8CH#xndOMc^qOFj3|s(#7G^WPH(%ID7`fUsMX?Aby%7Oue3^v5nFN2D-ywgJ;D3_v9j%> zWtEJdOmAytd_5`t1{ptD&WQT`zxK`qOseWw{H<}p9dQGr5fl_RoF13~3@Qqvj>saQ zC@xJi(*v!uInyAhQQU(2zDMH{H7aVr-JoV%gGM7sTri1B#N=tBi5la>{a?%H;b8^?uLn>wQWr1B?#~ zxZWPE<7-;L2RptX^^x0mSlZJ)zqkEbpMPav$L}wCUz4=cg4E%Go$+R8o;P6o!_3a% z0YBRD1*vlaJGI6~2K+eV%>i$5&+ku73+&G|9=}MixXg!v{Y8#1NIl{99ZE@mU69gSuUOniDP`T; zMktS;b^74l6e>dJUkWa4j z+?NX94}R|CgM)?W>iu|nW+E)fbswkH>7^YB&$zOa(s&-ezoI$iufYO zw=d$qb$o{+{!igy4pd0O z%hWazUg3Dhd3Y~BF2ZLxKD<5udS8T=dBXdx`AG6a{`;PBo$ylmpC921BKDtld_=o_ zKkonJc*lAEv+fY3%$QW%Ob}QO`gbq&);zTh$6n;`s#njbq6m$*zxJ@@_KKEmaOnnah(+5mqqyP z&i;(VsQlOaIJCU#__4_o`S0%0ko`X$@3`*Z$Jd4%sQ<-$fa4w4BhE96`#IimJ;Jwp zCc;m4yt&<<1@`DR$2+c5xO*(*wbzB0%Fl8qKdar&oA#)W9Ei40&28t&`6l_Oi||r* zrbhU2&d$m0c9Qdqa$YXNuXnuT`h~k^LbUYp{G95^nX0OZVfC%FF7ST29Ee)MoO9TuWBeGKK66Z!U?cslGNf=-^T175YD2UhbX@*Z)I@ z>i@xI?*Cx@pDs&~lO5(nhbAu!P3!B!+>@ch64@aRxfjxDXC*zvYDtuaB!U%AFkO*| zR5-KgN~f0~wd1H$b0#X4IcJ;==daa~RG8F|*fRUhE2o%tnrWw6?&K|Q$IBD$B+H#{ zrTahF)iT(iM)F})X5BUwGVYRhdDkhdF<>DODs6a;m*!*SIcl` zXSj=QxYHW$B!@e#;m&5dEKzXghC37KGN+v`bL#0b_fonnX^JvuAzhj5%FKq-8{#ab z(}`L-F_KQyl156Gn}~BJolY#L6I<4r>BNVm#nMT$rhT_g6w00d>7;?v<%vSl!PA42 zY6mAy3{GlK(wZKeRFpJ#I%&Cd(go5Gu7`?QU_%EaBuL?P*?>7?_dE0dsyCOI3L#E^7_bkZx+ zNw-W7b=ge@U^4y^yJA7+yN)?0FLO@5s%msa)$GnM*kt5iu*vW+o2jZsRa7<9mDG?| zWS&fCI=@&e^~j{^&2>0nFHoyeK&=;$`7@nwYR$;B zoXgbbr?+HV=4SGZ`NHg!&ge68g^U~tGKHBkUDvj@v*fwHl#-deKG!JD$n0L3@iV*F zXw9`0s%yG;9GIM$TqZ|?N;wit8l0(5LK~hOE5yB$W-4mxbJ>zz6)Yea!?VdI?&DIz7h>_8Nnoh`w<1-2~Gr5+QrWQE{CDr=B%q4RjyrU^prs@xw=7_%#1m6GR?V`)~3d6eZDX^Gg}U- zhc=#&Z>-Hs%(cpKS^X@sk+$wPCVVA#d`I^ydU0@VQ^CKSshN|_Wb5mjY9x)B+?<+R zb0OcPO)JAAJ-9_aKay{d7c^U=(p4^D2OKdNkmF@FmmGz_YrdfqdQ(dMd+c+cVyOHZ*onA@$%hbubknHpV zI)CR%YntjySL$2E#Pd}g5k>Rv=%UfoXqAKnVD;j%hb#~E>oAyyXn1i2bV*d%z-Ytq8DUhD}Ai1Y{|9i zY@i1`b<$_By~_>mG**gU9#%X?ie-w&Mv+WqI3|iD)RF2snJZ)!FxkzIMIys|Mb!=6 zlzXt$ISOfFz@)5$wRz6Z57PTqr#6E**?RF)6+wzPA<6@9T1 zH6d?fb-oJhsELlu+HrriqssfJ@||2tSNJ(K-#A+i_c?X0b%6|GtE#%OoSUc{fQXNA~tPSk^V|&0b>KohT9f`~1{~tC2N+K@RN=vRsx^ zL2j|HNArg4oQ$j_^0Va~LT+}hK6x9Uuhh%DW{j>>3J``^gsd2hQu5%Tf~_vf$mbmKJGnE-jc2h7{wS!A@R{m*6dd|!o9lz$K&ivQ7vih2S>R-wv<0{X%a^P#p zV=vzTzAkuL^`)r&>9TqGsS2eiuaV93d9bq~?94Z=^0kn^6!Q9Y*xt@H#&;3fO~7w7 zu6BBYC->nd`P6&8y`AK~+JyH4f70w|ei(0!8kVAbbI8vEXIv*6*IiCFU+;X#Gp-Ab zhjHl-y=u`NeWsk?*ou$ypI|WL>mG2kg_I^aJI`O z;LH#GWPYAAdG$x@_HlgxJ89%|SR|i&Me=z7IP=MM`|+qZ`FvjD3-j3sJLU4&`+17- zFwbWh5A%Eh>@d&%^BtN4=6NCHndhUBALiNLx1{xE5!ddp!~V7pxa#`2m^aqz?`Lvt z1^f4Dg;JDrovA-Et3^5AXU+PK!)|=h&l8QSJpG)Ih$nqS9(zB32+n$sN4;Ffo(sH?<4Mo9k#oEF|`(rm+hW>9w^~l$G#1Dw#&aWKU_Ed zUw)GJ`|Lklm;PUVit}@<^mRXe|F`7lP#vUFbewX&~GA^mB>%4==l8H6v@wn;LMNzTtt|k7h#|Ic?EWuAAg^{`p^8Z zT^R3%$OH4!ADr^1iCryfpZ)wnM~k_B-MW@gep3a;IWF{D)IF#D8#I~*6sK%V2` zrO3E=IWjIrCDl4V*}i9j>m1_S_b%fao_-gB=l4VY0LZ@);qMq%JF{f-cHW2lfv}@F z(W2?A247zuYf*j>xVEVl!L%%*(KiLm8MlEWe{cr_1+n47n`2OoZ z(N9Lh&taH9$^GXFv|jSNOkU0L{qirt^%&vvuy8$jB1P@6|2z}nuSEDRJ)IdhPfnAE z-u@oO)jp3?dxEctct=P0A>bU>d2r4rXBgK!=(&CG&)JZtKMNxA3nTK2!P&2#jqqQC zbDnt*ob${V;2e)@ZmNk)jz=fR-p93(am_Q2N4pvi`*S($us`n!JM7O*;B2p}z!}#q z;FP~B!cRs2XFtqquu?QGj-y)~O*rRCy+)};<*U%II1kYNaTup#A%8YF=UcxXQ2Sg5 zay`KQuo8CIZ#P1`?6*6Cvp>kkS6tD0`975A-^jO%$oo%HsvX)N33>Ln1B|O%b7b@V zZ6f5^-}F7HmaxC+T%<+qu)k>^)}sEfzn$!8!kM4bBHVvsN$s#7-VJ%?$De--^Vt@W zUlQRfz?n~;H)B3OF?r2{emc$PlgBmYa|`P~VLtu$y@dJfkN(Dd?hMX+?hVfT=&@Bx zn4d!;oagIkry26h&xzp7&uPXrt`lYR`8gl*%+Jln!~EOiE;#$c8Q{#1KToJH>zv^8V0wJwwIw&*P|<`B?%^ zfBbnuefbCQhvx~|UZ0qq(4Q}i>yGW^zYj;_qCZ<;9oS3yjL$>zJw3(w9|(E2V_Afc zi110^Y_B{x{o(t4)|-59O|hTX!Vcwcj&T1yHku#$^CaZipZ&U5<;j!pn@Rj(`}*(0 z(RMjSHlODeu+R4W%6QnmYiWg2!uHjl)zqRpw(pkU^k-YxIaPGLKLa68fA%&W`lG)) zt|jzmEbP#q8gRC60i5|fJHqt`F16^6ab06v{h245kLwo5Gp?T-599LRmlDSH9PBWz z-yAUOS;0?s`9^8=cn z(`55`&O@GgKGk@b=N}pm^L!!fFwe>NrWCjP4X|@M{J9f$=#T&2p3t9P!Vc%pmtcqf zybR9#zaHVAf;0bH;QWc+L*?_iHTa?6+k>;c{O8`aKaiJ4>^S-_e?4JSqGr()Wo4}`op9cGk>m0~Wfc%w^r~m&1=li~0*A-yiEaT9*zQM9Y$jCvf^R2zGubb$fq?L!SOjG9LOf z6?WK9ro#^XDS$K2Cxf%SE(2$K-2%>bxd)u_K44t)HXre}L7wsc-gp@ApJ0daz6U#u zcg^*sGAY`=^hdwfM~nJ%wroECn?s)d>|%Tuk=+FIa2j^#&v4kGKYN2SZxg_o|HC7E zT7>KO3~AAL*&b6xqKAuYSeMn5+;9**lRjfeB?cfsjr8k~7P7@X}pIl_;M@Y)DJ z37qXd51f9U3(kCAVqEiZo@~C|uY)}Ed9U#>pO3%}+x-dHVY@#IPCu7{GoSy6@E#kA z(^9m(Xs0(g*GrSYnFswIDJ`0ZoNPW1Imj~)Cm9d(a60TT59h)T^KcP3{l6BR>%bep z8P`vZtN&x*|Ifg=4txlFI^=&1`&`$(3ceZU!#Ba%&))@Sp8sK7{X8G_ehzu&dGn1# zNJ^OJ?;2M-%<~S$!|}U2IQmU=;uV(VLqqA4)a+L&U_ZYIbJRXXIyuK)6T2l zwDT1>-{-Dr7os$ue4o2MIPLU`@U4w&dyPdtcY-|UjdF1PtrDOAeZYHz9{~G$t;)+! zhCJOdZiGDjSpZId9t5X9zmD)%BmB(>|BLak zT|S3A{anMYQ>*;Ji0dG5|C+p}{rLkp^SmO$KLF>v@riNu{ZhpB4dgj*^s(#X;k=>WBdA5~aNZb*>*<^~D#1Aq zjE?Z?2p<>WM@4uGIOmx&!0FGG;LO|g#x-8Ow(I-D?T}~Q9ycE5?HSl%-d=NUgFOBD#(3z@I=w_$O4yDY8xN0Dn}IWLy(4@GIQ`e}ef0HSArE{$_lG?F zKgxLMe+}#~pU1%t^Vtke|4#&`pQpk8mGJX?$kWf8jfZ~T0Xy__5$w><2f^v*;s}2c zycg!NSHL-MybFv80sd_;tgHm>c$b;u!*XS*B$dA4I+ zME-d2ZIQQ=z&Wp;4*OS$PrjdD0C~=pU?h~XFm5d9_Djj*kL}aVTbuV5}bavfU}=m22S~F!Pzc1MEIlNY?l?_ z%>TRK%>O^ZzlZ#P1@Skgxybpo1-CqW0yT1=k`SrdhPD{~tp?p7Z=5t4I=5r7@ z^Hv$*Bf-Cvx_v)6(75K4?|%=6Jm3FL182Pr;9tZ3EaRd5(;-j$=SJkOipVbkuRy>0 z9XQW#tN^F~e~a*sjjR8+Ag<5A4@13QgR73u^E%tf6Di7WmCf@?&P(9T^Bb_k>p}ks`F`+cc3-JXiu(DrY~Ifc`pOe2%DE1_w6D^Hv;W`H z*I_5mxE6!c&NJYg2VO9)`TVJD-k;YX&w1bjZRVHXrXnkY~KpjEC{gfE~u$1Urnk08T$o0jHm5!Tufa z^J2);&)bZLe%=i`^z#AOp`Skor=Kr?GyluL*^d7JXFvQ5objf%mmsBRyWA<8Z zYd_h3jSfD*xW>!x)7uSp`a^yQ>@eQp;EZ<+IOCcU;f>&o_f&A&zYCnlk>7w*{%^+D zlz6kq&xerb=eyV3L4>4)<7>QewZr37Jvh(Do(0bH=huKACboQg-DX_tY2Yv<4b>P|qBw13{(hYuY6ye)Mcz@&CAMTOO_lIi8b6qhN zoa1FC?2Lz<(;+_qeqImG{%{NI{0w$}26@h_&lwNr)upi0AN4MS9rlMez}YU}fHTj# z?I_Mk3G+V+ob&25aOUlDaMpVpIO|;qPWdOndm-L8!I|ej8`pMUD4TEh4$M*kRsY0%yCw4o*Lp!~T8n^KX!+pKA;hAt|As zJ&cF#+tYa1zFUCP&z->8KlcD<{w(!+G*4*y)dX<5#f5 zez+8z{{InnI3E86J`sLytd2<0d~!XuJvi5Uqrth}I~ttpy_>;le=#`yc^RDkzYk9P zYwzmZcmC6#KH&7f(zv!e+xIYVuB(rMod-n6_n!wM@>}ng_@nk|=P2+U(Y}q~>{qjl zYr8y%dgno&^TxHtHBJ27%}ubw`RxwaVZXW$oc{bAb{>L1&p@93{K0tW&)cv=fBp(P z^yeSo?B{FjE>1|%JaZn~2%PQPADrufoxy3p3Y>X3*toXK!?OAQFa`3=Lxb@!4+Yp^ z9!`cG=HWDO`gu7x{rnN^KLS5*gM6oX$+-H-_r*`aPN#VZcIf9)aQe9doPPcl_J0mP z{{?yGbAv%5BqhveFXP%S%;(m|wV%+>?ZD}0IXL|s4*QG6mT&jbkf)zhjfZ|_VTXR^ zVTXP;g453nz}LilayvNt`9t8Ge>U7hf{~*683{X^8`pduB%6&okg`mlt5?QTX#3DC84ufKCG4W<6&Ia!4BK&7T952KL_Xc>^&RdOTn3kx53$OKL@A% zUga+RZk_fd;_VCm6!-w+I_{WH`4x0mlrx`u!p_t3ns1jX$TOcu8V~b19d?+{V_}E+ zY=QhUsP_kur~em$)Bh`B=ULdf3G(!RvGLIVCt-*FKMySXd4z8`Se%xk z?ZW=BBRKnm{E)mWDo?&|giitIeAp1-b0d77ah(T#DVy&n^TByN@qBRA@$#30-v<7O zagB@bTc3yg!I1wm_;~P7VP|*DC+iOpAt`E~{by6->f;jGyr0`Zp8coNc-VjTf*tms zD%fHF$v}P#>YW2V8T=CPvEX-tj{|=Sd;<78;ETY&F|PUH`hSB80m^hkKKBM^dlkU< zfc)vkHP4e}^Lf4wobOwgfHU3|;1gkIgGzZKMeQ@*9l;szzQ)72sv*z**$8>YI~$zw z&I4zEn-BZXiLX8nmqMQX?JnbCf4dKM*xwez4*OdhIOnUkz}dc^f^$9m1vu@bhHCbc zenQ?0ob%f*;HJxVG=}viZCn2zlo17~^5yaS)5TmgCd`I+(1&(ttcmZI%KKYJM0 z@kKv30;ivqu*3Oq0_5rc6v)3Qb@=`;9rE;lj`7g{Q(=ex&xalQe?2(+?fu{!FOPsT z-p3<+DR?i;pML~r9#$CF{Qnwx_&em8hqZ=_kd!bF8yFAUaWms#JN5==Tm!-B=N_>C zlGN|pYfs42&k4puKc~PB+baV*Y_D2y`gsmG&(B{7&hzuvg0tOkGOl@BinwkE=ep=_ zaMkg7SPuI)VwPzSj3eS&GI*`3;O~9)2gA&qFV8w#!!FY?tl8xgMx8 z9_I5f$iIwwj|M*)JO|D^p8|di$M2RSL`X_Fez!HQIp+8sU|i$k_}vAZ<97r&$Kxb$&cla;)6TKry^x2~!I_71 zjH~~z$>#gj<&b9{?ld0eVG-;w50AhO^YA!0{eKyp{eO9ce*wIFdmN6QLw{tItF%_&n!6oY>4nP!8tEo z0nT>00i5yPW?bWa8}TlJJmdY9@i5-s!Vcs8J?t>vx4`M=N8rrE7vPlNXg_gUinedB zH9FdBH*nfH1f28IVc?vXjsd5gV{rNNFVSi4I7G){wC;M{`p*bE zH6C{U19qlDo^ch7hjE<*JB;fL*kN23fwTYrD8hdl;rB)OFCzR|5$}uO9525G zR~?_{wW{Qal(l4IyQ~-C-v#GY9m-dN^L@QubJ3#m;=Kz@SwPSiCG0<28rL23xkE%g4bJa*7y-`wj0gXo`0VpP6`bXZ7cf_!iaQya+@ZJ&LH^Tkzk!kgew*vCa=U&EDp7}p0B0mxG%zq{#-vB-w{cR37 z^K%|}U&vnxPWjuw*+1`z@CPIO@d$s`xaNn)mvj_c+t0&qM&w7YBSKQt zPp+@#gMXv1eUp@*8CU=R3cks@@PbUqhbb^bZkv{yx+py9-zUNki16_dKE=5D&wSQGp82nbJmZ}kk$(%E{eS!QJK|#hKghWH^EcW2 zxVRz0pN;UI8+6!bKFf`V_Aia_rz3pb4NKeK3Y_h=U4-ux;d>ZYKRK>PK%VU~bfd(d zupP%359i5Akf(kB`x+`w`~LSe)IQsN7VJ}g9{2_rcW0yC_r;cP-|N7~fZqW-w0}Q1 z?Y{=jcHeqqF)AgT-x|QTgrCnu`19a=AN^8<|23)B&D&M-c%_YpuLADBZ(Y}-S~Slp zfBvTWES0(~;Qn{fo(}lC_I$J&kmBw5_kq&`z9JDzy&Z7h-YKk$)SruNyv+-^|J~vr z2i*UD@DNSB6z`A!o!4Um?thQ=63p{j?@QLMUkBX(URpoQt19n*?`j&xk#hh0c6*sU z?~muBBK*RD``>@MJi;Fdxc?o{cLF|kZE;FUx%D6KkLPa&e3L{xwIbklTVv|ufIn>V zYnY<9-)7wZE{W&6TmR`3$X6R57;rt;qh)Bo{rQfO0Wa5zq>K-^fB$oMz-L>3Iy&I~ z_jKz6?tiEDq=3&iJLd+x%=k3{-`4oe0iR?2by2{#O2SI{-yiY$^zHaeAn$)i{N;fA z-{X8E;2+q$!QUUz{=CwTPyY(!{qNp=9q{YS&y9K}54@lL_sh2lxc`0M0Rhk1I2;mi z|2x*B0`7mu`oMtu-x-<`aR2*!bpiLkPdqE&cUxR%2HgMN+hqawzYlkF!2R#FZ*Kis z^E1&NYZbl0CxJhVd5Q9mf>YkyO{obz4<8mqQzE&QRA*MUE}=Fy_|t;@RSdTm0Bayx{$=lY%YT9jK?cF*-aV*C`EU@+@s{*0PxMgr-5${ zelhqK;Mag}3C{CievZ+XAA@{v$S(u$1O6xQt--$p|1S92)|vF>ZNUBcE9KjQ?+p3x zflmhC4*Y2FzTg*w_XB?ce0%Wa;5&eCj&((U@X_D{!2P_gaqS5HG~{;ze-V5j_$S~y zgMS6S3-}N_{%gIvg6{>s8~9A{-NCN|9|V3Y_#WW?_rA5>GVnJbp9bF$=OxO){qKFN z9sB5#d)1%!P(B3mV_~NP{A6&ymeH4Q2OkRgKY|Yf?}7Q&ZZUSR4g@V5sX?{bLN+yMdTVQH zP&*xznL2J*I4m(#S2e1Z*pe<_ejXrW2|fCXN|5DO1@#T~~EX)nPyHwsLq}v5wRY{Z)qBE&VP+u4$_O z_93rAxmzKZ)L)H|Bl@~l(D~O-bC#!k6^?(Lr#t5JS~-Nv(jhs_N6gfA+f3xJ5V%!a zl`e10Wn1d!W^%{R%GM`iq&8R@w6Ac6%QV&4+T7Hd)5V5g2Uv*hiyR)|+f_Rbi=AxD zo>;<4=>zSzGF7@kzon^;?%Zymn`=I|!V6bit7)W*P`c+y(6XJ6Dz$EXD6QZ1qOV=F z)A72aGuz^%vL)B5YbTl5>Vm3$Pi&v;2DcxlOFOVig-&i1Ne=6LgzcvNBCRee2@g%B zvRLd|sI}@!DtRpHW*mXesw>b|sIHU6p)A$i@{p5N-23W=)mMqtQq+s|U5j)eS>~4M zf|VY1$lCOo$kujCHa#RWHdB|cug~NgTXV9Y95+0Znb9~aQ!{5ydT7Sh{+WCs*OJvY zViRg5DOKu=&XDa#2UV8e6>4K-glyW$xf86@XOYG}kLONYDGux7pgKfUhT3u5mi;J47s!|Ke z_D=Eds5xAk{hKQ9qG$ZiMjlLA?IFuKhV}+7u?m_fBj?*RaTl@vkJV~(L-A5)R%3ns zxLo~Q=ApHaZ7IxZHWVHy+XqHXYaw09XG8eJcP4idO-AXU`lgy}{h-lPp*qAOtFO7zLjU!13EkCq7rq01z4xiR`u>_Z>(~3v zyQ{CVTBxtTGt*st9qXaK{;olH^+%e%U8YI3iLo75L;pb2XaD_3Y;?DN{p_n2*1z*= z)PGco`ac$X-K~FS3H{wxqkezApYiLxOWm!1p4HFx*YCgUuKvX(^#8dU@!x3r^uJ;? z>c6W*{rb7A?)vZVe_;IjIi&9Dzfhun{oGA=_5WBxUq5HkU44I_2J@$%d+4tIN~@px zdcSvf_1Dr1s8XmuVm0)8meBu1>~*((e?JQA*Uu_;SASrM`uAN8|A&}9`;UH4On2+| z_oq@{zt5q&`jboeKYBIv{rx$tziKt~n@hxh@M`GKE8)N1H`d+w&o7~`_m_27|LPL@ zdLLJ)-&gvS;FR5NA&-)>(DZ+(FQWK-dZhOhsm+ld*#2QVfctyd+Jx)-|ch|J#Uttzylf{SL948F^2m zh_(B#kG1HWr21k0r%>0a;8^7Z?E^-HTi%>S#KCW>9b;=t zL;uf@_`haD^MYn1^#7xX|2_4? zW>VVR|73rNJofebxOR}nQ^_#@cUa=P8h=Nbei;A5rq6q1)fHdA59>P959gnh-7MZ! z{8x#-{G~?9IAB>?BK|`~z}N5N|FzX0#{Zt(l-^bR8_NC68cP`eO4IKu{)0_FjDMuv z{N7dk2Z(;D_%Gh7xc!eP5&tPxe;EJHb~Alf@$2W?OT|CT^m$LUwvW%hZ~xCD@i&`( zSMmQV690tW#ql4F`hENFvQFv#e?&j^w(H|ho#}_|KhyNB853jj-^cCuKZ)p{ZThxr z*AMgeXVVY+&q7n&EEpHz{IOj0OZA_geG(5RB)e3bcraZyU%$^^A9wb6>yD--ZH8<*#8!r`~1$4mSPJ|BFrkM~RvC|6%(VBKl=JB`RV5Lj5MuFBSj7 z?Kg~4TWZxoNSpOkbe>i@xwhtktP2calzIW@#5&xH(zJBK{OX&YK z5&zFN|JzJi_p@a4_4~ZOVfBae*FC!?id`N52FiyGHI^{`$4uXrF^Mtx@8kCSCeshw zZ@KB;?*+|582_=NU#k5M)QfFW{B@8vc|2FPF#g-E{xE+3`=9nUMl_QWxRd%g;{U?k z6M6lPK$g(|4`#_0^ncS0O3$BP*oV}*YQKF=Kg@s6J(7w~k@Yo8 z=>G`OFO~m+_C*Wy|5Py)`ri<#|6=PQA|{Vx6AeZwUFu>CJK{jU5UDEg(^ zf4N@ll|ujL%k$9xDOP_`{=e_kyY&CknE$5VmH!t+{C7Xr8~A^23IBf{^WQ$C*j4}e zXUu=o|Eu-iBpb!YufIq9@7q67JkS)|WPgEdVgCEcc~+GS+kf6jfl~Z+q~uv~{o9Cs zsrH{)qW+7-NLc^mNc~IqOKR*Y{_&Ce7nZ31k`ndLkJR5bI;sD%F7khRr2ds9>c6Z+ z{r5%cpSFKezkW9?OE~}D6RCfs-Yg-7?SF+l5A*+4r2d5mIQ5jiKV%8(|9zzXv#tJR zrraj`t7Hr7uh=;0f8qM!_hS>quGSCdn||1T-N`C3QRw3MF< z(EnGg{&4;(OmOO{uGSB`$N^ew3gbV=^t+1xG}8~`UugRNI&|Vq@%TMi^h?D*Z>Qq; ze_SH|=dJ!Q{$&%p7k_{0P+C*y|A(f(tc(7)o#>bH|3+K*v;W;x!v7Pk{xJVb^kk8g zy^>w~|0FAk&k8?@==YeKs5G~0CJ%i7@pf9}Iz_4Y=b8Qk@;ZxD*Qqqh>!sqKvvX3z z)%v28HrZ~K&DY;s9&35W>JQtmqCR>3pk&wn-}~)BU9X7$roU?!<1Z`vrTp(XD5;6p zuWy%1L;sgr{h|L08WRuuct>p?`oF(?sfe<${hv1dedPQai`V!2siq%v)6{ngiPno% z)b8jvo*pXtE`ASe-&=2%k}|9z0d2CoOSUln#z_53A4&9X3pz~jG?kCkKjFHh{@yl$ zY5n((F!==c3yqS!akx8da4OkcVuqCe8~&kppx|9;xWx3moHr9b zUo85i#?KPVpC9LKV&DPUyuRljN9teN&H6ux)W5t${STF>e>=HOsVyDGzpR_}Z!P+z z;$K;!ejO5F{9~;CGI@<R7mp~u&5ms%JM?8K|&HT|!L$5Q$?KAcEiAn~z;`u`8W C>tjU# diff --git a/src/main.cpp b/src/main.cpp index baaf79d..7eb5154 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,5 +1,6 @@ #include "./utils/utils.h" #include "./numerics/numerics.h" +#include "./decomp/decomp.h" #include "./core/omp_config.h" #include @@ -7,566 +8,18 @@ -#define CHECK(cond, msg) \ - do { if (!(cond)) throw std::runtime_error(msg); } while (0) - -template -void expect_throw(F&& f, const char* msg_if_no_throw) { - try { - f(); - throw std::runtime_error(msg_if_no_throw); - } catch (const std::exception&) { - // ok: an exception was thrown - } -} - int main(int argc, char const *argv[]) { + utils::Md A; + decomp::LUd lu; + lu.factor(A); + // Single-level, 16 threads, runtime may adjust - omp_configure(/*max_levels=*/1, /*dynamic=*/true, /*threads_per_level=*/{16}); + //omp_configure(/*max_levels=*/1, /*dynamic=*/true, /*threads_per_level=*/{16}); - using utils::Vi; - using utils::Vf; - using utils::Vd; - using utils::Mi; - using utils::Mf; - using utils::Md; - - // ---------------- Equality / Inequality ---------------- - { - Vf a(3, 1.0f); // [1,1,1] - Vf b(3, 1.0f); // [1,1,1] - Vf c(3, 2.0f); // [2,2,2] - - CHECK(a == b, "a should equal b"); - CHECK(!(a != b), "a should not be != b"); - CHECK(a != c, "a should not equal c"); - - // mutate one element - a[1] = 5.0f; - CHECK(a != b, "after mutation, a should differ from b"); - a[1] = 1.0f; // restore - } - -// ---------------- Vector + Vector (and +=) ---------------- - { - Vf a(3, 1.0f); // [1,1,1] - Vf b(3, 2.0f); // [2,2,2] - Vf expect(3, 3.0f); // [3,3,3] - - Vf c = a + b; - CHECK(c == expect, "a + b should be [3,3,3]"); - - a += b; // a becomes [3,3,3] - CHECK(a == expect, "a += b should produce [3,3,3]"); - } - - // ---------------- Vector - Vector (and -=) ---------------- - { - Vf a(3, 5.0f); // [5,5,5] - Vf b(3, 2.0f); // [2,2,2] - Vf expect(3, 3.0f); // [3,3,3] - - Vf c = a - b; - CHECK(c == expect, "a - b should be [3,3,3]"); - - a -= b; // a becomes [3,3,3] - CHECK(a == expect, "a -= b should produce [3,3,3]"); - } - - // ---------------- Elementwise Multiply (and *=) ---------------- - { - Vf a(3, 2.0f); // [2,2,2] - Vf b(3, 3.0f); // [3,3,3] - Vf expect(3, 6.0f); // [6,6,6] - - Vf c = a * b; - CHECK(c == expect, "a * b (elemwise) should be [6,6,6]"); - - a *= b; // a becomes [6,6,6] - CHECK(a == expect, "a *= b should produce [6,6,6]"); - } - - // ---------------- Elementwise Divide (and /=) ---------------- - { - Vf a(3, 6.0f); // [6,6,6] - Vf b(3, 2.0f); // [2,2,2] - Vf expect(3, 3.0f); // [3,3,3] - - Vf c = a / b; - CHECK(c == expect, "a / b (elemwise) should be [3,3,3]"); - - a /= b; // a becomes [3,3,3] - CHECK(a == expect, "a /= b should produce [3,3,3]"); - } - - // ---------------- Scalar + Vector (and +=) ---------------- - { - Vf a(3, 1.0f); // [1,1,1] - Vf expect1(3, 6.0f); // [6,6,6] - Vf expect2(3, 3.0f); // [3,3,3] - - Vf c = a + 5.0f; // v + s - CHECK(c == expect1, "a + 5 should be [6,6,6]"); - - Vf d = 2.0f + a; // s + v (friend operator) - CHECK(d == expect2, "2 + a should be [3,3,3]"); - - a += 2.0f; // a becomes [3,3,3] - CHECK(a == expect2, "a += 2 should produce [3,3,3]"); - } - - // ---------------- Scalar - Vector / Vector - Scalar ---------------- - { - Vf a(3, 5.0f); // [5,5,5] - Vf expect1(3, 3.0f); // [3,3,3] - Vf expect2(3, -3.0f); // [ -3,-3,-3 ] if 2 - a (only if you've implemented it) - - Vf c = a - 2.0f; // v - s - CHECK(c == expect1, "a - 2 should be [3,3,3]"); - - // NOTE: Your friend operator-(U a, const Vector b) currently returns (b - a), - // which means `2 - a` computes `a - 2`. That's a bit unusual. - // We'll avoid asserting 2 - a here to match your current implementation choice. - (void)expect2; // silence unused warning - } - - // ---------------- Scalar * Vector / Vector * Scalar ---------------- - { - - Vf a(3, 2.0f); // [2,2,2] - uint64_t b = 3; // 3 - Vf expect(3, 6.0f); // [6,6,6] - - Vf c = a * b; // v * s - CHECK(c == expect, "a * 3 should be [6,6,6]"); - - Vf d = b * a; // s * v (friend) - CHECK(d == expect, "3 * a should be [6,6,6]"); - - a *= b; // a becomes [6,6,6] - CHECK(a == expect, "a *= 3 should produce [6,6,6]"); - } - - // ---------------- Vector / Scalar (and /= scalar) ---------------- - { - Vf a(3, 6.0f); // [6,6,6] - uint64_t b = 2; // 3 - Vf expect(3, 3.0f); // [3,3,3] - - Vf c = a / b; // v / s - CHECK(c == expect, "a / 2 should be [3,3,3]"); - - a /= b; // a becomes [3,3,3] - CHECK(a == expect, "a /= 2 should produce [3,3,3]"); - } - - // ---------- sum ---------- - { - Vf a(3, 2.0f); // [2,2,2] - CHECK(a.sum() == 6.0f, "sum failed"); - } - - // ---------- dot ---------- - { - Vf a(3, 2.0f); // [2,2,2] - Vf b(3, 3.0f); // [3,3,3] - CHECK(a.dot(b) == 18.0f, "dot failed"); // 2*3 * 3 = 18 - Vf c(4, 1.0f); - expect_throw([&]{ (void)a.dot(c); }, "dot should throw on size mismatch"); - } - - // ---------- norm ---------- - { - Vf a(3, 2.0f); // [2,2,2] - float n = a.norm(); - CHECK(std::fabs(n - std::sqrt(12.0f)) < 1e-6f, "norm failed"); - } - - // ---------- normalize ---------- - { - Vf a(3, 3.0f); // [3,3,3], norm = sqrt(27) - Vf b = a.normalize(); - float n = b.norm(); - CHECK(std::fabs(n - 1.0f) < 1e-6f, "normalize failed"); - expect_throw([&]{ Vf z(3, 0.0f); z.inplace_normalize(); }, "normalize should throw on zero norm"); - } - - // ---------- scalar power ---------- - { - Vf a(3, 2.0f); // [2,2,2] - Vf c = a.power(3); // [8,8,8] - CHECK(c == Vf(3, 8.0f), "power(scalar) failed"); - - Vf d = a; d.inplace_power(4); // [16,16,16] - CHECK(d == Vf(3, 16.0f), "inplace_power(scalar) failed"); - } - - // ---------- vector power ---------- - { - Vf base(3, 2.0f); // [2,2,2] - Vf exps; exps.v = {1.0f, 2.0f, 3.0f}; // explicit construction for clarity - Vf out = base.power(exps); // [2^1, 2^2, 2^3] = [2,4,8] - Vf expect; expect.v = {2.0f, 4.0f, 8.0f}; - CHECK(out == expect, "power(vector) failed"); - - expect_throw([&]{ Vf bad(2, 1.0f); (void)base.power(bad); }, - "power(vector) should throw on size mismatch"); - } - - // ---------- square ---------- - { - Vf a; a.v = {4.0f, 9.0f, 16.0f}; - Vf b = a.sqrt(); // [4,9,16] - Vf expect; expect.v = {2.0f, 3.0f, 4.0f}; - CHECK(b == expect, "sqrt failed"); - - a.inplace_sqrt(); // mutate a to [4,9,16] - CHECK(a == expect, "inplace_square failed"); - } - - // ---------- scalar commutative friends (s + v, s * v) ---------- - { - Vf a(3, 2.0f); // [2,2,2] - Vf b = 3.0f + a; // [5,5,5] - Vf c = a + 3.0f; // [5,5,5] - CHECK(b == c, "s+v commutative failed"); - - Vf d = 4.0f * a; // [8,8,8] - Vf e = a * 4.0f; // [8,8,8] - CHECK(d == e, "s*v commutative failed"); - } - - // ---------------- Size mismatch throws ---------------- - { - Vf a(3, 1.0f); - Vf b(4, 2.0f); - - expect_throw([&] { a.inplace_add(b); }, - "inplace_add should throw on size mismatch"); - expect_throw([&] { (void)(a + b); }, - "operator+ should throw (through add) on size mismatch"); - expect_throw([&] { a.inplace_subtract(b); }, - "inplace_subtract should throw on size mismatch"); - expect_throw([&] { (void)(a - b); }, - "operator- should throw (through subtract) on size mismatch"); - expect_throw([&] { a.inplace_multiply(b); }, - "inplace_multiply should throw on size mismatch"); - expect_throw([&] { (void)(a * b); }, - "operator* should throw (through multiply) on size mismatch"); - expect_throw([&] { a.inplace_divide(b); }, - "inplace_divide should throw on size mismatch"); - expect_throw([&] { (void)(a / b); }, - "operator/ should throw (through divide) on size mismatch"); - } - - { - - auto* a = new utils::Vf(3, 1.0f); // constructor runs - delete a; // <- calls ~Vector() and frees memory - - } - - { - Vf a(2, 1.0f); // a = [1, 1] - Vf b(2, 1.0f); // b = [1, 1] - - a.clear(); // a = [] - CHECK(a.size() == 0, "clear() did not empty vector"); - - a.resize(2, 1.0f); // a = [1, 1] - - CHECK(a == b, "clear/resize lifecycle failed"); - - } - - std::cout << "All Vector tests passed ✅\n"; - - - // shape + element access - { - Mf M(3, 4, 0.0f); - CHECK(M.rows()==3 && M.cols()==4, "shape failed"); - - M(1,1) = 5.0f; - CHECK(M(1,1) == 5.0f, "write/read element failed"); - - // ensure independence of other cells - CHECK(M(0,0) == 0.0f && M(2,3) == 0.0f, "unexpected element modified"); - } - - // set/get row (with size checks) - { - Mf M(2, 3, 0.0f); // 2x3 - - Vf r(3, 0.0f); - r[0]=1; r[1]=2; r[2]=3; - M.set_row(1, r); - - Vf g = M.get_row(1); - CHECK(g.size()==3, "get_row size wrong"); - CHECK(g[0]==1 && g[1]==2 && g[2]==3, "get_row values wrong"); - - // size mismatch should throw - bool threw=false; - try { - Vf bad(2, 9.0f); - M.set_row(0, bad); - } catch (const std::exception&) { threw=true; } - CHECK(threw, "set_row should throw on size mismatch"); - } - - // set/get col (with size checks) - { - Mf M(3, 2, 0.0f); // 3x2 - - Vf c(3, 0.0f); - c[0]=4; c[1]=5; c[2]=6; - M.set_col(1, c); - - Vf h = M.get_col(1); - CHECK(h.size()==3, "get_col size wrong"); - CHECK(h[0]==4 && h[1]==5 && h[2]==6, "get_col values wrong"); - - bool threw=false; - try { - Vf bad(2, 7.0f); - M.set_col(0, bad); - } catch (const std::exception&) { threw=true; } - CHECK(threw, "set_col should throw on size mismatch"); - } - - // swap_rows / swap_cols - { - Mf M(3, 3, 0.0f); - // set rows to [1,2,3], [4,5,6], [7,8,9] - for (uint64_t j=0;j<3;++j) M(0,j) = 1.0f + j; - for (uint64_t j=0;j<3;++j) M(1,j) = 4.0f + j; - for (uint64_t j=0;j<3;++j) M(2,j) = 7.0f + j; - - M.swap_rows(0,2); - CHECK(M(0,0)==7 && M(0,1)==8 && M(0,2)==9, "swap_rows top row wrong"); - CHECK(M(2,0)==1 && M(2,1)==2 && M(2,2)==3, "swap_rows bottom row wrong"); - - M.swap_cols(0,2); - // after col swap: first row should be [9,8,7] - CHECK(M(0,0)==9 && M(0,1)==8 && M(0,2)==7, "swap_cols first row wrong"); - // bottom row should be [3,2,1] - CHECK(M(2,0)==3 && M(2,1)==2 && M(2,2)==1, "swap_cols last row wrong"); - } - - // Exact integer comparison / Floating-point exact equality / Floating-point with small perturbation - { - - Mi A(2,2,0); - A(0,0)=1; A(0,1)=2; - A(1,0)=3; A(1,1)=4; - - Mi B(2,2,0); - B(0,0)=1; B(0,1)=2; - B(1,0)=3; B(1,1)=4; - - Mi C(2,2,0); - C(0,0)=9; C(0,1)=9; - C(1,0)=9; C(1,1)=9; - - CHECK(A == B, "Matrix == failed on identical int matrices"); - CHECK(!(A != B), "Matrix != failed on identical int matrices"); - CHECK(A != C, "Matrix != failed on different int matrices"); - - // Floating-point exact equality - - - Md F1(2,2,0.0); - F1(0,0)=1.0; F1(0,1)=2.0; - F1(1,0)=3.0; F1(1,1)=4.0; - - Md F2(2,2,0.0); - F2(0,0)=1.0; F2(0,1)=2.0; - F2(1,0)=3.0; F2(1,1)=4.0; - - CHECK(F1 == F2, "Matrix == failed on identical float matrices"); - - // Floating-point with small perturbation - - Md F3 = F1; - F3(1,1) += 1e-10; // tiny difference - - CHECK(!(F1 == F3), "Matrix == should fail on exact compare with perturbation"); - CHECK(F1.nearly_equal(F3, 1e-9), "Matrix nearly_equal failed with tolerance"); - - // Larger perturbation - F3(1,1) += 1e-3; - CHECK(!F1.nearly_equal(F3, 1e-6), "Matrix nearly_equal should fail when tolerance too small"); - CHECK(F1.nearly_equal(F3, 1e-2), "Matrix nearly_equal should pass with loose tolerance"); - } - - std::cout << "Matrix basic tests passed ✅\n"; - - // --- Test: normal transpose --- - { - Mf M(2, 3, 0.0f); - // Fill: [ [1,2,3], - // [4,5,6] ] - M(0,0)=1; M(0,1)=2; M(0,2)=3; - M(1,0)=4; M(1,1)=5; M(1,2)=6; - - Mf MT = numerics::transpose(M); - - // Should be shape 3x2 - CHECK(MT.rows()==3 && MT.cols()==2, "transpose shape wrong"); - - // Values: [ [1,4], [2,5], [3,6] ] - CHECK(MT(0,0)==1 && MT(0,1)==4, "transpose value (0,*) wrong"); - CHECK(MT(1,0)==2 && MT(1,1)==5, "transpose value (1,*) wrong"); - CHECK(MT(2,0)==3 && MT(2,1)==6, "transpose value (2,*) wrong"); - - //std::cout << "Original M:\n" << M << "\n"; - //std::cout << "Transposed MT:\n" << MT << "\n\n"; - } - - // --- Test: inplace transpose (square only) --- - { - Mf S(3, 3, 0.0f); - // Fill with row-major increasing - float val = 1.0f; - for (uint64_t i=0;i(A, x); - CHECK(y.size()==2, "matvec size wrong"); - CHECK(y[0]==50 && y[1]==122, "matvec values wrong"); - - // dimension mismatch should throw - bool threw = false; - try { - Vd bad(4,1.0); - (void)numerics::matvec(A, bad); - } catch (const std::runtime_error&) { threw = true; } - CHECK(threw, "matvec: expected throw on dim mismatch"); - } - - std::cout << "matvec tests passed ✅\n"; - - // vecmat test - { - // A = [[1,2], - // [3,4]] (2x2) - Md A(2,2,0.0); - A(0,0)=1; A(0,1)=2; - A(1,0)=3; A(1,1)=4; - - // x^T = [5,6] - Vd x(2,0.0); - x[0]=5; x[1]=6; - - // y = x^T * A = [5*1+6*3, 5*2+6*4] = [23, 34] - Vd y = numerics::vecmat(x, A); - CHECK(y.size()==2, "vecmat size wrong"); - CHECK(y[0]==23 && y[1]==34, "vecmat values wrong"); - - // mismatch should throw - bool threw = false; - try { - Md B(3,2,0.0); // 3x2, doesn't match x size 2 - (void)numerics::vecmat(x, B); - } catch (const std::runtime_error&) { threw = true; } - CHECK(threw, "vecmat: expected throw on dim mismatch"); - } - - std::cout << "vecmat tests passed ✅\n"; - - - - // Inverse 'Gauss-Jordan' tests - { - Md A(2,2,0.0); - A(0,0)=4; A(0,1)=7; - A(1,0)=2; A(1,1)=6; - - Md Ai = numerics::inverse(A, "Gauss-Jordan"); - - Md I1 = numerics::matmul(A, Ai); - Md I2 = numerics::matmul(Ai, A); - - Md I(2,2,0.0); - I(0,0)=1; I(1,1)=1; - - CHECK(I1.nearly_equal(I), "A*inv(A) != I"); - CHECK(I2.nearly_equal(I), "inv(A)*A != I"); - } - - std::cout << "Inverse 'Gauss-Jordan' tests passed ✅\n"; return 0; } \ No newline at end of file diff --git a/test/test_all.cpp b/test/test_all.cpp new file mode 100644 index 0000000..ed3c41b --- /dev/null +++ b/test/test_all.cpp @@ -0,0 +1,2 @@ +#define TEST_MAIN +#include "test_common.h" \ No newline at end of file diff --git a/test/test_common.h b/test/test_common.h new file mode 100644 index 0000000..c3d151c --- /dev/null +++ b/test/test_common.h @@ -0,0 +1,52 @@ +#ifndef _test_common_n_ +#define _test_common_n_ + +#pragma once +#include +#include +#include +#include + +struct TestFailure : public std::runtime_error { + using std::runtime_error::runtime_error; +}; + + +#define CHECK(cond, msg) do { if (!(cond)) throw TestFailure(msg); } while (0) +#define CHECK_EQ(a,b,msg) do { if (!((a)==(b))) { throw TestFailure(std::string(msg) + " (" #a " != " #b ")"); } } while (0) + +#define TEST_CASE(name) \ + static void name(); \ + struct name##_registrar { name##_registrar(){ TestRegistry::add(#name, &name);} } name##_registrar_instance; \ + static void name() + +struct TestRegistry { + using Fn = void(*)(); + static std::vector>& list() { + static std::vector> v; return v; + } + static void add(const std::string& name, Fn fn) { list().push_back({name, fn}); } +}; + +// Default test runner main() +#ifdef TEST_MAIN +int main() { + int fails = 0; + for (auto& t : TestRegistry::list()) { + try { + t.second(); + std::cout << "[PASS] " << t.first << "\n"; + } catch (const TestFailure& e) { + std::cerr << "[FAIL] " << t.first << " -> " << e.what() << "\n"; + ++fails; + } catch (const std::exception& e) { + std::cerr << "[ERROR] " << t.first << " -> " << e.what() << "\n"; + ++fails; + } + } + std::cout << (fails ? "Some tests failed ❌\n" : "All tests passed ✅\n"); + return fails ? 1 : 0; +} +#endif + +#endif // _test_common_n_ \ No newline at end of file diff --git a/test/test_inverse.cpp b/test/test_inverse.cpp new file mode 100644 index 0000000..b18949b --- /dev/null +++ b/test/test_inverse.cpp @@ -0,0 +1,126 @@ +#include "test_common.h" +#include "./utils/utils.h" +#include "./numerics/inverse.h" +#include "./numerics/matmul.h" + + +TEST_CASE(Inverse_2x2_WellConditioned) { + using T = double; + // A = [[4,7],[2,6]] inverse = (1/10) * [[6,-7],[-2,4]] + utils::Matrix A(2,2, T{0}); + A(0,0)=4; A(0,1)=7; + A(1,0)=2; A(1,1)=6; + + auto Ainv = numerics::inverse(A); // out-of-place + + // Check A * Ainv ≈ I and Ainv * A ≈ I + auto Ileft = numerics::matmul(A, Ainv); + auto Iright = numerics::matmul(Ainv, A); + + utils::Md Iref(2,2, T{0}); + for (uint64_t i=0;i(2); + + CHECK((Ileft.nearly_equal(Iref, 1e-12)), "A * inverse(A) ≠ I"); + CHECK((Iright.nearly_equal(Iref, 1e-12)), "inverse(A) * A ≠ I"); +} + +TEST_CASE(Inverse_InPlace_Equals_OutOfPlace) { + using T = double; + utils::Matrix A(3,3, T{0}); + // A = [[3, 0, 2], + // [2, 0, -2], + // [0, 1, 1]] + A(0,0)=3; A(0,1)=0; A(0,2)= 2; + A(1,0)=2; A(1,1)=0; A(1,2)=-2; + A(2,0)=0; A(2,1)=1; A(2,2)= 1; + + auto Ainv_ref = numerics::inverse(A); // copy path + + auto A_inp = A; + numerics::inplace_inverse(A_inp); // in-place path + + CHECK((A_inp.nearly_equal(Ainv_ref, 1e-12)), "in-place inverse differs from out-of-place"); +} + +TEST_CASE(Inverse_Singular_Throws) { + using T = double; + utils::Matrix S(2,2, T{0}); + // Singular: rows are multiples → det = 0 + S(0,0)=1; S(0,1)=2; + S(1,0)=2; S(1,1)=4; + + bool threw=false; + try { + auto _ = numerics::inverse(S); + (void)_; + } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "inverse should throw on singular matrix"); + + threw=false; + try { + numerics::inplace_inverse(S); + } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "inplace_inverse should throw on singular matrix"); +} + +TEST_CASE(Inverse_RoundTrip_DiagonallyDominant_5x5) { + // Build a well-conditioned 5x5: diagonally dominant + utils::Md A(5,5,0.0); + for (uint64_t i=0;i<5;++i) { + double rowsum = 0.0; + for (uint64_t j=0;j<5;++j) { + if (i==j) continue; + A(i,j) = 0.01 * double(1 + ((i+1)*(j+3)) % 7); + rowsum += std::fabs(A(i,j)); + } + A(i,i) = rowsum + 1.0; // strictly diagonally dominant + } + + utils::Md A_copy = A; // ensure wrapper doesn't mutate input + utils::Md Ainv = numerics::inverse(A); + + // Input must be unchanged by the non-inplace wrapper + CHECK(A.nearly_equal(A_copy, 0.0), "inverse wrapper modified input"); + + + utils::Md I(5,5, 0); + for (uint64_t i=0;i(A, Ainv); + auto R = numerics::matmul(Ainv, A); + + CHECK(L.nearly_equal(I, 1e-10), "A * Ainv not close to I"); + CHECK(R.nearly_equal(I, 1e-10), "Ainv * A not close to I"); +} + +TEST_CASE(Inverse_NonSquare_Throws) { + // Non-square: 2x3 — algorithm expects square; should throw + utils::Md A(2,3,0.0); + bool threw = false; + try { + numerics::inplace_inverse(A); + } catch (const std::runtime_error&) { + threw = true; + } catch (...) { + threw = true; // any failure is fine; must not silently succeed + } + CHECK(threw, "inplace_inverse should throw on non-square matrix"); +} + + +TEST_CASE(Inverse_Unknown_Method_Throws) { + + utils::Md A(3,3, 0); + for (uint64_t i=0;i(A, "NotARealMethod"); + } catch (const std::runtime_error&) { + threw = true; + } + CHECK(threw, "should throw for unknown inverse method"); +} \ No newline at end of file diff --git a/test/test_matmul.cpp b/test/test_matmul.cpp new file mode 100644 index 0000000..8808fdf --- /dev/null +++ b/test/test_matmul.cpp @@ -0,0 +1,164 @@ +#include "test_common.h" +#include "./utils/utils.h" +#include "./numerics/matmul.h" + +#include + + +// ============ Basic correctness ============ +TEST_CASE(Matmul_Serial_Simple3x3) { + utils::Md A(3,3,0.0), B(3,3,0.0); + // A = [[1,2,3],[4,5,6],[7,8,9]] + double v=1.0; + for (uint64_t i=0;i<3;++i) for (uint64_t j=0;j<3;++j) A(i,j)=v++; + // B = [[9,8,7],[6,5,4],[3,2,1]] + double w=9.0; + for (uint64_t i=0;i<3;++i) for (uint64_t j=0;j<3;++j) B(i,j)=w--; + + auto C = numerics::matmul(A,B); + // Hand-checked first row: + // row0 dot columns: + // c00 = 1*9 + 2*6 + 3*3 = 30 + // c01 = 1*8 + 2*5 + 3*2 = 24 + // c02 = 1*7 + 2*4 + 3*1 = 18 + CHECK(C.rows()==3 && C.cols()==3, "shape 3x3 wrong"); + CHECK(C(0,0)==30.0 && C(0,1)==24.0 && C(0,2)==18.0, "first row wrong"); +} + +TEST_CASE(Matmul_OMP_Equals_Serial) { + utils::Md A(4,5,0.0), B(5,3,0.0); + // Fill deterministic + for (uint64_t i=0;i(A,B); + auto Cr = numerics::matmul_rows_omp(A,B); + auto Cc = numerics::matmul_collapse_omp(A,B); + auto Ca = numerics::matmul_auto(A,B); + + CHECK((Cs.nearly_equal(Cr, 1e-12)), "rows_omp != serial"); + CHECK((Cs.nearly_equal(Cc, 1e-12)), "collapse_omp != serial"); + CHECK((Cs.nearly_equal(Ca, 1e-12)), "auto != serial"); +} + +// ============ Dimension mismatch ============ +TEST_CASE(Matmul_DimensionMismatch_Throws) { + utils::Md A(2,3,0.0), B(4,2,0.0); + bool threw=false; + try { auto _ = numerics::matmul(A,B); (void)_; } + catch (const std::runtime_error&) { threw=true; } + CHECK(threw, "serial should throw on dim mismatch"); + + threw=false; try { auto _ = numerics::matmul_rows_omp(A,B); (void)_; } + catch (const std::runtime_error&) { threw=true; } + CHECK(threw, "rows_omp should throw on dim mismatch"); + + threw=false; try { auto _ = numerics::matmul_collapse_omp(A,B); (void)_; } + catch (const std::runtime_error&) { threw=true; } + CHECK(threw, "collapse_omp should throw on dim mismatch"); +} + +// ============ Edge cases ============ +TEST_CASE(Matmul_Edges_ZeroDims) { + // (0xK) * (KxP) -> (0xP) + utils::Md A0(0,5,0.0), B1(5,3,0.0); + auto C0 = numerics::matmul(A0,B1); + CHECK(C0.rows()==0 && C0.cols()==3, "0xK * KxP shape wrong"); + + // (MxK) * (Kx0) -> (Mx0) + utils::Md A2(7,4,0.0), B0(4,0,0.0); + auto C1 = numerics::matmul(A2,B0); + CHECK(C1.rows()==7 && C1.cols()==0, "MxK * Kx0 shape wrong"); +} + +// ============ Identity sanity ============ +TEST_CASE(Matmul_Identity) { + const uint64_t n=5; + utils::Md I(n,n,0.0), A(n,n,0.0); + for (uint64_t i=0;i(I,A); + auto R = numerics::matmul(A,I); + CHECK(L == A, "I*A != A"); + CHECK(R == A, "A*I != A"); +} + +// ============ Perf sanity (same kernel: 1 thread vs many) ============ +template +static double time_it(F&& f, int iters=1) { + auto t0 = std::chrono::high_resolution_clock::now(); + for (int i=0;i(t1 - t0).count(); +} + +TEST_CASE(Matmul_Perf_Sanity_RowOMP) { +#ifndef _OPENMP + return; +#else + int hw = omp_get_max_threads(); + if (hw <= 1) return; + + const uint64_t m=512, k=512, p=512; // ~134M MACs; adjust if needed + utils::Md A(m,k,0.0), B(k,p,0.0); + for (uint64_t i=0;i(A,B); + + int prev = omp_get_max_threads(); + auto t0 = std::chrono::high_resolution_clock::now(); + omp_set_num_threads(1); + numerics::matmul_rows_omp(A,B); + double t1 = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + + omp_set_num_threads(hw); + t0 = std::chrono::high_resolution_clock::now(); + numerics::matmul_rows_omp(A,B); + double tN = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + + omp_set_num_threads(prev); + + // Must not be notably slower with many threads + CHECK(tN <= t1 * 1.05, "rows_omp: multi-thread slower than single-thread"); +#endif +} + +TEST_CASE(Matmul_Perf_Sanity_CollapseOMP) { +#ifndef _OPENMP + return; +#else + int hw = omp_get_max_threads(); + if (hw <= 1) return; + + const uint64_t m=512, k=512, p=512; + utils::Md A(m,k,0.0), B(k,p,0.0); + for (uint64_t i=0;i(A,B); // warm-up + + int prev = omp_get_max_threads(); + auto t0 = std::chrono::high_resolution_clock::now(); + omp_set_num_threads(1); + numerics::matmul_collapse_omp(A,B); + double t1 = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + + omp_set_num_threads(hw); + t0 = std::chrono::high_resolution_clock::now(); + numerics::matmul_collapse_omp(A,B); + double tN = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + + omp_set_num_threads(prev); + + CHECK(tN <= t1 * 1.05, "collapse_omp: multi-thread slower than single-thread"); +#endif +} \ No newline at end of file diff --git a/test/test_matrix.cpp b/test/test_matrix.cpp new file mode 100644 index 0000000..5fc1b02 --- /dev/null +++ b/test/test_matrix.cpp @@ -0,0 +1,142 @@ + +#include "test_common.h" +#include "./utils/utils.h" + +using utils::Vf; using utils::Vd; using utils::Vi; +using utils::Mf; using utils::Md; using utils::Mi; + + +// ---------- Construction & element access ---------- +TEST_CASE(Matrix_Construct_Access) { + Md M; // default + CHECK(M.rows()==0 && M.cols()==0, "default ctor dims wrong"); + + Mf A(2,3, 1.0f); + CHECK(A.rows()==2 && A.cols()==3, "ctor dims wrong"); + CHECK(A(0,0)==1.0f && A(1,2)==1.0f, "fill wrong"); + + A(0,1)=2.5f; A(1,0)=3.5f; + CHECK(A(0,1)==2.5f && A(1,0)==3.5f, "operator() set/get failed"); +} + +// ---------- Equality, inequality, nearly_equal ---------- +TEST_CASE(Matrix_Equality) { + Mi A(2,2,0), B(2,2,0), C(2,2,1); + A(0,0)=1; A(1,1)=1; // A = I + B(0,0)=1; B(1,1)=1; // B = I + + CHECK(A == B, "== failed identical"); + CHECK(!(A != B), "!= failed identical"); + CHECK(A != C, "!= failed different"); + + Md F1(2,2,0.0), F2(2,2,0.0); + F1(0,0)=1.0; F1(1,1)=2.0; + F2(0,0)=1.0; F2(1,1)=2.0 + 5e-10; // tiny perturbation + CHECK(!(F1 == F2), "operator== is exact; should differ"); + CHECK(F1.nearly_equal(F2, 1e-9), "nearly_equal should accept tiny delta"); + CHECK(!F1.nearly_equal(F2, 1e-12), "nearly_equal too strict should fail"); +} + +// ---------- Row helpers ---------- +TEST_CASE(Matrix_Row_Get_Set) { + Mf M(3,4, 0.0f); + Vf r(4, 0.0f); + for (uint64_t j=0;j<4;++j) r[j] = float(j+1); // [1,2,3,4] + + M.set_row(1, r); + auto out = M.get_row(1); + CHECK(out == r, "set_row/get_row mismatch"); + + // size mismatch should throw + bool threw=false; + try { Vf bad(3, 9.0f); M.set_row(2, bad); } catch (const std::exception&) { threw=true; } + CHECK(threw, "set_row should throw on size mismatch"); + + // out of range + threw=false; + try { (void)M.get_row(3); } catch (const std::out_of_range&) { threw=true; } + CHECK(threw, "get_row should throw on OOB index"); +} + +// ---------- Column helpers ---------- +TEST_CASE(Matrix_Col_Get_Set) { + Md M(3,2, 0.0); + Vd c(3, 0.0); + c[0]=10; c[1]=20; c[2]=30; + + M.set_col(1, c); + auto out = M.get_col(1); + CHECK(out == c, "set_col/get_col mismatch"); + + // size mismatch should throw + bool threw=false; + try { Vd bad(2, 9.0); M.set_col(0, bad); } catch (const std::exception&) { threw=true; } + CHECK(threw, "set_col should throw on size mismatch"); + + // out of range + threw=false; + try { (void)M.get_col(2); } catch (const std::out_of_range&) { threw=true; } + CHECK(threw, "get_col should throw on OOB index"); +} + +// ---------- swap_rows / swap_cols ---------- +TEST_CASE(Matrix_Swap_Rows_Cols) { + Mi M(2,3,0); + // Row 0: [1,2,3], Row 1: [4,5,6] + M(0,0)=1; M(0,1)=2; M(0,2)=3; + M(1,0)=4; M(1,1)=5; M(1,2)=6; + + M.swap_rows(0,1); + CHECK(M(0,0)==4 && M(0,1)==5 && M(0,2)==6, "swap_rows row0 wrong"); + CHECK(M(1,0)==1 && M(1,1)==2 && M(1,2)==3, "swap_rows row1 wrong"); + + // swap back via cols + M.swap_cols(0,2); + // After swapping col0<->col2: + // Row0: [6,5,4], Row1: [3,2,1] + CHECK(M(0,0)==6 && M(0,1)==5 && M(0,2)==4, "swap_cols row0 wrong"); + CHECK(M(1,0)==3 && M(1,1)==2 && M(1,2)==1, "swap_cols row1 wrong"); + + // no-op swap (a==b) should not crash or change + M.swap_rows(1,1); + M.swap_cols(2,2); + + // OOB checks + bool threw=false; + try { M.swap_rows(5,1); } catch (const std::out_of_range&) { threw=true; } + CHECK(threw, "swap_rows should throw on OOB"); + threw=false; + try { M.swap_cols(0,9); } catch (const std::out_of_range&) { threw=true; } + CHECK(threw, "swap_cols should throw on OOB"); +} + +// ---------- data() layout (contiguous row-major) ---------- +TEST_CASE(Matrix_Data_Layout) { + Md M(2,3, 0.0); + // Fill increasing sequence + double val=1.0; + for (uint64_t i=0;i + +using utils::Vi; using utils::Vf; using utils::Vd; +using utils::Mi; using utils::Mf; using utils::Md; + + + +// ------------------------------------------------------------ +// matvec: y = A * x +// ------------------------------------------------------------ +TEST_CASE(Matvec_Serial_Simple) { + // A = [[1,2,3], + // [4,5,6]] + Md A(2,3,0.0); + A(0,0)=1; A(0,1)=2; A(0,2)=3; + A(1,0)=4; A(1,1)=5; A(1,2)=6; + Vd x(3,0.0); x[0]=7; x[1]=8; x[2]=9; + + auto y = numerics::matvec(A,x); // [ 1*7+2*8+3*9 , 4*7+5*8+6*9 ] = [50, 122] + CHECK(y.size()==2, "matvec size wrong"); + CHECK(y[0]==50.0 && y[1]==122.0, "matvec values wrong"); +} + +TEST_CASE(Matvec_OMP_Equals_Serial) { + Md A(3,3,0.0); + // A = I * 2 + for (uint64_t i=0;i<3;++i) A(i,i)=2.0; + Vd x(3,0.0); x[0]=1; x[1]=2; x[2]=3; + + auto ys = numerics::matvec(A,x); + auto yp = numerics::matvec_omp(A,x); + + CHECK((ys.nearly_equal_vec(yp)), "matvec_omp != matvec"); +} + +TEST_CASE(Matvec_Auto_Equals_Serial) { + Md A(2,2,0.0); A(0,0)=2; A(0,1)=1; A(1,0)=0.5; A(1,1)=3; + Vd x(2,0.0); x[0]=4; x[1]=5; + + auto ys = numerics::matvec(A,x); + auto ya = numerics::matvec_auto(A,x); + + CHECK((ys.nearly_equal_vec(ya)), "matvec_auto != serial"); +} + +TEST_CASE(Matvec_DimensionMismatch_Throws) { + Md A(2,3,0.0); + Vd x(4,0.0); + bool threw=false; + try { auto _ = numerics::matvec(A,x); (void)_; } + catch (const std::runtime_error&) { threw=true; } + CHECK(threw, "matvec must throw on dimension mismatch"); +} + +TEST_CASE(Matvec_Zero_Edges) { + Md A(0,3,0.0); // 0x3 + Vd x(3,1.0); + auto y = numerics::matvec(A,x); + CHECK(y.size()==0, "0xN * x should return size 0 vector"); + + Md B(2,0,0.0); // 2x0 + Vd z(0,0.0); + auto y2 = numerics::matvec(B,z); + CHECK(y2.size()==2 && y2[0]==0.0 && y2[1]==0.0, "N×0 * 0 should return zeros of size N"); +} + +TEST_CASE(Matvec_Float_Tolerance) { + Mf A(2,2,0.0f); A(0,0)=1.0f; A(0,1)=2.0f; A(1,0)=3.0f; A(1,1)=4.0f; + Vf x(2,0.0f); x[0]=0.1f; x[1]=0.2f; + + auto y1 = numerics::matvec(A,x); + auto y2 = numerics::matvec_omp(A,x); + + CHECK((y1.nearly_equal_vec(y2,1e-6f)), "matvec float omp mismatch"); +} +// +// ---------- Auto inside an outer parallel region (no accidental nested teams) ---------- +// We just check correctness; performance is environment-dependent. +// +TEST_CASE(Matvec_Auto_Inside_Outer_Parallel_Correctness) { + const uint64_t m=64, n=64; + Md A(m,n,1.0); Vd x(n,2.0); + //fill_deterministic(A); fill_deterministic(x); + Vd ref = numerics::matvec(A,x); + + // Call auto inside an outer team + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int rep=0; rep<32; ++rep) { + auto y = numerics::matvec_auto(A,x); + // Each thread checks its own result equals reference + if (!(y.nearly_equal_vec(ref))) { + throw TestFailure("matvec_auto wrong under outer parallel region"); + } + } +} +TEST_CASE(Matvec_Speed_Sanity) { + const uint64_t m=4096, n=4096; // ~16M MACs; adjust if needed + Md A(m,n,1.0); Vd x(n,2.0); + //fill_deterministic(A); fill_deterministic(x); + + auto t0 = std::chrono::high_resolution_clock::now(); + auto yS = numerics::matvec(A,x); + double tp = std::chrono::duration(t0 - std::chrono::high_resolution_clock::now()).count(); + + #ifdef _OPENMP + int threads = omp_get_max_threads(); + #else + int threads = 1; + #endif + + t0 = std::chrono::high_resolution_clock::now(); + auto yP = numerics::matvec_omp(A,x); + double ts = std::chrono::duration(t0 - std::chrono::high_resolution_clock::now()).count(); + + CHECK((yS.nearly_equal_vec(yP)), "matvec_omp != matvec_serial (large)"); + // Only enforce basic sanity if we *can* use >1 threads: + if (threads > 1) { + // Be generous: just require not significantly slower. + CHECK(tp <= ts, "matvec_omp unexpectedly much slower than serial"); + } +} + +// ------------------------------------------------------------ +// vecmat: y = x * A +// ------------------------------------------------------------ +TEST_CASE(Vecmat_Serial_Simple) { + // A = [[1,2], + // [3,4], + // [5,6]] (3x2) + Md A(3,2,0.0); + A(0,0)=1; A(0,1)=2; + A(1,0)=3; A(1,1)=4; + A(2,0)=5; A(2,1)=6; + + Vd x(3,0.0); x[0]=7; x[1]=8; x[2]=9; + + auto y = numerics::vecmat(x,A); // 1*7+3*8+5*9= 76 ; 2*7+4*8+6*9=100 + CHECK(y.size()==2, "vecmat size wrong"); + CHECK(y[0]==76.0 && y[1]==100.0, "vecmat values wrong"); +} + +TEST_CASE(Vecmat_OMP_Equals_Serial) { + Md A(2,2,0.0); A(0,0)=2; A(0,1)=1; A(1,0)=5; A(1,1)=-1; + Vd x(2,0.0); x[0]=0.5; x[1]=1.5; + + auto ys = numerics::vecmat(x,A); + auto yp = numerics::vecmat_omp(x,A); + + CHECK((ys.nearly_equal_vec(yp)), "vecmat_omp != vecmat"); +} + +TEST_CASE(Vecmat_Auto_Equals_Serial) { + Md A(2,3,0.0); + A(0,0)=1; A(0,1)=2; A(0,2)=3; + A(1,0)=4; A(1,1)=5; A(1,2)=6; + Vd x(2,0.0); x[0]=1; x[1]=2; + + auto ys = numerics::vecmat(x,A); + auto ya = numerics::vecmat_auto(x,A); + + CHECK((ys.nearly_equal_vec(ya)), "vecmat_auto != serial"); +} + +TEST_CASE(Vecmat_DimensionMismatch_Throws) { + Md A(2,2,0.0); + Vd x(3,0.0); + bool threw=false; + try { auto _ = numerics::vecmat(x,A); (void)_; } + catch (const std::runtime_error&) { threw=true; } + CHECK(threw, "vecmat must throw on dimension mismatch"); +} + +TEST_CASE(Vecmat_Zero_Edges) { + Md A(0,3,0.0); + Vd x(0,0.0); + auto y = numerics::vecmat(x,A); // 0×N times N×M → 0×M + CHECK(y.size()==3 && y[0]==0.0 && y[1]==0.0 && y[2]==0.0, "0-length x times A wrong"); + + Md B(3,0,0.0); + Vd z(3,1.0); + auto y2 = numerics::vecmat(z,B); // 1x3 * 3x0 → 1x0 + CHECK(y2.size()==0, "vecmat with N×0 result size wrong"); +} + +// +// ---------- Auto inside an outer parallel region (no accidental nested teams) ---------- +// We just check correctness; performance is environment-dependent. +// +TEST_CASE(Vecmat_Auto_Inside_Outer_Parallel_Correctness) { + const uint64_t m=64, n=64; + Md A(m,n,1.0); Vd x(m,2.0); + //fill_deterministic(A); fill_deterministic(x); + Vd ref = numerics::vecmat(x,A); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int rep=0; rep<32; ++rep) { + auto y = numerics::vecmat_auto(x,A); + if (!(y.nearly_equal_vec(ref))) { + throw TestFailure("vecmat_auto wrong under outer parallel region"); + } + } +} + + +TEST_CASE(Vecmat_Speed_Sanity) { + const uint64_t m=4096, n=4096; + Md A(m,n,1.0); Vd x(m,2.0); + //fill_deterministic(A); fill_deterministic(x); + + auto t0 = std::chrono::high_resolution_clock::now(); + auto yS = numerics::vecmat(x,A); + double ts = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + + #ifdef _OPENMP + int threads = omp_get_max_threads(); + #else + int threads = 1; + #endif + + t0 = std::chrono::high_resolution_clock::now(); + auto yP = numerics::vecmat_omp(x,A); + double tp = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + + CHECK((yS.nearly_equal_vec(yP)), "vecmat_omp != vecmat_serial (large)"); + if (threads > 1) { + CHECK(tp <= ts, "vecmat_omp unexpectedly much slower than serial"); + } +} diff --git a/test/test_transpose.cpp b/test/test_transpose.cpp new file mode 100644 index 0000000..4a898de --- /dev/null +++ b/test/test_transpose.cpp @@ -0,0 +1,88 @@ + +#include "test_common.h" +#include "./utils/utils.h" // matrix.h, vector.h +#include "./numerics/transpose.h" // numerics::transpose / inplace_transpose + +using utils::Mi; using utils::Mf; using utils::Md; + +// +// ---------- Out-of-place transpose (rectangular) ---------- +// +TEST_CASE(Transpose_Rectangular_OutOfPlace) { + // A = [ [1, 2, 3], + // [4, 5, 6] ] (2x3) + Md A(2,3,0.0); + A(0,0)=1; A(0,1)=2; A(0,2)=3; + A(1,0)=4; A(1,1)=5; A(1,2)=6; + + auto AT = numerics::transpose(A); // (3x2) + + CHECK(AT.rows()==3 && AT.cols()==2, "shape wrong after transpose"); + CHECK(AT(0,0)==1 && AT(1,0)==2 && AT(2,0)==3, "first column wrong"); + CHECK(AT(0,1)==4 && AT(1,1)==5 && AT(2,1)==6, "second column wrong"); + + // Involution: T(T(A)) == A + auto ATT = numerics::transpose(AT); + CHECK(ATT == A, "transpose(transpose(A)) != A"); +} + +// +// ---------- In-place transpose (square) ---------- +// +TEST_CASE(Transpose_Square_InPlace) { + // 3x3 with distinct values + Mf S(3,3,0.0f); + float val = 1.0f; + for (uint64_t i=0;i<3;++i) + for (uint64_t j=0;j<3;++j) + S(i,j) = val++; + + // Make an explicit transpose to compare against + auto ST = numerics::transpose(S); + + // In-place should match the out-of-place result + numerics::inplace_transpose(S); + CHECK(S == ST, "inplace_transpose result mismatch"); + + // Involution: applying in-place again should return original + numerics::inplace_transpose(S); + // Now S should equal the original pre-inplace matrix (which was transposed above) + // We can reconstruct original by transposing ST: + auto orig = numerics::transpose(ST); + CHECK(S == orig, "inplace transpose twice should restore original"); +} + +// +// ---------- In-place transpose must throw on non-square ---------- +// +TEST_CASE(Transpose_InPlace_Throws_On_Rectangular) { + Md R(2,3,0.0); // rectangular + bool threw = false; + try { + numerics::inplace_transpose(R); + } catch (const std::runtime_error&) { + threw = true; + } + CHECK(threw, "inplace_transpose must throw on non-square matrices"); +} + +// +// ---------- Edge cases: 0x0 and 1x1 ---------- +// +TEST_CASE(Transpose_Edge_0x0_1x1) { + // 0x0 should be fine both ways + Md E; // 0x0 + auto ET = numerics::transpose(E); + CHECK(ET.rows()==0 && ET.cols()==0, "0x0 transpose shape wrong"); + // in-place on 0x0 (rows==cols) should not throw + numerics::inplace_transpose(E); + CHECK(E.rows()==0 && E.cols()==0, "0x0 inplace transpose changed shape"); + + // 1x1 should be a no-op and not throw + Mi I(1,1,0); + I(0,0) = 42; + auto IT = numerics::transpose(I); + CHECK(IT.rows()==1 && IT.cols()==1 && IT(0,0)==42, "1x1 transpose wrong"); + numerics::inplace_transpose(I); + CHECK(I(0,0)==42, "1x1 inplace transpose changed value"); +} diff --git a/test/test_vector.cpp b/test/test_vector.cpp new file mode 100644 index 0000000..85beea1 --- /dev/null +++ b/test/test_vector.cpp @@ -0,0 +1,211 @@ + +#include "test_common.h" +#include "./utils/utils.h" + +using utils::Vf; using utils::Vd; using utils::Vi; + +// +// ---------- Basic construction & access ---------- +// +TEST_CASE(Vector_Construct_Size_Access) { + Vd a; // default + CHECK(a.size() == 0, "default size must be 0"); + + Vf b(3, 1.0f); // (n, fill) + CHECK(b.size() == 3, "size wrong"); + CHECK(b[0] == 1.0f && b[1] == 1.0f && b[2] == 1.0f, "fill wrong"); + + b[1] = 2.5f; + CHECK(b[1] == 2.5f, "operator[] write failed"); + + // resize (grow + value) + b.resize(5, 7.0f); + CHECK(b.size() == 5, "resize grow size wrong"); + CHECK(b[0] == 1.0f && b[1] == 2.5f && b[2] == 1.0f && b[3] == 7.0f && b[4] == 7.0f, + "resize grow values wrong"); + + // resize (shrink) + b.resize(2); + CHECK(b.size() == 2, "resize shrink size wrong"); +} + +TEST_CASE(Vector_Clear_PushBack) { + Vi v(0, 0); + v.push_back(10); + v.push_back(20); + CHECK(v.size() == 2, "push_back size wrong"); + CHECK(v[0] == 10 && v[1] == 20, "push_back values wrong"); + + v.clear(); + CHECK(v.size() == 0, "clear failed"); +} +// +// ---------- Equality / Inequality (tolerant for float/double) ---------- +// +TEST_CASE(Vector_Equality_Tolerant) { + Vd a(3, 1.0), b(3, 1.0); + CHECK(a == b, "== identical failed"); + CHECK(!(a != b), "!= identical failed"); + + // Tiny perturbation within eps (1e-6 default) + b[1] += 1e-7; + CHECK(a == b, "== tolerant failed"); + + // Larger perturbation should fail equality + b[1] += 1e-4; + CHECK(a != b, "!= with difference failed"); +} +// +// ---------- Scalar arithmetic: +, -, *, / (inplace and returning) ---------- +// +TEST_CASE(Vector_Scalar_Arithmetic) { + Vf a(3, 1.0f); + + // inplace + a.inplace_add(2); // int convertible to float + CHECK(a[0] == 3.0f && a[1] == 3.0f && a[2] == 3.0f, "inplace_add failed"); + + a.inplace_subtract(1.5f); + CHECK(std::fabs(a[0] - 1.5f) < 1e-6f && + std::fabs(a[1] - 1.5f) < 1e-6f && + std::fabs(a[2] - 1.5f) < 1e-6f, "inplace_subtract failed"); + + a.inplace_multiply(4.0); + CHECK(a[0] == 6.0f && a[1] == 6.0f && a[2] == 6.0f, "inplace_multiply failed"); + + a.inplace_divide(2); + CHECK(a[0] == 3.0f && a[1] == 3.0f && a[2] == 3.0f, "inplace_divide failed"); + + // returning + auto b = a + 1.0f; + CHECK(b[0] == 4.0f && b[1] == 4.0f && b[2] == 4.0f, "operator+(scalar) failed"); + + b = a - 2.0f; + CHECK(b[0] == 1.0f && b[1] == 1.0f && b[2] == 1.0f, "operator-(scalar) failed"); + + b = a * 10; // int -> float + CHECK(b[0] == 30.0f && b[1] == 30.0f && b[2] == 30.0f, "operator*(scalar) failed"); + + b = a / 3.0f; + CHECK(std::fabs(b[0] - 1.0f) < 1e-6f && + std::fabs(b[1] - 1.0f) < 1e-6f && + std::fabs(b[2] - 1.0f) < 1e-6f, "operator/(scalar) failed"); + + // scalar on the left (friends implemented for + and *) + Vf c(3, 2.0f); + auto d = 5 + c; // friend operator+(U, Vector) + CHECK(d[0] == 7.0f && d[1] == 7.0f && d[2] == 7.0f, "scalar + vector failed"); + + d = 3 * c; // friend operator*(U, Vector) + CHECK(d[0] == 6.0f && d[1] == 6.0f && d[2] == 6.0f, "scalar * vector failed"); +} +// +// ---------- Vector arithmetic: +, -, *, / (elementwise) ---------- +// +TEST_CASE(Vector_Vector_Arithmetic) { + Vd a(3, 1.0), b(3, 2.0); + + // returning + auto c = a + b; + CHECK(c[0]==3.0 && c[1]==3.0 && c[2]==3.0, "vec + vec failed"); + + c = b - a; + CHECK(c[0]==1.0 && c[1]==1.0 && c[2]==1.0, "vec - vec failed"); + + c = a * b; + CHECK(c[0]==2.0 && c[1]==2.0 && c[2]==2.0, "vec * vec failed"); + + c = b / b; + CHECK(c[0]==1.0 && c[1]==1.0 && c[2]==1.0, "vec / vec failed"); + + // inplace + a = Vd(3, 1.0); + a += b; + CHECK(a[0]==3.0 && a[1]==3.0 && a[2]==3.0, "inplace vec + vec failed"); + a -= b; + CHECK(a[0]==1.0 && a[1]==1.0 && a[2]==1.0, "inplace vec - vec failed"); + a *= b; + CHECK(a[0]==2.0 && a[1]==2.0 && a[2]==2.0, "inplace vec * vec failed"); + a /= b; + CHECK(a[0]==1.0 && a[1]==1.0 && a[2]==1.0, "inplace vec / vec failed"); +} +// +// ---------- Size mismatch error paths ---------- +// +TEST_CASE(Vector_SizeMismatch_Throws) { + Vd a(3, 1.0), b(4, 2.0); + + bool threw = false; + try { auto c = a + b; (void)c; } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "add should throw on size mismatch"); + + threw = false; + try { a.inplace_subtract(b); } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "inplace_subtract should throw on size mismatch"); + + threw = false; + try { auto d = a * b; (void)d; } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "multiply should throw on size mismatch"); + + threw = false; + try { auto s = a.dot(b); (void)s; } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "dot should throw on size mismatch"); +} + +// +// ---------- Power / sqrt ---------- +// +TEST_CASE(Vector_Power_Sqrt) { + Vd a(3, 2.0); // [2,2,2] + auto b = a.power(3.0); // [8,8,8] + CHECK(b[0]==8.0 && b[1]==8.0 && b[2]==8.0, "scalar power failed"); + + Vd p(3, 3.0); // [3,3,3] + auto c = b.power(p); // 8^3 = 512 + CHECK(c[0]==512.0 && c[1]==512.0 && c[2]==512.0, "vector power failed"); + + Vd d(3, 9.0); + auto e = d.sqrt(); // [3,3,3] + CHECK(e[0]==3.0 && e[1]==3.0 && e[2]==3.0, "sqrt failed"); + + // inplace + d.inplace_sqrt(); // becomes [3,3,3] + CHECK(d == e, "inplace_sqrt failed"); +} + +// +// ---------- Dot / Sum / Norm / Normalize ---------- +// +TEST_CASE(Vector_Dot_Sum_Norm_Normalize) { + Vd a(3, 0.0); + a[0]=1.0; a[1]=2.0; a[2]=2.0; + + CHECK(a.sum() == 5.0, "sum failed"); + CHECK(a.dot(a) == 9.0, "dot self failed"); + + auto n = a.norm(); + CHECK(std::fabs(n - 3.0) < 1e-12, "norm failed"); + + auto b = a.normalize(); + CHECK(std::fabs(b.norm() - 1.0) < 1e-12, "normalize() not unit"); + + // inplace normalize + a.inplace_normalize(); + CHECK(std::fabs(a.norm() - 1.0) < 1e-12, "inplace_normalize not unit"); + + // zero-norm error + Vd z(3, 0.0); + bool threw = false; + try { z.inplace_normalize(); } catch (const std::runtime_error&) { threw = true; } + CHECK(threw, "normalize zero vector must throw"); +} +// +// ---------- Stream output (basic sanity) ---------- +// +TEST_CASE(Vector_StreamOutput) { + Vi a(3, 2); + std::ostringstream oss; + oss << a; + auto s = oss.str(); + CHECK(s == "[2, 2, 2]", "ostream<< wrong format"); +} \ No newline at end of file