summaryrefslogtreecommitdiff
path: root/doc/software_manual.org
diff options
context:
space:
mode:
Diffstat (limited to 'doc/software_manual.org')
-rw-r--r--doc/software_manual.org64
1 files changed, 33 insertions, 31 deletions
diff --git a/doc/software_manual.org b/doc/software_manual.org
index 08a554e..6d6be25 100644
--- a/doc/software_manual.org
+++ b/doc/software_manual.org
@@ -373,49 +373,51 @@ void format_vector_into(Array_double *v, char *s) {
decomposed
#+BEGIN_SRC c
-Matrix_double **lu_decomp(Matrix_double *m) {
- assert(m->cols == m->rows);
+ Matrix_double **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 = copy_matrix(m);
+ Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+ Matrix_double *l = put_identity_diagonal(l_empt);
+ free(l_empt);
- 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);
+ 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 (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);
- }
+ if (denom == 0) {
+ printf("ERROR: non-factorable matrix\n");
+ exit(-1);
+ }
- double factor = -(u->data[y]->data[x] / denom);
+ 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]);
+ 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->data[y] = added;
+ l->data[y]->data[x] = -factor;
+ }
}
}
- }
- u_l[0] = u;
- u_l[1] = l;
- return u_l;
-}
+ u_l[0] = u;
+ u_l[1] = l;
+ return u_l;
+ }
#+END_SRC
*** ~bsubst~
+ Author: Elizabeth Hunt