util/X86MAC64/include/gsl/gsl_linalg.h

   1 /* linalg/gsl_linalg.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007, 2019 Gerard Jungman, Brian Gough, 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 #ifndef __GSL_LINALG_H__
  21 #define __GSL_LINALG_H__
  22
  23 #include <stdlib.h>
  24 #include <gsl/gsl_mode.h>
  25 #include <gsl/gsl_permutation.h>
  26 #include <gsl/gsl_vector.h>
  27 #include <gsl/gsl_matrix.h>
  28 #include <gsl/gsl_math.h>
  29 #include <gsl/gsl_inline.h>
  30 #include <gsl/gsl_blas.h>
  31
  32 #undef __BEGIN_DECLS
  33 #undef __END_DECLS
  34 #ifdef __cplusplus
  35 #define __BEGIN_DECLS extern "C" {
  36 #define __END_DECLS }
  37 #else
  38 #define __BEGIN_DECLS           /* empty */
  39 #define __END_DECLS             /* empty */
  40 #endif
  41
  42 __BEGIN_DECLS
  43
  44 typedef enum
  45   {
  46     GSL_LINALG_MOD_NONE = 0,
  47     GSL_LINALG_MOD_TRANSPOSE = 1,
  48     GSL_LINALG_MOD_CONJUGATE = 2
  49   }
  50 gsl_linalg_matrix_mod_t;
  51
  52 /* Note: You can now use the gsl_blas_dgemm function instead of matmult */
  53
  54 /* Simple implementation of matrix multiply.
  55  * Calculates C = A.B
  56  *
  57  * exceptions: GSL_EBADLEN
  58  */
  59 int gsl_linalg_matmult (const gsl_matrix * A,
  60                         const gsl_matrix * B,
  61                         gsl_matrix * C);
  62
  63
  64 /* Simple implementation of matrix multiply.
  65  * Allows transposition of either matrix, so it
  66  * can compute A.B or Trans(A).B or A.Trans(B) or Trans(A).Trans(B)
  67  *
  68  * exceptions: GSL_EBADLEN
  69  */
  70 int gsl_linalg_matmult_mod (const gsl_matrix * A,
  71                             gsl_linalg_matrix_mod_t modA,
  72                             const gsl_matrix * B,
  73                             gsl_linalg_matrix_mod_t modB,
  74                             gsl_matrix * C);
  75
  76 /* Calculate the matrix exponential by the scaling and
  77  * squaring method described in Moler + Van Loan,
  78  * SIAM Rev 20, 801 (1978). The mode argument allows
  79  * choosing an optimal strategy, from the table
  80  * given in the paper, for a given precision.
  81  *
  82  * exceptions: GSL_ENOTSQR, GSL_EBADLEN
  83  */
  84 int gsl_linalg_exponential_ss(
  85   const gsl_matrix * A,
  86   gsl_matrix * eA,
  87   gsl_mode_t mode
  88   );
  89
  90
  91 /* Householder Transformations */
  92
  93 double gsl_linalg_householder_transform (gsl_vector * v);
  94 double gsl_linalg_householder_transform2 (double * alpha, gsl_vector * v);
  95 gsl_complex gsl_linalg_complex_householder_transform (gsl_vector_complex * v);
  96
  97 int gsl_linalg_householder_hm (double tau,
  98                                const gsl_vector * v,
  99                                gsl_matrix * A);
 100
 101 int gsl_linalg_householder_mh (double tau,
 102                                const gsl_vector * v,
 103                                gsl_matrix * A);
 104
 105 int gsl_linalg_householder_hv (double tau,
 106                                const gsl_vector * v,
 107                                gsl_vector * w);
 108
 109 int gsl_linalg_householder_left(const double tau,
 110                                 const gsl_vector * v,
 111                                 gsl_matrix * A,
 112                                 gsl_vector * work);
 113
 114 int gsl_linalg_householder_right(const double tau,
 115                                  const gsl_vector * v,
 116                                  gsl_matrix * A,
 117                                  gsl_vector * work);
 118
 119 int gsl_linalg_householder_hm1 (double tau,
 120                                 gsl_matrix * A);
 121
 122 int gsl_linalg_complex_householder_hm (gsl_complex tau,
 123                                        const gsl_vector_complex * v,
 124                                        gsl_matrix_complex * A);
 125
 126 int gsl_linalg_complex_householder_mh (gsl_complex tau,
 127                                        const gsl_vector_complex * v,
 128                                        gsl_matrix_complex * A);
 129
 130 int gsl_linalg_complex_householder_hv (gsl_complex tau,
 131                                        const gsl_vector_complex * v,
 132                                        gsl_vector_complex * w);
 133
 134 int gsl_linalg_complex_householder_left (const gsl_complex tau,
 135                                          const gsl_vector_complex * v,
 136                                          gsl_matrix_complex * A,
 137                                          gsl_vector_complex * work);
 138
 139 /* Hessenberg reduction */
 140
 141 int gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau);
 142 int gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
 143                                  gsl_matrix * U);
 144 int gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
 145                                        gsl_matrix * U);
 146 int gsl_linalg_hessenberg_set_zero(gsl_matrix * H);
 147 int gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A,
 148                                     size_t top, gsl_vector *tau);
 149
 150 /* Hessenberg-Triangular reduction */
 151
 152 int gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B,
 153                               gsl_matrix * U, gsl_matrix * V,
 154                               gsl_vector * work);
 155
 156 /* Singular Value Decomposition
 157
 158  * exceptions:
 159  */
 160
 161 int
 162 gsl_linalg_SV_decomp (gsl_matrix * A,
 163                       gsl_matrix * V,
 164                       gsl_vector * S,
 165                       gsl_vector * work);
 166
 167 int
 168 gsl_linalg_SV_decomp_mod (gsl_matrix * A,
 169                           gsl_matrix * X,
 170                           gsl_matrix * V,
 171                           gsl_vector * S,
 172                           gsl_vector * work);
 173
 174 int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A,
 175                                  gsl_matrix * Q,
 176                                  gsl_vector * S);
 177
 178 int
 179 gsl_linalg_SV_solve (const gsl_matrix * U,
 180                      const gsl_matrix * Q,
 181                      const gsl_vector * S,
 182                      const gsl_vector * b,
 183                      gsl_vector * x);
 184
 185 int gsl_linalg_SV_leverage(const gsl_matrix *U, gsl_vector *h);
 186
 187
 188 /* LU Decomposition, Gaussian elimination with partial pivoting
 189  */
 190
 191 int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum);
 192
 193 int gsl_linalg_LU_solve (const gsl_matrix * LU,
 194                          const gsl_permutation * p,
 195                          const gsl_vector * b,
 196                          gsl_vector * x);
 197
 198 int gsl_linalg_LU_svx (const gsl_matrix * LU,
 199                        const gsl_permutation * p,
 200                        gsl_vector * x);
 201
 202 int gsl_linalg_LU_refine (const gsl_matrix * A,
 203                           const gsl_matrix * LU,
 204                           const gsl_permutation * p,
 205                           const gsl_vector * b,
 206                           gsl_vector * x,
 207                           gsl_vector * work);
 208
 209 int gsl_linalg_LU_invert (const gsl_matrix * LU,
 210                           const gsl_permutation * p,
 211                           gsl_matrix * inverse);
 212 int gsl_linalg_LU_invx (gsl_matrix * LU, const gsl_permutation * p);
 213
 214 double gsl_linalg_LU_det (gsl_matrix * LU, int signum);
 215 double gsl_linalg_LU_lndet (gsl_matrix * LU);
 216 int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum);
 217
 218 /* Banded LU decomposition */
 219
 220 int gsl_linalg_LU_band_decomp (const size_t M, const size_t lb, const size_t ub, gsl_matrix * AB, gsl_vector_uint * piv);
 221
 222 int gsl_linalg_LU_band_solve (const size_t lb, const size_t ub, const gsl_matrix * LUB,
 223                               const gsl_vector_uint * piv, const gsl_vector * b, gsl_vector * x);
 224
 225 int gsl_linalg_LU_band_svx (const size_t lb, const size_t ub, const gsl_matrix * LUB,
 226                             const gsl_vector_uint * piv, gsl_vector * x);
 227
 228 int gsl_linalg_LU_band_unpack (const size_t M, const size_t lb, const size_t ub, const gsl_matrix * LUB,
 229                                const gsl_vector_uint * piv, gsl_matrix * L, gsl_matrix * U);
 230
 231 /* Complex LU Decomposition */
 232
 233 int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A,
 234                                   gsl_permutation * p,
 235                                   int *signum);
 236
 237 int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU,
 238                                  const gsl_permutation * p,
 239                                  const gsl_vector_complex * b,
 240                                  gsl_vector_complex * x);
 241
 242 int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU,
 243                                const gsl_permutation * p,
 244                                gsl_vector_complex * x);
 245
 246 int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A,
 247                                   const gsl_matrix_complex * LU,
 248                                   const gsl_permutation * p,
 249                                   const gsl_vector_complex * b,
 250                                   gsl_vector_complex * x,
 251                                   gsl_vector_complex * work);
 252
 253 int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU,
 254                                   const gsl_permutation * p,
 255                                   gsl_matrix_complex * inverse);
 256 int gsl_linalg_complex_LU_invx (gsl_matrix_complex * LU, const gsl_permutation * p);
 257
 258 gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU,
 259                                        int signum);
 260
 261 double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU);
 262
 263 gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU,
 264                                           int signum);
 265
 266 /* QR decomposition */
 267
 268 int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau);
 269
 270 int gsl_linalg_QR_decomp_old (gsl_matrix * A, gsl_vector * tau);
 271
 272 int gsl_linalg_QR_decomp_r (gsl_matrix * A, gsl_matrix * T);
 273
 274 int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x);
 275
 276 int gsl_linalg_QR_solve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x);
 277
 278 int gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x);
 279
 280 int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b,
 281                            gsl_vector * x, gsl_vector * residual);
 282
 283 int gsl_linalg_QR_lssolve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b,
 284                              gsl_vector * x, gsl_vector * work);
 285
 286 int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x);
 287
 288 int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x);
 289
 290 int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x);
 291
 292 int gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R, gsl_vector * w, const gsl_vector * v);
 293
 294 int gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);
 295
 296 int gsl_linalg_QR_QTvec_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
 297
 298 int gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);
 299
 300 int gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);
 301
 302 int gsl_linalg_QR_QTmat_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * B, gsl_matrix * work);
 303
 304 int gsl_linalg_QR_matQ (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);
 305
 306 int gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R);
 307
 308 int gsl_linalg_QR_unpack_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q, gsl_matrix * R);
 309
 310 int gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x);
 311
 312 int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x);
 313
 314 int gsl_linalg_QR_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work);
 315
 316 /* complex QR decomposition */
 317
 318 int gsl_linalg_complex_QR_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau);
 319
 320 int gsl_linalg_complex_QR_decomp_r (gsl_matrix_complex * A, gsl_matrix_complex * T);
 321
 322 int gsl_linalg_complex_QR_solve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
 323                                  const gsl_vector_complex * b, gsl_vector_complex * x);
 324
 325 int gsl_linalg_complex_QR_solve_r (const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
 326                                    const gsl_vector_complex * b, gsl_vector_complex * x);
 327
 328 int gsl_linalg_complex_QR_svx (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * x);
 329
 330 int gsl_linalg_complex_QR_lssolve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
 331                                    const gsl_vector_complex * b, gsl_vector_complex * x,
 332                                    gsl_vector_complex * residual);
 333
 334 int gsl_linalg_complex_QR_lssolve_r (const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
 335                                      const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * work);
 336
 337 int gsl_linalg_complex_QR_QHvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v);
 338
 339 int gsl_linalg_complex_QR_QHvec_r(const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
 340                                   gsl_vector_complex * b, gsl_vector_complex * work);
 341
 342 int gsl_linalg_complex_QR_Qvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v);
 343
 344 int gsl_linalg_complex_QR_unpack (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
 345                                   gsl_matrix_complex * Q, gsl_matrix_complex * R);
 346
 347 int gsl_linalg_complex_QR_unpack_r(const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
 348                                    gsl_matrix_complex * Q, gsl_matrix_complex * R);
 349
 350 /* banded QR decomposition */
 351
 352 int gsl_linalg_QR_band_decomp_L2 (const size_t M, const size_t p, const size_t q,
 353                                   gsl_matrix * AB, gsl_vector * tau);
 354
 355 int gsl_linalg_QR_band_unpack_L2 (const size_t p, const size_t q, const gsl_matrix * QRB,
 356                                   const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R);
 357
 358 /* Q R P^T decomposition */
 359
 360 int gsl_linalg_QRPT_decomp (gsl_matrix * A,
 361                             gsl_vector * tau,
 362                             gsl_permutation * p,
 363                             int *signum,
 364                             gsl_vector * norm);
 365
 366 int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
 367                              gsl_matrix * q, gsl_matrix * r,
 368                              gsl_vector * tau,
 369                              gsl_permutation * p,
 370                              int *signum,
 371                              gsl_vector * norm);
 372
 373 int gsl_linalg_QRPT_solve (const gsl_matrix * QR,
 374                            const gsl_vector * tau,
 375                            const gsl_permutation * p,
 376                            const gsl_vector * b,
 377                            gsl_vector * x);
 378
 379 int gsl_linalg_QRPT_lssolve (const gsl_matrix * QR,
 380                              const gsl_vector * tau,
 381                              const gsl_permutation * p,
 382                              const gsl_vector * b,
 383                              gsl_vector * x,
 384                              gsl_vector * residual);
 385
 386 int gsl_linalg_QRPT_lssolve2 (const gsl_matrix * QR,
 387                               const gsl_vector * tau,
 388                               const gsl_permutation * p,
 389                               const gsl_vector * b,
 390                               const size_t rank,
 391                               gsl_vector * x,
 392                               gsl_vector * residual);
 393
 394 int gsl_linalg_QRPT_svx (const gsl_matrix * QR,
 395                          const gsl_vector * tau,
 396                          const gsl_permutation * p,
 397                          gsl_vector * x);
 398
 399 int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,
 400                              const gsl_matrix * R,
 401                              const gsl_permutation * p,
 402                              const gsl_vector * b,
 403                              gsl_vector * x);
 404
 405 int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
 406                              const gsl_permutation * p,
 407                              const gsl_vector * b,
 408                              gsl_vector * x);
 409
 410 int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
 411                            const gsl_permutation * p,
 412                            gsl_vector * x);
 413
 414 int gsl_linalg_QRPT_update (gsl_matrix * Q,
 415                             gsl_matrix * R,
 416                             const gsl_permutation * p,
 417                             gsl_vector * u,
 418                             const gsl_vector * v);
 419
 420 size_t gsl_linalg_QRPT_rank (const gsl_matrix * QR, const double tol);
 421
 422 int gsl_linalg_QRPT_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work);
 423
 424 /* triangle on top of diagonal QR decomposition */
 425
 426 int gsl_linalg_QR_UD_decomp (gsl_matrix * U, const gsl_vector * D, gsl_matrix * Y, gsl_matrix * T);
 427
 428 int gsl_linalg_QR_UD_lssolve (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
 429                               const gsl_vector * b, gsl_vector * x, gsl_vector * work);
 430
 431 /* triangle on top of rectangle QR decomposition */
 432
 433 int gsl_linalg_QR_UR_decomp (gsl_matrix * S, gsl_matrix * A, gsl_matrix * T);
 434
 435 /* triangle on top of triangle QR decomposition */
 436
 437 int gsl_linalg_QR_UU_decomp (gsl_matrix * U, gsl_matrix * S, gsl_matrix * T);
 438
 439 int gsl_linalg_QR_UU_lssolve (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
 440                               const gsl_vector * b, gsl_vector * x, gsl_vector * work);
 441
 442 int gsl_linalg_QR_UU_QTvec(const gsl_matrix * Y, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
 443
 444 /* triangle on top of trapezoidal QR decomposition */
 445
 446 int gsl_linalg_QR_UZ_decomp (gsl_matrix * S, gsl_matrix * A, gsl_matrix * T);
 447
 448 /* QL decomposition */
 449
 450 int gsl_linalg_QL_decomp (gsl_matrix * A, gsl_vector * tau);
 451
 452 int gsl_linalg_QL_unpack (const gsl_matrix * QL, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L);
 453
 454 /* COD decomposition */
 455
 456 int gsl_linalg_COD_decomp(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
 457                           gsl_permutation * p, size_t * rank, gsl_vector * work);
 458
 459 int gsl_linalg_COD_decomp_e(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
 460                             gsl_permutation * p, double tol, size_t * rank, gsl_vector * work);
 461
 462 int gsl_linalg_COD_lssolve (const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
 463                             const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
 464                             gsl_vector * x, gsl_vector * residual);
 465
 466 int
 467 gsl_linalg_COD_lssolve2 (const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
 468                          const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
 469                          gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work);
 470
 471 int gsl_linalg_COD_unpack(const gsl_matrix * QRZT, const gsl_vector * tau_Q,
 472                           const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q,
 473                           gsl_matrix * R, gsl_matrix * Z);
 474
 475 int gsl_linalg_COD_matZ(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
 476                         gsl_matrix * A, gsl_vector * work);
 477
 478 /* LQ decomposition */
 479
 480 int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);
 481
 482 int gsl_linalg_LQ_lssolve(const gsl_matrix * LQ, const gsl_vector * tau,
 483                           const gsl_vector * b, gsl_vector * x, gsl_vector * residual);
 484
 485 int gsl_linalg_LQ_QTvec(const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v);
 486
 487 int gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau,
 488                          const gsl_vector * b, gsl_vector * x);
 489
 490 int gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau,
 491                          gsl_vector * x);
 492
 493 int gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau,
 494                            const gsl_vector * b, gsl_vector * x,
 495                            gsl_vector * residual);
 496
 497 int gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b,
 498                           gsl_vector * x);
 499
 500 int gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x);
 501
 502 int gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b,
 503                         gsl_vector * x);
 504
 505 int gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau,
 506                         gsl_vector * v);
 507
 508 int gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau,
 509                          gsl_vector * v);
 510
 511 int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
 512                           gsl_matrix * Q, gsl_matrix * L);
 513
 514 int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
 515                           const gsl_vector * v, gsl_vector * w);
 516 int gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L,
 517                            const gsl_vector * b, gsl_vector * x);
 518
 519 /* P^T L Q decomposition */
 520
 521 int gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau,
 522                             gsl_permutation * p, int *signum,
 523                             gsl_vector * norm);
 524
 525 int gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q,
 526                              gsl_matrix * r, gsl_vector * tau,
 527                              gsl_permutation * p, int *signum,
 528                              gsl_vector * norm);
 529
 530 int gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
 531                            const gsl_vector * tau,
 532                            const gsl_permutation * p,
 533                            const gsl_vector * b,
 534                            gsl_vector * x);
 535
 536 int gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
 537                            const gsl_vector * tau,
 538                            const gsl_permutation * p,
 539                            gsl_vector * x);
 540
 541 int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
 542                              const gsl_permutation * p,
 543                              const gsl_vector * b,
 544                              gsl_vector * x);
 545
 546 int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
 547                             const gsl_permutation * p,
 548                             const gsl_vector * b,
 549                             gsl_vector * x);
 550
 551 int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
 552                           const gsl_permutation * p,
 553                           gsl_vector * x);
 554
 555 int gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
 556                             const gsl_permutation * p,
 557                             const gsl_vector * v, gsl_vector * w);
 558
 559 /* Cholesky Decomposition */
 560
 561 int gsl_linalg_cholesky_decomp (gsl_matrix * A);
 562 int gsl_linalg_cholesky_decomp1 (gsl_matrix * A);
 563
 564 int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky,
 565                                const gsl_vector * b,
 566                                gsl_vector * x);
 567 int gsl_linalg_cholesky_solve_mat (const gsl_matrix * cholesky,
 568                                    const gsl_matrix * B,
 569                                    gsl_matrix * X);
 570
 571 int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky,
 572                              gsl_vector * x);
 573 int gsl_linalg_cholesky_svx_mat (const gsl_matrix * cholesky,
 574                                  gsl_matrix * X);
 575
 576 int gsl_linalg_cholesky_invert(gsl_matrix * cholesky);
 577
 578 /* Cholesky decomposition with unit-diagonal triangular parts.
 579  *   A = L D L^T, where diag(L) = (1,1,...,1).
 580  *   Upon exit, A contains L and L^T as for Cholesky, and
 581  *   the diagonal of A is (1,1,...,1). The vector Dis set
 582  *   to the diagonal elements of the diagonal matrix D.
 583  */
 584 int gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D);
 585
 586 int gsl_linalg_cholesky_scale(const gsl_matrix * A, gsl_vector * S);
 587
 588 int gsl_linalg_cholesky_scale_apply(gsl_matrix * A, const gsl_vector * S);
 589
 590 int gsl_linalg_cholesky_decomp2(gsl_matrix * A, gsl_vector * S);
 591
 592 int gsl_linalg_cholesky_svx2 (const gsl_matrix * LLT,
 593                               const gsl_vector * S,
 594                               gsl_vector * x);
 595
 596 int gsl_linalg_cholesky_solve2 (const gsl_matrix * LLT,
 597                                 const gsl_vector * S,
 598                                 const gsl_vector * b,
 599                                 gsl_vector * x);
 600
 601 int gsl_linalg_cholesky_rcond (const gsl_matrix * LLT, double * rcond,
 602                                gsl_vector * work);
 603
 604 /* Complex Cholesky Decomposition */
 605
 606 int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A);
 607
 608 int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
 609                                        const gsl_vector_complex * b,
 610                                        gsl_vector_complex * x);
 611
 612 int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
 613                                      gsl_vector_complex * x);
 614
 615 int gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * cholesky);
 616
 617 /* Pivoted Cholesky LDLT decomposition */
 618
 619 int gsl_linalg_pcholesky_decomp (gsl_matrix * A, gsl_permutation * p);
 620
 621 int gsl_linalg_pcholesky_solve(const gsl_matrix * LDLT,
 622                                const gsl_permutation * p,
 623                                const gsl_vector * b,
 624                                gsl_vector * x);
 625
 626 int gsl_linalg_pcholesky_svx(const gsl_matrix * LDLT,
 627                              const gsl_permutation * p,
 628                              gsl_vector * x);
 629
 630 int gsl_linalg_pcholesky_decomp2(gsl_matrix * A, gsl_permutation * p,
 631                                  gsl_vector * S);
 632
 633 int gsl_linalg_pcholesky_solve2(const gsl_matrix * LDLT,
 634                                 const gsl_permutation * p,
 635                                 const gsl_vector * S,
 636                                 const gsl_vector * b,
 637                                 gsl_vector * x);
 638
 639 int gsl_linalg_pcholesky_svx2(const gsl_matrix * LDLT,
 640                               const gsl_permutation * p,
 641                               const gsl_vector * S,
 642                               gsl_vector * x);
 643
 644 int gsl_linalg_pcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
 645                                 gsl_matrix * Ainv);
 646
 647 int gsl_linalg_pcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
 648                                 double * rcond, gsl_vector * work);
 649
 650 /* Modified Cholesky decomposition */
 651
 652 int gsl_linalg_mcholesky_decomp (gsl_matrix * A, gsl_permutation * p,
 653                                  gsl_vector * E);
 654
 655 int gsl_linalg_mcholesky_solve(const gsl_matrix * LDLT,
 656                                const gsl_permutation * p,
 657                                const gsl_vector * b,
 658                                gsl_vector * x);
 659
 660 int gsl_linalg_mcholesky_svx(const gsl_matrix * LDLT,
 661                              const gsl_permutation * p,
 662                              gsl_vector * x);
 663
 664 int gsl_linalg_mcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
 665                                 double * rcond, gsl_vector * work);
 666
 667 int gsl_linalg_mcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
 668                                 gsl_matrix * Ainv);
 669
 670 /* Banded Cholesky decomposition */
 671
 672 int gsl_linalg_cholesky_band_decomp(gsl_matrix * A);
 673
 674 int gsl_linalg_cholesky_band_solve (const gsl_matrix * LLT, const gsl_vector * b,
 675                                     gsl_vector * x);
 676
 677 int gsl_linalg_cholesky_band_svx (const gsl_matrix * LLT, gsl_vector * x);
 678
 679 int gsl_linalg_cholesky_band_solvem (const gsl_matrix * LLT, const gsl_matrix * B,
 680                                      gsl_matrix * X);
 681
 682 int gsl_linalg_cholesky_band_svxm (const gsl_matrix * LLT, gsl_matrix * X);
 683
 684 int gsl_linalg_cholesky_band_invert (const gsl_matrix * LLT, gsl_matrix * Ainv);
 685
 686 int gsl_linalg_cholesky_band_unpack (const gsl_matrix * LLT, gsl_matrix * L);
 687
 688 int gsl_linalg_cholesky_band_scale(const gsl_matrix * A, gsl_vector * S);
 689
 690 int gsl_linalg_cholesky_band_scale_apply(gsl_matrix * A, const gsl_vector * S);
 691
 692 int gsl_linalg_cholesky_band_rcond (const gsl_matrix * LLT, double * rcond, gsl_vector * work);
 693
 694 /* L D L^T decomposition */
 695
 696 int gsl_linalg_ldlt_decomp (gsl_matrix * A);
 697
 698 int gsl_linalg_ldlt_solve (const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x);
 699
 700 int gsl_linalg_ldlt_svx (const gsl_matrix * LDLT, gsl_vector * x);
 701
 702 int gsl_linalg_ldlt_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work);
 703
 704 /* Banded L D L^T decomposition */
 705
 706 int gsl_linalg_ldlt_band_decomp (gsl_matrix * A);
 707
 708 int gsl_linalg_ldlt_band_solve (const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x);
 709
 710 int gsl_linalg_ldlt_band_svx (const gsl_matrix * LDLT, gsl_vector * x);
 711
 712 int gsl_linalg_ldlt_band_unpack (const gsl_matrix * LDLT, gsl_matrix * L, gsl_vector * D);
 713
 714 int gsl_linalg_ldlt_band_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work);
 715
 716 /* Symmetric to symmetric tridiagonal decomposition */
 717
 718 int gsl_linalg_symmtd_decomp (gsl_matrix * A,
 719                               gsl_vector * tau);
 720
 721 int gsl_linalg_symmtd_unpack (const gsl_matrix * A,
 722                               const gsl_vector * tau,
 723                               gsl_matrix * Q,
 724                               gsl_vector * diag,
 725                               gsl_vector * subdiag);
 726
 727 int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
 728                                 gsl_vector * diag,
 729                                 gsl_vector * subdiag);
 730
 731 /* Hermitian to symmetric tridiagonal decomposition */
 732
 733 int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A,
 734                               gsl_vector_complex * tau);
 735
 736 int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
 737                               const gsl_vector_complex * tau,
 738                               gsl_matrix_complex * U,
 739                               gsl_vector * diag,
 740                               gsl_vector * sudiag);
 741
 742 int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
 743                                 gsl_vector * diag,
 744                                 gsl_vector * subdiag);
 745
 746 /* Linear Solve Using Householder Transformations
 747
 748  * exceptions:
 749  */
 750
 751 int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x);
 752 int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x);
 753
 754 /* Linear solve for a symmetric tridiagonal system.
 755
 756  * The input vectors represent the NxN matrix as follows:
 757  *
 758  *     diag[0]  offdiag[0]             0    ...
 759  *  offdiag[0]     diag[1]    offdiag[1]    ...
 760  *           0  offdiag[1]       diag[2]    ...
 761  *           0           0    offdiag[2]    ...
 762  *         ...         ...           ...    ...
 763  */
 764 int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag,
 765                                    const gsl_vector * offdiag,
 766                                    const gsl_vector * b,
 767                                    gsl_vector * x);
 768
 769 /* Linear solve for a nonsymmetric tridiagonal system.
 770
 771  * The input vectors represent the NxN matrix as follows:
 772  *
 773  *       diag[0]  abovediag[0]              0    ...
 774  *  belowdiag[0]       diag[1]   abovediag[1]    ...
 775  *             0  belowdiag[1]        diag[2]    ...
 776  *             0             0   belowdiag[2]    ...
 777  *           ...           ...            ...    ...
 778  */
 779 int gsl_linalg_solve_tridiag (const gsl_vector * diag,
 780                                    const gsl_vector * abovediag,
 781                                    const gsl_vector * belowdiag,
 782                                    const gsl_vector * b,
 783                                    gsl_vector * x);
 784
 785
 786 /* Linear solve for a symmetric cyclic tridiagonal system.
 787
 788  * The input vectors represent the NxN matrix as follows:
 789  *
 790  *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
 791  *   offdiag[0]     diag[1]    offdiag[1]   .....
 792  *            0  offdiag[1]       diag[2]   .....
 793  *            0           0    offdiag[2]   .....
 794  *          ...         ...
 795  * offdiag[N-1]         ...
 796  */
 797 int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag,
 798                                        const gsl_vector * offdiag,
 799                                        const gsl_vector * b,
 800                                        gsl_vector * x);
 801
 802 /* Linear solve for a nonsymmetric cyclic tridiagonal system.
 803
 804  * The input vectors represent the NxN matrix as follows:
 805  *
 806  *        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
 807  *   belowdiag[0]       diag[1]  abovediag[1]   .....
 808  *              0  belowdiag[1]       diag[2]
 809  *              0             0  belowdiag[2]   .....
 810  *            ...           ...
 811  * abovediag[N-1]           ...
 812  */
 813 int gsl_linalg_solve_cyc_tridiag (const gsl_vector * diag,
 814                                   const gsl_vector * abovediag,
 815                                   const gsl_vector * belowdiag,
 816                                   const gsl_vector * b,
 817                                   gsl_vector * x);
 818
 819
 820 /* Bidiagonal decomposition */
 821
 822 int gsl_linalg_bidiag_decomp (gsl_matrix * A,
 823                               gsl_vector * tau_U,
 824                               gsl_vector * tau_V);
 825
 826 int gsl_linalg_bidiag_unpack (const gsl_matrix * A,
 827                               const gsl_vector * tau_U,
 828                               gsl_matrix * U,
 829                               const gsl_vector * tau_V,
 830                               gsl_matrix * V,
 831                               gsl_vector * diag,
 832                               gsl_vector * superdiag);
 833
 834 int gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
 835                                gsl_vector * tau_U,
 836                                gsl_vector * tau_V,
 837                                gsl_matrix * V);
 838
 839 int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
 840                                 gsl_vector * diag,
 841                                 gsl_vector * superdiag);
 842
 843 /* Balancing */
 844
 845 int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D);
 846 int gsl_linalg_balance_accum (gsl_matrix * A, gsl_vector * D);
 847 int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D);
 848
 849 /* condition estimation */
 850
 851 int gsl_linalg_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work);
 852 int gsl_linalg_tri_upper_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work);
 853 int gsl_linalg_tri_lower_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work);
 854 int gsl_linalg_invnorm1(const size_t N,
 855                         int (* Ainvx)(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params),
 856                         void * params, double * Ainvnorm, gsl_vector * work);
 857
 858 /* triangular matrices */
 859
 860 int gsl_linalg_tri_upper_invert(gsl_matrix * T);
 861 int gsl_linalg_tri_lower_invert(gsl_matrix * T);
 862 int gsl_linalg_tri_upper_unit_invert(gsl_matrix * T);
 863 int gsl_linalg_tri_lower_unit_invert(gsl_matrix * T);
 864
 865 int gsl_linalg_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T);
 866 int gsl_linalg_complex_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T);
 867
 868 int gsl_linalg_tri_LTL(gsl_matrix * L);
 869 int gsl_linalg_tri_UL(gsl_matrix * LU);
 870 int gsl_linalg_complex_tri_LHL(gsl_matrix_complex * L);
 871 int gsl_linalg_complex_tri_UL(gsl_matrix_complex * LU);
 872
 873 INLINE_DECL void gsl_linalg_givens (const double a, const double b,
 874                                     double *c, double *s);
 875 INLINE_DECL void gsl_linalg_givens_gv (gsl_vector * v, const size_t i,
 876                                        const size_t j, const double c,
 877                                        const double s);
 878
 879 #ifdef HAVE_INLINE
 880
 881 /* Generate a Givens rotation (cos,sin) which takes v=(x,y) to (|v|,0)
 882    From Golub and Van Loan, "Matrix Computations", Section 5.1.8 */
 883 INLINE_FUN
 884 void
 885 gsl_linalg_givens (const double a, const double b, double *c, double *s)
 886 {
 887   if (b == 0)
 888     {
 889       *c = 1;
 890       *s = 0;
 891     }
 892   else if (fabs (b) > fabs (a))
 893     {
 894       double t = -a / b;
 895       double s1 = 1.0 / sqrt (1 + t * t);
 896       *s = s1;
 897       *c = s1 * t;
 898     }
 899   else
 900     {
 901       double t = -b / a;
 902       double c1 = 1.0 / sqrt (1 + t * t);
 903       *c = c1;
 904       *s = c1 * t;
 905     }
 906 } /* gsl_linalg_givens() */
 907
 908 INLINE_FUN
 909 void
 910 gsl_linalg_givens_gv (gsl_vector * v, const size_t i, const size_t j,
 911                       const double c, const double s)
 912 {
 913   /* Apply rotation to vector v' = G^T v */
 914
 915   double vi = gsl_vector_get (v, i);
 916   double vj = gsl_vector_get (v, j);
 917   gsl_vector_set (v, i, c * vi - s * vj);
 918   gsl_vector_set (v, j, s * vi + c * vj);
 919 }
 920
 921 #endif /* HAVE_INLINE */
 922
 923 __END_DECLS
 924
 925 #endif /* __GSL_LINALG_H__ */