3 * Copyright (C) 2001, 2007 Brian Gough
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.
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.
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.
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_matrix.h>
26 #include <gsl/gsl_complex.h>
27 #include <gsl/gsl_complex_math.h>
28 #include <gsl/gsl_permute_vector.h>
29 #include <gsl/gsl_blas.h>
30 #include <gsl/gsl_complex_math.h>
32 #include <gsl/gsl_linalg.h>
34 /* Factorise a general N x N complex matrix A into,
38 * where P is a permutation matrix, L is unit lower triangular and U
39 * is upper triangular.
41 * L is stored in the strict lower triangular part of the input
42 * matrix. The diagonal elements of L are unity and are not stored.
44 * U is stored in the diagonal and upper triangular part of the
47 * P is stored in the permutation p. Column j of P is column k of the
48 * identity matrix, where k = permutation->data[j]
50 * signum gives the sign of the permutation, (-1)^n, where n is the
51 * number of interchanges in the permutation.
53 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
54 * Elimination with Partial Pivoting).
58 gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)
60 if (A->size1 != A->size2)
62 GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
64 else if (p->size != A->size1)
66 GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
70 const size_t N = A->size1;
74 gsl_permutation_init (p);
76 for (j = 0; j < N - 1; j++)
78 /* Find maximum in the j-th column */
80 gsl_complex ajj = gsl_matrix_complex_get (A, j, j);
81 double max = gsl_complex_abs (ajj);
84 for (i = j + 1; i < N; i++)
86 gsl_complex aij = gsl_matrix_complex_get (A, i, j);
87 double ai = gsl_complex_abs (aij);
98 gsl_matrix_complex_swap_rows (A, j, i_pivot);
99 gsl_permutation_swap (p, j, i_pivot);
100 *signum = -(*signum);
103 ajj = gsl_matrix_complex_get (A, j, j);
105 if (!(GSL_REAL(ajj) == 0.0 && GSL_IMAG(ajj) == 0.0))
107 for (i = j + 1; i < N; i++)
109 gsl_complex aij_orig = gsl_matrix_complex_get (A, i, j);
110 gsl_complex aij = gsl_complex_div (aij_orig, ajj);
111 gsl_matrix_complex_set (A, i, j, aij);
113 for (k = j + 1; k < N; k++)
115 gsl_complex aik = gsl_matrix_complex_get (A, i, k);
116 gsl_complex ajk = gsl_matrix_complex_get (A, j, k);
118 /* aik = aik - aij * ajk */
120 gsl_complex aijajk = gsl_complex_mul (aij, ajk);
121 gsl_complex aik_new = gsl_complex_sub (aik, aijajk);
123 gsl_matrix_complex_set (A, i, k, aik_new);
134 gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
136 if (LU->size1 != LU->size2)
138 GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
140 else if (LU->size1 != p->size)
142 GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
144 else if (LU->size1 != b->size)
146 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
148 else if (LU->size2 != x->size)
150 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
156 gsl_vector_complex_memcpy (x, b);
160 gsl_linalg_complex_LU_svx (LU, p, x);
168 gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
170 if (LU->size1 != LU->size2)
172 GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
174 else if (LU->size1 != p->size)
176 GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
178 else if (LU->size1 != x->size)
180 GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
184 /* Apply permutation to RHS */
186 gsl_permute_vector_complex (p, x);
188 /* Solve for c using forward-substitution, L c = P b */
190 gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);
192 /* Perform back-substitution, U x = c */
194 gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);
202 gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * residual)
204 if (A->size1 != A->size2)
206 GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
208 if (LU->size1 != LU->size2)
210 GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
212 else if (A->size1 != LU->size2)
214 GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
216 else if (LU->size1 != p->size)
218 GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
220 else if (LU->size1 != b->size)
222 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
224 else if (LU->size1 != x->size)
226 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
230 /* Compute residual, residual = (A * x - b) */
232 gsl_vector_complex_memcpy (residual, b);
235 gsl_complex one = GSL_COMPLEX_ONE;
236 gsl_complex negone = GSL_COMPLEX_NEGONE;
237 gsl_blas_zgemv (CblasNoTrans, one, A, x, negone, residual);
240 /* Find correction, delta = - (A^-1) * residual, and apply it */
242 gsl_linalg_complex_LU_svx (LU, p, residual);
245 gsl_complex negone= GSL_COMPLEX_NEGONE;
246 gsl_blas_zaxpy (negone, residual, x);
254 gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
256 size_t i, n = LU->size1;
258 int status = GSL_SUCCESS;
260 gsl_matrix_complex_set_identity (inverse);
262 for (i = 0; i < n; i++)
264 gsl_vector_complex_view c = gsl_matrix_complex_column (inverse, i);
265 int status_i = gsl_linalg_complex_LU_svx (LU, p, &(c.vector));
275 gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
277 size_t i, n = LU->size1;
279 gsl_complex det = gsl_complex_rect((double) signum, 0.0);
281 for (i = 0; i < n; i++)
283 gsl_complex zi = gsl_matrix_complex_get (LU, i, i);
284 det = gsl_complex_mul (det, zi);
292 gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
294 size_t i, n = LU->size1;
298 for (i = 0; i < n; i++)
300 gsl_complex z = gsl_matrix_complex_get (LU, i, i);
301 lndet += log (gsl_complex_abs (z));
309 gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
311 size_t i, n = LU->size1;
313 gsl_complex phase = gsl_complex_rect((double) signum, 0.0);
315 for (i = 0; i < n; i++)
317 gsl_complex z = gsl_matrix_complex_get (LU, i, i);
319 double r = gsl_complex_abs(z);
323 phase = gsl_complex_rect(0.0, 0.0);
328 z = gsl_complex_div_real(z, r);
329 phase = gsl_complex_mul(phase, z);