| 409 | */ |
| 410 | template <typename MatrixType, typename OrderingType> |
| 411 | void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) |
| 412 | { |
| 413 | |
| 414 | //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. |
| 415 | |
| 416 | // Firstly, copy the whole input matrix. |
| 417 | m_mat = mat; |
| 418 | |
| 419 | // Compute fill-in ordering |
| 420 | OrderingType ord; |
| 421 | ord(m_mat,m_perm_c); |
| 422 | |
| 423 | // Apply the permutation to the column of the input matrix |
| 424 | if (m_perm_c.size()) |
| 425 | { |
| 426 | 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. |
| 427 | // Then, permute only the column pointers |
| 428 | ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?const_cast<StorageIndex*>(mat.outerIndexPtr()):0); |
| 429 | |
| 430 | // If the input matrix 'mat' is uncompressed, then the outer-indices do not match the ones of m_mat, and a copy is thus needed. |
| 431 | if(!mat.isCompressed()) |
| 432 | IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1); |
| 433 | |
| 434 | // Apply the permutation and compute the nnz per column. |
| 435 | for (Index i = 0; i < mat.cols(); i++) |
| 436 | { |
| 437 | m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; |
| 438 | m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; |
| 439 | } |
| 440 | } |
| 441 | |
| 442 | // Compute the column elimination tree of the permuted matrix |
| 443 | IndexVector firstRowElt; |
| 444 | internal::coletree(m_mat, m_etree,firstRowElt); |
| 445 | |
| 446 | // In symmetric mode, do not do postorder here |
| 447 | if (!m_symmetricmode) { |
| 448 | IndexVector post, iwork; |
| 449 | // Post order etree |
| 450 | internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post); |
| 451 | |
| 452 | |
| 453 | // Renumber etree in postorder |
| 454 | Index m = m_mat.cols(); |
| 455 | iwork.resize(m+1); |
| 456 | for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i)); |
| 457 | m_etree = iwork; |
| 458 | |
| 459 | // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree |
| 460 | PermutationType post_perm(m); |
| 461 | for (Index i = 0; i < m; i++) |
| 462 | post_perm.indices()(i) = post(i); |
| 463 | |
| 464 | // Combine the two permutations : postorder the permutation for future use |
| 465 | if(m_perm_c.size()) { |
| 466 | m_perm_c = post_perm * m_perm_c; |
| 467 | } |
| 468 | |