我正在编写自己的矩阵类,其中包含您希望伴随矩阵类的方法.
除了我的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得到了单元测试的数字.
我对分解方法有什么看法？
解决方法
在哪里选择结果列的最大值：
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));
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 [email protected] 举报，一经查实，本站将立刻删除。
