diff --git a/hws/hw3/README.md b/hws/hw3/README.md new file mode 100644 index 0000000..07d18fe --- /dev/null +++ b/hws/hw3/README.md @@ -0,0 +1,27 @@ +# Numerical Analysis Homework 3 + +## Compile Alert + +This source file contain UNIX `fork` function!!! please compile and run in Ubuntu. + +when compiling source, must be included below c file + +```makefile +nr/ansiother/nrutil.c nr/ansirecipes/gaussj.c nr/ansirecipes/ludcmp.c nr/ansirecipes/svdcmp.c nr/ansirecipes/mprove.c nr/ansirecipes/lubksb.c nr/ansirecipes/pythag.c main.c matutil.c +``` + +and include dir `nr/ansi/other` and `$(workspace)`(this folder) + +cflags are: +```makefile +CFLAGS := -Wno-unused-but-set-variable -Wno-unused-parameter -Wno-unused-variable -Wall -g -std=c90 +LDFLAGS := -lm +``` + +**executable should be located in the same folder as the *.dat file.** + +## Results & Report + +Compile, build and run results are included in hw3.txt + +The report is the output of a compilation executable. diff --git a/hws/hw3/hw3.txt b/hws/hw3/hw3.txt new file mode 100644 index 0000000..2d0674a --- /dev/null +++ b/hws/hw3/hw3.txt @@ -0,0 +1,294 @@ +================================= +data from lineq1.dat: +----------------------------- +matrix data A: + 4.0000000 2.0000000 3.0000000 -1.0000000 +-2.0000000 -1.0000000 -2.0000000 2.0000000 + 5.0000000 3.0000000 4.0000000 -1.0000000 +11.0000000 4.0000000 6.0000000 1.0000000 +matrix data b: + 4.0000000 +-3.0000000 + 4.0000000 +11.0000000 +----------------------------- +created process for gauss-jordan method +----------------------------- +gauss-jordan method: +Numerical Recipes run-time error... +gaussj: Singular Matrix +...now exiting to system... +created process for LU Decomposition +----------------------------- +LU Decomposition: +11.0000000 4.0000000 6.0000000 1.0000000 + 0.4545455 1.1818181 1.2727273 -1.4545455 +-0.1818182 -0.2307692 -0.6153846 1.8461540 + 0.3636364 0.4615384 -0.3750000 0.0000000 +index: 4 3 3 4 +solution x: 1.0000000 -3.0000000 2.0000000 0.0000000 +----------------------------- +created process for Singular Value Decomposition +----------------------------- +Singular Value Decompoisition: +U: +-0.3346812 -0.3418892 0.0040041 0.8781140 + 0.1847021 0.6724079 0.6366467 0.3292924 +-0.4362502 -0.4101597 0.7300825 -0.3292934 +-0.8145915 0.5125896 -0.2482825 -0.1097641 +w: 16.0844917 2.9560738 0.7420990 0.0000003 +S: +16.0844917 0.0000000 0.0000000 0.0000000 + 0.0000000 2.9560738 0.0000000 0.0000000 + 0.0000000 0.0000000 0.7420990 0.0000000 + 0.0000000 0.0000000 0.0000000 0.0000003 +V: +-0.7988990 0.2961076 -0.4554271 -0.2581990 +-0.3370440 -0.1814251 0.7660415 -0.5163977 +-0.4977463 -0.3164953 0.2282085 0.7745966 + 0.0202520 0.8827433 0.3920297 0.2581990 + +Orig = U * S * V^T + 4.0000000 2.0000005 3.0000007 -0.9999996 +-1.9999998 -0.9999998 -1.9999995 1.9999996 + 4.9999981 2.9999988 3.9999988 -0.9999996 +10.9999971 4.0000000 5.9999990 0.9999993 +we want to get x by using x = V * (S^-1 * (U^T * b)): +x: + 1.2935313 +-2.4129376 + 1.1194062 +-0.2935313 +----------------------------- +created process for LUDecomp with mprove +----------------------------- +mprove method: +x_init: 1.0000000 -3.0000000 2.0000000 0.0000000 +x using mprove: 1.0000000 -3.0000000 2.0000000 0.0000000 +----------------------------- +last we want to get det and inv of A: +created process for det and inv +----------------------------- +det: 0.000000 +There is no inverse matrix, because det ~ 0.0f +----------------------------- +last we want to get det and inv of A: +end for lineq1.dat +================================= +data from lineq2.dat: +----------------------------- +matrix data A: + 2.0000000 -4.0000000 -5.0000000 5.0000000 0.0000000 +-1.0000000 1.0000000 2.0000000 0.0000000 4.0000000 +-1.0000000 6.0000000 0.0000000 3.0000000 2.0000000 + 0.0000000 1.0000000 3.0000000 7.0000000 5.0000000 + 5.0000000 0.0000000 8.0000000 7.0000000 -2.0000000 +matrix data b: +-5.0000000 + 2.0000000 + 0.0000000 + 4.0000000 +-1.0000000 +----------------------------- +created process for gauss-jordan method +----------------------------- +gauss-jordan method: +-2.8735671 +-0.6123569 + 0.9762775 + 0.6358188 +-0.5534414 +----------------------------- +created process for LU Decomposition +----------------------------- +LU Decomposition: + 5.0000000 0.0000000 8.0000000 7.0000000 -2.0000000 +-0.2000000 6.0000000 1.6000000 4.4000001 1.6000000 + 0.4000000 -0.6666667 -7.1333332 5.1333332 1.8666668 + 0.0000000 0.1666667 -0.3831776 8.2336445 5.4485979 +-0.2000000 0.1666667 -0.4672897 0.3723042 2.1770713 +index: 5 3 5 4 5 +solution x: -2.8735664 -0.6123567 0.9762774 0.6358187 -0.5534411 +----------------------------- +created process for Singular Value Decomposition +----------------------------- +Singular Value Decompoisition: +U: +-0.0712607 -0.3099797 -0.1001519 0.6951792 0.6368123 +-0.1167440 -0.7382753 0.5349323 -0.3748873 0.1209441 +-0.2206404 -0.1961647 -0.7456251 -0.5069756 0.3160007 +-0.5691139 0.5318492 0.3635619 -0.1996495 0.4703279 +-0.7802050 -0.1936958 -0.1252318 0.2816048 -0.5087020 +w: 13.8655672 0.7772490 4.6621399 9.1695461 8.3262081 +S: +13.8655672 0.0000000 0.0000000 0.0000000 0.0000000 + 0.0000000 0.7772490 0.0000000 0.0000000 0.0000000 + 0.0000000 0.0000000 4.6621399 0.0000000 0.0000000 + 0.0000000 0.0000000 0.0000000 9.1695461 0.0000000 + 0.0000000 0.0000000 0.0000000 0.0000000 8.3262081 +V: +-0.2672927 -0.8414277 -0.1320789 0.4013552 -0.2049950 +-0.1243843 -0.1846177 -0.6809420 -0.6976467 -0.0072029 +-0.5644318 0.1535274 0.3559430 -0.2804697 -0.6726719 +-0.7546358 0.2942181 -0.2293633 0.2757668 0.4640102 +-0.1581916 -0.3844222 0.5827267 -0.4444011 0.5386392 + +Orig = U * S * V^T + 1.9999996 -3.9999971 -4.9999990 5.0000005 0.0000007 +-0.9999999 0.9999996 1.9999998 -0.0000004 4.0000000 +-0.9999992 5.9999981 -0.0000005 3.0000002 1.9999995 + 0.0000012 1.0000011 2.9999993 6.9999962 5.0000000 + 5.0000024 0.0000014 8.0000000 6.9999962 -1.9999998 +we want to get x by using x = V * (S^-1 * (U^T * b)): +x: +-2.8735659 +-0.6123567 + 0.9762776 + 0.6358188 +-0.5534413 +----------------------------- +created process for LUDecomp with mprove +----------------------------- +mprove method: +x_init: -2.8735664 -0.6123567 0.9762774 0.6358187 -0.5534411 +x using mprove: -2.8735662 -0.6123565 0.9762774 0.6358185 -0.5534410 +----------------------------- +last we want to get det and inv of A: +created process for det and inv +----------------------------- +det: 3836.000000 +inv: + 0.3545361 0.7669450 0.2077686 -0.5954121 0.2531283 + 0.0354536 0.1266945 0.1957769 -0.1595412 0.0503128 +-0.1386861 -0.0985402 -0.0967153 0.1240876 0.0164233 +-0.0521377 -0.3039626 -0.0232013 0.2346195 -0.0445777 + 0.1491137 0.4593328 0.0513556 -0.1710115 0.0424922 +----------------------------- +last we want to get det and inv of A: +end for lineq2.dat +================================= +data from lineq3.dat: +----------------------------- +matrix data A: + 0.4000000 8.1999998 6.6999998 1.9000000 2.2000000 5.3000002 + 7.8000002 8.3000002 7.6999998 3.3000000 1.9000000 4.8000002 + 5.5000000 8.8000002 3.0000000 1.0000000 5.0999999 6.4000001 + 5.0999999 5.0999999 3.5999999 5.8000002 5.6999998 4.9000001 + 3.5000000 2.7000000 5.6999998 8.1999998 9.6000004 2.9000001 + 3.0000000 5.3000002 5.5999999 3.5000000 6.8000002 5.6999998 +matrix data b: +-2.9000001 +-8.1999998 + 7.6999998 +-1.0000000 + 5.6999998 + 3.0000000 +----------------------------- +created process for gauss-jordan method +----------------------------- +gauss-jordan method: +-0.3266082 + 1.5322930 +-1.0448254 +-1.5874470 + 2.9284797 +-2.2189310 +----------------------------- +created process for LU Decomposition +----------------------------- +LU Decomposition: + 7.8000002 8.3000002 7.6999998 3.3000000 1.9000000 4.8000002 + 0.0512820 7.7743587 6.3051281 1.7307692 2.1025641 5.0538464 + 0.7051282 0.3791227 -4.8199039 -1.9830968 2.9631264 1.0993568 + 0.6538461 -0.0420514 0.2426345 4.1962576 3.8271508 1.7073183 + 0.3846154 0.2711082 -0.1927610 0.3286928 4.8124266 2.1344366 + 0.4487179 -0.1317612 -0.6381130 1.3540252 1.1913373 -2.7410173 +index: 2 2 3 4 6 6 +solution x: -0.3266079 1.5322925 -1.0448257 -1.5874470 2.9284797 -2.2189300 +----------------------------- +created process for Singular Value Decomposition +----------------------------- +Singular Value Decompoisition: +U: +-0.3537291 0.3737144 0.7512395 0.2092905 0.3221614 0.1525169 +-0.4576060 0.3859336 -0.4780837 0.6014348 -0.2266329 0.0012832 +-0.4140943 0.3215108 -0.2003601 -0.7149711 -0.0887913 0.4073634 +-0.3948520 -0.2070790 -0.2807029 -0.1138602 0.7322581 -0.4162332 +-0.4182500 -0.7414790 0.0732381 0.1737806 -0.1087733 0.4773684 +-0.4039270 -0.1238982 0.2877024 -0.2003150 -0.5375242 -0.6400477 +w: 30.6345253 9.9541349 5.3502665 4.8673272 1.8196821 1.1195914 +S: +30.6345253 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 + 0.0000000 9.9541349 0.0000000 0.0000000 0.0000000 0.0000000 + 0.0000000 0.0000000 5.3502665 0.0000000 0.0000000 0.0000000 + 0.0000000 0.0000000 0.0000000 4.8673272 0.0000000 0.0000000 + 0.0000000 0.0000000 0.0000000 0.0000000 1.8196821 0.0000000 + 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.1195914 +V: +-0.3485524 0.0909273 -0.9051291 0.0553008 -0.2121177 -0.0541523 +-0.5100967 0.5407045 0.1345508 -0.1554864 0.3139328 0.5537202 +-0.4309946 0.0777932 0.3306456 0.6877020 -0.4654545 -0.0963503 +-0.3176094 -0.5434607 -0.0693887 0.3556214 0.6865202 -0.0344047 +-0.4169183 -0.5973302 0.1461581 -0.4902155 -0.3847928 0.2442064 +-0.3973150 0.2028947 0.1647232 -0.3647657 0.1429302 -0.7876277 + +Orig = U * S * V^T + 0.3999988 8.1999960 6.6999979 1.8999976 2.1999958 5.2999964 + 7.7999992 8.2999983 7.6999993 3.3000000 1.9000000 4.7999992 + 5.4999986 8.7999973 2.9999979 0.9999974 5.0999985 6.3999977 + 5.0999994 5.0999970 3.5999978 5.8000007 5.6999979 4.8999968 + 3.4999995 2.6999977 5.6999984 8.1999989 9.5999975 2.8999982 + 2.9999990 5.2999978 5.5999980 3.4999983 6.7999973 5.6999974 +we want to get x by using x = V * (S^-1 * (U^T * b)): +x: +-0.3266090 + 1.5322922 +-1.0448247 +-1.5874470 + 2.9284794 +-2.2189293 +----------------------------- +created process for LUDecomp with mprove +----------------------------- +mprove method: +x_init: -0.3266079 1.5322925 -1.0448257 -1.5874470 2.9284797 -2.2189300 +x using mprove: -0.3266079 1.5322924 -1.0448252 -1.5874476 2.9284801 -2.2189305 +----------------------------- +last we want to get det and inv of A: +created process for det and inv +----------------------------- +det: 16178.396484 +inv: +-0.1622052 0.1228010 0.0240679 -0.0164308 -0.0228398 0.0461325 + 0.1694071 -0.0411167 0.2283134 -0.0876241 0.1803061 -0.3956549 +-0.0116364 0.1227448 -0.1174069 -0.1809809 0.0159104 0.1867664 + 0.1056689 -0.0517256 -0.1089164 0.2997742 0.0008588 -0.1905406 +-0.0530261 -0.0423615 0.1605081 -0.2240342 0.1618107 0.0150243 +-0.0623407 -0.0646943 -0.2342164 0.3511257 -0.3648281 0.4346332 +----------------------------- +last we want to get det and inv of A: +end for lineq3.dat +======================================= +===============< Report >============== +======================================= +# 1. discuss three method. +In terms of the accuracy of the solution, +there is no noticable difference between three methods.(at 1e-6) +But `gaussj` method is unstable when there are many solutions. +In terms of stability, svdcmp is the most stable +because the solution can always come out (because it is, by definition, expressed as matmul). +I could not feel it in this example +but I think svd is the slowest for that; matmul is too expensive. +so I think that LUdcmp is balanced: +Not as unstable as gaussj, not as expensive as svdcmp. +# 2. mprove +`mprove` is a method of increasing the accuracy of the solution +by using an iterative method using the alud and initial solution obtained by `ludcmp`. +The accuracy of the LU can be compensated using this method, +so the LUdcmp is getting more useful. +# 3. det and inv +Without first example, we can get inverse matrix with gaussj. +We can also get determinants using cofactor expansion with recursive method. +======================================= +=============< End Report >============ +======================================= diff --git a/hws/hw3/main.c b/hws/hw3/main.c index 4bff324..5657978 100644 --- a/hws/hw3/main.c +++ b/hws/hw3/main.c @@ -4,52 +4,160 @@ #include #include #include +#include #include "matutil.h" -void try_gaussj(JMatrixData* data) { +void try_gaussj(JMatrixData *data) { printf("gauss-jordan method:\n"); + fflush(stdout); + fflush(stderr); gaussj(data->a, data->m, data->b, 1); print_matrix(data->m, 1, data->b); } -void try_lu(JMatrixData* data) { +void try_lu(JMatrixData *data) { printf("LU Decomposition:\n"); - - int* indx = calloc(data->m + 1, sizeof(int)); - float* d = calloc(data->m + 1, sizeof(float)); + + int *indx = calloc(data->m + 1, sizeof(int)); + float *d = calloc(data->m + 1, sizeof(float)); ludcmp(data->a, data->m, indx, d); print_matrix(data->m, data->n, data->a); - + printf("index: "); print_vector_int(data->m, indx); - printf("d: "); - print_vector_float(data->m, d); - printf("solution:\n"); - float * b = calloc(data->m + 1, sizeof(float)); + printf("solution x: "); + float *b = calloc(data->m + 1, sizeof(float)); int i; for (i = 1; i <= data->m; i++) { b[i] = data->b[i][1]; } lubksb(data->a, data->m, indx, b); - for (i = 1; i <= data->m; i++) { - printf("%6.2f ", b[i]); - } + print_vector_float(data->m, b); free(indx); free(d); - + free(b); } -void processMatrix(JMatrixData* data) { +void try_svd(JMatrixData *data) { + int i; + printf("Singular Value Decompoisition:\n"); + fflush(stdout); + fflush(stderr); + float *w = calloc(data->m + 1, sizeof(float)); + float **v = new_matrix(data->m, data->n); + svdcmp(data->a, data->m, data->n, w, v); + printf("U:\n"); + print_matrix(data->m, data->n, data->a); + printf("w: "); + print_vector_float(data->m, w); + printf("S:\n"); + float **s = new_diagonal(data->m, w); + print_matrix(data->m, data->n, s); + printf("V:\n"); + print_matrix(data->m, data->n, v); + printf("\n"); + + float **t1 = matmul(data->m, data->m, data->m, data->a, s); + float **vt = mattranspose(data->m, data->n, v); + float **t2 = matmul(data->m, data->n, data->m, t1, vt); + + printf("Orig = U * S * V^T\n"); + print_matrix(data->m, data->n, t2); + + printf("we want to get x by using x = V * (S^-1 * (U^T * b)):\n"); + printf("x:\n"); + + float **invS = new_matrix(data->m, data->m); + for (i = 1; i <= data->m; i++) { + invS[i][i] = 1 / w[i]; + } + + float **UT = mattranspose(data->m, data->n, data->a); + float **l1 = matmul(data->m, data->m, 1, UT, data->b); + float **l2 = matmul(data->m, data->m, 1, invS, l1); + float **l3 = matmul(data->m, data->m, 1, v, l2); + + print_matrix(data->m, 1, l3); + + simple_free_matrix(data->m, data->n, v); + simple_free_matrix(data->m, data->m, vt); + simple_free_matrix(data->m, data->m, s); + simple_free_matrix(data->m, data->n, t1); + simple_free_matrix(data->m, data->n, t2); + simple_free_matrix(data->m, data->m, invS); + simple_free_matrix(data->m, data->m, UT); + simple_free_matrix(data->m, 1, l1); + simple_free_matrix(data->m, 1, l2); + simple_free_matrix(data->m, 1, l3); + + free(w); +} + +void try_det_inv(JMatrixData *data) { + int m = data->m; + + float det = matdet(m, data->a); + printf("det: %f\n", det); + + if (det < 1e-6) { + printf("There is no inverse matrix, because det ~ 0.0f\n"); + } else { + float **inv = matinv(m, data->a); + printf("inv:\n"); + print_matrix(m, m, inv); + simple_free_matrix(m, m, inv); + } +} + +void try_mprove(JMatrixData *data) { + printf("mprove method:\n"); + fflush(stdout); + fflush(stderr); + float **a = new_matrix(data->m, data->m); + copy_matrix(data->m, data->m, data->a, a); + + int *indx = calloc(data->m + 1, sizeof(int)); + float *d = malloc(sizeof(float)); + ludcmp(a, data->m, indx, d); + float *x_init = calloc(data->m + 1, sizeof(float)); + float *b = calloc(data->m + 1, sizeof(float)); + int i; + for (i = 1; i <= data->m; i++) { + x_init[i] = data->b[i][1]; + b[i] = data->b[i][1]; + } + + lubksb(a, data->m, indx, x_init); + float *x = calloc(data->m + 1, sizeof(float)); + for (i = 1; i <= data->m; i++) { + x[i] = x_init[i]; + } + printf("x_init: "); + print_vector_float(data->m, x_init); + + mprove(data->a, a, data->m, indx, b, x); + printf("x using mprove: "); + print_vector_float(data->m, x); + + free(x); + free(x_init); + free(indx); + free(d); + simple_free_matrix(data->m, data->m, a); +} + +void processMatrix(JMatrixData *data) { printf("-----------------------------\n"); printf("matrix data A:\n"); print_matrix(data->m, data->n, data->a); printf("matrix data b:\n"); print_matrix(data->m, 1, data->b); printf("-----------------------------\n"); - + fflush(stdout); + fflush(stderr); pid_t pid = fork(); if (pid < 0) { printf("process 실행 실패\n"); @@ -57,7 +165,7 @@ void processMatrix(JMatrixData* data) { } else if (pid == 0) { printf("created process for gauss-jordan method\n"); printf("-----------------------------\n"); - JMatrixData* copied = new_jmatdata(data->m, data->n); + JMatrixData *copied = new_jmatdata(data->m, data->n); copy_jmatdata(data, copied); try_gaussj(copied); free_jmatdata(copied); @@ -66,6 +174,9 @@ void processMatrix(JMatrixData* data) { } wait(NULL); + + fflush(stdout); + fflush(stderr); pid = fork(); if (pid < 0) { printf("process 실행 실패\n"); @@ -73,18 +184,88 @@ void processMatrix(JMatrixData* data) { } else if (pid == 0) { printf("created process for LU Decomposition\n"); printf("-----------------------------\n"); - JMatrixData* copied = new_jmatdata(data->m, data->n); + fflush(stdout); + fflush(stderr); + JMatrixData *copied = new_jmatdata(data->m, data->n); copy_jmatdata(data, copied); try_lu(copied); free_jmatdata(copied); printf("-----------------------------\n"); exit(0); - } wait(NULL); + + fflush(stdout); + fflush(stderr); + + pid = fork(); + if (pid < 0) { + printf("process 실행 실패\n"); + return; + } else if (pid == 0) { + printf("created process for Singular Value Decomposition\n"); + printf("-----------------------------\n"); + fflush(stdout); + fflush(stderr); + JMatrixData *copied = new_jmatdata(data->m, data->n); + copy_jmatdata(data, copied); + try_svd(copied); + free_jmatdata(copied); + printf("-----------------------------\n"); + exit(0); + } + + wait(NULL); + + fflush(stdout); + fflush(stderr); + + pid = fork(); + if (pid < 0) { + printf("process 실행 실패\n"); + return; + } else if (pid == 0) { + printf("created process for LUDecomp with mprove\n"); + printf("-----------------------------\n"); + fflush(stdout); + fflush(stderr); + JMatrixData *copied = new_jmatdata(data->m, data->n); + copy_jmatdata(data, copied); + try_mprove(copied); + free_jmatdata(copied); + printf("-----------------------------\n"); + exit(0); + } + + wait(NULL); + + fflush(stdout); + fflush(stderr); + + printf("last we want to get det and inv of A:\n"); + pid = fork(); + if (pid < 0) { + printf("process 실행 실패\n"); + return; + } else if (pid == 0) { + printf("created process for det and inv\n"); + printf("-----------------------------\n"); + fflush(stdout); + fflush(stderr); + JMatrixData *copied = new_jmatdata(data->m, data->n); + copy_jmatdata(data, copied); + try_det_inv(copied); + free_jmatdata(copied); + printf("-----------------------------\n"); + exit(0); + } + + wait(NULL); + fflush(stdout); + fflush(stderr); } -JMatrixData*readDataFrom(char *filename) { +JMatrixData *readDataFrom(char *filename) { FILE *fp = fopen(filename, "r"); if (fp == NULL) { @@ -96,7 +277,7 @@ JMatrixData*readDataFrom(char *filename) { fscanf(fp, "%d", &m); fscanf(fp, "%d", &n); - JMatrixData*data = new_jmatdata(m, n); + JMatrixData *data = new_jmatdata(m, n); int i, j; for (i = 1; i <= data->m; i++) { @@ -117,19 +298,20 @@ int main() { int m, n; int i, j; - const char *filenames[] = { "lineq1.dat", "lineq2.dat", "lineq3.dat"}; pid_t pid; for (i = 0; i < 3; i++) { pid = fork(); + fflush(stdout); + fflush(stderr); if (pid < 0) { printf("process 실행 실패\n"); return 2; } else if (pid == 0) { printf("=================================\n"); printf("data from %s:\n", filenames[i]); - JMatrixData* data = readDataFrom(filenames[i]); + JMatrixData *data = readDataFrom(filenames[i]); processMatrix(data); printf("end for %s\n", filenames[i]); exit(0); @@ -137,9 +319,27 @@ int main() { wait(NULL); } - - - printf("end\n"); + fflush(stdout); + fflush(stderr); + printf("=======================================\n"); + printf("===============< Report >==============\n"); + printf("=======================================\n"); + printf("# 1. discuss three method.\n"); + printf("In terms of the accuracy of the solution,\nthere is no noticable difference between three methods.(at 1e-6)\n"); + printf("But `gaussj` method is unstable when there are many solutions.\n"); + printf("In terms of stability, svdcmp is the most stable\nbecause the solution can always come out (because it is, by definition, expressed as matmul).\n"); + printf("I could not feel it in this example\nbut I think svd is the slowest for that; matmul is too expensive.\n"); + printf("so I think that LUdcmp is balanced:\nNot as unstable as gaussj, not as expensive as svdcmp.\n"); + printf("# 2. mprove\n"); + printf("`mprove` is a method of increasing the accuracy of the solution\nby using an iterative method using the alud and initial solution obtained by `ludcmp`.\n"); + printf("The accuracy of the LU can be compensated using this method,\nso the LUdcmp is getting more useful.\n"); + printf("# 3. det and inv\n"); + printf("Without first example, we can get inverse matrix with gaussj.\n"); + printf("We can also get determinants using cofactor expansion with recursive method.\n"); + printf("=======================================\n"); + printf("=============< End Report >============\n"); + printf("=======================================\n"); + return 0; } \ No newline at end of file diff --git a/hws/hw3/matutil.c b/hws/hw3/matutil.c index 9480e79..e039a0e 100644 --- a/hws/hw3/matutil.c +++ b/hws/hw3/matutil.c @@ -1,23 +1,42 @@ #include "matutil.h" +#define FLOAT_FORMAT "%10.7f " + void print_vector_float(int size, float *v) { int i; - for (i = 0; i < size; i++) { - printf("%6.2f ", v[i]); + for (i = 1; i <= size; i++) { + printf(FLOAT_FORMAT, v[i]); } printf("\n"); } void print_vector_int(int size, int *v) { int i; - for (i = 0; i < size; i++) { + for (i = 1; i <= size; i++) { printf("%d ", v[i]); } printf("\n"); } float **new_matrix(int m, int n) { - return matrix(1, m, 1, n); + float **mat = matrix(1, m, 1, n); + int i, j; + for (i = 1; i <= m; i++) { + for (j = 1; j <= n; j++) { + mat[i][j] = 0; + } + } + return mat; +} + +float **new_diagonal(int n, float *v) { + float **mat = new_matrix(n, n); + int i; + for (i = 1; i <= n; i++) { + mat[i][i] = v[i]; + } + + return mat; } void simple_free_matrix(int m, int n, float **mat) { @@ -28,7 +47,7 @@ void print_matrix(int m, int n, float **mat) { int i, j; for (i = 1; i <= m; i++) { for (j = 1; j <= n; j++) { - printf("%6.2f ", mat[i][j]); + printf(FLOAT_FORMAT, mat[i][j]); } printf("\n"); } @@ -43,6 +62,92 @@ void copy_matrix(int m, int n, float **src, float **dst) { } } +float **matmul(int m, int n, int p, float **a, float **b) { + float **c = new_matrix(m, p); + int i, j, k; + + for (i = 1; i <= m; i++) { + + for (j = 1; j <= p; j++) { + c[i][j] = 0; + + for (k = 1; k <= n; k++) { + + c[i][j] += a[i][k] * b[k][j]; + } + } + } + + return c; +} + +float **mattranspose(int m, int n, float **mat) { + float **t = new_matrix(n, m); + int i, j; + for (i = 1; i <= m; i++) { + for (j = 1; j <= n; j++) { + t[j][i] = mat[i][j]; + } + } + return t; +} + +float **matinv(int m, float **mat) { + float **inv = new_matrix(m, m); + int i, j; + for (i = 1; i <= m; i++) { + for (j = 1; j <= m; j++) { + if (i == j) { + inv[i][j] = 1; + } else { + inv[i][j] = 0; + } + } + } + gaussj(mat, m, inv, m); + return inv; +} + +void create_submatrix(float **matrix, float **sub_matrix, int n, int exclude_col) { + int sub_i = 1; + int i, j; + for (i = 2; i <= n; i++) { + int sub_j = 1; + for (j = 1; j <= n; j++) { + if (j == exclude_col) { + continue; + } + sub_matrix[sub_i][sub_j] = matrix[i][j]; + sub_j++; + } + sub_i++; + } +} + +float matdet(int m, float **mat) { + if (m == 1) { + return mat[1][1]; + } else if (m == 2) { + return mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]; + } + + float det = .0f; + + float **submat = new_matrix(m - 1, m - 1); + int i, j; + for(i = 1; i <= m; i ++) { + create_submatrix(mat, submat, m, i); + float sign = (i % 2 == 0) ? -1 : 1; + + float el = mat[1][i]; + + det += el * sign * matdet(m-1, submat); + } + + simple_free_matrix(m - 1, m - 1, submat); + return det; +} + /* for JMatrixData */ void copy_jmatdata(JMatrixData *src, JMatrixData *dst) { diff --git a/hws/hw3/matutil.h b/hws/hw3/matutil.h index 1fde058..9ea8101 100644 --- a/hws/hw3/matutil.h +++ b/hws/hw3/matutil.h @@ -11,15 +11,24 @@ void print_vector_float(int size, float* v); void print_vector_int(int size, int* v); - float **new_matrix(int m, int n); +float **new_diagonal(int n, float *v); + void print_matrix(int m, int n, float **mat); void copy_matrix(int m, int n, float **src, float **dst); void simple_free_matrix(int m, int n, float **mat); +float **matmul(int m, int n, int p, float **a, float **b); + +float **mattranspose(int m, int n, float **mat); + +float **matinv(int m, float **mat); + +float matdet(int m, float **mat); + /** * a struct for get solution of ax = b */