代码的内容取自here
除了我的Decompose方法之外,一切看起来都很完美(有关LU分解矩阵意义的更多信息,请查看here)：
public Matrix Decompose(out int[] permutationArray,out int toggle) { // Doolittle LUP decomposition with partial pivoting. // rerturns: result is L (with 1s on diagonal) and U; perm holds row permutations; toggle is +1 or -1 (even or odd) if (!this.IsSquare()) { throw new Exception("LU-Decomposition can only be found for a square matrix"); } Matrix decomposedMatrix = new Matrix(this); // copy this matrix before messing with it permutationArray = new int[NumberOfRows]; // set up row permutation result for (int i = 0; i < NumberOfRows; ++i) { permutationArray[i] = i; } toggle = 1; // toggle tracks row swaps. +1 -> even,-1 -> odd. used by MatrixDeterminant for (int columnIndex = 0; columnIndex < NumberOfRows - 1; columnIndex++) // each column { // find largest value in col j double maxOfColumn = Math.Abs(decomposedMatrix.GetElement(columnIndex,columnIndex)); int pivotRowIndex = columnIndex; for (int rowIndex = columnIndex + 1; rowIndex < NumberOfRows; rowIndex++) { if (Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)) > maxOfColumn) { maxOfColumn = decomposedMatrix.GetElement(rowIndex,columnIndex); pivotRowIndex = rowIndex; } } if (pivotRowIndex != columnIndex) // if largest value not on pivot,swap rows { double[] rowPtr = decomposedMatrix.GetRow(pivotRowIndex); decomposedMatrix.SetRowOfMatrix(pivotRowIndex,decomposedMatrix.GetRow(columnIndex)); decomposedMatrix.SetRowOfMatrix(columnIndex,rowPtr); int tmp = permutationArray[pivotRowIndex]; // and swap permutation info permutationArray[pivotRowIndex] = permutationArray[columnIndex]; permutationArray[columnIndex] = tmp; toggle = -toggle; // adjust the row-swap toggle } if (decomposedMatrix.GetElement(columnIndex,columnIndex) == 0.0) { // find a good row to swap int goodRow = -1; for (int row = columnIndex + 1; row < decomposedMatrix.NumberOfRows; row++) { if (decomposedMatrix.GetElement(row,columnIndex) != 0.0) { goodRow = row; } } if (goodRow == -1) { throw new Exception("Cannot use Doolittle's method on this matrix"); } // swap rows so 0.0 no longer on diagonal double[] rowPtr = decomposedMatrix.GetRow(goodRow); decomposedMatrix.SetRowOfMatrix(goodRow,rowPtr); int tmp = permutationArray[goodRow]; // and swap perm info permutationArray[goodRow] = permutationArray[columnIndex]; permutationArray[columnIndex] = tmp; toggle = -toggle; // adjust the row-swap toggle } for (int rowIndex = columnIndex + 1; rowIndex < this.NumberOfRows; ++rowIndex) { double valuetoStore = decomposedMatrix.GetElement(rowIndex,columnIndex) / decomposedMatrix.GetElement(columnIndex,columnIndex); decomposedMatrix.SetElement(rowIndex,columnIndex,valuetoStore); for (int nextColumnIndex = columnIndex + 1; nextColumnIndex < NumberOfRows; ++nextColumnIndex) { double valuetoStore2 = decomposedMatrix.GetElement(rowIndex,nextColumnIndex) - decomposedMatrix.GetElement(rowIndex,columnIndex) * decomposedMatrix.GetElement(columnIndex,nextColumnIndex); decomposedMatrix.SetElement(rowIndex,nextColumnIndex,valuetoStore2); } } } // main j column loop return decomposedMatrix; } // MatrixDecompose 这是它失败的单元测试:(请注意！如果开始矩阵使用正数而不是负数,则完全相同的单元测试通过)
[TestMethod()] public voID Matrix_DecomposeTest3_DifferentSigns() { Matrix matrixA = new Matrix(9); //Set up matrix A double[] row1OfMatrixA = { 3,-7,0 }; double[] row2OfMatrixA = { 0,1,8,-4,2,0 }; double[] row3OfMatrixA = { 0,-9,3,0 }; double[] row4OfMatrixA = { -5,4,7,-1,3 }; double[] row5OfMatrixA = { 0,-3,-8,0 }; double[] row6OfMatrixA = { 0,5,-2,4 }; double[] row7OfMatrixA = { 0,7 }; double[] row8OfMatrixA = { 0,0 }; double[] row9OfMatrixA = { 0,6 }; matrixA.SetRowOfMatrix(0,row1OfMatrixA); matrixA.SetRowOfMatrix(1,row2OfMatrixA); matrixA.SetRowOfMatrix(2,row3OfMatrixA); matrixA.SetRowOfMatrix(3,row4OfMatrixA); matrixA.SetRowOfMatrix(4,row5OfMatrixA); matrixA.SetRowOfMatrix(5,row6OfMatrixA); matrixA.SetRowOfMatrix(6,row7OfMatrixA); matrixA.SetRowOfMatrix(7,row8OfMatrixA); matrixA.SetRowOfMatrix(8,row9OfMatrixA); Console.Out.Write(matrixA); //The LUP Decomposition //Correct L Part Matrix correctLPartOfLUPDecomposition = new Matrix(9); double[] row1OfLMatrix = { 1,0 }; double[] row2OfLMatrix = { 0,0 }; double[] row3OfLMatrix = { 0,.143,0 }; double[] row4OfLMatrix = { -.6,0 }; double[] row5OfLMatrix = { 0,0 }; double[] row6OfLMatrix = { 0,.435,0 }; double[] row7OfLMatrix = { 0,-.217,-.5,0 }; double[] row8OfLMatrix = { 0,-.429,-.796,-.342,-.319,.282,-.226,0 }; double[] row9OfLMatrix = { 0.000,0.286,0.963,0.161,0.523,0.065,0.013,0.767,1.000 }; correctLPartOfLUPDecomposition.SetRowOfMatrix(0,row1OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(1,row2OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(2,row3OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(3,row4OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(4,row5OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(5,row6OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(6,row7OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(7,row8OfLMatrix); correctLPartOfLUPDecomposition.SetRowOfMatrix(8,row9OfLMatrix); //Correct U Part Matrix correctUPartOfLUPDecomposition = new Matrix(9); double[] row1OfUMatrix = { -5.000,0.000,4.000,7.000,-1.000,3.000 }; double[] row2OfUMatrix = { 0.000,2.000,5.000,-3.000,-2.000,4.000 }; double[] row3OfUMatrix = { 0.000,7.714,-0.714,-3.571,1.000,-0.571 }; double[] row4OfUMatrix = { 0.000,-4.600,4.200,-0.600,1.800 }; double[] row5OfUMatrix = { 0.000,-9.000,3.000,0.000 }; double[] row6OfUMatrix = { 0.000,-9.826,4.261,6.217 }; double[] row7OfUMatrix = { 0.000,9.000,9.500 }; double[] row8OfUMatrix = { 0.000,-2.043,2.272 }; double[] row9OfUMatrix = { 0.000,-3.153 }; correctUPartOfLUPDecomposition.SetRowOfMatrix(0,row1OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(1,row2OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(2,row3OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(3,row4OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(4,row5OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(5,row6OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(6,row7OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(7,row8OfUMatrix); correctUPartOfLUPDecomposition.SetRowOfMatrix(8,row9OfUMatrix); //Calculate values for the above int[] permutationArray; int toggleValue; Matrix decomposeResult = matrixA.Decompose(out permutationArray,out toggleValue); Matrix calculatedLPartOfLUPDecomposition = decomposeResult.ExtractLower(); Matrix calculatedUPartOfLUPDecomposition = decomposeResult.ExtractUpper(); //Compare the two matrices double tolerance = .001; Matrix LDifferenceMatrix = calculatedLPartOfLUPDecomposition - correctLPartOfLUPDecomposition; Matrix UDifferenceMatrix = calculatedUPartOfLUPDecomposition - correctUPartOfLUPDecomposition; for (int i = 0; i < LDifferenceMatrix.NumberOfRows; i++) { for (int j = 0; j < LDifferenceMatrix.NumberOfColumns; j++) { (Math.Abs(LDifferenceMatrix.GetElement(i,j)) < tolerance).Should().BeTrue(); } } for (int i = 0; i < UDifferenceMatrix.NumberOfRows; i++) { for (int j = 0; j < UDifferenceMatrix.NumberOfColumns; j++) { (Math.Abs(UDifferenceMatrix.GetElement(i,j)) < tolerance).Should().BeTrue(); } } } 我使用找到的工具here得到了单元测试的数字.
我有单元测试Matrix的所有其他方法,他们都通过.
我对分解方法有什么看法？
解决方法 在哪里选择结果列的最大值：maxOfColumn = decomposedMatrix.GetElement(rowIndex,columnIndex);
你应该真的选择那里找到的最大绝对值(你在if语句中做到这一点,之后似乎已经忘记了).我会通过添加注释来明确指出,以便下一个阅读此代码的人可以更清楚地看到它
//if a value with a larger absolute value has been found replace the currently selectedmax if (Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)) > maxOfColumn) { maxOfColumn = Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)); 总结 以上是内存溢出为你收集整理的c# – 在负面情况下几乎正确的矩阵分解代码失败全部内容，希望文章能够帮你解决c# – 在负面情况下几乎正确的矩阵分解代码失败所遇到的程序开发问题。
如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）