From 78bd8c4a95049fcd738156a6244635c822044915 Mon Sep 17 00:00:00 2001 From: Elizabeth Hunt Date: Wed, 11 Oct 2023 15:56:20 -0600 Subject: [PATCH] lu factorization --- inc/lizfcm.h | 16 +++++--- lib/lizfcm.a | Bin 11512 -> 14496 bytes src/matrix.c | 62 +++++++++++++++++++++++++++++- src/vector.c | 23 ++++++++++- test/main.c | 106 +++++++++++++++++++++++++++++---------------------- 5 files changed, 153 insertions(+), 54 deletions(-) diff --git a/inc/lizfcm.h b/inc/lizfcm.h index 2dc9e42..99acbf2 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -12,22 +12,26 @@ extern double forward_derivative_at(double (*f)(double), double a, double h); extern double backward_derivative_at(double (*f)(double), double a, double h); extern double sum_v(Array_double *v); +extern Array_double *scale_v(Array_double *v1, double m); extern Array_double *add_v(Array_double *v1, Array_double *v2); extern Array_double *minus_v(Array_double *v1, Array_double *v2); extern double dot_v(Array_double *v1, Array_double *v2); extern double l2_norm(Array_double *v); extern double l1_norm(Array_double *v); extern double linf_norm(Array_double *v); - extern double l2_distance(Array_double *v1, Array_double *v2); extern double l1_distance(Array_double *v1, Array_double *v2); extern double linf_distance(Array_double *v1, Array_double *v2); - +extern Array_double *copy_vector(Array_double *v1); +extern void free_vector(Array_double *v); extern void format_vector_into(Array_double *v, char *s); - -extern void format_matrix_into(Matrix_double *m, char *s); -extern void put_identity_diagonal(Matrix_double *m); - extern Line *least_squares_lin_reg(Array_double *x, Array_double *y); +extern void put_identity_diagonal(Matrix_double *m); +extern Matrix_double *copy_matrix(Matrix_double *m); +extern void free_matrix(Matrix_double *m); +extern Matrix_double **put_lu_decomp(Matrix_double *m); + +extern void format_matrix_into(Matrix_double *m, char *s); + #endif // LIZFCM_H diff --git a/lib/lizfcm.a b/lib/lizfcm.a index 7a21e875b0aa821425c37f48edfa9bfdb51bf990..969c89475c1ecb7220497fad381d5da52f74143a 100644 GIT binary patch literal 14496 zcmd5@U2I&(bslm>awXXo5wnqNOVpK}z;vpWNk*3FR9D%=QXm7HN>CI)1uWfNFPG*f z%O$hBl10lkh|muznif)^2O}upmD(bt;|C*X00poRA{G)o7^&T=fgS>MmALKFb<)@g zP_pHg{l1wqmv=8OMbb`z40!juXMWE7oik_V+~LmsnNoK6z@~@#9_{aSmO|g&$M^N_ z?fdk;y^ix(Zyyp&a$j%n-aUH|_c#)CNN@7I-Cz3CFCIMf1?NjgjvYUAP$9g9=MiFQ zJfzf*g#Ks26TtI#D^(NvfbjjY;H<^1 zOhtLwT(MHh6vEL%qot_agCm)2Zmg`lvV9rJJ%1 zubj;kz#K29M+k*nrd;vLUpt*C<;q?mU-U}3XQZiUSSf$jZ=)fDZL}H2PFKAAAciVm z8Att@XGV*eg7O$$IE}H(j*g6Bw3VVf4v2&ie0hYBM#*E}040?ny(u#`RvJAUn2Ft^ znkRq$;lF?RzjrBxnKVBefLBR5co6;&&wx(VpjY0+|`M8 zqj#I}eiRQwTpsdgyd7G+t~o3|(eifV+=l5!tmR^La>raXwqxP)WcNaK^3+1yNq5iw z_~rZSv5l|9Taq8tx4wM;;TE+I^>m?v3wSQ#{lZS*6*X}T?YNM|JK`@fPc`O&*Qo{7 zlj*MAs~SE%0D6mROig-owR@8d)$-Y&);0^hMd%LIsBKO*Zg+r2{@cj+bENGYc)DBN zm+5YtJ2;KDUz~|~KmA?x;lb9m_2Kui?N?&aZU6rq_vod2`6`lC!_EPxS4FSsD;S25 z1Yb9FNIR91c=h@U=@Sq{hLzH9A|1Mucnb*mO}%N*y$*y9re+MfuLGfqDXw`y_dX!Y zbpct9wi)Qgfh@Pdpi8Ahqf-lL1JLDC#&XnQmP;9Q4+2>(Y0%vdMA6iML6`BM+Z{oG zrfHk+L4pljGwuN+yRLml=wpaucoPrPhlFN3NEd|ubD_U3H1_~Oelw)e@T#PL4fGb& zKOt=(khbA*q3t|<7jz&Q8-r*mDtMkZbh?GK*1Fxvf%B8 zDS?IZGUalvRPjocqAFBI#`^VdpE4PC;{}nuW&~i5`}&Q6)xW+y#?~;(9uM?llc+;_ zliA`JBOGOaxvx@=sfYHrZeAmQ9U1iVf|j=SSEh5#_O~Xrke@B9qn6v>JK9UjCn@}f z!3^T1rhvk?iT!H!@+xve&(s{y0YcY&dr3z2VV(~Y7N1aWtIplDU)O4z@oY(6yS+Ji zP4(>Tu63wu7>=%52jUQ~ZeD&N?l|3b*ig)swwthN*IGIq*f8XeDF6QA;#THCJ;t_8 zvc96Wm&o5heGSyt&^8q|QtE49kW| zW0lhJAT>J6(;cLtP1vm%q+ABqr+Z0b5^6vO924|D#?OwaeGB_3GEY zy+QrT7wr}|O24?)D+dq4YWrK;e(pyZEN%6RGR?5o{p}u;-QV9H#+Ydy%k}Gy`&-LL z?~TFpE-FAEgh;gS&imV=Qq~ddk#?nl()A{mtLZECox*Yhh@-6v>MM1W5Bb98TK;%j zdw^VHf7>$MfIh{sN3CyG-~YiIJFbGm1w4AM3Vr*2K&QNJ|#3|&T z8_JICq?k4Z8P{dYh=W{cYhoJ^?M~4i16=^I{Xt+*AJ{kQgPkXq3(_AajfuD8&OLw7 zzTN|xJp5p-lsOlV;pkHtYPM3^IIUN2_Gs?~4O4oPIZCf=m*VJ;a!iY;H`eBunC9OZ z8Cl8^gQczB*a2Rv-dyDXA&?(g)^T|1W5hDODda0TI3BpcLI|G=Egxa0YG|j`Vdl5( z;yc4?=gb*ll+?M=u8>4Y_a#;HyQAGy8%@?3nfx9o4%;^dfuKwvG_hPw4{Ip21BKq| zRO&Ez2KA69QaQkMt9R6+v_p#>Y97^e)jN{&vs;q$kF+J{<1HKi zqrMMk%vNJ7>xjp>bAA?i)Oda&NVkCgIqhKlb34e-@`|^qi*@jR*Vxg9n#puqySC{hb+QdDGilguvA+@Z{!;1<*NN%BpjwhR)84u7!?EbTYn$(6y6fM5`TqC{ZL7sp z19Q{G_VJx#Qj3{7!{c*1+|H06;Bsuo9M9)*%*%Alo3ryX+p$J}K*V`ewNCsa1mzOW zr1O^d?7hkPo;EdK2eK~p{-6C5^N_~EHopj&(QSt;Ft2R~TmH$bx; zzD!__e3?K$dngm^Z%`&Uzmy5gp`KsNV|eZGWx|dX@-)c=Xx%>W&@zF1Trb$KDad@& zIAE@WGGWKTmkIQhYl1Hmnikm$jRSmNnExlH7vX2qKC94t+(vu+{-cey-nHu6*V`lf zKWm&{f4$|{$lOxCI3||%=Ni7esJYhbZFsJK$9=rgvAmDkp5PqRG;bmKwfzgv8`lBb zKW$IwYotAqb5%Vshtj`BlPpK{k2(cCnYVg{v!Hp$dc5?W{Qc%LxwL1-SZLjVo=cw$ z2KUR5=_Z|~?V!B3V87j}CccS%iSHM-G9GVH@7G~VdSq=|DL?f*XamA))0C{!EKfUS z=lET&72{u7z&dnsV~d){rX1_j^7e%7#J^+RY|_(3hW}i&pLb;Mn^kA}5(mop=W>qg zI%oRLyE31l%!enH8yTyNt524lPY&&ZIGFGvA>!Pui1HJXexB(OTfDe%#JokOck+OS@Q#6H_5W7VYb&JR)E)!uJy%HI#sdMu zO6h~Eq`xZZy@t)Q_WexKf3`yUrMOZ!n#mKhoeY zwE+G=_f6n-@VPGdlHfOhTM<72#F;U1RN@Z;;WJ4LNPJx47cgEp(k|ba%MSxte$b%%6p-al7<7*VQ8abbpqoTu zH_B}jdJV8#408qM-Tz)}Z?yko@Zg-FJXsk*XPV#gCp+pux(eoa3*3 zLNG0u6m$fYU=8z$*Y+8~s^ElRS}-Z-2m*Au$5A+VUjGznG{l|{KO3SSL^+<WH^L} z{sXJ8gF;*VMv_jq`uc*MnTm7vH zZPz>c513xX!>}I@X}exMCbV7Oz6n~g;!UY*A)4Pmg6aPgqV2VTz~=zDM_?^LAjL-6 zPyPO3TnD+<=&)A*(m;S8mI0Lp1Gt#Dh)ob|Sj`_~`S6Yo zgOy)H>9x+jYFy1PJ1p&Tq*M~^V&7))BS~zoA6D}R-9S0|fOvml0vPlMVg08K(0nb&uzO9O_3c@=8pGWN>=V)8OSP@$41C$ryrVGb&0WT$ z({V?E_aM}*%nSU!awT?Z;cxJKKA!DX=QG;xq^wT^`()l_;C+Xk$j`jL*qUtEyxM<* zdkZ{MVc;4yq*=AeZTiO=%5|_ zW6~&S^k{ADt>SY8Hu)SulkXyY*5pecGV|IoZt_F3Ed6KnT?hTppyTw3f`|4wB3~W;PWm@> z8&T#Jw6|ycQ-0f-e+qs13C4_l^ZnNJQDVMzTQL6GHy`Mu^@sRAN#uc#Qrk1@8bLx1)SWxA zzqpGy=gS1f!j}n-RZu2Kw_z@UO)>#m%f#e%Gp>^-;Ul`UOjKjr1K&MP3g46ZJBTl9 zv6I4goA7-Kd`~v{im}hs-zKII&+(m}d}=Nz6H}n|98YeWqf9X0RBRg_@JUv@xF>h}AFoqy!nA|>osKHU9^ERp!uRs z|69|tw#jz~#&a$(Uf;u9=rOZni1GeAiDSL94}mw$+7>hGAosQ$gC?7>USG17>5t$`4(P;`e0x=ne5b@5`fH~4C+B_N zKG(!o@bDYLEPTxH3+ZoGQ9fon*Ln=VN6!z||9`<;1bxlXm;7ZnmsKN)1EtYt%XB#7 zH9J};tLFRv5-S;)R~syx8<74}oJkV&o6--O0)lnwutE1l^pNS#1DVcq zH&~}C23`7B*-rj{LMPKPZlsS3I)d+ETp3>jVxO6~3B*F3cngTCQg0e`uLJRynlb3U z4g~9z+@Ew0f~J4D7s!4MLU_qHDM-H}>020o;t3$@k#k|{C}=dvJq!d3{ZA4=dyhG_ zAMs$}P8e)Y3nm2}K>+i$pTJM>{P<(U^L%K}0lecw`dthv!x82O{Trb(h$G!D`+W9; zbeGV5LO(1t)=bgNdbzXtint*||Kio;J8fPJ1E=uj1p5zRrZ)S}RuTUs)_m5F tUYtXLcOb2zDv8LD>tG zx-JS$GZ$zF=Y?;I-b0OAFJ0Cu>4w%p^;y-$NBf6f>=};X8R~r{*89?MKRpO~D(d`O zo3P)1#lLNPQ0F29CO*%4BiCjQpUgX$EyMQ<0@+?e`#)9k>mmB9sIyo{IxG#?FN^M= zXPvvs&gGhv)okg?SE`Xdb_Oa }uc)TaoLZH#F!*F{CeKIf>oM8xGwN^II)JtNv~ z6oVx#6R)mODxcmh4itD$vMN~S)*>JJwb*ULs>(y(<_9PAlGn3owzNaDQEl08S_9>| zg0#N8UAsmJPszkc`F@SIy9=qRayo~jQ>yZzHb9?ka629v3Z>?epxZHHW3y6i%rp<} ztM+Q${C12#W?F8uCEks6%dzR_fO$$qyG#r#(J#UqJBX%D49@`Be;SC|M5jy)p9fYz zKLRWRz6x{!qoR)ho#4%)&oTA!#xxQwf#Dv|m;givH?{%si#}yycqb4$jD}1MHv^Gl zR5vla2F6OvOW`Xq-Uf1UMIYvS-g_>aQph5sZxE!;|AfjJ6b_32LVLL_kpi?Cp5aI2jT3OC#48KrKB!h@o* zB=B&6a|ToDJK`AD|vDipQG80Zkk`$ zLN&Esx*4pcQ?(!9yU_PVF1ak8sM}#PuGGMPUl#sgh^S$I@c^c=%*lEKS@dH|*8vd5UF)Fkbz-(9Usn22cY@aN5?fmwVJIPB<4zPIf(m{a1!DdUr7 z#O4dZMnB{S=*;Hpn{Bg;{iVsJ(`Ct}R$BrcF>2Q$8Ggpk)}-evRpyRcW#-+H3=Q~8 zt5TIPI@m$|naaIX(V(dyaxOI-#{A`bOH_t=sw9$GfbS;5*YQUx2;Xee*NA}y7T}{! z^``IKRf#)=DucYvrwSv+V!yj7uc@}x3WePyQV4x_dD)#jdpNx^PrHWj?n@QNl^LGH zES=t;Tdata[y]->data[y] = 1.0; } +Matrix_double *copy_matrix(Matrix_double *m) { + Matrix_double *copy = InitMatrixWithSize(double, m->rows, m->cols, 0.0); + for (size_t y = 0; y < copy->rows; y++) { + free_vector(copy->data[y]); + copy->data[y] = copy_vector(m->data[y]); + } + return copy; +} + +Matrix_double **put_lu_decomp(Matrix_double *m) { + assert(m->cols == m->rows); + + Matrix_double *u = copy_matrix(m); + Matrix_double *l = InitMatrixWithSize(double, m->rows, m->cols, 0.0); + put_identity_diagonal(l); + + Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2); + + for (size_t y = 0; y < m->rows; y++) { + if (u->data[y]->data[y] == 0) { + printf("ERROR: a pivot is zero in given matrix\n"); + exit(-1); + } + } + + if (u && l) { + for (size_t x = 0; x < m->cols; x++) { + for (size_t y = x + 1; y < m->rows; y++) { + double denom = u->data[x]->data[x]; + + if (denom == 0) { + printf("ERROR: non-factorable matrix\n"); + exit(-1); + } + + double factor = -(u->data[y]->data[x] / denom); + + Array_double *scaled = scale_v(u->data[x], factor); + Array_double *added = add_v(scaled, u->data[y]); + free_vector(scaled); + free_vector(u->data[y]); + + u->data[y] = added; + l->data[y]->data[x] = -factor; + } + } + } + + u_l[0] = u; + u_l[1] = l; + return u_l; +} + +void free_matrix(Matrix_double *m) { + for (size_t y = 0; y < m->rows; ++y) + free_vector(m->data[y]); + free(m); +} + void format_matrix_into(Matrix_double *m, char *s) { sprintf(s, ""); if (m->rows == 0) sprintf(s, "empty"); for (size_t y = 0; y < m->rows; ++y) { - char row_s[256]; + char *row_s = malloc(sizeof(char) * 256); format_vector_into(m->data[y], row_s); sprintf(s, "%s %s \n", s, row_s); + free(row_s); } } diff --git a/src/vector.c b/src/vector.c index 61692e1..e6e2544 100644 --- a/src/vector.c +++ b/src/vector.c @@ -41,12 +41,19 @@ double sum_v(Array_double *v) { return sum; } +Array_double *scale_v(Array_double *v, double m) { + Array_double *copy = copy_vector(v); + for (size_t i = 0; i < v->size; i++) + copy->data[i] *= m; + return copy; +} + Array_double *add_v(Array_double *v1, Array_double *v2) { assert(v1->size == v2->size); - Array_double *sum = InitArrayWithSize(double, v1->size, 0); + Array_double *sum = copy_vector(v1); for (size_t i = 0; i < v1->size; i++) - sum->data[i] = v1->data[i] + v2->data[i]; + sum->data[i] += v2->data[i]; return sum; } @@ -80,6 +87,18 @@ double linf_distance(Array_double *v1, Array_double *v2) { return dist; } +Array_double *copy_vector(Array_double *v) { + Array_double *copy = InitArrayWithSize(double, v->size, 0.0); + for (size_t i = 0; i < copy->size; ++i) + copy->data[i] = v->data[i]; + return copy; +} + +void free_vector(Array_double *v) { + free(v->data); + free(v); +} + void format_vector_into(Array_double *v, char *s) { sprintf(s, ""); if (v->size == 0) diff --git a/test/main.c b/test/main.c index 798bfbf..3551af2 100644 --- a/test/main.c +++ b/test/main.c @@ -5,56 +5,72 @@ double f(double x) { return (x - 1) / (x + 1); } int main() { - printf("smaceps(): %.10e\n", smaceps()); - printf("dmaceps(): %.10e\n", dmaceps()); - - Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3}); - Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5}); - - char v_s[256]; - char w_s[256]; - format_vector_into(v, v_s); - format_vector_into(w, w_s); - - printf("v: %s\n", v_s); - printf("w: %s\n", w_s); - printf("l1_norm(v): %f\n", l1_norm(v)); - printf("l2_norm(v): %f\n", l2_norm(v)); - printf("linf_norm(v): %f\n", linf_norm(v)); - - printf("l1_dist(v, w): %f\n", l1_distance(v, w)); - printf("l2_dist(v, w): %f\n", l2_distance(v, w)); - printf("linf_dist(v, w): %f\n", linf_distance(v, w)); - - double h = 0.001; - printf("approx f'(1) w/ c.d.: %f\n", central_derivative_at(&f, 1, h)); - printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h)); - printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h)); - - v = InitArray(double, {1, 2, 3, 4, 5}); - w = InitArray(double, {2, 3, 4, 5, 6}); - format_vector_into(v, v_s); - format_vector_into(w, w_s); - printf("v: %s\n", v_s); - printf("w: %s\n", w_s); - - Line *line = least_squares_lin_reg(v, w); - printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a); - - v = InitArray(double, {1, 2, 3, 4, 5, 6, 7}); - w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5}); - format_vector_into(v, v_s); - format_vector_into(w, w_s); - printf("v: %s\n", v_s); - printf("w: %s\n", w_s); - line = least_squares_lin_reg(v, w); - printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a); + // printf("smaceps(): %.10e\n", smaceps()); + // printf("dmaceps(): %.10e\n", dmaceps()); + // + // Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3}); + // Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5}); + // + // char v_s[256]; + // char w_s[256]; + // format_vector_into(v, v_s); + // format_vector_into(w, w_s); + // + // printf("v: %s\n", v_s); + // printf("w: %s\n", w_s); + // printf("l1_norm(v): %f\n", l1_norm(v)); + // printf("l2_norm(v): %f\n", l2_norm(v)); + // printf("linf_norm(v): %f\n", linf_norm(v)); + // + // printf("l1_dist(v, w): %f\n", l1_distance(v, w)); + // printf("l2_dist(v, w): %f\n", l2_distance(v, w)); + // printf("linf_dist(v, w): %f\n", linf_distance(v, w)); + // + // double h = 0.001; + // printf("approx f'(1) w/ c.d.: %f\n", central_derivative_at(&f, 1, h)); + // printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h)); + // printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h)); + // + // v = InitArray(double, {1, 2, 3, 4, 5}); + // w = InitArray(double, {2, 3, 4, 5, 6}); + // format_vector_into(v, v_s); + // format_vector_into(w, w_s); + // printf("v: %s\n", v_s); + // printf("w: %s\n", w_s); + // + // Line *line = least_squares_lin_reg(v, w); + // printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a); + // + // v = InitArray(double, {1, 2, 3, 4, 5, 6, 7}); + // w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5}); + // + // format_vector_into(v, v_s); + // format_vector_into(w, w_s); + // printf("v: %s\n", v_s); + // printf("w: %s\n", w_s); + // line = least_squares_lin_reg(v, w); + // printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a); char m_s[256]; Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0); put_identity_diagonal(m); + for (int i = 0; i < 8; i++) { + for (int j = 0; j < 8; j++) { + m->data[i]->data[j] = (i + 1.0) + j * 3 + (rand() % 12); + } + } + format_matrix_into(m, m_s); - printf("%s", m_s); + printf("%s\n", m_s); + + Matrix_double **u_l = put_lu_decomp(m); + Matrix_double *u = u_l[0]; + Matrix_double *l = u_l[1]; + + format_matrix_into(u, m_s); + printf("%s\n", m_s); + format_matrix_into(l, m_s); + printf("%s\n", m_s); return 0; }