计算LU的逆错误。

2023207/input.txt

@@ -0,0 +1,9 @@
+3
+3 2
+44 62 
+44 43 30 
+3 34 
+27 63 53 
+14 52 19
+
+3 3 2 44 62 44 43 30 3 34 27 63 53 14 52 19

2023207/main.cpp

@@ -1,5 +1,5 @@
-#include <stdio.h>
 #include <math.h>
+#include <stdio.h>
 int p = 0;
 int m = 0, n = 0;
 double A[5][10000], Z[10000];
@@ -17,7 +17,7 @@ double get_content(int row, int column) {
     return A[(row - column + p / 2)][column];
 }
 
-int get_content(int row, int column, double val) {
+int set_content(int row, int column, double val) {
     if (abs(row - column) > p / 2) {
         return 0;
     }
@@ -31,30 +31,114 @@ int get_content(int row, int column, double val) {
     return 0;
 }
 
+int output() {
+    for (int i = 0; i < p; i++) {
+        for (int j = 0; j < n; j++) {
+            printf("%lf ", get_content(i, j));
+        }
+        printf("\n");
+    }
+    printf("\n");
+    return 0;
+}
+
+int outputZ() {
+    for (int i = 0; i < n; i++) {
+        printf("%lf ", Z[i]);
+    }
+    printf("\n");
+    return 0;
+}
+
 int main() {
     scanf("%d", &p);
     scanf("%d %d", &n, &m);
-    for (int i = 0; i < p; i++) {
-        for (int j = 0; j < n - abs(i - p / 2); j++) {
+    for (int i = 0; i < p / 2; i++) {
+        for (int j = p / 2 - i; j < n; j++) {
+            scanf("%lf", &A[i][j]);
+        }
+    }
+    for (int i = p / 2; i < p; i++) {
+        for (int j = 0; j < n - (i - p / 2); j++) {
             scanf("%lf", &A[i][j]);
         }
     }
 
+    printf("Readin\n");
+    output();
+
     for (int col = 0; col < n; col++) {
+        double diag_elem = get_content(col, col);
+        // make everything below the diag to 0 row by row
         for (int row = col + 1; row < n; row++) {
-            double diag_elem = get_content(row, row);
             double base = get_content(row, col);
             if (base != 0) {
-
+                double coefficient = base / diag_elem;
+                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));
+                }
             }
         }
     }
 
-    // for (int i = 0; i < p; i++) {
-    //     for (int j = 0; j < n; j++) {
-    //         printf("%lf ", get_content(i, j));
-    //     }
-    //     printf("\n");
-    // }
+    printf("L ready\n");
+    output();
+
+    for (int row = 0; row < n; row++) {
+        set_content(row, row, 1 / get_content(row, row));
+        for (int col = row + 1; col < n; col++) {
+            set_content(row, col, -get_content(row, col));
+        }
+    }
+
+    printf("All ready\n");
+    output();
+
+    for (int i = 0; i < m; i++) {
+        // read in one col of z
+        for (int j = 0; j < n; j++) {
+            scanf("%lf", &Z[j]);
+        }
+
+        printf("Readin\n");
+        outputZ();
+
+        // calculate Z = L-1 * Z
+        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);
+            }
+        }
+
+        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();
+
+        // calculate Z =  U-1 * Z
+        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 = j + 1; k < n; k++) {
+                Z[j] += Z[k] * get_content(j, k);
+            }
+        }
+
+        printf("U-1 D-1 L-1 * Z\n");
+        outputZ();
+
+        // output
+        for (int j = 0; j < n; j++) {
+            printf("%lf ", Z[j]);
+        }
+        printf("\n");
+    }
     return 0;
 }