include/flint/arb_mat.h

0001 /*
0002     Copyright (C) 2012 Fredrik Johansson
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef ARB_MAT_H
0013 #define ARB_MAT_H
0014
0015 #ifdef ARB_MAT_INLINES_C
0016 #define ARB_MAT_INLINE
0017 #else
0018 #define ARB_MAT_INLINE static inline
0019 #endif
0020
0021 #include "fmpq_types.h"
0022 #include "arb_types.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 #define arb_mat_entry(mat,i,j) ((mat)->entries + (i) * (mat)->stride + (j))
0029 #define arb_mat_nrows(mat) ((mat)->r)
0030 #define arb_mat_ncols(mat) ((mat)->c)
0031
0032 ARB_MAT_INLINE arb_ptr
0033 arb_mat_entry_ptr(arb_mat_t mat, slong i, slong j)
0034 {
0035     return arb_mat_entry(mat, i, j);
0036 }
0037
0038 /* Memory management */
0039
0040 void arb_mat_init(arb_mat_t mat, slong r, slong c);
0041
0042 void arb_mat_clear(arb_mat_t mat);
0043
0044 ARB_MAT_INLINE void
0045 arb_mat_swap(arb_mat_t mat1, arb_mat_t mat2)
0046 {
0047     FLINT_SWAP(arb_mat_struct, *mat1, *mat2);
0048 }
0049
0050 void arb_mat_swap_entrywise(arb_mat_t mat1, arb_mat_t mat2);
0051
0052 /* Window matrices */
0053
0054 void arb_mat_window_init(arb_mat_t window, const arb_mat_t mat, slong r1, slong c1, slong r2, slong c2);
0055
0056 ARB_MAT_INLINE void
0057 arb_mat_window_clear(arb_mat_t FLINT_UNUSED(window))
0058 {
0059 }
0060
0061 /* Conversions */
0062
0063 void arb_mat_set(arb_mat_t dest, const arb_mat_t src);
0064
0065 void arb_mat_set_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src);
0066
0067 void arb_mat_set_round_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src, slong prec);
0068
0069 void arb_mat_set_fmpq_mat(arb_mat_t dest, const fmpq_mat_t src, slong prec);
0070
0071 /* Random generation */
0072
0073 void arb_mat_randtest(arb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits);
0074
0075 void arb_mat_randtest_cho(arb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits);
0076
0077 void arb_mat_randtest_spd(arb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits);
0078
0079 /* I/O */
0080
0081 #ifdef FLINT_HAVE_FILE
0082 void arb_mat_fprintd(FILE * file, const arb_mat_t mat, slong digits);
0083 #endif
0084
0085 void arb_mat_printd(const arb_mat_t mat, slong digits);
0086
0087 /* Comparisons */
0088
0089 int arb_mat_eq(const arb_mat_t mat1, const arb_mat_t mat2);
0090
0091 int arb_mat_ne(const arb_mat_t mat1, const arb_mat_t mat2);
0092
0093 int arb_mat_equal(const arb_mat_t mat1, const arb_mat_t mat2);
0094
0095 int arb_mat_overlaps(const arb_mat_t mat1, const arb_mat_t mat2);
0096
0097 int arb_mat_contains(const arb_mat_t mat1, const arb_mat_t mat2);
0098
0099 int arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2);
0100
0101 int arb_mat_contains_fmpz_mat(const arb_mat_t mat1, const fmpz_mat_t mat2);
0102
0103 ARB_MAT_INLINE int
0104 arb_mat_is_empty(const arb_mat_t mat)
0105 {
0106     return (mat->r == 0) || (mat->c == 0);
0107 }
0108
0109 ARB_MAT_INLINE int
0110 arb_mat_is_square(const arb_mat_t mat)
0111 {
0112     return (mat->r == mat->c);
0113 }
0114
0115 int arb_mat_is_exact(const arb_mat_t A);
0116
0117 int arb_mat_is_zero(const arb_mat_t mat);
0118 int arb_mat_is_finite(const arb_mat_t mat);
0119 int arb_mat_is_triu(const arb_mat_t mat);
0120 int arb_mat_is_tril(const arb_mat_t mat);
0121
0122 ARB_MAT_INLINE int
0123 arb_mat_is_diag(const arb_mat_t mat)
0124 {
0125     return arb_mat_is_tril(mat) && arb_mat_is_triu(mat);
0126 }
0127
0128 /* Radius and interval operations */
0129
0130 void arb_mat_get_mid(arb_mat_t B, const arb_mat_t A);
0131
0132 void arb_mat_add_error_mag(arb_mat_t mat, const mag_t err);
0133
0134 /* Special matrices */
0135
0136 void arb_mat_zero(arb_mat_t mat);
0137
0138 void arb_mat_one(arb_mat_t mat);
0139
0140 void arb_mat_ones(arb_mat_t mat);
0141
0142 void arb_mat_indeterminate(arb_mat_t mat);
0143
0144 void arb_mat_hilbert(arb_mat_t mat, slong prec);
0145
0146 void arb_mat_pascal(arb_mat_t mat, int triangular, slong prec);
0147
0148 void arb_mat_stirling(arb_mat_t mat, int kind, slong prec);
0149
0150 void arb_mat_dct(arb_mat_t mat, int type, slong prec);
0151
0152 void arb_mat_transpose(arb_mat_t mat1, const arb_mat_t mat2);
0153
0154 /* Norms */
0155
0156 void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A);
0157
0158 void arb_mat_frobenius_norm(arb_t res, const arb_mat_t A, slong prec);
0159
0160 void arb_mat_bound_frobenius_norm(mag_t b, const arb_mat_t A);
0161
0162 /* Arithmetic */
0163
0164 void arb_mat_neg(arb_mat_t dest, const arb_mat_t src);
0165
0166 void arb_mat_add(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec);
0167
0168 void arb_mat_sub(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec);
0169
0170 void arb_mat_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec);
0171
0172 void arb_mat_mul_classical(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec);
0173
0174 void arb_mat_mul_threaded(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec);
0175
0176 void _arb_mat_addmul_rad_mag_fast(arb_mat_t C, mag_srcptr A, mag_srcptr B, slong ar, slong ac, slong bc);
0177
0178 void arb_mat_mul_block(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec);
0179
0180 void arb_mat_mul_entrywise(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec);
0181
0182 void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec);
0183 void arb_mat_sqr(arb_mat_t B, const arb_mat_t A, slong prec);
0184
0185 void arb_mat_pow_ui(arb_mat_t B, const arb_mat_t A, ulong exp, slong prec);
0186
0187 /* Scalar arithmetic */
0188
0189 void arb_mat_scalar_mul_2exp_si(arb_mat_t B, const arb_mat_t A, slong c);
0190
0191 void arb_mat_scalar_mul_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec);
0192 void arb_mat_scalar_mul_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec);
0193 void arb_mat_scalar_mul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec);
0194
0195 void arb_mat_scalar_addmul_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec);
0196 void arb_mat_scalar_addmul_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec);
0197 void arb_mat_scalar_addmul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec);
0198
0199 void arb_mat_scalar_div_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec);
0200 void arb_mat_scalar_div_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec);
0201 void arb_mat_scalar_div_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec);
0202
0203 /* Vector arithmetic */
0204
0205 void _arb_mat_vector_mul_row(arb_ptr res, arb_srcptr v, const arb_mat_t A, slong prec);
0206
0207 void _arb_mat_vector_mul_col(arb_ptr res, const arb_mat_t A, arb_srcptr v, slong prec);
0208
0209 void arb_mat_vector_mul_row(arb_ptr res, arb_srcptr v, const arb_mat_t A, slong prec);
0210
0211 void arb_mat_vector_mul_col(arb_ptr res, const arb_mat_t A, arb_srcptr v, slong prec);
0212
0213 /* Solving */
0214
0215 void arb_mat_swap_rows(arb_mat_t mat, slong * perm, slong r, slong s);
0216
0217 slong arb_mat_find_pivot_partial(const arb_mat_t mat,
0218                                     slong start_row, slong end_row, slong c);
0219
0220 void arb_mat_solve_tril_classical(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec);
0221 void arb_mat_solve_tril_recursive(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec);
0222 void arb_mat_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec);
0223
0224 void arb_mat_solve_triu_classical(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec);
0225 void arb_mat_solve_triu_recursive(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec);
0226 void arb_mat_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec);
0227
0228 int arb_mat_lu_classical(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec);
0229 int arb_mat_lu_recursive(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec);
0230 int arb_mat_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec);
0231
0232 void arb_mat_solve_lu_precomp(arb_mat_t X, const slong * perm,
0233     const arb_mat_t A, const arb_mat_t B, slong prec);
0234
0235 int arb_mat_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec);
0236
0237 int arb_mat_solve_lu(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec);
0238
0239 int arb_mat_solve_precond(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec);
0240
0241 int arb_mat_solve_preapprox(arb_mat_t X, const arb_mat_t A,
0242     const arb_mat_t B, const arb_mat_t R, const arb_mat_t T, slong prec);
0243
0244 void arb_mat_approx_mul(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec);
0245 void arb_mat_approx_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec);
0246 void arb_mat_approx_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec);
0247 int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec);
0248 void arb_mat_approx_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t A, const arb_mat_t B, slong prec);
0249 int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec);
0250 int arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec);
0251
0252 int arb_mat_inv(arb_mat_t X, const arb_mat_t A, slong prec);
0253
0254 void arb_mat_det_lu(arb_t det, const arb_mat_t A, slong prec);
0255 void arb_mat_det_precond(arb_t det, const arb_mat_t A, slong prec);
0256 void arb_mat_det(arb_t det, const arb_mat_t A, slong prec);
0257
0258 int _arb_mat_cholesky_banachiewicz(arb_mat_t A, slong prec);
0259
0260 int arb_mat_cho(arb_mat_t L, const arb_mat_t A, slong prec);
0261
0262 void arb_mat_solve_cho_precomp(arb_mat_t X,
0263     const arb_mat_t L, const arb_mat_t B, slong prec);
0264
0265 void arb_mat_inv_cho_precomp(arb_mat_t X, const arb_mat_t L, slong prec);
0266
0267 int arb_mat_spd_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec);
0268
0269 int arb_mat_spd_inv(arb_mat_t X, const arb_mat_t A, slong prec);
0270
0271 int _arb_mat_ldl_inplace(arb_mat_t A, slong prec);
0272
0273 int _arb_mat_ldl_golub_and_van_loan(arb_mat_t A, slong prec);
0274
0275 int arb_mat_ldl(arb_mat_t L, const arb_mat_t A, slong prec);
0276
0277 void arb_mat_solve_ldl_precomp(arb_mat_t X,
0278     const arb_mat_t L, const arb_mat_t B, slong prec);
0279
0280 void arb_mat_inv_ldl_precomp(arb_mat_t X, const arb_mat_t L, slong prec);
0281
0282 /* Special functions */
0283
0284 void arb_mat_exp_taylor_sum(arb_mat_t S, const arb_mat_t A, slong N, slong prec);
0285
0286 void arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec);
0287
0288 void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec);
0289 void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec);
0290 void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec);
0291 void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec);
0292
0293 void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec);
0294
0295 void _arb_mat_diag_prod(arb_t res, const arb_mat_t A, slong a, slong b, slong prec);
0296 void arb_mat_diag_prod(arb_t res, const arb_mat_t A, slong prec);
0297
0298 /* Sparsity structure */
0299
0300 slong arb_mat_count_is_zero(const arb_mat_t mat);
0301
0302 ARB_MAT_INLINE slong
0303 arb_mat_count_not_is_zero(const arb_mat_t mat)
0304 {
0305     slong size;
0306     size = arb_mat_nrows(mat) * arb_mat_ncols(mat);
0307     return size - arb_mat_count_is_zero(mat);
0308 }
0309
0310 slong arb_mat_allocated_bytes(const arb_mat_t x);
0311
0312 /* LLL reduction */
0313
0314 int arb_mat_spd_get_fmpz_mat(fmpz_mat_t B, const arb_mat_t A, slong prec);
0315
0316 void arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, double delta, double eta, slong prec);
0317
0318 int arb_mat_spd_is_lll_reduced(const arb_mat_t A, double delta, double eta, slong prec);
0319
0320 #ifdef __cplusplus
0321 }
0322 #endif
0323
0324 #endif