99 using Real =
typename MatrixType::ElemType;
101 if (A.getNRow() != A.getNCol() || A.getNRow() != AI.getNRow() ||
102 A.getNCol() != AI.getNCol()) {
106 int size = A.getNRow();
107 int *index = NULL, iScratch[10];
108 Real *column = NULL, dScratch[10];
116 index =
new int[size];
117 column =
new Real[size];
143 typename MatrixType::ElemPoinerType tmp2Size) {
144 if (A.getNRow() != A.getNCol() || A.getNRow() != AI.getNRow() ||
145 A.getNCol() != AI.getNCol() || A.getNRow() != size) {
160 for (j = 0; j < size; j++) {
161 for (i = 0; i < size; i++) {
168 for (i = 0; i < size; i++) {
169 AI(i, j) = tmp2Size[i];
188 typename MatrixType::ElemPoinerType tmpSize) {
189 using Real =
typename MatrixType::ElemType;
193 Real largest, temp1, temp2, sum;
198 for (i = 0; i < size; i++) {
199 for (largest = 0.0, j = 0; j < size; j++) {
200 if ((temp2 = abs(A(i, j))) > largest) { largest = temp2; }
203 if (largest == 0.0) {
204 std::cerr <<
"Unable to factor linear system";
207 tmpSize[i] = 1.0 / largest;
212 for (j = 0; j < size; j++) {
213 for (i = 0; i < j; i++) {
215 for (k = 0; k < i; k++) {
216 sum -= A(i, k) * A(k, j);
223 for (largest = 0.0, i = j; i < size; i++) {
225 for (k = 0; k < j; k++) {
226 sum -= A(i, k) * A(k, j);
230 if ((temp1 = tmpSize[i] * abs(sum)) >= largest) {
239 for (k = 0; k < size; k++) {
241 A(maxI, k) = A(j, k);
244 tmpSize[maxI] = tmpSize[j];
251 if (abs(A(j, j)) <= std::numeric_limits<Real>::epsilon()) {
252 std::cerr <<
"Unable to factor linear system";
256 if (j != (size - 1)) {
257 temp1 = 1.0 / A(j, j);
258 for (i = j + 1; i < size; i++) {
280 typename MatrixType::ElemPoinerType x,
int size) {
281 using Real =
typename MatrixType::ElemType;
289 for (ii = -1, i = 0; i < size; i++) {
295 for (j = ii; j <= (i - 1); j++) {
296 sum -= A(i, j) * x[j];
307 for (i = size - 1; i >= 0; i--) {
309 for (j = i + 1; j < size; j++) {
310 sum -= A(i, j) * x[j];
312 x[i] = sum / A(i, i);
void LUSolveLinearSystem(MatrixType &A, int *index, typename MatrixType::ElemPoinerType x, int size)
Solve linear equations Ax = b using LU decompostion A = LU where L is lower triangular matrix and U i...
int LUFactorLinearSystem(MatrixType &A, int *index, int size, typename MatrixType::ElemPoinerType tmpSize)
Factor linear equations Ax = b using LU decompostion A = LU where L is lower triangular matrix and U ...