This commit is contained in:
unlockable
2023-12-10 19:48:11 +08:00
parent 9e9e66d90b
commit 172a84a600

View File

@@ -64,8 +64,8 @@ int main() {
}
}
printf("Readin\n");
output();
// printf("Readin\n");
// output();
for (int col = 0; col < n; col++) {
double diag_elem = get_content(col, col);
@@ -74,7 +74,7 @@ int main() {
double base = get_content(row, col);
if (base != 0) {
double coefficient = base / diag_elem;
set_content(row, col, coefficient);
set_content(row, col, -coefficient);
for (int l = col + 1; l < n; l++) {
set_content(row, l, get_content(row, l) - coefficient * get_content(col, l));
}
@@ -82,18 +82,19 @@ int main() {
}
}
printf("L ready\n");
output();
// printf("L ready\n");
// output();
for (int row = 0; row < n; row++) {
set_content(row, row, 1 / get_content(row, row));
double diag_elem = get_content(row, row);
set_content(row, row, 1 / diag_elem);
for (int col = row + 1; col < n; col++) {
set_content(row, col, -get_content(row, col));
set_content(row, col, -get_content(row, col) / diag_elem);
}
}
printf("All ready\n");
output();
// printf("All ready\n");
// output();
for (int i = 0; i < m; i++) {
// read in one col of z
@@ -101,42 +102,46 @@ int main() {
scanf("%lf", &Z[j]);
}
printf("Readin\n");
outputZ();
// printf("Readin\n");
// outputZ();
// calculate Z = L-1 * Z
// do this column by column
for (int j = 0; j < n; j++) {
// diag of U is always 1, z[j] * U[j][j] = z[j], so k doesn't have to equal to j
for (int k = 0; k < j; k++) {
Z[j] += Z[k] * get_content(j, k);
for (int k = j + 1; k < n; k++) {
Z[k] += Z[j] * get_content(k, j);
}
}
printf("L-1 * Z\n");
outputZ();
// printf("L-1 * Z\n");
// outputZ();
// calculate Z = D-1 * Z
for (int j = 0; j < n; j++) {
Z[j] = Z[j] * get_content(j, j);
}
printf("D-1 L-1 * Z\n");
outputZ();
// printf("D-1 L-1 * Z\n");
// outputZ();
// calculate Z = U-1 * Z
for (int j = 0; j < n; j++) {
// do this column by column
// printf("Starting U-1\n");
for (int j = n - 1; j >= 0; j--) {
// diag of U is always 1, z[j] * U[j][j] = z[j], so k doesn't have to equal to j
for (int k = j + 1; k < n; k++) {
Z[j] += Z[k] * get_content(j, k);
for (int k = 0; k < j; k++) {
Z[k] += Z[j] * get_content(k, j);
}
// outputZ();
}
printf("U-1 D-1 L-1 * Z\n");
outputZ();
// printf("U-1 D-1 L-1 * Z\n");
// outputZ();
// output
for (int j = 0; j < n; j++) {
printf("%lf ", Z[j]);
printf("%.4lf ", Z[j]);
}
printf("\n");
}