43 size_t i = 0, j = 0, k = 0;
44 size_t rowP = 0, colP = Cols - 1;
46 for (i = 0; i < Rows; i++)
49 for (i = 0; i < Cols; i++)
58 hProcess = cb->
addItem(
"LU decomposition...",
67 for (i = 0; i < Dim; i++)
69 Completion += Dim - i;
79 for (j = i, pTmp1 = &A(i, i); j < Rows; j++, pTmp1 += Cols)
80 if (fabs(*pTmp1) > tmp)
89 if (tmp < std::numeric_limits< C_FLOAT64 >::epsilon())
99 for (j = i, pTmp1 = &A(i, colP); j < Rows; j++, pTmp1 += Cols)
100 if (fabs(*pTmp1) > tmp)
106 if (tmp >= std::numeric_limits< C_FLOAT64 >::epsilon())
break;
124 for (k = 0; k < Rows; k++, pTmp1 += Cols, pTmp2 += Cols)
131 std::swap(col[i], col[colP]);
139 for (k = 0; k < Cols; k++, pTmp1++, pTmp2++)
146 std::swap(row[i], row[rowP]);
155 pTmp1 = &A(i + 1, i);
157 for (k = i + 1; k < Rows; k++, pTmp1 += Cols)
171 pTmp1 = &A(i + 1, i);
172 pTmp3 = &A(i + 1, i + 1);
174 for (ii = i + 1; ii < Rows; ii++, pTmp1 += Cols, pTmp3 += i + 1)
176 pTmp2 = &A(i, i + 1);
178 for (jj = i + 1; jj < Cols; jj++, pTmp2++, pTmp3++)
180 *pTmp3 -= *pTmp1 * *pTmp2;
187 DebugFile << rowP <<
"\tx\t" << colP + 1 << std::endl;
188 DebugFile <<
"Row\t" << row << std::endl;
189 DebugFile <<
"Col\t" << col << std::endl;
190 DebugFile << A << std::endl;
214 if (lRow.size() != Rows) lRow.resize(Rows);
218 size_t i = 0, j = 0, k = 0;
219 size_t ip, jp = Cols - 1;
221 for (i = 0; i < Rows; i++)
227 for (i = 0; i < Cols; i++)
233 if (Rows == 0 || Cols == 0)
238 for (i = 0; i < Cols; i++)
239 for (j = 0; j < Rows; j++)
253 C_INT subRows = Rows;
254 C_INT subCols = Cols;
258 size_t LastCols = Cols;
312 dgetrf_(&subRows, &subCols, &T(lCol[iTop], iTop), &Rows, RowPivot.array(), &Info);
314 if (Info < 0)
return false;
317 DebugFile <<
"After dgetrf_\n";
318 DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
321 for (i = 0; i < RowPivot.size(); i++)
322 if (i != RowPivot[i] - 1)
324 std::swap(row[i + iTop], row[RowPivot[i] - 1 + iTop]);
326 pTmp1 = &T(0, i + iTop);
327 pTmp2 = &T(0, RowPivot[i] - 1 + iTop);
329 for (k = 0; k < iTop; k++, pTmp1 += Rows, pTmp2 += Rows)
338 DebugFile << row << std::endl;
339 DebugFile <<
"After row swap\n";
340 DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
343 if (Info == 0)
break;
348 for (i = iTop + 1, found =
false; i < imax; i++)
349 if (T(lCol[i], i) != 0.0)
358 for (i = Rows - 1; jTop <= i; i--)
360 for (j = Cols - 1; lCol[jTop] <= j; j--)
367 for (k = iTop; k < i && k < j; k++, pTmp1 += Rows, pTmp2++)
368 tmp += *pTmp1 * *pTmp2;
370 if (i <= j) tmp += *pTmp2;
376 subRows = Rows - iTop;
377 subCols = Cols - lCol[iTop];
379 pTmp1 = &T(lCol[iTop], 0);
380 pTmp2 = &T(lCol[jTop], 0);
382 for (k = 0; k < Rows; k++, pTmp1++, pTmp2++)
389 std::swap(col[iTop], col[jTop]);
392 DebugFile <<
"After partial inversion: " << iTop <<
" , " << jTop << std::endl;
393 DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
394 DebugFile <<
"After column swap\n";
395 DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
400 for (j = iTop + 1, found =
false; j < Cols; j++)
402 for (i = iTop; i < j && i < Rows; i++)
403 if (T(lCol[j], i) != 0)
416 subRows = Rows - iTop;
417 subCols = Cols - lCol[jTop];
422 for (j = jTop; j < LastCols; j += Delta)
423 for (i = j - Delta; i < j && i + Delta < LastCols; i++)
425 std::swap(col[i], col[i + Delta]);
426 std::swap(lCol[i], lCol[i + Delta]);
432 RowPivot.resize(
std::min(subRows, subCols));
435 DebugFile << lCol << std::endl;
436 DebugFile << col << std::endl;
440 for (i = 0; i < Cols; i++)
441 for (j = 0; j < Rows; j++)
442 A(j, i) = T(lCol[i], j);
int dgetrf_(integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info)
virtual size_t numRows() const
void resize(size_t size, const bool ©=false)
virtual bool progressItem(const size_t &handle)
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
bool LUfactor(CMatrix< C_FLOAT64 > &A, CVector< size_t > &row, CVector< size_t > &col, CProcessReport *cb)
virtual bool finishItem(const size_t &handle)
virtual size_t numCols() const