| 367 | */ |
| 368 | template <typename MatrixType, typename OrderingType> |
| 369 | void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) |
| 370 | { |
| 371 | |
| 372 | //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. |
| 373 | |
| 374 | OrderingType ord; |
| 375 | ord(mat,m_perm_c); |
| 376 | |
| 377 | // Apply the permutation to the column of the input matrix |
| 378 | //First copy the whole input matrix. |
| 379 | m_mat = mat; |
| 380 | if (m_perm_c.size()) { |
| 381 | m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. |
| 382 | //Then, permute only the column pointers |
| 383 | const Index * outerIndexPtr; |
| 384 | if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr(); |
| 385 | else |
| 386 | { |
| 387 | Index *outerIndexPtr_t = new Index[mat.cols()+1]; |
| 388 | for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; |
| 389 | outerIndexPtr = outerIndexPtr_t; |
| 390 | } |
| 391 | for (Index i = 0; i < mat.cols(); i++) |
| 392 | { |
| 393 | m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; |
| 394 | m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; |
| 395 | } |
| 396 | if(!mat.isCompressed()) delete[] outerIndexPtr; |
| 397 | } |
| 398 | // Compute the column elimination tree of the permuted matrix |
| 399 | IndexVector firstRowElt; |
| 400 | internal::coletree(m_mat, m_etree,firstRowElt); |
| 401 | |
| 402 | // In symmetric mode, do not do postorder here |
| 403 | if (!m_symmetricmode) { |
| 404 | IndexVector post, iwork; |
| 405 | // Post order etree |
| 406 | internal::treePostorder(m_mat.cols(), m_etree, post); |
| 407 | |
| 408 | |
| 409 | // Renumber etree in postorder |
| 410 | Index m = m_mat.cols(); |
| 411 | iwork.resize(m+1); |
| 412 | for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i)); |
| 413 | m_etree = iwork; |
| 414 | |
| 415 | // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree |
| 416 | PermutationType post_perm(m); |
| 417 | for (Index i = 0; i < m; i++) |
| 418 | post_perm.indices()(i) = post(i); |
| 419 | |
| 420 | // Combine the two permutations : postorder the permutation for future use |
| 421 | if(m_perm_c.size()) { |
| 422 | m_perm_c = post_perm * m_perm_c; |
| 423 | } |
| 424 | |
| 425 | } // end postordering |
| 426 |
nothing calls this directly
no test coverage detected