Line data Source code
1 : /*
2 : * File: CholeskyFactorization.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on August 4, 2015, 8:14 PM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include "CholeskyFactorization.h"
22 :
23 8 : CholeskyFactorization::CholeskyFactorization(Matrix& matrix) :
24 8 : FactoredSolver(matrix) {
25 8 : m_L = NULL;
26 8 : m_factor = NULL;
27 8 : if (matrix.getNrows() != matrix.getNcols()){
28 0 : throw std::invalid_argument("CholeskyFactorization factorization can only be applied to square matrices");
29 : }
30 8 : if (matrix.getType() != Matrix::MATRIX_SPARSE) {
31 6 : this->m_L = new double[matrix.length()]();
32 : }
33 8 : }
34 :
35 21 : CholeskyFactorization::~CholeskyFactorization() {
36 7 : if (m_L != NULL) {
37 5 : delete[] m_L;
38 5 : m_L = NULL;
39 : }
40 7 : if (m_factor != NULL) {
41 2 : cholmod_free_factor(&m_factor, Matrix::cholmod_handle());
42 2 : m_factor = NULL;
43 : }
44 14 : }
45 :
46 8 : int CholeskyFactorization::factorize() {
47 8 : if (m_matrix_type == Matrix::MATRIX_SPARSE) {
48 : /* Cholesky decomposition of a SPARSE matrix: */
49 2 : if (m_matrix->m_sparse == NULL) {
50 2 : m_matrix->_createSparse();
51 : }
52 : /* analyze */
53 2 : m_factor = cholmod_analyze(m_matrix->m_sparse, Matrix::cholmod_handle());
54 : /* factorize */
55 2 : cholmod_factorize(m_matrix->m_sparse, m_factor, Matrix::cholmod_handle());
56 : /* Success: status = 0, else 1*/
57 2 : return (m_factor->minor == m_matrix->m_nrows) ? ForBESUtils::STATUS_OK : ForBESUtils::STATUS_NUMERICAL_PROBLEMS;
58 : } else { /* If this is any non-sparse matrix: */
59 6 : memcpy(m_L, m_matrix->getData(), m_matrix->length() * sizeof (double)); /* m_L := m_matrix.m_data */
60 6 : int info = ForBESUtils::STATUS_OK;
61 6 : if (m_matrix_type == Matrix::MATRIX_DENSE) { /* This is a dense matrix */
62 4 : info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, m_L, m_matrix_nrows);
63 : #ifdef SET_L_OFFDIAG_TO_ZERO
64 : for (size_t i = 0; i < m_matrix_nrows; i++) {
65 : for (size_t j = i + 1; j < m_matrix_nrows; j++) {
66 : L.set(i, j, 0.0);
67 : }
68 : }
69 : #endif
70 2 : } else if (m_matrix_type == Matrix::MATRIX_SYMMETRIC) { /* This is a symmetric matrix */
71 2 : info = LAPACKE_dpptrf(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, m_L);
72 : }
73 6 : return info;
74 : }
75 : }
76 :
77 98 : int CholeskyFactorization::solve(Matrix& rhs, Matrix& solution) {
78 98 : if (m_matrix_type == Matrix::MATRIX_SPARSE) {
79 : cholmod_dense *x;
80 :
81 : /* cast the RHS as a cholmod_dense b = rhs */
82 2 : if (rhs.m_type == Matrix::MATRIX_DENSE) {
83 :
84 : cholmod_dense *b;
85 2 : b = cholmod_allocate_dense(rhs.m_nrows, rhs.m_ncols, rhs.m_nrows, CHOLMOD_REAL, Matrix::cholmod_handle());
86 2 : b->x = rhs.m_data;
87 :
88 : /* Solve - rhs is dense*/
89 2 : x = cholmod_solve(CHOLMOD_A, m_factor, b, Matrix::cholmod_handle());
90 2 : solution = Matrix(rhs.m_nrows, rhs.m_ncols);
91 2 : solution.m_delete_data = false;
92 2 : memcpy(solution.m_data, static_cast<double*> (x->x), rhs.m_nrows * rhs.m_ncols * sizeof (double));
93 2 : cholmod_free_dense(&x, Matrix::cholmod_handle());
94 :
95 0 : } else if (rhs.m_type == Matrix::MATRIX_SPARSE) {
96 : // still untested!
97 : cholmod_sparse * rhs_sparse;
98 0 : if (rhs.m_sparse == NULL) {
99 0 : const_cast<Matrix&> (rhs)._createSparse();
100 : }
101 0 : rhs_sparse = rhs.m_sparse;
102 0 : cholmod_sparse * result = cholmod_spsolve(CHOLMOD_LDLt, m_factor, rhs_sparse, Matrix::cholmod_handle());
103 0 : solution = Matrix(rhs.m_nrows, rhs.m_ncols, Matrix::MATRIX_SPARSE);
104 0 : solution.m_sparse = result;
105 0 : solution._createTriplet();
106 : } else {
107 0 : throw std::logic_error("Not supported");
108 : }
109 2 : return ForBESUtils::STATUS_OK;
110 : } else { /* the matrix to be factorized is not sparse */
111 96 : int info = ForBESUtils::STATUS_UNDEFINED_FUNCTION;
112 96 : solution = Matrix(rhs);
113 96 : if (m_matrix_type == Matrix::MATRIX_DENSE) {
114 85 : info = LAPACKE_dpotrs(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, rhs.m_ncols, m_L, m_matrix_nrows, solution.m_data, m_matrix_nrows);
115 11 : } else if (m_matrix_type == Matrix::MATRIX_SYMMETRIC) {
116 11 : info = LAPACKE_dpptrs(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, rhs.m_ncols, m_L, solution.m_data, m_matrix_nrows);
117 : } else {
118 0 : throw std::invalid_argument("This matrix type is not supported - only DENSE, SPARSE and SYMMETRIC are supported");
119 : }
120 96 : return info;
121 : }
122 : }
123 :
124 :
|
