Line data Source code
1 : /*
2 : * File: LDLFactorization_AAt.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on November 5, 2015, 1:12 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 "S_LDLFactorization.h"
22 :
23 2 : Matrix S_LDLFactorization::multiply_AAtr_betaI(Matrix& A, double beta) {
24 2 : size_t n = A.getNrows();
25 2 : size_t m = A.getNcols();
26 2 : Matrix result(n, n, Matrix::MATRIX_SYMMETRIC);
27 :
28 7 : for (size_t i = 0; i < n; i++) {
29 14 : for (size_t j = 0; j <= i; j++) {
30 9 : result.set(i, j, 0.0);
31 81 : for (size_t k = 0; k < m; k++) {
32 72 : result.set(i, j, result.get(i, j) + A.get(i, k) * A.get(j, k));
33 : }
34 9 : if (i == j) {
35 5 : result.set(i, j, result.get(i, j) + beta);
36 : }
37 : }
38 : }
39 2 : return result;
40 : }
41 :
42 3 : S_LDLFactorization::S_LDLFactorization(Matrix& matrix, double beta) : FactoredSolver(matrix), m_beta(beta) {
43 3 : m_factor = NULL;
44 3 : m_delegated_solver = NULL;
45 3 : }
46 :
47 3 : int S_LDLFactorization::factorize() {
48 3 : if (m_matrix_type == Matrix::MATRIX_SPARSE) {
49 : double beta_temp[2];
50 1 : beta_temp[0] = m_beta;
51 1 : beta_temp[1] = 0.0;
52 :
53 1 : if (m_matrix->m_sparse == NULL) {
54 1 : m_matrix->_createSparse();
55 : }
56 :
57 1 : m_matrix->m_sparse->stype = 0;
58 1 : m_factor = cholmod_analyze(m_matrix->m_sparse, Matrix::cholmod_handle());
59 1 : cholmod_factorize_p(m_matrix->m_sparse, beta_temp, NULL, 0, m_factor, Matrix::cholmod_handle());
60 1 : return (m_factor->minor == m_matrix->m_nrows) ? ForBESUtils::STATUS_OK : ForBESUtils::STATUS_NUMERICAL_PROBLEMS;
61 2 : } else if (m_matrix_type == Matrix::MATRIX_DENSE) {
62 : /*
63 : * We here need to factorize a dense matrix
64 : * We distinguish between two cases
65 : * 1. A is short (more columns than rows)
66 : * 2. A is tall (more rows than columns)
67 : */
68 2 : if (m_matrix->getNrows() <= m_matrix->getNcols()) {
69 : /* this is a ###SHORT### matrix */
70 : /*
71 : * Note: here we're creating matrix F on the fly and we're passing
72 : * its REFERENCE to LDLFactorization. Eventually, m_delegated_solver will
73 : * lose track of the original matrix.
74 : */
75 1 : Matrix F = multiply_AAtr_betaI(*m_matrix, m_beta);
76 1 : m_delegated_solver = new LDLFactorization(F);
77 1 : int status = m_delegated_solver->factorize();
78 1 : return status;
79 : } else {
80 : /* this is a ~~~TALL~~~ matrix */
81 1 : m_matrix->transpose();
82 : /*
83 : * Note: here we're creating matrix F_tilde on the fly and we're passing
84 : * its REFERENCE to LDLFactorization. Eventually, m_delegated_solver will
85 : * lose track of the original matrix.
86 : */
87 1 : Matrix F_tilde = multiply_AAtr_betaI(*m_matrix, m_beta);
88 1 : m_matrix->transpose();
89 1 : m_delegated_solver = new LDLFactorization(F_tilde);
90 1 : int status = m_delegated_solver->factorize();
91 1 : return status;
92 : }
93 : } else {
94 0 : throw std::invalid_argument("[uoe] Unsupported operation");
95 : }
96 : }
97 :
98 3 : int S_LDLFactorization::solve(Matrix& rhs, Matrix& solution) {
99 3 : solution = Matrix(rhs.m_nrows, rhs.m_ncols);
100 3 : if (m_matrix_type == Matrix::MATRIX_SPARSE) {
101 1 : if (m_factor == NULL) {
102 0 : throw std::invalid_argument(__FCT_MISS_EXCPT);
103 : }
104 : cholmod_dense *b;
105 : cholmod_dense *x;
106 1 : b = cholmod_allocate_dense(rhs.m_nrows, rhs.m_ncols, rhs.m_nrows, CHOLMOD_REAL, Matrix::cholmod_handle());
107 1 : b->x = rhs.m_data;
108 1 : x = cholmod_solve(CHOLMOD_A, m_factor, b, Matrix::cholmod_handle());
109 1 : solution.m_delete_data = false;
110 1 : memcpy(solution.m_data, static_cast<double*> (x->x), rhs.m_nrows * rhs.m_ncols * sizeof (double));
111 1 : cholmod_free_dense(&x, Matrix::cholmod_handle());
112 1 : return ForBESUtils::STATUS_OK;
113 2 : } else if (m_matrix_type == Matrix::MATRIX_DENSE) {
114 2 : if (m_delegated_solver == NULL) {
115 0 : throw std::invalid_argument(__FCT_MISS_EXCPT);
116 : }
117 2 : if (m_matrix_nrows <= m_matrix_ncols) {
118 : /* m_matrix is ###SHORT### and dense */
119 1 : int status = m_delegated_solver->solve(rhs, solution);
120 1 : return status;
121 : } else {
122 : /* m_matrix is ~~~TALL~~~ and dense */
123 1 : m_matrix->transpose();
124 1 : Matrix temp = (*m_matrix) * const_cast<Matrix&>(rhs);
125 1 : m_matrix->transpose();
126 2 : Matrix c;
127 1 : int status = m_delegated_solver->solve(temp, c);
128 1 : if (status != ForBESUtils::STATUS_OK){
129 0 : return status;
130 : }
131 1 : solution = *m_matrix * c - rhs;
132 1 : double beta_inv = -1.0/m_beta;
133 1 : solution *= beta_inv;
134 2 : return ForBESUtils::STATUS_OK;
135 : }
136 :
137 : } else {
138 0 : throw std::invalid_argument("[uoe] Unsupported operation");
139 : }
140 : }
141 :
142 9 : S_LDLFactorization::~S_LDLFactorization() {
143 3 : if (m_factor != NULL) {
144 1 : cholmod_free_factor(&m_factor, Matrix::cholmod_handle());
145 1 : m_factor = NULL;
146 : }
147 3 : if (m_delegated_solver != NULL) {
148 2 : delete m_delegated_solver;
149 : }
150 6 : }
151 :
|
