| 526 | */ |
| 527 | template <typename MatrixType, typename OrderingType> |
| 528 | void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) { |
| 529 | // TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. |
| 530 | |
| 531 | // Firstly, copy the whole input matrix. |
| 532 | m_mat = mat; |
| 533 | |
| 534 | // Compute fill-in ordering |
| 535 | OrderingType ord; |
| 536 | ord(m_mat, m_perm_c); |
| 537 | |
| 538 | // Apply the permutation to the column of the input matrix |
| 539 | if (m_perm_c.size()) { |
| 540 | m_mat.uncompress(); // NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This |
| 541 | // vector is filled but not subsequently used. |
| 542 | // Then, permute only the column pointers |
| 543 | ei_declare_aligned_stack_constructed_variable( |
| 544 | StorageIndex, outerIndexPtr, mat.cols() + 1, |
| 545 | mat.isCompressed() ? const_cast<StorageIndex*>(mat.outerIndexPtr()) : 0); |
| 546 | |
| 547 | // If the input matrix 'mat' is uncompressed, then the outer-indices do not match the ones of m_mat, and a copy is |
| 548 | // thus needed. |
| 549 | if (!mat.isCompressed()) |
| 550 | IndexVector::Map(outerIndexPtr, mat.cols() + 1) = IndexVector::Map(m_mat.outerIndexPtr(), mat.cols() + 1); |
| 551 | |
| 552 | // Apply the permutation and compute the nnz per column. |
| 553 | for (Index i = 0; i < mat.cols(); i++) { |
| 554 | m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; |
| 555 | m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i + 1] - outerIndexPtr[i]; |
| 556 | } |
| 557 | } |
| 558 | |
| 559 | // Compute the column elimination tree of the permuted matrix |
| 560 | IndexVector firstRowElt; |
| 561 | internal::coletree(m_mat, m_etree, firstRowElt); |
| 562 | |
| 563 | // In symmetric mode, do not do postorder here |
| 564 | if (!m_symmetricmode) { |
| 565 | IndexVector post, iwork; |
| 566 | // Post order etree |
| 567 | internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post); |
| 568 | |
| 569 | // Renumber etree in postorder |
| 570 | Index m = m_mat.cols(); |
| 571 | iwork.resize(m + 1); |
| 572 | for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i)); |
| 573 | m_etree = iwork; |
| 574 | |
| 575 | // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree |
| 576 | PermutationType post_perm(m); |
| 577 | for (Index i = 0; i < m; i++) post_perm.indices()(i) = post(i); |
| 578 | |
| 579 | // Combine the two permutations : postorder the permutation for future use |
| 580 | if (m_perm_c.size()) { |
| 581 | m_perm_c = post_perm * m_perm_c; |
| 582 | } |
| 583 | |
| 584 | } // end postordering |
| 585 |
nothing calls this directly
no test coverage detected