include/gsl/gsl_linalg.h

0001 /* linalg/gsl_linalg.h
0002  *
0003  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007, 2019 Gerard Jungman, Brian Gough, Patrick Alken
0004  *
0005  * This program is free software; you can redistribute it and/or modify
0006  * it under the terms of the GNU General Public License as published by
0007  * the Free Software Foundation; either version 3 of the License, or (at
0008  * your option) any later version.
0009  *
0010  * This program is distributed in the hope that it will be useful, but
0011  * WITHOUT ANY WARRANTY; without even the implied warranty of
0012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0013  * General Public License for more details.
0014  *
0015  * You should have received a copy of the GNU General Public License
0016  * along with this program; if not, write to the Free Software
0017  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0018  */
0019
0020 #ifndef __GSL_LINALG_H__
0021 #define __GSL_LINALG_H__
0022
0023 #include <stdlib.h>
0024 #include <gsl/gsl_mode.h>
0025 #include <gsl/gsl_permutation.h>
0026 #include <gsl/gsl_vector.h>
0027 #include <gsl/gsl_matrix.h>
0028 #include <gsl/gsl_math.h>
0029 #include <gsl/gsl_inline.h>
0030 #include <gsl/gsl_blas.h>
0031
0032 #undef __BEGIN_DECLS
0033 #undef __END_DECLS
0034 #ifdef __cplusplus
0035 #define __BEGIN_DECLS extern "C" {
0036 #define __END_DECLS }
0037 #else
0038 #define __BEGIN_DECLS           /* empty */
0039 #define __END_DECLS             /* empty */
0040 #endif
0041
0042 __BEGIN_DECLS
0043
0044 typedef enum
0045   {
0046     GSL_LINALG_MOD_NONE = 0,
0047     GSL_LINALG_MOD_TRANSPOSE = 1,
0048     GSL_LINALG_MOD_CONJUGATE = 2
0049   }
0050 gsl_linalg_matrix_mod_t;
0051
0052 /* Note: You can now use the gsl_blas_dgemm function instead of matmult */
0053
0054 /* Simple implementation of matrix multiply.
0055  * Calculates C = A.B
0056  *
0057  * exceptions: GSL_EBADLEN
0058  */
0059 int gsl_linalg_matmult (const gsl_matrix * A,
0060                         const gsl_matrix * B,
0061                         gsl_matrix * C);
0062
0063
0064 /* Simple implementation of matrix multiply.
0065  * Allows transposition of either matrix, so it
0066  * can compute A.B or Trans(A).B or A.Trans(B) or Trans(A).Trans(B)
0067  *
0068  * exceptions: GSL_EBADLEN
0069  */
0070 int gsl_linalg_matmult_mod (const gsl_matrix * A,
0071                             gsl_linalg_matrix_mod_t modA,
0072                             const gsl_matrix * B,
0073                             gsl_linalg_matrix_mod_t modB,
0074                             gsl_matrix * C);
0075
0076 /* Calculate the matrix exponential by the scaling and
0077  * squaring method described in Moler + Van Loan,
0078  * SIAM Rev 20, 801 (1978). The mode argument allows
0079  * choosing an optimal strategy, from the table
0080  * given in the paper, for a given precision.
0081  *
0082  * exceptions: GSL_ENOTSQR, GSL_EBADLEN
0083  */
0084 int gsl_linalg_exponential_ss(
0085   const gsl_matrix * A,
0086   gsl_matrix * eA,
0087   gsl_mode_t mode
0088   );
0089
0090
0091 /* Householder Transformations */
0092
0093 double gsl_linalg_householder_transform (gsl_vector * v);
0094 double gsl_linalg_householder_transform2 (double * alpha, gsl_vector * v);
0095 gsl_complex gsl_linalg_complex_householder_transform (gsl_vector_complex * v);
0096
0097 int gsl_linalg_householder_hm (double tau,
0098                                const gsl_vector * v,
0099                                gsl_matrix * A);
0100
0101 int gsl_linalg_householder_mh (double tau,
0102                                const gsl_vector * v,
0103                                gsl_matrix * A);
0104
0105 int gsl_linalg_householder_hv (double tau,
0106                                const gsl_vector * v,
0107                                gsl_vector * w);
0108
0109 int gsl_linalg_householder_left(const double tau,
0110                                 const gsl_vector * v,
0111                                 gsl_matrix * A,
0112                                 gsl_vector * work);
0113
0114 int gsl_linalg_householder_right(const double tau,
0115                                  const gsl_vector * v,
0116                                  gsl_matrix * A,
0117                                  gsl_vector * work);
0118
0119 int gsl_linalg_householder_hm1 (double tau,
0120                                 gsl_matrix * A);
0121
0122 int gsl_linalg_complex_householder_hm (gsl_complex tau,
0123                                        const gsl_vector_complex * v,
0124                                        gsl_matrix_complex * A);
0125
0126 int gsl_linalg_complex_householder_mh (gsl_complex tau,
0127                                        const gsl_vector_complex * v,
0128                                        gsl_matrix_complex * A);
0129
0130 int gsl_linalg_complex_householder_hv (gsl_complex tau,
0131                                        const gsl_vector_complex * v,
0132                                        gsl_vector_complex * w);
0133
0134 int gsl_linalg_complex_householder_left (const gsl_complex tau,
0135                                          const gsl_vector_complex * v,
0136                                          gsl_matrix_complex * A,
0137                                          gsl_vector_complex * work);
0138
0139 /* Hessenberg reduction */
0140
0141 int gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau);
0142 int gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
0143                                  gsl_matrix * U);
0144 int gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
0145                                        gsl_matrix * U);
0146 int gsl_linalg_hessenberg_set_zero(gsl_matrix * H);
0147 int gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A,
0148                                     size_t top, gsl_vector *tau);
0149
0150 /* Hessenberg-Triangular reduction */
0151
0152 int gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B,
0153                               gsl_matrix * U, gsl_matrix * V,
0154                               gsl_vector * work);
0155
0156 /* Singular Value Decomposition
0157
0158  * exceptions:
0159  */
0160
0161 int
0162 gsl_linalg_SV_decomp (gsl_matrix * A,
0163                       gsl_matrix * V,
0164                       gsl_vector * S,
0165                       gsl_vector * work);
0166
0167 int
0168 gsl_linalg_SV_decomp_mod (gsl_matrix * A,
0169                           gsl_matrix * X,
0170                           gsl_matrix * V,
0171                           gsl_vector * S,
0172                           gsl_vector * work);
0173
0174 int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A,
0175                                  gsl_matrix * Q,
0176                                  gsl_vector * S);
0177
0178 int
0179 gsl_linalg_SV_solve (const gsl_matrix * U,
0180                      const gsl_matrix * Q,
0181                      const gsl_vector * S,
0182                      const gsl_vector * b,
0183                      gsl_vector * x);
0184
0185 int
0186 gsl_linalg_SV_solve2 (const double tol,
0187                       const gsl_matrix * U,
0188                       const gsl_matrix * V,
0189                       const gsl_vector * S,
0190                       const gsl_vector * b,
0191                       gsl_vector * x,
0192                       gsl_vector * work);
0193
0194 int
0195 gsl_linalg_SV_lssolve (const double lambda,
0196                        const gsl_matrix * U,
0197                        const gsl_matrix * V,
0198                        const gsl_vector * S,
0199                        const gsl_vector * b,
0200                        gsl_vector * x,
0201                        double * rnorm,
0202                        gsl_vector * work);
0203
0204 int gsl_linalg_SV_leverage(const gsl_matrix *U, gsl_vector *h);
0205
0206
0207 /* LU Decomposition, Gaussian elimination with partial pivoting
0208  */
0209
0210 int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum);
0211
0212 int gsl_linalg_LU_solve (const gsl_matrix * LU,
0213                          const gsl_permutation * p,
0214                          const gsl_vector * b,
0215                          gsl_vector * x);
0216
0217 int gsl_linalg_LU_svx (const gsl_matrix * LU,
0218                        const gsl_permutation * p,
0219                        gsl_vector * x);
0220
0221 int gsl_linalg_LU_refine (const gsl_matrix * A,
0222                           const gsl_matrix * LU,
0223                           const gsl_permutation * p,
0224                           const gsl_vector * b,
0225                           gsl_vector * x,
0226                           gsl_vector * work);
0227
0228 int gsl_linalg_LU_invert (const gsl_matrix * LU,
0229                           const gsl_permutation * p,
0230                           gsl_matrix * inverse);
0231 int gsl_linalg_LU_invx (gsl_matrix * LU, const gsl_permutation * p);
0232
0233 double gsl_linalg_LU_det (gsl_matrix * LU, int signum);
0234 double gsl_linalg_LU_lndet (gsl_matrix * LU);
0235 int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum);
0236
0237 /* Banded LU decomposition */
0238
0239 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);
0240
0241 int gsl_linalg_LU_band_solve (const size_t lb, const size_t ub, const gsl_matrix * LUB,
0242                               const gsl_vector_uint * piv, const gsl_vector * b, gsl_vector * x);
0243
0244 int gsl_linalg_LU_band_svx (const size_t lb, const size_t ub, const gsl_matrix * LUB,
0245                             const gsl_vector_uint * piv, gsl_vector * x);
0246
0247 int gsl_linalg_LU_band_unpack (const size_t M, const size_t lb, const size_t ub, const gsl_matrix * LUB,
0248                                const gsl_vector_uint * piv, gsl_matrix * L, gsl_matrix * U);
0249
0250 /* Complex LU Decomposition */
0251
0252 int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A,
0253                                   gsl_permutation * p,
0254                                   int *signum);
0255
0256 int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU,
0257                                  const gsl_permutation * p,
0258                                  const gsl_vector_complex * b,
0259                                  gsl_vector_complex * x);
0260
0261 int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU,
0262                                const gsl_permutation * p,
0263                                gsl_vector_complex * x);
0264
0265 int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A,
0266                                   const gsl_matrix_complex * LU,
0267                                   const gsl_permutation * p,
0268                                   const gsl_vector_complex * b,
0269                                   gsl_vector_complex * x,
0270                                   gsl_vector_complex * work);
0271
0272 int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU,
0273                                   const gsl_permutation * p,
0274                                   gsl_matrix_complex * inverse);
0275 int gsl_linalg_complex_LU_invx (gsl_matrix_complex * LU, const gsl_permutation * p);
0276
0277 gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU,
0278                                        int signum);
0279
0280 double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU);
0281
0282 gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU,
0283                                           int signum);
0284
0285 /* QR decomposition */
0286
0287 int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau);
0288
0289 int gsl_linalg_QR_decomp_old (gsl_matrix * A, gsl_vector * tau);
0290
0291 int gsl_linalg_QR_decomp_r (gsl_matrix * A, gsl_matrix * T);
0292
0293 int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x);
0294
0295 int gsl_linalg_QR_solve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x);
0296
0297 int gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x);
0298
0299 int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b,
0300                            gsl_vector * x, gsl_vector * residual);
0301
0302 int gsl_linalg_QR_lssolve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b,
0303                              gsl_vector * x, gsl_vector * work);
0304
0305 int gsl_linalg_QR_lssolvem_r (const gsl_matrix * QR, const gsl_matrix * T,
0306                               const gsl_matrix * B, gsl_matrix * X, gsl_matrix * work);
0307
0308 int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x);
0309
0310 int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x);
0311
0312 int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x);
0313
0314 int gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R, gsl_vector * w, const gsl_vector * v);
0315
0316 int gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);
0317
0318 int gsl_linalg_QR_QTvec_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
0319
0320 int gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);
0321
0322 int gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);
0323
0324 int gsl_linalg_QR_QTmat_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * B, gsl_matrix * work);
0325
0326 int gsl_linalg_QR_matQ (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);
0327
0328 int gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R);
0329
0330 int gsl_linalg_QR_unpack_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q, gsl_matrix * R);
0331
0332 int gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x);
0333
0334 int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x);
0335
0336 int gsl_linalg_QR_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work);
0337
0338 /* complex QR decomposition */
0339
0340 int gsl_linalg_complex_QR_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau);
0341
0342 int gsl_linalg_complex_QR_decomp_r (gsl_matrix_complex * A, gsl_matrix_complex * T);
0343
0344 int gsl_linalg_complex_QR_solve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
0345                                  const gsl_vector_complex * b, gsl_vector_complex * x);
0346
0347 int gsl_linalg_complex_QR_solve_r (const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
0348                                    const gsl_vector_complex * b, gsl_vector_complex * x);
0349
0350 int gsl_linalg_complex_QR_svx (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * x);
0351
0352 int gsl_linalg_complex_QR_lssolve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
0353                                    const gsl_vector_complex * b, gsl_vector_complex * x,
0354                                    gsl_vector_complex * residual);
0355
0356 int gsl_linalg_complex_QR_lssolve_r (const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
0357                                      const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * work);
0358
0359 int gsl_linalg_complex_QR_lssolvem_r (const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
0360                                       const gsl_matrix_complex * B, gsl_matrix_complex * X, gsl_matrix_complex * work);
0361
0362 int gsl_linalg_complex_QR_QHvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v);
0363
0364 int gsl_linalg_complex_QR_QHvec_r(const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
0365                                   gsl_vector_complex * b, gsl_vector_complex * work);
0366
0367 int gsl_linalg_complex_QR_QHmat_r(const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
0368                                   gsl_matrix_complex * B, gsl_matrix_complex * work);
0369
0370 int gsl_linalg_complex_QR_Qvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v);
0371
0372 int gsl_linalg_complex_QR_unpack (const gsl_matrix_complex * QR, const gsl_vector_complex * tau,
0373                                   gsl_matrix_complex * Q, gsl_matrix_complex * R);
0374
0375 int gsl_linalg_complex_QR_unpack_r(const gsl_matrix_complex * QR, const gsl_matrix_complex * T,
0376                                    gsl_matrix_complex * Q, gsl_matrix_complex * R);
0377
0378 /* banded QR decomposition */
0379
0380 int gsl_linalg_QR_band_decomp_L2 (const size_t M, const size_t p, const size_t q,
0381                                   gsl_matrix * AB, gsl_vector * tau);
0382
0383 int gsl_linalg_QR_band_unpack_L2 (const size_t p, const size_t q, const gsl_matrix * QRB,
0384                                   const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R);
0385
0386 /* Q R P^T decomposition */
0387
0388 int gsl_linalg_QRPT_decomp (gsl_matrix * A,
0389                             gsl_vector * tau,
0390                             gsl_permutation * p,
0391                             int *signum,
0392                             gsl_vector * norm);
0393
0394 int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
0395                              gsl_matrix * q, gsl_matrix * r,
0396                              gsl_vector * tau,
0397                              gsl_permutation * p,
0398                              int *signum,
0399                              gsl_vector * norm);
0400
0401 int gsl_linalg_QRPT_solve (const gsl_matrix * QR,
0402                            const gsl_vector * tau,
0403                            const gsl_permutation * p,
0404                            const gsl_vector * b,
0405                            gsl_vector * x);
0406
0407 int gsl_linalg_QRPT_lssolve (const gsl_matrix * QR,
0408                              const gsl_vector * tau,
0409                              const gsl_permutation * p,
0410                              const gsl_vector * b,
0411                              gsl_vector * x,
0412                              gsl_vector * residual);
0413
0414 int gsl_linalg_QRPT_lssolve2 (const gsl_matrix * QR,
0415                               const gsl_vector * tau,
0416                               const gsl_permutation * p,
0417                               const gsl_vector * b,
0418                               const size_t rank,
0419                               gsl_vector * x,
0420                               gsl_vector * residual);
0421
0422 int gsl_linalg_QRPT_svx (const gsl_matrix * QR,
0423                          const gsl_vector * tau,
0424                          const gsl_permutation * p,
0425                          gsl_vector * x);
0426
0427 int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,
0428                              const gsl_matrix * R,
0429                              const gsl_permutation * p,
0430                              const gsl_vector * b,
0431                              gsl_vector * x);
0432
0433 int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
0434                              const gsl_permutation * p,
0435                              const gsl_vector * b,
0436                              gsl_vector * x);
0437
0438 int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
0439                            const gsl_permutation * p,
0440                            gsl_vector * x);
0441
0442 int gsl_linalg_QRPT_update (gsl_matrix * Q,
0443                             gsl_matrix * R,
0444                             const gsl_permutation * p,
0445                             gsl_vector * u,
0446                             const gsl_vector * v);
0447
0448 size_t gsl_linalg_QRPT_rank (const gsl_matrix * QR, const double tol);
0449
0450 int gsl_linalg_QRPT_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work);
0451
0452 /* triangle on top of diagonal QR decomposition */
0453
0454 int gsl_linalg_QR_UD_decomp (gsl_matrix * U, const gsl_vector * D, gsl_matrix * Y, gsl_matrix * T);
0455
0456 int gsl_linalg_QR_UD_lssolve (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
0457                               const gsl_vector * b, gsl_vector * x, gsl_vector * work);
0458
0459 int gsl_linalg_QR_UD_lssvx (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
0460                             gsl_vector * x, gsl_vector * work);
0461
0462 int gsl_linalg_QR_UD_QTvec(const gsl_matrix * Y, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
0463
0464 /* triangle on top of rectangle QR decomposition */
0465
0466 int gsl_linalg_QR_UR_decomp (gsl_matrix * S, gsl_matrix * A, gsl_matrix * T);
0467
0468 int gsl_linalg_QR_UR_lssolve (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
0469                               const gsl_vector * b, gsl_vector * x, gsl_vector * work);
0470
0471 int gsl_linalg_QR_UR_lssvx (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
0472                             gsl_vector * x, gsl_vector * work);
0473
0474 int gsl_linalg_QR_UR_QTvec(const gsl_matrix * Y, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
0475
0476 /* triangle on top of triangle QR decomposition */
0477
0478 int gsl_linalg_QR_UU_decomp (gsl_matrix * U, gsl_matrix * S, gsl_matrix * T);
0479
0480 int gsl_linalg_QR_UU_lssolve (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
0481                               const gsl_vector * b, gsl_vector * x, gsl_vector * work);
0482
0483 int gsl_linalg_QR_UU_lssvx (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
0484                             gsl_vector * x, gsl_vector * work);
0485
0486 int gsl_linalg_QR_UU_QTvec(const gsl_matrix * Y, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
0487
0488 /* triangle on top of trapezoidal QR decomposition */
0489
0490 int gsl_linalg_QR_UZ_decomp (gsl_matrix * S, gsl_matrix * A, gsl_matrix * T);
0491
0492 /* QL decomposition */
0493
0494 int gsl_linalg_QL_decomp (gsl_matrix * A, gsl_vector * tau);
0495
0496 int gsl_linalg_QL_unpack (const gsl_matrix * QL, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L);
0497
0498 /* COD decomposition */
0499
0500 int gsl_linalg_COD_decomp(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
0501                           gsl_permutation * p, size_t * rank, gsl_vector * work);
0502
0503 int gsl_linalg_COD_decomp_e(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
0504                             gsl_permutation * p, double tol, size_t * rank, gsl_vector * work);
0505
0506 int gsl_linalg_COD_lssolve (const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
0507                             const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
0508                             gsl_vector * x, gsl_vector * residual);
0509
0510 int
0511 gsl_linalg_COD_lssolve2 (const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
0512                          const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
0513                          gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work);
0514
0515 int gsl_linalg_COD_unpack(const gsl_matrix * QRZT, const gsl_vector * tau_Q,
0516                           const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q,
0517                           gsl_matrix * R, gsl_matrix * Z);
0518
0519 int gsl_linalg_COD_matZ(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
0520                         gsl_matrix * A, gsl_vector * work);
0521
0522 /* LQ decomposition */
0523
0524 int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);
0525
0526 int gsl_linalg_LQ_lssolve(const gsl_matrix * LQ, const gsl_vector * tau,
0527                           const gsl_vector * b, gsl_vector * x, gsl_vector * residual);
0528
0529 int gsl_linalg_LQ_QTvec(const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v);
0530
0531 int gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau,
0532              const gsl_vector * b, gsl_vector * x);
0533
0534 int gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau,
0535                          gsl_vector * x);
0536
0537 int gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau,
0538                const gsl_vector * b, gsl_vector * x,
0539                gsl_vector * residual);
0540
0541 int gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b,
0542               gsl_vector * x);
0543
0544 int gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x);
0545
0546 int gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b,
0547             gsl_vector * x);
0548
0549 int gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau,
0550             gsl_vector * v);
0551
0552 int gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau,
0553              gsl_vector * v);
0554
0555 int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
0556               gsl_matrix * Q, gsl_matrix * L);
0557
0558 int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
0559               const gsl_vector * v, gsl_vector * w);
0560 int gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L,
0561                const gsl_vector * b, gsl_vector * x);
0562
0563 /* P^T L Q decomposition */
0564
0565 int gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau,
0566                 gsl_permutation * p, int *signum,
0567                 gsl_vector * norm);
0568
0569 int gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q,
0570                  gsl_matrix * r, gsl_vector * tau,
0571                  gsl_permutation * p, int *signum,
0572                  gsl_vector * norm);
0573
0574 int gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
0575                const gsl_vector * tau,
0576                const gsl_permutation * p,
0577                const gsl_vector * b,
0578                gsl_vector * x);
0579
0580 int gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
0581                            const gsl_vector * tau,
0582                            const gsl_permutation * p,
0583                            gsl_vector * x);
0584
0585 int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
0586                  const gsl_permutation * p,
0587                  const gsl_vector * b,
0588                  gsl_vector * x);
0589
0590 int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
0591                 const gsl_permutation * p,
0592                 const gsl_vector * b,
0593                 gsl_vector * x);
0594
0595 int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
0596               const gsl_permutation * p,
0597               gsl_vector * x);
0598
0599 int gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
0600                 const gsl_permutation * p,
0601                 const gsl_vector * v, gsl_vector * w);
0602
0603 /* Cholesky Decomposition */
0604
0605 int gsl_linalg_cholesky_decomp (gsl_matrix * A);
0606 int gsl_linalg_cholesky_decomp1 (gsl_matrix * A);
0607
0608 int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky,
0609                                const gsl_vector * b,
0610                                gsl_vector * x);
0611 int gsl_linalg_cholesky_solve_mat (const gsl_matrix * cholesky,
0612                                    const gsl_matrix * B,
0613                                    gsl_matrix * X);
0614
0615 int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky,
0616                              gsl_vector * x);
0617 int gsl_linalg_cholesky_svx_mat (const gsl_matrix * cholesky,
0618                                  gsl_matrix * X);
0619
0620 int gsl_linalg_cholesky_invert(gsl_matrix * cholesky);
0621
0622 /* Cholesky decomposition with unit-diagonal triangular parts.
0623  *   A = L D L^T, where diag(L) = (1,1,...,1).
0624  *   Upon exit, A contains L and L^T as for Cholesky, and
0625  *   the diagonal of A is (1,1,...,1). The vector Dis set
0626  *   to the diagonal elements of the diagonal matrix D.
0627  */
0628 int gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D);
0629
0630 int gsl_linalg_cholesky_scale(const gsl_matrix * A, gsl_vector * S);
0631
0632 int gsl_linalg_cholesky_scale_apply(gsl_matrix * A, const gsl_vector * S);
0633
0634 int gsl_linalg_cholesky_decomp2(gsl_matrix * A, gsl_vector * S);
0635
0636 int gsl_linalg_cholesky_svx2 (const gsl_matrix * LLT,
0637                               const gsl_vector * S,
0638                               gsl_vector * x);
0639
0640 int gsl_linalg_cholesky_solve2 (const gsl_matrix * LLT,
0641                                 const gsl_vector * S,
0642                                 const gsl_vector * b,
0643                                 gsl_vector * x);
0644
0645 int gsl_linalg_cholesky_rcond (const gsl_matrix * LLT, double * rcond,
0646                                gsl_vector * work);
0647
0648 /* Complex Cholesky Decomposition */
0649
0650 int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A);
0651
0652 int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
0653                                        const gsl_vector_complex * b,
0654                                        gsl_vector_complex * x);
0655
0656 int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
0657                                      gsl_vector_complex * x);
0658
0659 int gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * cholesky);
0660
0661 int gsl_linalg_complex_cholesky_scale(const gsl_matrix_complex * A, gsl_vector * S);
0662
0663 int gsl_linalg_complex_cholesky_scale_apply(gsl_matrix_complex * A, const gsl_vector * S);
0664
0665 int gsl_linalg_complex_cholesky_decomp2(gsl_matrix_complex * A, gsl_vector * S);
0666
0667 int gsl_linalg_complex_cholesky_svx2 (const gsl_matrix_complex * LLT, const gsl_vector * S, gsl_vector_complex * x);
0668
0669 int gsl_linalg_complex_cholesky_solve2 (const gsl_matrix_complex * LLT, const gsl_vector * S,
0670                                         const gsl_vector_complex * b, gsl_vector_complex * x);
0671
0672 /* Pivoted Cholesky LDLT decomposition */
0673
0674 int gsl_linalg_pcholesky_decomp (gsl_matrix * A, gsl_permutation * p);
0675
0676 int gsl_linalg_pcholesky_solve(const gsl_matrix * LDLT,
0677                                const gsl_permutation * p,
0678                                const gsl_vector * b,
0679                                gsl_vector * x);
0680
0681 int gsl_linalg_pcholesky_svx(const gsl_matrix * LDLT,
0682                              const gsl_permutation * p,
0683                              gsl_vector * x);
0684
0685 int gsl_linalg_pcholesky_decomp2(gsl_matrix * A, gsl_permutation * p,
0686                                  gsl_vector * S);
0687
0688 int gsl_linalg_pcholesky_solve2(const gsl_matrix * LDLT,
0689                                 const gsl_permutation * p,
0690                                 const gsl_vector * S,
0691                                 const gsl_vector * b,
0692                                 gsl_vector * x);
0693
0694 int gsl_linalg_pcholesky_svx2(const gsl_matrix * LDLT,
0695                               const gsl_permutation * p,
0696                               const gsl_vector * S,
0697                               gsl_vector * x);
0698
0699 int gsl_linalg_pcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
0700                                 gsl_matrix * Ainv);
0701
0702 int gsl_linalg_pcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
0703                                 double * rcond, gsl_vector * work);
0704
0705 /* Modified Cholesky decomposition */
0706
0707 int gsl_linalg_mcholesky_decomp (gsl_matrix * A, gsl_permutation * p,
0708                                  gsl_vector * E);
0709
0710 int gsl_linalg_mcholesky_solve(const gsl_matrix * LDLT,
0711                                const gsl_permutation * p,
0712                                const gsl_vector * b,
0713                                gsl_vector * x);
0714
0715 int gsl_linalg_mcholesky_svx(const gsl_matrix * LDLT,
0716                              const gsl_permutation * p,
0717                              gsl_vector * x);
0718
0719 int gsl_linalg_mcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
0720                                 double * rcond, gsl_vector * work);
0721
0722 int gsl_linalg_mcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
0723                                 gsl_matrix * Ainv);
0724
0725 /* Banded Cholesky decomposition */
0726
0727 int gsl_linalg_cholesky_band_decomp(gsl_matrix * A);
0728
0729 int gsl_linalg_cholesky_band_solve (const gsl_matrix * LLT, const gsl_vector * b,
0730                                     gsl_vector * x);
0731
0732 int gsl_linalg_cholesky_band_svx (const gsl_matrix * LLT, gsl_vector * x);
0733
0734 int gsl_linalg_cholesky_band_solvem (const gsl_matrix * LLT, const gsl_matrix * B,
0735                                      gsl_matrix * X);
0736
0737 int gsl_linalg_cholesky_band_svxm (const gsl_matrix * LLT, gsl_matrix * X);
0738
0739 int gsl_linalg_cholesky_band_invert (const gsl_matrix * LLT, gsl_matrix * Ainv);
0740
0741 int gsl_linalg_cholesky_band_unpack (const gsl_matrix * LLT, gsl_matrix * L);
0742
0743 int gsl_linalg_cholesky_band_scale(const gsl_matrix * A, gsl_vector * S);
0744
0745 int gsl_linalg_cholesky_band_scale_apply(gsl_matrix * A, const gsl_vector * S);
0746
0747 int gsl_linalg_cholesky_band_rcond (const gsl_matrix * LLT, double * rcond, gsl_vector * work);
0748
0749 /* L D L^T decomposition */
0750
0751 int gsl_linalg_ldlt_decomp (gsl_matrix * A);
0752
0753 int gsl_linalg_ldlt_solve (const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x);
0754
0755 int gsl_linalg_ldlt_svx (const gsl_matrix * LDLT, gsl_vector * x);
0756
0757 int gsl_linalg_ldlt_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work);
0758
0759 /* Banded L D L^T decomposition */
0760
0761 int gsl_linalg_ldlt_band_decomp (gsl_matrix * A);
0762
0763 int gsl_linalg_ldlt_band_solve (const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x);
0764
0765 int gsl_linalg_ldlt_band_svx (const gsl_matrix * LDLT, gsl_vector * x);
0766
0767 int gsl_linalg_ldlt_band_unpack (const gsl_matrix * LDLT, gsl_matrix * L, gsl_vector * D);
0768
0769 int gsl_linalg_ldlt_band_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work);
0770
0771 /* Symmetric to symmetric tridiagonal decomposition */
0772
0773 int gsl_linalg_symmtd_decomp (gsl_matrix * A,
0774                               gsl_vector * tau);
0775
0776 int gsl_linalg_symmtd_unpack (const gsl_matrix * A,
0777                               const gsl_vector * tau,
0778                               gsl_matrix * Q,
0779                               gsl_vector * diag,
0780                               gsl_vector * subdiag);
0781
0782 int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
0783                                 gsl_vector * diag,
0784                                 gsl_vector * subdiag);
0785
0786 /* Hermitian to symmetric tridiagonal decomposition */
0787
0788 int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A,
0789                               gsl_vector_complex * tau);
0790
0791 int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
0792                               const gsl_vector_complex * tau,
0793                               gsl_matrix_complex * U,
0794                               gsl_vector * diag,
0795                               gsl_vector * sudiag);
0796
0797 int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
0798                                 gsl_vector * diag,
0799                                 gsl_vector * subdiag);
0800
0801 /* Linear Solve Using Householder Transformations
0802
0803  * exceptions:
0804  */
0805
0806 int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x);
0807 int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x);
0808
0809 /* Linear solve for a symmetric tridiagonal system.
0810
0811  * The input vectors represent the NxN matrix as follows:
0812  *
0813  *     diag[0]  offdiag[0]             0    ...
0814  *  offdiag[0]     diag[1]    offdiag[1]    ...
0815  *           0  offdiag[1]       diag[2]    ...
0816  *           0           0    offdiag[2]    ...
0817  *         ...         ...           ...    ...
0818  */
0819 int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag,
0820                                    const gsl_vector * offdiag,
0821                                    const gsl_vector * b,
0822                                    gsl_vector * x);
0823
0824 /* Linear solve for a nonsymmetric tridiagonal system.
0825
0826  * The input vectors represent the NxN matrix as follows:
0827  *
0828  *       diag[0]  abovediag[0]              0    ...
0829  *  belowdiag[0]       diag[1]   abovediag[1]    ...
0830  *             0  belowdiag[1]        diag[2]    ...
0831  *             0             0   belowdiag[2]    ...
0832  *           ...           ...            ...    ...
0833  */
0834 int gsl_linalg_solve_tridiag (const gsl_vector * diag,
0835                                    const gsl_vector * abovediag,
0836                                    const gsl_vector * belowdiag,
0837                                    const gsl_vector * b,
0838                                    gsl_vector * x);
0839
0840
0841 /* Linear solve for a symmetric cyclic tridiagonal system.
0842
0843  * The input vectors represent the NxN matrix as follows:
0844  *
0845  *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
0846  *   offdiag[0]     diag[1]    offdiag[1]   .....
0847  *            0  offdiag[1]       diag[2]   .....
0848  *            0           0    offdiag[2]   .....
0849  *          ...         ...
0850  * offdiag[N-1]         ...
0851  */
0852 int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag,
0853                                        const gsl_vector * offdiag,
0854                                        const gsl_vector * b,
0855                                        gsl_vector * x);
0856
0857 /* Linear solve for a nonsymmetric cyclic tridiagonal system.
0858
0859  * The input vectors represent the NxN matrix as follows:
0860  *
0861  *        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
0862  *   belowdiag[0]       diag[1]  abovediag[1]   .....
0863  *              0  belowdiag[1]       diag[2]
0864  *              0             0  belowdiag[2]   .....
0865  *            ...           ...
0866  * abovediag[N-1]           ...
0867  */
0868 int gsl_linalg_solve_cyc_tridiag (const gsl_vector * diag,
0869                                   const gsl_vector * abovediag,
0870                                   const gsl_vector * belowdiag,
0871                                   const gsl_vector * b,
0872                                   gsl_vector * x);
0873
0874
0875 /* Bidiagonal decomposition */
0876
0877 int gsl_linalg_bidiag_decomp (gsl_matrix * A,
0878                               gsl_vector * tau_U,
0879                               gsl_vector * tau_V);
0880
0881 int gsl_linalg_bidiag_unpack (const gsl_matrix * A,
0882                               const gsl_vector * tau_U,
0883                               gsl_matrix * U,
0884                               const gsl_vector * tau_V,
0885                               gsl_matrix * V,
0886                               gsl_vector * diag,
0887                               gsl_vector * superdiag);
0888
0889 int gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
0890                                gsl_vector * tau_U,
0891                                gsl_vector * tau_V,
0892                                gsl_matrix * V);
0893
0894 int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
0895                                 gsl_vector * diag,
0896                                 gsl_vector * superdiag);
0897
0898 /* Balancing */
0899
0900 int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D);
0901 int gsl_linalg_balance_accum (gsl_matrix * A, gsl_vector * D);
0902 int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D);
0903
0904 /* condition estimation */
0905
0906 int gsl_linalg_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work);
0907 int gsl_linalg_tri_upper_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work);
0908 int gsl_linalg_tri_lower_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work);
0909 int gsl_linalg_invnorm1(const size_t N,
0910                         int (* Ainvx)(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params),
0911                         void * params, double * Ainvnorm, gsl_vector * work);
0912
0913 /* triangular matrices */
0914
0915 int gsl_linalg_tri_upper_invert(gsl_matrix * T);
0916 int gsl_linalg_tri_lower_invert(gsl_matrix * T);
0917 int gsl_linalg_tri_upper_unit_invert(gsl_matrix * T);
0918 int gsl_linalg_tri_lower_unit_invert(gsl_matrix * T);
0919
0920 int gsl_linalg_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T);
0921 int gsl_linalg_complex_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T);
0922
0923 int gsl_linalg_tri_LTL(gsl_matrix * L);
0924 int gsl_linalg_tri_UL(gsl_matrix * LU);
0925 int gsl_linalg_complex_tri_LHL(gsl_matrix_complex * L);
0926 int gsl_linalg_complex_tri_UL(gsl_matrix_complex * LU);
0927
0928 INLINE_DECL void gsl_linalg_givens (const double a, const double b,
0929                                     double *c, double *s);
0930 INLINE_DECL void gsl_linalg_givens_gv (gsl_vector * v, const size_t i,
0931                                        const size_t j, const double c,
0932                                        const double s);
0933
0934 #ifdef HAVE_INLINE
0935
0936 /* Generate a Givens rotation (cos,sin) which takes v=(x,y) to (|v|,0)
0937    From Golub and Van Loan, "Matrix Computations", Section 5.1.8 */
0938 INLINE_FUN
0939 void
0940 gsl_linalg_givens (const double a, const double b, double *c, double *s)
0941 {
0942   if (b == 0)
0943     {
0944       *c = 1;
0945       *s = 0;
0946     }
0947   else if (fabs (b) > fabs (a))
0948     {
0949       double t = -a / b;
0950       double s1 = 1.0 / sqrt (1 + t * t);
0951       *s = s1;
0952       *c = s1 * t;
0953     }
0954   else
0955     {
0956       double t = -b / a;
0957       double c1 = 1.0 / sqrt (1 + t * t);
0958       *c = c1;
0959       *s = c1 * t;
0960     }
0961 } /* gsl_linalg_givens() */
0962
0963 INLINE_FUN
0964 void
0965 gsl_linalg_givens_gv (gsl_vector * v, const size_t i, const size_t j,
0966                       const double c, const double s)
0967 {
0968   /* Apply rotation to vector v' = G^T v */
0969
0970   double vi = gsl_vector_get (v, i);
0971   double vj = gsl_vector_get (v, j);
0972   gsl_vector_set (v, i, c * vi - s * vj);
0973   gsl_vector_set (v, j, s * vi + c * vj);
0974 }
0975
0976 #endif /* HAVE_INLINE */
0977
0978 __END_DECLS
0979
0980 #endif /* __GSL_LINALG_H__ */