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 int gsl_linalg_SV_leverage(const gsl_matrix *U, gsl_vector *h); 171 172 173 /* LU Decomposition, Gaussian elimination with partial pivoting 174 */ 175 176 int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum); 177 178 int gsl_linalg_LU_solve (const gsl_matrix * LU, 179 const gsl_permutation * p, 180 const gsl_vector * b, 181 gsl_vector * x); 182 183 int gsl_linalg_LU_svx (const gsl_matrix * LU, 184 const gsl_permutation * p, 185 gsl_vector * x); 186 187 int gsl_linalg_LU_refine (const gsl_matrix * A, 188 const gsl_matrix * LU, 189 const gsl_permutation * p, 190 const gsl_vector * b, 191 gsl_vector * x, 192 gsl_vector * residual); 193 194 int gsl_linalg_LU_invert (const gsl_matrix * LU, 195 const gsl_permutation * p, 196 gsl_matrix * inverse); 197 198 double gsl_linalg_LU_det (gsl_matrix * LU, int signum); 199 double gsl_linalg_LU_lndet (gsl_matrix * LU); 200 int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum); 201 202 /* Complex LU Decomposition */ 203 204 int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, 205 gsl_permutation * p, 206 int *signum); 207 208 int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, 209 const gsl_permutation * p, 210 const gsl_vector_complex * b, 211 gsl_vector_complex * x); 212 213 int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, 214 const gsl_permutation * p, 215 gsl_vector_complex * x); 216 217 int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, 218 const gsl_matrix_complex * LU, 219 const gsl_permutation * p, 220 const gsl_vector_complex * b, 221 gsl_vector_complex * x, 222 gsl_vector_complex * residual); 223 224 int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, 225 const gsl_permutation * p, 226 gsl_matrix_complex * inverse); 227 228 gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, 229 int signum); 230 231 double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU); 232 233 gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, 234 int signum); 235 236 /* QR decomposition */ 237 238 int gsl_linalg_QR_decomp (gsl_matrix * A, 239 gsl_vector * tau); 240 241 int gsl_linalg_QR_solve (const gsl_matrix * QR, 242 const gsl_vector * tau, 243 const gsl_vector * b, 244 gsl_vector * x); 245 246 int gsl_linalg_QR_svx (const gsl_matrix * QR, 247 const gsl_vector * tau, 248 gsl_vector * x); 249 250 int gsl_linalg_QR_lssolve (const gsl_matrix * QR, 251 const gsl_vector * tau, 252 const gsl_vector * b, 253 gsl_vector * x, 254 gsl_vector * residual); 255 256 257 int gsl_linalg_QR_QRsolve (gsl_matrix * Q, 258 gsl_matrix * R, 259 const gsl_vector * b, 260 gsl_vector * x); 261 262 int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, 263 const gsl_vector * b, 264 gsl_vector * x); 265 266 int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, 267 gsl_vector * x); 268 269 int gsl_linalg_QR_update (gsl_matrix * Q, 270 gsl_matrix * R, 271 gsl_vector * w, 272 const gsl_vector * v); 273 274 int gsl_linalg_QR_QTvec (const gsl_matrix * QR, 275 const gsl_vector * tau, 276 gsl_vector * v); 277 278 int gsl_linalg_QR_Qvec (const gsl_matrix * QR, 279 const gsl_vector * tau, 280 gsl_vector * v); 281 282 int gsl_linalg_QR_QTmat (const gsl_matrix * QR, 283 const gsl_vector * tau, 284 gsl_matrix * A); 285 286 int gsl_linalg_QR_unpack (const gsl_matrix * QR, 287 const gsl_vector * tau, 288 gsl_matrix * Q, 289 gsl_matrix * R); 290 291 int gsl_linalg_R_solve (const gsl_matrix * R, 292 const gsl_vector * b, 293 gsl_vector * x); 294 295 int gsl_linalg_R_svx (const gsl_matrix * R, 296 gsl_vector * x); 297 298 299 /* Q R P^T decomposition */ 300 301 int gsl_linalg_QRPT_decomp (gsl_matrix * A, 302 gsl_vector * tau, 303 gsl_permutation * p, 304 int *signum, 305 gsl_vector * norm); 306 307 int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A, 308 gsl_matrix * q, gsl_matrix * r, 309 gsl_vector * tau, 310 gsl_permutation * p, 311 int *signum, 312 gsl_vector * norm); 313 314 int gsl_linalg_QRPT_solve (const gsl_matrix * QR, 315 const gsl_vector * tau, 316 const gsl_permutation * p, 317 const gsl_vector * b, 318 gsl_vector * x); 319 320 321 int gsl_linalg_QRPT_svx (const gsl_matrix * QR, 322 const gsl_vector * tau, 323 const gsl_permutation * p, 324 gsl_vector * x); 325 326 int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q, 327 const gsl_matrix * R, 328 const gsl_permutation * p, 329 const gsl_vector * b, 330 gsl_vector * x); 331 332 int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR, 333 const gsl_permutation * p, 334 const gsl_vector * b, 335 gsl_vector * x); 336 337 int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR, 338 const gsl_permutation * p, 339 gsl_vector * x); 340 341 int gsl_linalg_QRPT_update (gsl_matrix * Q, 342 gsl_matrix * R, 343 const gsl_permutation * p, 344 gsl_vector * u, 345 const gsl_vector * v); 346 347 /* LQ decomposition */ 348 349 int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau); 350 351 int gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau, 352 const gsl_vector * b, gsl_vector * x); 353 354 int gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau, 355 gsl_vector * x); 356 357 int gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau, 358 const gsl_vector * b, gsl_vector * x, 359 gsl_vector * residual); 360 361 int gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b, 362 gsl_vector * x); 363 364 int gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x); 365 366 int gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b, 367 gsl_vector * x); 368 369 int gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau, 370 gsl_vector * v); 371 372 int gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau, 373 gsl_vector * v); 374 375 int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau, 376 gsl_matrix * Q, gsl_matrix * L); 377 378 int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R, 379 const gsl_vector * v, gsl_vector * w); 380 int gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L, 381 const gsl_vector * b, gsl_vector * x); 382 383 /* P^T L Q decomposition */ 384 385 int gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau, 386 gsl_permutation * p, int *signum, 387 gsl_vector * norm); 388 389 int gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q, 390 gsl_matrix * r, gsl_vector * tau, 391 gsl_permutation * p, int *signum, 392 gsl_vector * norm); 393 394 int gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR, 395 const gsl_vector * tau, 396 const gsl_permutation * p, 397 const gsl_vector * b, 398 gsl_vector * x); 399 400 int gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ, 401 const gsl_vector * tau, 402 const gsl_permutation * p, 403 gsl_vector * x); 404 405 int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L, 406 const gsl_permutation * p, 407 const gsl_vector * b, 408 gsl_vector * x); 409 410 int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ, 411 const gsl_permutation * p, 412 const gsl_vector * b, 413 gsl_vector * x); 414 415 int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ, 416 const gsl_permutation * p, 417 gsl_vector * x); 418 419 int gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L, 420 const gsl_permutation * p, 421 const gsl_vector * v, gsl_vector * w); 422 423 /* Cholesky Decomposition */ 424 425 int gsl_linalg_cholesky_decomp (gsl_matrix * A); 426 427 int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky, 428 const gsl_vector * b, 429 gsl_vector * x); 430 431 int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky, 432 gsl_vector * x); 433 434 int gsl_linalg_cholesky_invert(gsl_matrix * cholesky); 435 436 /* Cholesky decomposition with unit-diagonal triangular parts. 437 * A = L D L^T, where diag(L) = (1,1,...,1). 438 * Upon exit, A contains L and L^T as for Cholesky, and 439 * the diagonal of A is (1,1,...,1). The vector Dis set 440 * to the diagonal elements of the diagonal matrix D. 441 */ 442 int gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D); 443 444 /* Complex Cholesky Decomposition */ 445 446 int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A); 447 448 int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky, 449 const gsl_vector_complex * b, 450 gsl_vector_complex * x); 451 452 int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky, 453 gsl_vector_complex * x); 454 455 int gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * cholesky); 456 457 458 /* Symmetric to symmetric tridiagonal decomposition */ 459 460 int gsl_linalg_symmtd_decomp (gsl_matrix * A, 461 gsl_vector * tau); 462 463 int gsl_linalg_symmtd_unpack (const gsl_matrix * A, 464 const gsl_vector * tau, 465 gsl_matrix * Q, 466 gsl_vector * diag, 467 gsl_vector * subdiag); 468 469 int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A, 470 gsl_vector * diag, 471 gsl_vector * subdiag); 472 473 /* Hermitian to symmetric tridiagonal decomposition */ 474 475 int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, 476 gsl_vector_complex * tau); 477 478 int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A, 479 const gsl_vector_complex * tau, 480 gsl_matrix_complex * U, 481 gsl_vector * diag, 482 gsl_vector * sudiag); 483 484 int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A, 485 gsl_vector * diag, 486 gsl_vector * subdiag); 487 488 /* Linear Solve Using Householder Transformations 489 490 * exceptions: 491 */ 492 493 int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x); 494 int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x); 495 496 /* Linear solve for a symmetric tridiagonal system. 497 498 * The input vectors represent the NxN matrix as follows: 499 * 500 * diag[0] offdiag[0] 0 ... 501 * offdiag[0] diag[1] offdiag[1] ... 502 * 0 offdiag[1] diag[2] ... 503 * 0 0 offdiag[2] ... 504 * ... ... ... ... 505 */ 506 int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag, 507 const gsl_vector * offdiag, 508 const gsl_vector * b, 509 gsl_vector * x); 510 511 /* Linear solve for a nonsymmetric tridiagonal system. 512 513 * The input vectors represent the NxN matrix as follows: 514 * 515 * diag[0] abovediag[0] 0 ... 516 * belowdiag[0] diag[1] abovediag[1] ... 517 * 0 belowdiag[1] diag[2] ... 518 * 0 0 belowdiag[2] ... 519 * ... ... ... ... 520 */ 521 int gsl_linalg_solve_tridiag (const gsl_vector * diag, 522 const gsl_vector * abovediag, 523 const gsl_vector * belowdiag, 524 const gsl_vector * b, 525 gsl_vector * x); 526 527 528 /* Linear solve for a symmetric cyclic tridiagonal system. 529 530 * The input vectors represent the NxN matrix as follows: 531 * 532 * diag[0] offdiag[0] 0 ..... offdiag[N-1] 533 * offdiag[0] diag[1] offdiag[1] ..... 534 * 0 offdiag[1] diag[2] ..... 535 * 0 0 offdiag[2] ..... 536 * ... ... 537 * offdiag[N-1] ... 538 */ 539 int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag, 540 const gsl_vector * offdiag, 541 const gsl_vector * b, 542 gsl_vector * x); 543 544 /* Linear solve for a nonsymmetric cyclic tridiagonal system. 545 546 * The input vectors represent the NxN matrix as follows: 547 * 548 * diag[0] abovediag[0] 0 ..... belowdiag[N-1] 549 * belowdiag[0] diag[1] abovediag[1] ..... 550 * 0 belowdiag[1] diag[2] 551 * 0 0 belowdiag[2] ..... 552 * ... ... 553 * abovediag[N-1] ... 554 */ 555 int gsl_linalg_solve_cyc_tridiag (const gsl_vector * diag, 556 const gsl_vector * abovediag, 557 const gsl_vector * belowdiag, 558 const gsl_vector * b, 559 gsl_vector * x); 560 561 562 /* Bidiagonal decomposition */ 563 564 int gsl_linalg_bidiag_decomp (gsl_matrix * A, 565 gsl_vector * tau_U, 566 gsl_vector * tau_V); 567 568 int gsl_linalg_bidiag_unpack (const gsl_matrix * A, 569 const gsl_vector * tau_U, 570 gsl_matrix * U, 571 const gsl_vector * tau_V, 572 gsl_matrix * V, 573 gsl_vector * diag, 574 gsl_vector * superdiag); 575 576 int gsl_linalg_bidiag_unpack2 (gsl_matrix * A, 577 gsl_vector * tau_U, 578 gsl_vector * tau_V, 579 gsl_matrix * V); 580 581 int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A, 582 gsl_vector * diag, 583 gsl_vector * superdiag); 584 585 /* Balancing */ 586 587 int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D); 588 int gsl_linalg_balance_accum (gsl_matrix * A, gsl_vector * D); 589 int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D); 590 591 592 __END_DECLS 593 594 #endif /* __GSL_LINALG_H__ */ 595