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

   1 /* Cholesky Decomposition
   2  *
   3  * Copyright (C) 2000  Thomas Walter
   4  *
   5  * 03 May 2000: Modified for GSL by Brian Gough
   6  * 29 Jul 2005: Additions by Gerard Jungman
   7  * Copyright (C) 2000,2001, 2002, 2003, 2005, 2007 Brian Gough, Gerard Jungman
   8  *
   9  * This is free software; you can redistribute it and/or modify it
  10  * under the terms of the GNU General Public License as published by the
  11  * Free Software Foundation; either version 3, or (at your option) any
  12  * later version.
  13  *
  14  * This source is distributed in the hope that it will be useful, but WITHOUT
  15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  17  * for more details.
  18  */
  19
  20 /*
  21  * Cholesky decomposition of a symmetrix positive definite matrix.
  22  * This is useful to solve the matrix arising in
  23  *    periodic cubic splines
  24  *    approximating splines
  25  *
  26  * This algorithm does:
  27  *   A = L * L'
  28  * with
  29  *   L  := lower left triangle matrix
  30  *   L' := the transposed form of L.
  31  *
  32  */
  33
  34 #include <config.h>
  35
  36 #include <gsl/gsl_math.h>
  37 #include <gsl/gsl_errno.h>
  38 #include <gsl/gsl_vector.h>
  39 #include <gsl/gsl_matrix.h>
  40 #include <gsl/gsl_blas.h>
  41 #include <gsl/gsl_linalg.h>
  42
  43 static inline
  44 double
  45 quiet_sqrt (double x)
  46      /* avoids runtime error, for checking matrix for positive definiteness */
  47 {
  48   return (x >= 0) ? sqrt(x) : GSL_NAN;
  49 }
  50
  51 int
  52 gsl_linalg_cholesky_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("cholesky decomposition requires square matrix", GSL_ENOTSQR);
  60     }
  61   else
  62     {
  63       size_t i,j,k;
  64       int status = 0;
  65
  66       /* Do the first 2 rows explicitly.  It is simple, and faster.  And
  67        * one can return if the matrix has only 1 or 2 rows.
  68        */
  69
  70       double A_00 = gsl_matrix_get (A, 0, 0);
  71
  72       double L_00 = quiet_sqrt(A_00);
  73
  74       if (A_00 <= 0)
  75         {
  76           status = GSL_EDOM ;
  77         }
  78
  79       gsl_matrix_set (A, 0, 0, L_00);
  80
  81       if (M > 1)
  82         {
  83           double A_10 = gsl_matrix_get (A, 1, 0);
  84           double A_11 = gsl_matrix_get (A, 1, 1);
  85
  86           double L_10 = A_10 / L_00;
  87           double diag = A_11 - L_10 * L_10;
  88           double L_11 = quiet_sqrt(diag);
  89
  90           if (diag <= 0)
  91             {
  92               status = GSL_EDOM;
  93             }
  94
  95           gsl_matrix_set (A, 1, 0, L_10);
  96           gsl_matrix_set (A, 1, 1, L_11);
  97         }
  98
  99       for (k = 2; k < M; k++)
 100         {
 101           double A_kk = gsl_matrix_get (A, k, k);
 102
 103           for (i = 0; i < k; i++)
 104             {
 105               double sum = 0;
 106
 107               double A_ki = gsl_matrix_get (A, k, i);
 108               double A_ii = gsl_matrix_get (A, i, i);
 109
 110               gsl_vector_view ci = gsl_matrix_row (A, i);
 111               gsl_vector_view ck = gsl_matrix_row (A, k);
 112
 113               if (i > 0) {
 114                 gsl_vector_view di = gsl_vector_subvector(&ci.vector, 0, i);
 115                 gsl_vector_view dk = gsl_vector_subvector(&ck.vector, 0, i);
 116
 117                 gsl_blas_ddot (&di.vector, &dk.vector, &sum);
 118               }
 119
 120               A_ki = (A_ki - sum) / A_ii;
 121               gsl_matrix_set (A, k, i, A_ki);
 122             }
 123
 124           {
 125             gsl_vector_view ck = gsl_matrix_row (A, k);
 126             gsl_vector_view dk = gsl_vector_subvector (&ck.vector, 0, k);
 127
 128             double sum = gsl_blas_dnrm2 (&dk.vector);
 129             double diag = A_kk - sum * sum;
 130
 131             double L_kk = quiet_sqrt(diag);
 132
 133             if (diag <= 0)
 134               {
 135                 status = GSL_EDOM;
 136               }
 137
 138             gsl_matrix_set (A, k, k, L_kk);
 139           }
 140         }
 141
 142       /* Now copy the transposed lower triangle to the upper triangle,
 143        * the diagonal is common.
 144        */
 145
 146       for (i = 1; i < M; i++)
 147         {
 148           for (j = 0; j < i; j++)
 149             {
 150               double A_ij = gsl_matrix_get (A, i, j);
 151               gsl_matrix_set (A, j, i, A_ij);
 152             }
 153         }
 154
 155       if (status == GSL_EDOM)
 156         {
 157           GSL_ERROR ("matrix must be positive definite", GSL_EDOM);
 158         }
 159
 160       return GSL_SUCCESS;
 161     }
 162 }
 163
 164
 165 int
 166 gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
 167                            const gsl_vector * b,
 168                            gsl_vector * x)
 169 {
 170   if (LLT->size1 != LLT->size2)
 171     {
 172       GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
 173     }
 174   else if (LLT->size1 != b->size)
 175     {
 176       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 177     }
 178   else if (LLT->size2 != x->size)
 179     {
 180       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 181     }
 182   else
 183     {
 184       /* Copy x <- b */
 185
 186       gsl_vector_memcpy (x, b);
 187
 188       /* Solve for c using forward-substitution, L c = b */
 189
 190       gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);
 191
 192       /* Perform back-substitution, U x = c */
 193
 194       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);
 195
 196
 197       return GSL_SUCCESS;
 198     }
 199 }
 200
 201 int
 202 gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
 203                          gsl_vector * x)
 204 {
 205   if (LLT->size1 != LLT->size2)
 206     {
 207       GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
 208     }
 209   else if (LLT->size2 != x->size)
 210     {
 211       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 212     }
 213   else
 214     {
 215       /* Solve for c using forward-substitution, L c = b */
 216
 217       gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);
 218
 219       /* Perform back-substitution, U x = c */
 220
 221       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LLT, x);
 222
 223       return GSL_SUCCESS;
 224     }
 225 }
 226
 227
 228 int
 229 gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
 230 {
 231   const size_t N = A->size1;
 232   size_t i, j;
 233
 234   /* initial Cholesky */
 235   int stat_chol = gsl_linalg_cholesky_decomp(A);
 236
 237   if(stat_chol == GSL_SUCCESS)
 238   {
 239     /* calculate D from diagonal part of initial Cholesky */
 240     for(i = 0; i < N; ++i)
 241     {
 242       const double C_ii = gsl_matrix_get(A, i, i);
 243       gsl_vector_set(D, i, C_ii*C_ii);
 244     }
 245
 246     /* multiply initial Cholesky by 1/sqrt(D) on the right */
 247     for(i = 0; i < N; ++i)
 248     {
 249       for(j = 0; j < N; ++j)
 250       {
 251         gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
 252       }
 253     }
 254
 255     /* Because the initial Cholesky contained both L and transpose(L),
 256        the result of the multiplication is not symmetric anymore;
 257        but the lower triangle _is_ correct. Therefore we reflect
 258        it to the upper triangle and declare victory.
 259        */
 260     for(i = 0; i < N; ++i)
 261       for(j = i + 1; j < N; ++j)
 262         gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
 263   }
 264
 265   return stat_chol;
 266 }