1 /* linalg/ldlt.c
2 *
3 * Copyright (C) 2018 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* L D L^T decomposition of a symmetric positive semi-definite matrix */
21
22 #include <config.h>
23
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_errno.h>
26 #include <gsl/gsl_vector.h>
27 #include <gsl/gsl_matrix.h>
28 #include <gsl/gsl_blas.h>
29 #include <gsl/gsl_linalg.h>
30
31 static double ldlt_norm1(const gsl_matrix * A);
32 static int ldlt_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
33
34 /*
35 gsl_linalg_ldlt_decomp()
36 Perform L D L^T decomposition of a symmetric positive
37 semi-definite matrix using lower triangle
38
39 Inputs: A - (input) symmetric, positive semi-definite matrix
40 (output) lower triangle contains L factor;
41 diagonal contains D
42
43 Return: success/error
44
45 Notes:
46 1) Based on algorithm 4.1.1 of Golub and Van Loan, Matrix Computations (4th ed).
47 2) The first subrow A(1, 2:end) is used as temporary workspace
48 3) The 1-norm ||A||_1 of the original matrix is stored in the upper right corner on output
49 */
50
51 int
gsl_linalg_ldlt_decomp(gsl_matrix * A)52 gsl_linalg_ldlt_decomp (gsl_matrix * A)
53 {
54 const size_t M = A->size1;
55 const size_t N = A->size2;
56
57 if (M != N)
58 {
59 GSL_ERROR ("LDLT decomposition requires square matrix", GSL_ENOTSQR);
60 }
61 else
62 {
63 size_t i, j;
64 double a00, anorm;
65 gsl_vector_view work, v;
66
67 /* check for quick return */
68 if (N == 1)
69 return GSL_SUCCESS;
70
71 /* compute ||A||_1 */
72 anorm = ldlt_norm1(A);
73
74 /* special case first column */
75 a00 = gsl_matrix_get(A, 0, 0);
76 if (a00 == 0.0)
77 {
78 GSL_ERROR ("matrix is singular", GSL_EDOM);
79 }
80
81 v = gsl_matrix_subcolumn(A, 0, 1, N - 1);
82 gsl_vector_scale(&v.vector, 1.0 / a00);
83
84 /* use first subrow A(1, 2:end) as temporary workspace */
85 work = gsl_matrix_subrow(A, 0, 1, N - 1);
86
87 for (j = 1; j < N; ++j)
88 {
89 gsl_vector_view w = gsl_vector_subvector(&work.vector, 0, j);
90 double ajj = gsl_matrix_get(A, j, j);
91 double dval;
92
93 for (i = 0; i < j; ++i)
94 {
95 double aii = gsl_matrix_get(A, i, i);
96 double aji = gsl_matrix_get(A, j, i);
97 gsl_vector_set(&w.vector, i, aji * aii);
98 }
99
100 v = gsl_matrix_subrow(A, j, 0, j); /* A(j,1:j-1) */
101 gsl_blas_ddot(&v.vector, &w.vector, &dval);
102 ajj -= dval;
103
104 if (ajj == 0.0)
105 {
106 GSL_ERROR ("matrix is singular", GSL_EDOM);
107 }
108
109 gsl_matrix_set(A, j, j, ajj);
110
111 if (j < N - 1)
112 {
113 double ajjinv = 1.0 / ajj;
114 gsl_matrix_view m = gsl_matrix_submatrix(A, j + 1, 0, N - j - 1, j); /* A(j+1:n, 1:j-1) */
115 v = gsl_matrix_subcolumn(A, j, j + 1, N - j - 1); /* A(j+1:n, j) */
116 gsl_blas_dgemv(CblasNoTrans, -ajjinv, &m.matrix, &w.vector, ajjinv, &v.vector);
117 }
118 }
119
120 /* save ||A||_1 in upper right corner */
121 gsl_matrix_set(A, 0, N - 1, anorm);
122
123 return GSL_SUCCESS;
124 }
125 }
126
127 int
gsl_linalg_ldlt_solve(const gsl_matrix * LDLT,const gsl_vector * b,gsl_vector * x)128 gsl_linalg_ldlt_solve (const gsl_matrix * LDLT,
129 const gsl_vector * b,
130 gsl_vector * x)
131 {
132 if (LDLT->size1 != LDLT->size2)
133 {
134 GSL_ERROR ("LDLT matrix must be square", GSL_ENOTSQR);
135 }
136 else if (LDLT->size1 != b->size)
137 {
138 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
139 }
140 else if (LDLT->size2 != x->size)
141 {
142 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
143 }
144 else
145 {
146 int status;
147
148 /* copy x <- b */
149 gsl_vector_memcpy (x, b);
150
151 status = gsl_linalg_ldlt_svx(LDLT, x);
152
153 return status;
154 }
155 }
156
157 int
gsl_linalg_ldlt_svx(const gsl_matrix * LDLT,gsl_vector * x)158 gsl_linalg_ldlt_svx (const gsl_matrix * LDLT,
159 gsl_vector * x)
160 {
161 if (LDLT->size1 != LDLT->size2)
162 {
163 GSL_ERROR ("LDLT matrix must be square", GSL_ENOTSQR);
164 }
165 else if (LDLT->size2 != x->size)
166 {
167 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
168 }
169 else
170 {
171 gsl_vector_const_view diag = gsl_matrix_const_diagonal(LDLT);
172
173 /* solve for z using forward-substitution, L z = b */
174 gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasUnit, LDLT, x);
175
176 /* solve for y, D y = z */
177 gsl_vector_div(x, &diag.vector);
178
179 /* perform back-substitution, L^T x = y */
180 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasUnit, LDLT, x);
181
182 return GSL_SUCCESS;
183 }
184 }
185
186 int
gsl_linalg_ldlt_rcond(const gsl_matrix * LDLT,double * rcond,gsl_vector * work)187 gsl_linalg_ldlt_rcond (const gsl_matrix * LDLT, double * rcond,
188 gsl_vector * work)
189 {
190 const size_t M = LDLT->size1;
191 const size_t N = LDLT->size2;
192
193 if (M != N)
194 {
195 GSL_ERROR ("LDLT matrix must be square", GSL_ENOTSQR);
196 }
197 else if (work->size != 3 * N)
198 {
199 GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
200 }
201 else
202 {
203 int status;
204 double Anorm; /* ||A||_1 */
205 double Ainvnorm; /* ||A^{-1}||_1 */
206
207 if (N == 1)
208 Anorm = fabs(gsl_matrix_get(LDLT, 0, 0));
209 else
210 Anorm = gsl_matrix_get(LDLT, 0, N - 1);
211
212 *rcond = 0.0;
213
214 /* don't continue if matrix is singular */
215 if (Anorm == 0.0)
216 return GSL_SUCCESS;
217
218 /* estimate ||A^{-1}||_1 */
219 status = gsl_linalg_invnorm1(N, ldlt_Ainv, (void *) LDLT, &Ainvnorm, work);
220
221 if (status)
222 return status;
223
224 if (Ainvnorm != 0.0)
225 *rcond = (1.0 / Anorm) / Ainvnorm;
226
227 return GSL_SUCCESS;
228 }
229 }
230
231 /* compute 1-norm of symmetric matrix stored in lower triangle */
232 static double
ldlt_norm1(const gsl_matrix * A)233 ldlt_norm1(const gsl_matrix * A)
234 {
235 const size_t N = A->size1;
236 double max = 0.0;
237 size_t i, j;
238
239 for (j = 0; j < N; ++j)
240 {
241 gsl_vector_const_view v = gsl_matrix_const_subcolumn(A, j, j, N - j);
242 double sum = gsl_blas_dasum(&v.vector);
243
244 /* now add symmetric elements from above diagonal */
245 for (i = 0; i < j; ++i)
246 {
247 double Aij = gsl_matrix_get(A, i, j);
248 sum += fabs(Aij);
249 }
250
251 if (sum > max)
252 max = sum;
253 }
254
255 return max;
256 }
257
258 /* x := A^{-1} x = A^{-t} x, A = L D L^T */
259 static int
ldlt_Ainv(CBLAS_TRANSPOSE_t TransA,gsl_vector * x,void * params)260 ldlt_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
261 {
262 int status;
263 gsl_matrix * A = (gsl_matrix * ) params;
264 gsl_vector_const_view diag = gsl_matrix_const_diagonal(A);
265
266 (void) TransA; /* unused parameter warning */
267
268 /* compute L^{-1} x */
269 status = gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasUnit, A, x);
270 if (status)
271 return status;
272
273 /* compute D^{-1} x */
274 gsl_vector_div(x, &diag.vector);
275
276 /* compute L^{-t} x */
277 status = gsl_blas_dtrsv(CblasLower, CblasTrans, CblasUnit, A, x);
278 if (status)
279 return status;
280
281 return GSL_SUCCESS;
282 }
283