Line data Source code
1 : /*
2 : * File: LDLFactorization.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on July 30, 2015, 3:02 AM
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 "LDLFactorization.h"
22 :
23 43 : LDLFactorization::LDLFactorization(Matrix& matr) : FactoredSolver(matr) {
24 43 : this->LDL = NULL;
25 43 : this->ipiv = NULL;
26 43 : this->m_sparse_ldl_factor = NULL;
27 43 : this->m_matrix_type = m_matrix->getType();
28 43 : this->m_matrix_nrows = m_matrix->getNrows();
29 43 : if (matr.isEmpty()){
30 0 : throw std::invalid_argument("LDL factorization cannot be applied to empty matrices");
31 : }
32 43 : if (m_matrix_type == Matrix::MATRIX_LOWERTR) {
33 0 : throw std::invalid_argument("LDL factorization cannot be applied to non-symmetric matrices (e.g., Lower/Upper triangular)");
34 : }
35 105 : if (m_matrix_type != Matrix::MATRIX_DENSE &&
36 45 : matr.getType() != Matrix::MATRIX_SYMMETRIC &&
37 2 : matr.getType() != Matrix::MATRIX_SPARSE) { // Only DENSE and SYMMETRIC and SPARSE are supported!
38 0 : throw std::logic_error("This matrix type is not supported by LDLFactorization");
39 : }
40 43 : if (matr.getNrows() != matr.getNcols()) {
41 0 : throw std::invalid_argument("Matrix not square");
42 : }
43 43 : if (m_matrix_type == Matrix::MATRIX_SPARSE) {
44 2 : m_sparse_ldl_factor = new sparse_ldl_factor;
45 45 : return;
46 : }
47 41 : this->LDL = new double[matr.length()];
48 41 : this->ipiv = new int[matr.getNrows()];
49 41 : memcpy(this->LDL, matr.getData(), matr.length() * sizeof (double));
50 : }
51 :
52 129 : LDLFactorization::~LDLFactorization() {
53 43 : if (this->LDL != NULL) {
54 41 : delete[] LDL;
55 : }
56 43 : if (this->ipiv != NULL) {
57 41 : delete[] this->ipiv;
58 : }
59 86 : }
60 :
61 43 : int LDLFactorization::factorize() {
62 43 : int status = ForBESUtils::STATUS_UNDEFINED_FUNCTION;
63 43 : if (this->m_matrix_type == Matrix::MATRIX_DENSE) {
64 24 : status = LAPACKE_dsytrf(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, LDL, m_matrix_nrows, ipiv);
65 19 : } else if (this->m_matrix_type == Matrix::MATRIX_SYMMETRIC) {
66 17 : status = LAPACKE_dsptrf(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, LDL, ipiv);
67 2 : } else if (this->m_matrix_type == Matrix::MATRIX_SPARSE) {
68 : // Factorize sparse matrix
69 2 : m_matrix->_createSparse();
70 2 : int * Parent = new int[m_matrix_nrows];
71 2 : int * Lnz = new int[m_matrix_nrows];
72 2 : int * Flag = new int[m_matrix_nrows];
73 2 : m_sparse_ldl_factor->Lp = new int[m_matrix_nrows + 1];
74 : ldl_symbolic(m_matrix_nrows,
75 : static_cast<int*>(m_matrix->m_sparse->p),
76 : static_cast<int*>(m_matrix->m_sparse->i),
77 : m_sparse_ldl_factor->Lp,
78 : Parent,
79 : Lnz,
80 : Flag,
81 2 : NULL, NULL);
82 2 : int lnz = m_sparse_ldl_factor->Lp[m_matrix_nrows];
83 : int d;
84 2 : m_sparse_ldl_factor->Li = new int[lnz];
85 2 : m_sparse_ldl_factor->Lx = new double[lnz];
86 2 : m_sparse_ldl_factor->D = new double[m_matrix_nrows];
87 :
88 2 : double *Y = new double[m_matrix_nrows];
89 2 : int *Pattern = new int[m_matrix_nrows];
90 :
91 : d = ldl_numeric(m_matrix_nrows,
92 : static_cast<int*>(m_matrix->m_sparse->p),
93 : static_cast<int*>(m_matrix->m_sparse->i),
94 : static_cast<double*>(m_matrix->m_sparse->x),
95 : m_sparse_ldl_factor->Lp,
96 : Parent,
97 : Lnz,
98 : m_sparse_ldl_factor->Li,
99 : m_sparse_ldl_factor->Lx,
100 : m_sparse_ldl_factor->D,
101 : Y,
102 : Pattern,
103 : Flag,
104 2 : NULL, NULL);
105 :
106 2 : delete[] Parent;
107 2 : delete[] Pattern;
108 2 : delete[] Y;
109 2 : delete[] Flag;
110 2 : delete[] Lnz;
111 :
112 2 : return d == m_matrix_nrows ? ForBESUtils::STATUS_OK : ForBESUtils::STATUS_NUMERICAL_PROBLEMS;
113 :
114 : } else {
115 0 : throw std::invalid_argument("This matrix type is not supported by LDLFactorization");
116 : }
117 41 : return status;
118 : }
119 :
120 43 : int LDLFactorization::solve(Matrix& rhs, Matrix& solution) {
121 43 : solution = Matrix(m_matrix_nrows, 1, Matrix::MATRIX_DENSE); // solution = rhs (DENSE)
122 1124 : for (size_t i = 0; i < m_matrix_nrows; i++) {
123 1081 : solution.set(i, 0, rhs.get(i, 0));
124 : }
125 43 : int status = ForBESUtils::STATUS_OK;
126 43 : if (Matrix::MATRIX_DENSE == this->m_matrix_type) {
127 23 : status = LAPACKE_dsytrs(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, 1, LDL, m_matrix_nrows, ipiv, solution.getData(), m_matrix_nrows);
128 20 : } else if (Matrix::MATRIX_SYMMETRIC == this->m_matrix_type) {
129 17 : status = LAPACKE_dsptrs(LAPACK_COL_MAJOR, 'L', m_matrix_nrows, 1, LDL, ipiv, solution.getData(), m_matrix_nrows);
130 3 : } else if (Matrix::MATRIX_SPARSE == this->m_matrix_type) {
131 3 : double * b = solution.getData();
132 : ldl_lsolve(m_matrix_nrows, b,
133 : m_sparse_ldl_factor->Lp,
134 : m_sparse_ldl_factor->Li,
135 3 : m_sparse_ldl_factor->Lx);
136 3 : ldl_dsolve(m_matrix_nrows, b, m_sparse_ldl_factor->D);
137 3 : ldl_ltsolve(m_matrix_nrows, b, m_sparse_ldl_factor->Lp, m_sparse_ldl_factor->Li, m_sparse_ldl_factor->Lx);
138 3 : status = ForBESUtils::STATUS_OK;
139 : }
140 43 : return status;
141 : }
142 :
143 0 : double* LDLFactorization::getLDL() const {
144 0 : return LDL;
145 : }
146 :
147 0 : int* LDLFactorization::getIpiv() const {
148 0 : return ipiv;
149 : }
|
