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

   1 /* linalg/choleskyc.c
   2  *
   3  * Copyright (C) 2007 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 #include <config.h>
  21 #include <gsl/gsl_matrix.h>
  22 #include <gsl/gsl_vector.h>
  23 #include <gsl/gsl_linalg.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_complex.h>
  26 #include <gsl/gsl_complex_math.h>
  27 #include <gsl/gsl_blas.h>
  28 #include <gsl/gsl_errno.h>
  29
  30 /*
  31  * This module contains routines related to the Cholesky decomposition
  32  * of a complex Hermitian positive definite matrix.
  33  */
  34
  35 static void cholesky_complex_conj_vector(gsl_vector_complex *v);
  36
  37 /*
  38 gsl_linalg_complex_cholesky_decomp()
  39   Perform the Cholesky decomposition on a Hermitian positive definite
  40 matrix. See Golub & Van Loan, "Matrix Computations" (3rd ed),
  41 algorithm 4.2.2.
  42
  43 Inputs: A - (input/output) complex postive definite matrix
  44
  45 Return: success or error
  46
  47 The lower triangle of A is overwritten with the Cholesky decomposition
  48 */
  49
  50 int
  51 gsl_linalg_complex_cholesky_decomp(gsl_matrix_complex *A)
  52 {
  53   const size_t N = A->size1;
  54
  55   if (N != A->size2)
  56     {
  57       GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
  58     }
  59   else
  60     {
  61       size_t i, j;
  62       gsl_complex z;
  63       double ajj;
  64
  65       for (j = 0; j < N; ++j)
  66         {
  67           z = gsl_matrix_complex_get(A, j, j);
  68           ajj = GSL_REAL(z);
  69
  70           if (j > 0)
  71             {
  72               gsl_vector_complex_const_view aj =
  73                 gsl_matrix_complex_const_subrow(A, j, 0, j);
  74
  75               gsl_blas_zdotc(&aj.vector, &aj.vector, &z);
  76               ajj -= GSL_REAL(z);
  77             }
  78
  79           if (ajj <= 0.0)
  80             {
  81               GSL_ERROR("matrix is not positive definite", GSL_EDOM);
  82             }
  83
  84           ajj = sqrt(ajj);
  85           GSL_SET_COMPLEX(&z, ajj, 0.0);
  86           gsl_matrix_complex_set(A, j, j, z);
  87
  88           if (j < N - 1)
  89             {
  90               gsl_vector_complex_view av =
  91                 gsl_matrix_complex_subcolumn(A, j, j + 1, N - j - 1);
  92
  93               if (j > 0)
  94                 {
  95                   gsl_vector_complex_view aj =
  96                     gsl_matrix_complex_subrow(A, j, 0, j);
  97                   gsl_matrix_complex_view am =
  98                     gsl_matrix_complex_submatrix(A, j + 1, 0, N - j - 1, j);
  99
 100                   cholesky_complex_conj_vector(&aj.vector);
 101
 102                   gsl_blas_zgemv(CblasNoTrans,
 103                                  GSL_COMPLEX_NEGONE,
 104                                  &am.matrix,
 105                                  &aj.vector,
 106                                  GSL_COMPLEX_ONE,
 107                                  &av.vector);
 108
 109                   cholesky_complex_conj_vector(&aj.vector);
 110                 }
 111
 112               gsl_blas_zdscal(1.0 / ajj, &av.vector);
 113             }
 114         }
 115
 116       /* Now store L^H in upper triangle */
 117       for (i = 1; i < N; ++i)
 118         {
 119           for (j = 0; j < i; ++j)
 120             {
 121               z = gsl_matrix_complex_get(A, i, j);
 122               gsl_matrix_complex_set(A, j, i, gsl_complex_conjugate(z));
 123             }
 124         }
 125
 126       return GSL_SUCCESS;
 127     }
 128 } /* gsl_linalg_complex_cholesky_decomp() */
 129
 130 /*
 131 gsl_linalg_complex_cholesky_solve()
 132   Solve A x = b where A is in cholesky form
 133 */
 134
 135 int
 136 gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
 137                                    const gsl_vector_complex * b,
 138                                    gsl_vector_complex * x)
 139 {
 140   if (cholesky->size1 != cholesky->size2)
 141     {
 142       GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
 143     }
 144   else if (cholesky->size1 != b->size)
 145     {
 146       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 147     }
 148   else if (cholesky->size2 != x->size)
 149     {
 150       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 151     }
 152   else
 153     {
 154       gsl_vector_complex_memcpy (x, b);
 155
 156       /* solve for y using forward-substitution, L y = b */
 157
 158       gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);
 159
 160       /* perform back-substitution, L^H x = y */
 161
 162       gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);
 163
 164       return GSL_SUCCESS;
 165     }
 166 } /* gsl_linalg_complex_cholesky_solve() */
 167
 168 /*
 169 gsl_linalg_complex_cholesky_svx()
 170   Solve A x = b in place where A is in cholesky form
 171 */
 172
 173 int
 174 gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
 175                                  gsl_vector_complex * x)
 176 {
 177   if (cholesky->size1 != cholesky->size2)
 178     {
 179       GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
 180     }
 181   else if (cholesky->size2 != x->size)
 182     {
 183       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 184     }
 185   else
 186     {
 187       /* solve for y using forward-substitution, L y = b */
 188
 189       gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);
 190
 191       /* perform back-substitution, L^H x = y */
 192
 193       gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);
 194
 195       return GSL_SUCCESS;
 196     }
 197 } /* gsl_linalg_complex_cholesky_svx() */
 198
 199 /********************************************
 200  *           INTERNAL ROUTINES              *
 201  ********************************************/
 202
 203 static void
 204 cholesky_complex_conj_vector(gsl_vector_complex *v)
 205 {
 206   size_t i;
 207
 208   for (i = 0; i < v->size; ++i)
 209     {
 210       gsl_complex z = gsl_vector_complex_get(v, i);
 211       gsl_vector_complex_set(v, i, gsl_complex_conjugate(z));
 212     }
 213 } /* cholesky_complex_conj_vector() */