htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/luc.c

   1 /* linalg/luc.c
   2  *
   3  * Copyright (C) 2001, 2007 Brian Gough
   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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <string.h>
  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>
  31
  32 #include <gsl/gsl_linalg.h>
  33
  34 /* Factorise a general N x N complex matrix A into,
  35  *
  36  *   P A = L U
  37  *
  38  * where P is a permutation matrix, L is unit lower triangular and U
  39  * is upper triangular.
  40  *
  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.
  43  *
  44  * U is stored in the diagonal and upper triangular part of the
  45  * input matrix.
  46  *
  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]
  49  *
  50  * signum gives the sign of the permutation, (-1)^n, where n is the
  51  * number of interchanges in the permutation.
  52  *
  53  * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
  54  * Elimination with Partial Pivoting).
  55  */
  56
  57 int
  58 gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)
  59 {
  60   if (A->size1 != A->size2)
  61     {
  62       GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
  63     }
  64   else if (p->size != A->size1)
  65     {
  66       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
  67     }
  68   else
  69     {
  70       const size_t N = A->size1;
  71       size_t i, j, k;
  72
  73       *signum = 1;
  74       gsl_permutation_init (p);
  75
  76       for (j = 0; j < N - 1; j++)
  77         {
  78           /* Find maximum in the j-th column */
  79
  80           gsl_complex ajj = gsl_matrix_complex_get (A, j, j);
  81           double max = gsl_complex_abs (ajj);
  82           size_t i_pivot = j;
  83
  84           for (i = j + 1; i < N; i++)
  85             {
  86               gsl_complex aij = gsl_matrix_complex_get (A, i, j);
  87               double ai = gsl_complex_abs (aij);
  88
  89               if (ai > max)
  90                 {
  91                   max = ai;
  92                   i_pivot = i;
  93                 }
  94             }
  95
  96           if (i_pivot != j)
  97             {
  98               gsl_matrix_complex_swap_rows (A, j, i_pivot);
  99               gsl_permutation_swap (p, j, i_pivot);
 100               *signum = -(*signum);
 101             }
 102
 103           ajj = gsl_matrix_complex_get (A, j, j);
 104
 105           if (!(GSL_REAL(ajj) == 0.0 && GSL_IMAG(ajj) == 0.0))
 106             {
 107               for (i = j + 1; i < N; i++)
 108                 {
 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);
 112
 113                   for (k = j + 1; k < N; k++)
 114                     {
 115                       gsl_complex aik = gsl_matrix_complex_get (A, i, k);
 116                       gsl_complex ajk = gsl_matrix_complex_get (A, j, k);
 117
 118                       /* aik = aik - aij * ajk */
 119
 120                       gsl_complex aijajk = gsl_complex_mul (aij, ajk);
 121                       gsl_complex aik_new = gsl_complex_sub (aik, aijajk);
 122
 123                       gsl_matrix_complex_set (A, i, k, aik_new);
 124                     }
 125                 }
 126             }
 127         }
 128
 129       return GSL_SUCCESS;
 130     }
 131 }
 132
 133 int
 134 gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
 135 {
 136   if (LU->size1 != LU->size2)
 137     {
 138       GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
 139     }
 140   else if (LU->size1 != p->size)
 141     {
 142       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
 143     }
 144   else if (LU->size1 != b->size)
 145     {
 146       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 147     }
 148   else if (LU->size2 != x->size)
 149     {
 150       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 151     }
 152   else
 153     {
 154       /* Copy x <- b */
 155
 156       gsl_vector_complex_memcpy (x, b);
 157
 158       /* Solve for x */
 159
 160       gsl_linalg_complex_LU_svx (LU, p, x);
 161
 162       return GSL_SUCCESS;
 163     }
 164 }
 165
 166
 167 int
 168 gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
 169 {
 170   if (LU->size1 != LU->size2)
 171     {
 172       GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
 173     }
 174   else if (LU->size1 != p->size)
 175     {
 176       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
 177     }
 178   else if (LU->size1 != x->size)
 179     {
 180       GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
 181     }
 182   else
 183     {
 184       /* Apply permutation to RHS */
 185
 186       gsl_permute_vector_complex (p, x);
 187
 188       /* Solve for c using forward-substitution, L c = P b */
 189
 190       gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);
 191
 192       /* Perform back-substitution, U x = c */
 193
 194       gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);
 195
 196       return GSL_SUCCESS;
 197     }
 198 }
 199
 200
 201 int
 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)
 203 {
 204   if (A->size1 != A->size2)
 205     {
 206       GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
 207     }
 208   if (LU->size1 != LU->size2)
 209     {
 210       GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
 211     }
 212   else if (A->size1 != LU->size2)
 213     {
 214       GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
 215     }
 216   else if (LU->size1 != p->size)
 217     {
 218       GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
 219     }
 220   else if (LU->size1 != b->size)
 221     {
 222       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 223     }
 224   else if (LU->size1 != x->size)
 225     {
 226       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 227     }
 228   else
 229     {
 230       /* Compute residual, residual = (A * x  - b) */
 231
 232       gsl_vector_complex_memcpy (residual, b);
 233
 234       {
 235         gsl_complex one = GSL_COMPLEX_ONE;
 236         gsl_complex negone = GSL_COMPLEX_NEGONE;
 237         gsl_blas_zgemv (CblasNoTrans, one, A, x, negone, residual);
 238       }
 239
 240       /* Find correction, delta = - (A^-1) * residual, and apply it */
 241
 242       gsl_linalg_complex_LU_svx (LU, p, residual);
 243
 244       {
 245         gsl_complex negone= GSL_COMPLEX_NEGONE;
 246         gsl_blas_zaxpy (negone, residual, x);
 247       }
 248
 249       return GSL_SUCCESS;
 250     }
 251 }
 252
 253 int
 254 gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
 255 {
 256   size_t i, n = LU->size1;
 257
 258   int status = GSL_SUCCESS;
 259
 260   gsl_matrix_complex_set_identity (inverse);
 261
 262   for (i = 0; i < n; i++)
 263     {
 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));
 266
 267       if (status_i)
 268         status = status_i;
 269     }
 270
 271   return status;
 272 }
 273
 274 gsl_complex
 275 gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
 276 {
 277   size_t i, n = LU->size1;
 278
 279   gsl_complex det = gsl_complex_rect((double) signum, 0.0);
 280
 281   for (i = 0; i < n; i++)
 282     {
 283       gsl_complex zi = gsl_matrix_complex_get (LU, i, i);
 284       det = gsl_complex_mul (det, zi);
 285     }
 286
 287   return det;
 288 }
 289
 290
 291 double
 292 gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
 293 {
 294   size_t i, n = LU->size1;
 295
 296   double lndet = 0.0;
 297
 298   for (i = 0; i < n; i++)
 299     {
 300       gsl_complex z = gsl_matrix_complex_get (LU, i, i);
 301       lndet += log (gsl_complex_abs (z));
 302     }
 303
 304   return lndet;
 305 }
 306
 307
 308 gsl_complex
 309 gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
 310 {
 311   size_t i, n = LU->size1;
 312
 313   gsl_complex phase = gsl_complex_rect((double) signum, 0.0);
 314
 315   for (i = 0; i < n; i++)
 316     {
 317       gsl_complex z = gsl_matrix_complex_get (LU, i, i);
 318
 319       double r = gsl_complex_abs(z);
 320
 321       if (r == 0)
 322         {
 323           phase = gsl_complex_rect(0.0, 0.0);
 324           break;
 325         }
 326       else
 327         {
 328           z = gsl_complex_div_real(z, r);
 329           phase = gsl_complex_mul(phase, z);
 330         }
 331     }
 332
 333   return phase;
 334 }