htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/gsl_linalg.h

   1 /* linalg/gsl_linalg.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 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 <gsl/gsl_mode.h>
  24 #include <gsl/gsl_permutation.h>
  25 #include <gsl/gsl_vector.h>
  26 #include <gsl/gsl_matrix.h>
  27
  28 #undef __BEGIN_DECLS
  29 #undef __END_DECLS
  30 #ifdef __cplusplus
  31 #define __BEGIN_DECLS extern "C" {
  32 #define __END_DECLS }
  33 #else
  34 #define __BEGIN_DECLS           /* empty */
  35 #define __END_DECLS             /* empty */
  36 #endif
  37
  38 __BEGIN_DECLS
  39
  40 typedef enum
  41   {
  42     GSL_LINALG_MOD_NONE = 0,
  43     GSL_LINALG_MOD_TRANSPOSE = 1,
  44     GSL_LINALG_MOD_CONJUGATE = 2
  45   }
  46 gsl_linalg_matrix_mod_t;
  47
  48
  49 /* Note: You can now use the gsl_blas_dgemm function instead of matmult */
  50
  51 /* Simple implementation of matrix multiply.
  52  * Calculates C = A.B
  53  *
  54  * exceptions: GSL_EBADLEN
  55  */
  56 int gsl_linalg_matmult (const gsl_matrix * A,
  57                         const gsl_matrix * B,
  58                         gsl_matrix * C);
  59
  60
  61 /* Simple implementation of matrix multiply.
  62  * Allows transposition of either matrix, so it
  63  * can compute A.B or Trans(A).B or A.Trans(B) or Trans(A).Trans(B)
  64  *
  65  * exceptions: GSL_EBADLEN
  66  */
  67 int gsl_linalg_matmult_mod (const gsl_matrix * A,
  68                             gsl_linalg_matrix_mod_t modA,
  69                             const gsl_matrix * B,
  70                             gsl_linalg_matrix_mod_t modB,
  71                             gsl_matrix * C);
  72
  73 /* Calculate the matrix exponential by the scaling and
  74  * squaring method described in Moler + Van Loan,
  75  * SIAM Rev 20, 801 (1978). The mode argument allows
  76  * choosing an optimal strategy, from the table
  77  * given in the paper, for a given precision.
  78  *
  79  * exceptions: GSL_ENOTSQR, GSL_EBADLEN
  80  */
  81 int gsl_linalg_exponential_ss(
  82   const gsl_matrix * A,
  83   gsl_matrix * eA,
  84   gsl_mode_t mode
  85   );
  86
  87
  88 /* Householder Transformations */
  89
  90 double gsl_linalg_householder_transform (gsl_vector * v);
  91 gsl_complex gsl_linalg_complex_householder_transform (gsl_vector_complex * v);
  92
  93 int gsl_linalg_householder_hm (double tau,
  94                                const gsl_vector * v,
  95                                gsl_matrix * A);
  96
  97 int gsl_linalg_householder_mh (double tau,
  98                                const gsl_vector * v,
  99                                gsl_matrix * A);
 100
 101 int gsl_linalg_householder_hv (double tau,
 102                                const gsl_vector * v,
 103                                gsl_vector * w);
 104
 105 int gsl_linalg_householder_hm1 (double tau,
 106                                 gsl_matrix * A);
 107
 108 int gsl_linalg_complex_householder_hm (gsl_complex tau,
 109                                        const gsl_vector_complex * v,
 110                                        gsl_matrix_complex * A);
 111
 112 int gsl_linalg_complex_householder_mh (gsl_complex tau,
 113                                        const gsl_vector_complex * v,
 114                                        gsl_matrix_complex * A);
 115
 116 int gsl_linalg_complex_householder_hv (gsl_complex tau,
 117                                        const gsl_vector_complex * v,
 118                                        gsl_vector_complex * w);
 119
 120 /* Hessenberg reduction */
 121
 122 int gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau);
 123 int gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
 124                                  gsl_matrix * U);
 125 int gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
 126                                        gsl_matrix * U);
 127 int gsl_linalg_hessenberg_set_zero(gsl_matrix * H);
 128 int gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A,
 129                                     size_t top, gsl_vector *tau);
 130
 131 /* To support gsl-1.9 interface: DEPRECATED */
 132 int gsl_linalg_hessenberg(gsl_matrix *A, gsl_vector *tau);
 133
 134
 135 /* Hessenberg-Triangular reduction */
 136
 137 int gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B,
 138                               gsl_matrix * U, gsl_matrix * V,
 139                               gsl_vector * work);
 140
 141 /* Singular Value Decomposition
 142
 143  * exceptions:
 144  */
 145
 146 int
 147 gsl_linalg_SV_decomp (gsl_matrix * A,
 148                       gsl_matrix * V,
 149                       gsl_vector * S,
 150                       gsl_vector * work);
 151
 152 int
 153 gsl_linalg_SV_decomp_mod (gsl_matrix * A,
 154                           gsl_matrix * X,
 155                           gsl_matrix * V,
 156                           gsl_vector * S,
 157                           gsl_vector * work);
 158
 159 int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A,
 160                                  gsl_matrix * Q,
 161                                  gsl_vector * S);
 162
 163 int
 164 gsl_linalg_SV_solve (const gsl_matrix * U,
 165                      const gsl_matrix * Q,
 166                      const gsl_vector * S,
 167                      const gsl_vector * b,
 168                      gsl_vector * x);
 169
 170
 171 /* LU Decomposition, Gaussian elimination with partial pivoting
 172  */
 173
 174 int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum);
 175
 176 int gsl_linalg_LU_solve (const gsl_matrix * LU,
 177                          const gsl_permutation * p,
 178                          const gsl_vector * b,
 179                          gsl_vector * x);
 180
 181 int gsl_linalg_LU_svx (const gsl_matrix * LU,
 182                        const gsl_permutation * p,
 183                        gsl_vector * x);
 184
 185 int gsl_linalg_LU_refine (const gsl_matrix * A,
 186                           const gsl_matrix * LU,
 187                           const gsl_permutation * p,
 188                           const gsl_vector * b,
 189                           gsl_vector * x,
 190                           gsl_vector * residual);
 191
 192 int gsl_linalg_LU_invert (const gsl_matrix * LU,
 193                           const gsl_permutation * p,
 194                           gsl_matrix * inverse);
 195
 196 double gsl_linalg_LU_det (gsl_matrix * LU, int signum);
 197 double gsl_linalg_LU_lndet (gsl_matrix * LU);
 198 int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum);
 199
 200 /* Complex LU Decomposition */
 201
 202 int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A,
 203                                   gsl_permutation * p,
 204                                   int *signum);
 205
 206 int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU,
 207                                  const gsl_permutation * p,
 208                                  const gsl_vector_complex * b,
 209                                  gsl_vector_complex * x);
 210
 211 int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU,
 212                                const gsl_permutation * p,
 213                                gsl_vector_complex * x);
 214
 215 int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A,
 216                                   const gsl_matrix_complex * LU,
 217                                   const gsl_permutation * p,
 218                                   const gsl_vector_complex * b,
 219                                   gsl_vector_complex * x,
 220                                   gsl_vector_complex * residual);
 221
 222 int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU,
 223                                   const gsl_permutation * p,
 224                                   gsl_matrix_complex * inverse);
 225
 226 gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU,
 227                                        int signum);
 228
 229 double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU);
 230
 231 gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU,
 232                                           int signum);
 233
 234 /* QR decomposition */
 235
 236 int gsl_linalg_QR_decomp (gsl_matrix * A,
 237                           gsl_vector * tau);
 238
 239 int gsl_linalg_QR_solve (const gsl_matrix * QR,
 240                          const gsl_vector * tau,
 241                          const gsl_vector * b,
 242                          gsl_vector * x);
 243
 244 int gsl_linalg_QR_svx (const gsl_matrix * QR,
 245                        const gsl_vector * tau,
 246                        gsl_vector * x);
 247
 248 int gsl_linalg_QR_lssolve (const gsl_matrix * QR,
 249                            const gsl_vector * tau,
 250                            const gsl_vector * b,
 251                            gsl_vector * x,
 252                            gsl_vector * residual);
 253
 254
 255 int gsl_linalg_QR_QRsolve (gsl_matrix * Q,
 256                            gsl_matrix * R,
 257                            const gsl_vector * b,
 258                            gsl_vector * x);
 259
 260 int gsl_linalg_QR_Rsolve (const gsl_matrix * QR,
 261                           const gsl_vector * b,
 262                           gsl_vector * x);
 263
 264 int gsl_linalg_QR_Rsvx (const gsl_matrix * QR,
 265                         gsl_vector * x);
 266
 267 int gsl_linalg_QR_update (gsl_matrix * Q,
 268                           gsl_matrix * R,
 269                           gsl_vector * w,
 270                           const gsl_vector * v);
 271
 272 int gsl_linalg_QR_QTvec (const gsl_matrix * QR,
 273                          const gsl_vector * tau,
 274                          gsl_vector * v);
 275
 276 int gsl_linalg_QR_Qvec (const gsl_matrix * QR,
 277                         const gsl_vector * tau,
 278                         gsl_vector * v);
 279
 280 int gsl_linalg_QR_QTmat (const gsl_matrix * QR,
 281                          const gsl_vector * tau,
 282                          gsl_matrix * A);
 283
 284 int gsl_linalg_QR_unpack (const gsl_matrix * QR,
 285                           const gsl_vector * tau,
 286                           gsl_matrix * Q,
 287                           gsl_matrix * R);
 288
 289 int gsl_linalg_R_solve (const gsl_matrix * R,
 290                         const gsl_vector * b,
 291                         gsl_vector * x);
 292
 293 int gsl_linalg_R_svx (const gsl_matrix * R,
 294                       gsl_vector * x);
 295
 296
 297 /* Q R P^T decomposition */
 298
 299 int gsl_linalg_QRPT_decomp (gsl_matrix * A,
 300                             gsl_vector * tau,
 301                             gsl_permutation * p,
 302                             int *signum,
 303                             gsl_vector * norm);
 304
 305 int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
 306                              gsl_matrix * q, gsl_matrix * r,
 307                              gsl_vector * tau,
 308                              gsl_permutation * p,
 309                              int *signum,
 310                              gsl_vector * norm);
 311
 312 int gsl_linalg_QRPT_solve (const gsl_matrix * QR,
 313                            const gsl_vector * tau,
 314                            const gsl_permutation * p,
 315                            const gsl_vector * b,
 316                            gsl_vector * x);
 317
 318
 319 int gsl_linalg_QRPT_svx (const gsl_matrix * QR,
 320                          const gsl_vector * tau,
 321                          const gsl_permutation * p,
 322                          gsl_vector * x);
 323
 324 int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,
 325                              const gsl_matrix * R,
 326                              const gsl_permutation * p,
 327                              const gsl_vector * b,
 328                              gsl_vector * x);
 329
 330 int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
 331                              const gsl_permutation * p,
 332                              const gsl_vector * b,
 333                              gsl_vector * x);
 334
 335 int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
 336                            const gsl_permutation * p,
 337                            gsl_vector * x);
 338
 339 int gsl_linalg_QRPT_update (gsl_matrix * Q,
 340                             gsl_matrix * R,
 341                             const gsl_permutation * p,
 342                             gsl_vector * u,
 343                             const gsl_vector * v);
 344
 345 /* LQ decomposition */
 346
 347 int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);
 348
 349 int gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau,
 350                          const gsl_vector * b, gsl_vector * x);
 351
 352 int gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau,
 353                          gsl_vector * x);
 354
 355 int gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau,
 356                            const gsl_vector * b, gsl_vector * x,
 357                            gsl_vector * residual);
 358
 359 int gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b,
 360                           gsl_vector * x);
 361
 362 int gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x);
 363
 364 int gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b,
 365                         gsl_vector * x);
 366
 367 int gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau,
 368                         gsl_vector * v);
 369
 370 int gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau,
 371                          gsl_vector * v);
 372
 373 int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
 374                           gsl_matrix * Q, gsl_matrix * L);
 375
 376 int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
 377                           const gsl_vector * v, gsl_vector * w);
 378 int gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L,
 379                            const gsl_vector * b, gsl_vector * x);
 380
 381 /* P^T L Q decomposition */
 382
 383 int gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau,
 384                             gsl_permutation * p, int *signum,
 385                             gsl_vector * norm);
 386
 387 int gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q,
 388                              gsl_matrix * r, gsl_vector * tau,
 389                              gsl_permutation * p, int *signum,
 390                              gsl_vector * norm);
 391
 392 int gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
 393                            const gsl_vector * tau,
 394                            const gsl_permutation * p,
 395                            const gsl_vector * b,
 396                            gsl_vector * x);
 397
 398 int gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
 399                            const gsl_vector * tau,
 400                            const gsl_permutation * p,
 401                            gsl_vector * x);
 402
 403 int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
 404                              const gsl_permutation * p,
 405                              const gsl_vector * b,
 406                              gsl_vector * x);
 407
 408 int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
 409                             const gsl_permutation * p,
 410                             const gsl_vector * b,
 411                             gsl_vector * x);
 412
 413 int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
 414                           const gsl_permutation * p,
 415                           gsl_vector * x);
 416
 417 int gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
 418                             const gsl_permutation * p,
 419                             const gsl_vector * v, gsl_vector * w);
 420
 421 /* Cholesky Decomposition */
 422
 423 int gsl_linalg_cholesky_decomp (gsl_matrix * A);
 424
 425 int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky,
 426                                const gsl_vector * b,
 427                                gsl_vector * x);
 428
 429 int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky,
 430                              gsl_vector * x);
 431
 432
 433 /* Cholesky decomposition with unit-diagonal triangular parts.
 434  *   A = L D L^T, where diag(L) = (1,1,...,1).
 435  *   Upon exit, A contains L and L^T as for Cholesky, and
 436  *   the diagonal of A is (1,1,...,1). The vector Dis set
 437  *   to the diagonal elements of the diagonal matrix D.
 438  */
 439 int gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D);
 440
 441 /* Complex Cholesky Decomposition */
 442
 443 int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A);
 444
 445 int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
 446                                        const gsl_vector_complex * b,
 447                                        gsl_vector_complex * x);
 448
 449 int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
 450                                      gsl_vector_complex * x);
 451
 452 /* Symmetric to symmetric tridiagonal decomposition */
 453
 454 int gsl_linalg_symmtd_decomp (gsl_matrix * A,
 455                               gsl_vector * tau);
 456
 457 int gsl_linalg_symmtd_unpack (const gsl_matrix * A,
 458                               const gsl_vector * tau,
 459                               gsl_matrix * Q,
 460                               gsl_vector * diag,
 461                               gsl_vector * subdiag);
 462
 463 int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
 464                                 gsl_vector * diag,
 465                                 gsl_vector * subdiag);
 466
 467 /* Hermitian to symmetric tridiagonal decomposition */
 468
 469 int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A,
 470                               gsl_vector_complex * tau);
 471
 472 int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
 473                               const gsl_vector_complex * tau,
 474                               gsl_matrix_complex * Q,
 475                               gsl_vector * diag,
 476                               gsl_vector * sudiag);
 477
 478 int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
 479                                 gsl_vector * diag,
 480                                 gsl_vector * subdiag);
 481
 482 /* Linear Solve Using Householder Transformations
 483
 484  * exceptions:
 485  */
 486
 487 int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x);
 488 int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x);
 489
 490 /* Linear solve for a symmetric tridiagonal system.
 491
 492  * The input vectors represent the NxN matrix as follows:
 493  *
 494  *     diag[0]  offdiag[0]             0    ...
 495  *  offdiag[0]     diag[1]    offdiag[1]    ...
 496  *           0  offdiag[1]       diag[2]    ...
 497  *           0           0    offdiag[2]    ...
 498  *         ...         ...           ...    ...
 499  */
 500 int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag,
 501                                    const gsl_vector * offdiag,
 502                                    const gsl_vector * b,
 503                                    gsl_vector * x);
 504
 505 /* Linear solve for a nonsymmetric tridiagonal system.
 506
 507  * The input vectors represent the NxN matrix as follows:
 508  *
 509  *       diag[0]  abovediag[0]              0    ...
 510  *  belowdiag[0]       diag[1]   abovediag[1]    ...
 511  *             0  belowdiag[1]        diag[2]    ...
 512  *             0             0   belowdiag[2]    ...
 513  *           ...           ...            ...    ...
 514  */
 515 int gsl_linalg_solve_tridiag (const gsl_vector * diag,
 516                                    const gsl_vector * abovediag,
 517                                    const gsl_vector * belowdiag,
 518                                    const gsl_vector * b,
 519                                    gsl_vector * x);
 520
 521
 522 /* Linear solve for a symmetric cyclic tridiagonal system.
 523
 524  * The input vectors represent the NxN matrix as follows:
 525  *
 526  *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
 527  *   offdiag[0]     diag[1]    offdiag[1]   .....
 528  *            0  offdiag[1]       diag[2]   .....
 529  *            0           0    offdiag[2]   .....
 530  *          ...         ...
 531  * offdiag[N-1]         ...
 532  */
 533 int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag,
 534                                        const gsl_vector * offdiag,
 535                                        const gsl_vector * b,
 536                                        gsl_vector * x);
 537
 538 /* Linear solve for a nonsymmetric cyclic tridiagonal system.
 539
 540  * The input vectors represent the NxN matrix as follows:
 541  *
 542  *        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
 543  *   belowdiag[0]       diag[1]  abovediag[1]   .....
 544  *              0  belowdiag[1]       diag[2]
 545  *              0             0  belowdiag[2]   .....
 546  *            ...           ...
 547  * abovediag[N-1]           ...
 548  */
 549 int gsl_linalg_solve_cyc_tridiag (const gsl_vector * diag,
 550                                   const gsl_vector * abovediag,
 551                                   const gsl_vector * belowdiag,
 552                                   const gsl_vector * b,
 553                                   gsl_vector * x);
 554
 555
 556 /* Bidiagonal decomposition */
 557
 558 int gsl_linalg_bidiag_decomp (gsl_matrix * A,
 559                               gsl_vector * tau_U,
 560                               gsl_vector * tau_V);
 561
 562 int gsl_linalg_bidiag_unpack (const gsl_matrix * A,
 563                               const gsl_vector * tau_U,
 564                               gsl_matrix * U,
 565                               const gsl_vector * tau_V,
 566                               gsl_matrix * V,
 567                               gsl_vector * diag,
 568                               gsl_vector * superdiag);
 569
 570 int gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
 571                                gsl_vector * tau_U,
 572                                gsl_vector * tau_V,
 573                                gsl_matrix * V);
 574
 575 int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
 576                                 gsl_vector * diag,
 577                                 gsl_vector * superdiag);
 578
 579 /* Balancing */
 580
 581 int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D);
 582 int gsl_linalg_balance_accum (gsl_matrix * A, gsl_vector * D);
 583 int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D);
 584
 585
 586 __END_DECLS
 587
 588 #endif /* __GSL_LINALG_H__ */