1.. _arb-mat: 2 3**arb_mat.h** -- matrices over the real numbers 4=============================================================================== 5 6An :type:`arb_mat_t` represents a dense matrix over the real numbers, 7implemented as an array of entries of type :type:`arb_struct`. 8The dimension (number of rows and columns) of a matrix is fixed at 9initialization, and the user must ensure that inputs and outputs to 10an operation have compatible dimensions. The number of rows or columns 11in a matrix can be zero. 12 13.. note:: 14 15 Methods prefixed with *arb_mat_approx* treat all input entries 16 as floating-point numbers (ignoring the radii of the balls) and 17 compute floating-point output (balls with zero radius) representing 18 approximate solutions *without error bounds*. 19 All other methods compute rigorous error bounds. 20 The *approx* methods are typically useful for computing initial 21 values or preconditioners for rigorous solvers. 22 Some users may also find *approx* methods useful 23 for doing ordinary numerical linear algebra in applications where 24 error bounds are not needed. 25 26Types, macros and constants 27------------------------------------------------------------------------------- 28 29.. type:: arb_mat_struct 30 31.. type:: arb_mat_t 32 33 Contains a pointer to a flat array of the entries (entries), an array of 34 pointers to the start of each row (rows), and the number of rows (r) 35 and columns (c). 36 37 An *arb_mat_t* is defined as an array of length one of type 38 *arb_mat_struct*, permitting an *arb_mat_t* to 39 be passed by reference. 40 41.. macro:: arb_mat_entry(mat, i, j) 42 43 Macro giving a pointer to the entry at row *i* and column *j*. 44 45.. macro:: arb_mat_nrows(mat) 46 47 Returns the number of rows of the matrix. 48 49.. macro:: arb_mat_ncols(mat) 50 51 Returns the number of columns of the matrix. 52 53 54Memory management 55------------------------------------------------------------------------------- 56 57.. function:: void arb_mat_init(arb_mat_t mat, slong r, slong c) 58 59 Initializes the matrix, setting it to the zero matrix with *r* rows 60 and *c* columns. 61 62.. function:: void arb_mat_clear(arb_mat_t mat) 63 64 Clears the matrix, deallocating all entries. 65 66.. function:: slong arb_mat_allocated_bytes(const arb_mat_t x) 67 68 Returns the total number of bytes heap-allocated internally by this object. 69 The count excludes the size of the structure itself. Add 70 ``sizeof(arb_mat_struct)`` to get the size of the object as a whole. 71 72.. function:: void arb_mat_window_init(arb_mat_t window, const arb_mat_t mat, slong r1, slong c1, slong r2, slong c2) 73 74 Initializes *window* to a window matrix into the submatrix of *mat* 75 starting at the corner at row *r1* and column *c1* (inclusive) and ending 76 at row *r2* and column *c2* (exclusive). 77 78.. function:: void arb_mat_window_clear(arb_mat_t window) 79 80 Frees the window matrix. 81 82Conversions 83------------------------------------------------------------------------------- 84 85.. function:: void arb_mat_set(arb_mat_t dest, const arb_mat_t src) 86 87.. function:: void arb_mat_set_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src) 88 89.. function:: void arb_mat_set_round_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src, slong prec) 90 91.. function:: void arb_mat_set_fmpq_mat(arb_mat_t dest, const fmpq_mat_t src, slong prec) 92 93 Sets *dest* to *src*. The operands must have identical dimensions. 94 95Random generation 96------------------------------------------------------------------------------- 97 98.. function:: void arb_mat_randtest(arb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits) 99 100 Sets *mat* to a random matrix with up to *prec* bits of precision 101 and with exponents of width up to *mag_bits*. 102 103Input and output 104------------------------------------------------------------------------------- 105 106.. function:: void arb_mat_printd(const arb_mat_t mat, slong digits) 107 108 Prints each entry in the matrix with the specified number of decimal digits. 109 110.. function:: void arb_mat_fprintd(FILE * file, const arb_mat_t mat, slong digits) 111 112 Prints each entry in the matrix with the specified number of decimal 113 digits to the stream *file*. 114 115Comparisons 116------------------------------------------------------------------------------- 117 118Predicate methods return 1 if the property certainly holds and 0 otherwise. 119 120.. function:: int arb_mat_equal(const arb_mat_t mat1, const arb_mat_t mat2) 121 122 Returns whether the matrices have the same dimensions 123 and identical intervals as entries. 124 125.. function:: int arb_mat_overlaps(const arb_mat_t mat1, const arb_mat_t mat2) 126 127 Returns whether the matrices have the same dimensions 128 and each entry in *mat1* overlaps with the corresponding entry in *mat2*. 129 130.. function:: int arb_mat_contains(const arb_mat_t mat1, const arb_mat_t mat2) 131 132.. function:: int arb_mat_contains_fmpz_mat(const arb_mat_t mat1, const fmpz_mat_t mat2) 133 134.. function:: int arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2) 135 136 Returns whether the matrices have the same dimensions and each entry 137 in *mat2* is contained in the corresponding entry in *mat1*. 138 139.. function:: int arb_mat_eq(const arb_mat_t mat1, const arb_mat_t mat2) 140 141 Returns whether *mat1* and *mat2* certainly represent the same matrix. 142 143.. function:: int arb_mat_ne(const arb_mat_t mat1, const arb_mat_t mat2) 144 145 Returns whether *mat1* and *mat2* certainly do not represent the same matrix. 146 147.. function:: int arb_mat_is_empty(const arb_mat_t mat) 148 149 Returns whether the number of rows or the number of columns in *mat* is zero. 150 151.. function:: int arb_mat_is_square(const arb_mat_t mat) 152 153 Returns whether the number of rows is equal to the number of columns in *mat*. 154 155.. function:: int arb_mat_is_exact(const arb_mat_t mat) 156 157 Returns whether all entries in *mat* have zero radius. 158 159.. function:: int arb_mat_is_zero(const arb_mat_t mat) 160 161 Returns whether all entries in *mat* are exactly zero. 162 163.. function:: int arb_mat_is_finite(const arb_mat_t mat) 164 165 Returns whether all entries in *mat* are finite. 166 167.. function:: int arb_mat_is_triu(const arb_mat_t mat) 168 169 Returns whether *mat* is upper triangular; that is, all entries 170 below the main diagonal are exactly zero. 171 172.. function:: int arb_mat_is_tril(const arb_mat_t mat) 173 174 Returns whether *mat* is lower triangular; that is, all entries 175 above the main diagonal are exactly zero. 176 177.. function:: int arb_mat_is_diag(const arb_mat_t mat) 178 179 Returns whether *mat* is a diagonal matrix; that is, all entries 180 off the main diagonal are exactly zero. 181 182Special matrices 183------------------------------------------------------------------------------- 184 185.. function:: void arb_mat_zero(arb_mat_t mat) 186 187 Sets all entries in mat to zero. 188 189.. function:: void arb_mat_one(arb_mat_t mat) 190 191 Sets the entries on the main diagonal to ones, 192 and all other entries to zero. 193 194.. function:: void arb_mat_ones(arb_mat_t mat) 195 196 Sets all entries in the matrix to ones. 197 198.. function:: void arb_mat_indeterminate(arb_mat_t mat) 199 200 Sets all entries in the matrix to indeterminate (NaN). 201 202.. function:: void arb_mat_hilbert(arb_mat_t mat, slong prec) 203 204 Sets *mat* to the Hilbert matrix, which has entries `A_{j,k} = 1/(j+k+1)`. 205 206.. function:: void arb_mat_pascal(arb_mat_t mat, int triangular, slong prec) 207 208 Sets *mat* to a Pascal matrix, whose entries are binomial coefficients. 209 If *triangular* is 0, constructs a full symmetric matrix 210 with the rows of Pascal's triangle as successive antidiagonals. 211 If *triangular* is 1, constructs the upper triangular matrix with 212 the rows of Pascal's triangle as columns, and if *triangular* is -1, 213 constructs the lower triangular matrix with the rows of Pascal's 214 triangle as rows. 215 216 The entries are computed using recurrence relations. 217 When the dimensions get large, some precision loss is possible; in that 218 case, the user may wish to create the matrix at slightly higher precision 219 and then round it to the final precision. 220 221.. function:: void arb_mat_stirling(arb_mat_t mat, int kind, slong prec) 222 223 Sets *mat* to a Stirling matrix, whose entries are Stirling numbers. 224 If *kind* is 0, the entries are set to the unsigned Stirling numbers 225 of the first kind. If *kind* is 1, the entries are set to the signed 226 Stirling numbers of the first kind. If *kind* is 2, the entries are 227 set to the Stirling numbers of the second kind. 228 229 The entries are computed using recurrence relations. 230 When the dimensions get large, some precision loss is possible; in that 231 case, the user may wish to create the matrix at slightly higher precision 232 and then round it to the final precision. 233 234.. function:: void arb_mat_dct(arb_mat_t mat, int type, slong prec) 235 236 Sets *mat* to the DCT (discrete cosine transform) matrix of order *n* 237 where *n* is the smallest dimension of *mat* (if *mat* is not square, 238 the matrix is extended periodically along the larger dimension). 239 There are many different conventions for defining DCT matrices; here, 240 we use the normalized "DCT-II" transform matrix 241 242 .. math :: 243 244 A_{j,k} = \sqrt{\frac{2}{n}} \cos\left(\frac{\pi j}{n} \left(k+\frac{1}{2}\right)\right) 245 246 which satisfies `A^{-1} = A^T`. 247 The *type* parameter is currently ignored and should be set to 0. 248 In the future, it might be used to select a different convention. 249 250Transpose 251------------------------------------------------------------------------------- 252 253.. function:: void arb_mat_transpose(arb_mat_t dest, const arb_mat_t src) 254 255 Sets *dest* to the exact transpose *src*. The operands must have 256 compatible dimensions. Aliasing is allowed. 257 258Norms 259------------------------------------------------------------------------------- 260 261.. function:: void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A) 262 263 Sets *b* to an upper bound for the infinity norm (i.e. the largest 264 absolute value row sum) of *A*. 265 266.. function:: void arb_mat_frobenius_norm(arb_t res, const arb_mat_t A, slong prec) 267 268 Sets *res* to the Frobenius norm (i.e. the square root of the sum 269 of squares of entries) of *A*. 270 271.. function:: void arb_mat_bound_frobenius_norm(mag_t res, const arb_mat_t A) 272 273 Sets *res* to an upper bound for the Frobenius norm of *A*. 274 275Arithmetic 276------------------------------------------------------------------------------- 277 278.. function:: void arb_mat_neg(arb_mat_t dest, const arb_mat_t src) 279 280 Sets *dest* to the exact negation of *src*. The operands must have 281 the same dimensions. 282 283.. function:: void arb_mat_add(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec) 284 285 Sets res to the sum of *mat1* and *mat2*. The operands must have the same dimensions. 286 287.. function:: void arb_mat_sub(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec) 288 289 Sets *res* to the difference of *mat1* and *mat2*. The operands must have 290 the same dimensions. 291 292.. function:: void arb_mat_mul_classical(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec) 293 294.. function:: void arb_mat_mul_threaded(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec) 295 296.. function:: void arb_mat_mul_block(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec) 297 298.. function:: void arb_mat_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec) 299 300 Sets *res* to the matrix product of *mat1* and *mat2*. The operands must have 301 compatible dimensions for matrix multiplication. 302 303 The *classical* version performs matrix multiplication in the trivial way. 304 305 The *block* version decomposes the input matrices into one or several 306 blocks of uniformly scaled matrices and multiplies 307 large blocks via *fmpz_mat_mul*. It also invokes 308 :func:`_arb_mat_addmul_rad_mag_fast` for the radius matrix multiplications. 309 310 The *threaded* version performs classical multiplication but splits the 311 computation over the number of threads returned by *flint_get_num_threads()*. 312 313 The default version chooses an algorithm automatically. 314 315.. function:: void arb_mat_mul_entrywise(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec) 316 317 Sets *C* to the entrywise product of *A* and *B*. 318 The operands must have the same dimensions. 319 320.. function:: void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec) 321 322.. function:: void arb_mat_sqr(arb_mat_t res, const arb_mat_t mat, slong prec) 323 324 Sets *res* to the matrix square of *mat*. The operands must both be square 325 with the same dimensions. 326 327.. function:: void arb_mat_pow_ui(arb_mat_t res, const arb_mat_t mat, ulong exp, slong prec) 328 329 Sets *res* to *mat* raised to the power *exp*. Requires that *mat* 330 is a square matrix. 331 332.. function:: void _arb_mat_addmul_rad_mag_fast(arb_mat_t C, mag_srcptr A, mag_srcptr B, slong ar, slong ac, slong bc) 333 334 Helper function for matrix multiplication. 335 Adds to the radii of *C* the matrix product of the matrices represented 336 by *A* and *B*, where *A* is a linear array of coefficients in row-major 337 order and *B* is a linear array of coefficients in column-major order. 338 This function assumes that all exponents are small and is unsafe 339 for general use. 340 341.. function:: void arb_mat_approx_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec) 342 343 Approximate matrix multiplication. The input radii are ignored and 344 the output matrix is set to an approximate floating-point result. 345 The radii in the output matrix will *not* necessarily be zeroed. 346 347Scalar arithmetic 348------------------------------------------------------------------------------- 349 350.. function:: void arb_mat_scalar_mul_2exp_si(arb_mat_t B, const arb_mat_t A, slong c) 351 352 Sets *B* to *A* multiplied by `2^c`. 353 354.. function:: void arb_mat_scalar_addmul_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec) 355 356.. function:: void arb_mat_scalar_addmul_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec) 357 358.. function:: void arb_mat_scalar_addmul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec) 359 360 Sets *B* to `B + A \times c`. 361 362.. function:: void arb_mat_scalar_mul_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec) 363 364.. function:: void arb_mat_scalar_mul_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec) 365 366.. function:: void arb_mat_scalar_mul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec) 367 368 Sets *B* to `A \times c`. 369 370.. function:: void arb_mat_scalar_div_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec) 371 372.. function:: void arb_mat_scalar_div_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec) 373 374.. function:: void arb_mat_scalar_div_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec) 375 376 Sets *B* to `A / c`. 377 378 379Gaussian elimination and solving 380------------------------------------------------------------------------------- 381 382.. function:: int arb_mat_lu_classical(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec) 383 384.. function:: int arb_mat_lu_recursive(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec) 385 386.. function:: int arb_mat_lu(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec) 387 388 Given an `n \times n` matrix `A`, computes an LU decomposition `PLU = A` 389 using Gaussian elimination with partial pivoting. 390 The input and output matrices can be the same, performing the 391 decomposition in-place. 392 393 Entry `i` in the permutation vector perm is set to the row index in 394 the input matrix corresponding to row `i` in the output matrix. 395 396 The algorithm succeeds and returns nonzero if it can find `n` invertible 397 (i.e. not containing zero) pivot entries. This guarantees that the matrix 398 is invertible. 399 400 The algorithm fails and returns zero, leaving the entries in `P` and `LU` 401 undefined, if it cannot find `n` invertible pivot elements. 402 In this case, either the matrix is singular, the input matrix was 403 computed to insufficient precision, or the LU decomposition was 404 attempted at insufficient precision. 405 406 The *classical* version uses Gaussian elimination directly while 407 the *recursive* version performs the computation in a block recursive 408 way to benefit from fast matrix multiplication. The default version 409 chooses an algorithm automatically. 410 411.. function:: void arb_mat_solve_tril_classical(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec) 412 413.. function:: void arb_mat_solve_tril_recursive(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec) 414 415.. function:: void arb_mat_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec) 416 417.. function:: void arb_mat_solve_triu_classical(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec) 418 419.. function:: void arb_mat_solve_triu_recursive(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec) 420 421.. function:: void arb_mat_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec) 422 423 Solves the lower triangular system `LX = B` or the upper triangular system 424 `UX = B`, respectively. If *unit* is set, the main diagonal of *L* or *U* 425 is taken to consist of all ones, and in that case the actual entries on 426 the diagonal are not read at all and can contain other data. 427 428 The *classical* versions perform the computations iteratively while the 429 *recursive* versions perform the computations in a block recursive 430 way to benefit from fast matrix multiplication. The default versions 431 choose an algorithm automatically. 432 433.. function:: void arb_mat_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t LU, const arb_mat_t B, slong prec) 434 435 Solves `AX = B` given the precomputed nonsingular LU decomposition `A = PLU`. 436 The matrices `X` and `B` are allowed to be aliased with each other, 437 but `X` is not allowed to be aliased with `LU`. 438 439.. function:: int arb_mat_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec) 440 441.. function:: int arb_mat_solve_lu(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec) 442 443.. function:: int arb_mat_solve_precond(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec) 444 445 Solves `AX = B` where `A` is a nonsingular `n \times n` matrix 446 and `X` and `B` are `n \times m` matrices. 447 448 If `m > 0` and `A` cannot be inverted numerically (indicating either that 449 `A` is singular or that the precision is insufficient), the values in the 450 output matrix are left undefined and zero is returned. A nonzero return 451 value guarantees that `A` is invertible and that the exact solution 452 matrix is contained in the output. 453 454 Three algorithms are provided: 455 456 * The *lu* version performs LU decomposition directly in ball arithmetic. 457 This is fast, but the bounds typically blow up exponentially with *n*, 458 even if the system is well-conditioned. This algorithm is usually 459 the best choice at very high precision. 460 * The *precond* version computes an approximate inverse to precondition 461 the system [HS1967]_. This is usually several times slower than direct LU 462 decomposition, but the bounds do not blow up with *n* if the system is 463 well-conditioned. This algorithm is usually 464 the best choice for large systems at low to moderate precision. 465 * The default version selects between *lu* and *precomp* automatically. 466 467 The automatic choice should be reasonable most of the time, but users 468 may benefit from trying either *lu* or *precond* in specific applications. 469 For example, the *lu* solver often performs better for ill-conditioned 470 systems where use of very high precision is unavoidable. 471 472.. function:: int arb_mat_solve_preapprox(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, const arb_mat_t R, const arb_mat_t T, slong prec) 473 474 Solves `AX = B` where `A` is a nonsingular `n \times n` matrix 475 and `X` and `B` are `n \times m` matrices, given an approximation 476 `R` of the matrix inverse of `A`, and given the approximation `T` 477 of the solution `X`. 478 479 If `m > 0` and `A` cannot be inverted numerically (indicating either that 480 `A` is singular or that the precision is insufficient, or that `R` is 481 not a close enough approximation of the inverse of `A`), the values in the 482 output matrix are left undefined and zero is returned. A nonzero return 483 value guarantees that `A` is invertible and that the exact solution 484 matrix is contained in the output. 485 486.. function:: int arb_mat_inv(arb_mat_t X, const arb_mat_t A, slong prec) 487 488 Sets `X = A^{-1}` where `A` is a square matrix, computed by solving 489 the system `AX = I`. 490 491 If `A` cannot be inverted numerically (indicating either that 492 `A` is singular or that the precision is insufficient), the values in the 493 output matrix are left undefined and zero is returned. 494 A nonzero return value guarantees that the matrix is invertible 495 and that the exact inverse is contained in the output. 496 497.. function:: void arb_mat_det_lu(arb_t det, const arb_mat_t A, slong prec) 498 499.. function:: void arb_mat_det_precond(arb_t det, const arb_mat_t A, slong prec) 500 501.. function:: void arb_mat_det(arb_t det, const arb_mat_t A, slong prec) 502 503 Sets *det* to the determinant of the matrix *A*. 504 505 The *lu* version uses Gaussian elimination with partial pivoting. If at 506 some point an invertible pivot element cannot be found, the elimination is 507 stopped and the magnitude of the determinant of the remaining submatrix 508 is bounded using Hadamard's inequality. 509 510 The *precond* version computes an approximate LU factorization of *A* 511 and multiplies by the inverse *L* and *U* martices as preconditioners 512 to obtain a matrix close to the identity matrix [Rum2010]_. An enclosure 513 for this determinant is computed using Gershgorin circles. This is about 514 four times slower than direct Gaussian elimination, but much more 515 numerically stable. 516 517 The default version automatically selects between the *lu* and *precond* 518 versions and additionally handles small or triangular matrices 519 by direct formulas. 520 521.. function:: void arb_mat_approx_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec) 522 523.. function:: void arb_mat_approx_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec) 524 525.. function:: int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec) 526 527.. function:: 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) 528 529.. function:: int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec) 530 531.. function:: int arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec) 532 533 These methods perform approximate solving *without any error control*. 534 The radii in the input matrices are ignored, the computations are done 535 numerically with floating-point arithmetic (using ordinary 536 Gaussian elimination and triangular solving, accelerated through 537 the use of block recursive strategies for large matrices), and the 538 output matrices are set to the approximate floating-point results with 539 zeroed error bounds. 540 541 Approximate solutions are useful for computing preconditioning matrices 542 for certified solutions. Some users may also find these methods useful 543 for doing ordinary numerical linear algebra in applications where 544 error bounds are not needed. 545 546Cholesky decomposition and solving 547------------------------------------------------------------------------------- 548 549.. function:: int _arb_mat_cholesky_banachiewicz(arb_mat_t A, slong prec) 550 551.. function:: int arb_mat_cho(arb_mat_t L, const arb_mat_t A, slong prec) 552 553 Computes the Cholesky decomposition of *A*, returning nonzero iff 554 the symmetric matrix defined by the lower triangular part of *A* 555 is certainly positive definite. 556 557 If a nonzero value is returned, then *L* is set to the lower triangular 558 matrix such that `A = L * L^T`. 559 560 If zero is returned, then either the matrix is not symmetric positive 561 definite, the input matrix was computed to insufficient precision, 562 or the decomposition was attempted at insufficient precision. 563 564 The underscore method computes *L* from *A* in-place, leaving the 565 strict upper triangular region undefined. 566 567.. function:: void arb_mat_solve_cho_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec) 568 569 Solves `AX = B` given the precomputed Cholesky decomposition `A = L L^T`. 570 The matrices *X* and *B* are allowed to be aliased with each other, 571 but *X* is not allowed to be aliased with *L*. 572 573.. function:: int arb_mat_spd_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec) 574 575 Solves `AX = B` where *A* is a symmetric positive definite matrix 576 and *X* and *B* are `n \times m` matrices, using Cholesky decomposition. 577 578 If `m > 0` and *A* cannot be factored using Cholesky decomposition 579 (indicating either that *A* is not symmetric positive definite or that 580 the precision is insufficient), the values in the 581 output matrix are left undefined and zero is returned. A nonzero return 582 value guarantees that the symmetric matrix defined through the lower 583 triangular part of *A* is invertible and that the exact solution matrix 584 is contained in the output. 585 586.. function:: void arb_mat_inv_cho_precomp(arb_mat_t X, const arb_mat_t L, slong prec) 587 588 Sets `X = A^{-1}` where `A` is a symmetric positive definite matrix 589 whose Cholesky decomposition *L* has been computed with 590 :func:`arb_mat_cho`. 591 The inverse is calculated using the method of [Kri2013]_ which is more 592 efficient than solving `AX = I` with :func:`arb_mat_solve_cho_precomp`. 593 594.. function:: int arb_mat_spd_inv(arb_mat_t X, const arb_mat_t A, slong prec) 595 596 Sets `X = A^{-1}` where *A* is a symmetric positive definite matrix. 597 It is calculated using the method of [Kri2013]_ which computes fewer 598 intermediate results than solving `AX = I` with :func:`arb_mat_spd_solve`. 599 600 If *A* cannot be factored using Cholesky decomposition 601 (indicating either that *A* is not symmetric positive definite or that 602 the precision is insufficient), the values in the 603 output matrix are left undefined and zero is returned. A nonzero return 604 value guarantees that the symmetric matrix defined through the lower 605 triangular part of *A* is invertible and that the exact inverse 606 is contained in the output. 607 608.. function:: int _arb_mat_ldl_inplace(arb_mat_t A, slong prec) 609 610.. function:: int _arb_mat_ldl_golub_and_van_loan(arb_mat_t A, slong prec) 611 612.. function:: int arb_mat_ldl(arb_mat_t res, const arb_mat_t A, slong prec) 613 614 Computes the `LDL^T` decomposition of *A*, returning nonzero iff 615 the symmetric matrix defined by the lower triangular part of *A* 616 is certainly positive definite. 617 618 If a nonzero value is returned, then *res* is set to a lower triangular 619 matrix that encodes the `L * D * L^T` decomposition of *A*. 620 In particular, `L` is a lower triangular matrix with ones on its diagonal 621 and whose strictly lower triangular region is the same as that of *res*. 622 `D` is a diagonal matrix with the same diagonal as that of *res*. 623 624 If zero is returned, then either the matrix is not symmetric positive 625 definite, the input matrix was computed to insufficient precision, 626 or the decomposition was attempted at insufficient precision. 627 628 The underscore methods compute *res* from *A* in-place, leaving the 629 strict upper triangular region undefined. 630 The default method uses algorithm 4.1.2 from [GVL1996]_. 631 632.. function:: void arb_mat_solve_ldl_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec) 633 634 Solves `AX = B` given the precomputed `A = LDL^T` decomposition 635 encoded by *L*. The matrices *X* and *B* are allowed to be aliased 636 with each other, but *X* is not allowed to be aliased with *L*. 637 638.. function:: void arb_mat_inv_ldl_precomp(arb_mat_t X, const arb_mat_t L, slong prec) 639 640 Sets `X = A^{-1}` where `A` is a symmetric positive definite matrix 641 whose `LDL^T` decomposition encoded by *L* has been computed with 642 :func:`arb_mat_ldl`. 643 The inverse is calculated using the method of [Kri2013]_ which is more 644 efficient than solving `AX = I` with :func:`arb_mat_solve_ldl_precomp`. 645 646Characteristic polynomial and companion matrix 647------------------------------------------------------------------------------- 648 649.. function:: void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec) 650 651.. function:: void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec) 652 653 Sets *poly* to the characteristic polynomial of *mat* which must be 654 a square matrix. If the matrix has *n* rows, the underscore method 655 requires space for `n + 1` output coefficients. 656 Employs a division-free algorithm using `O(n^4)` operations. 657 658.. function:: void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec) 659 660.. function:: void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec) 661 662 Sets the *n* by *n* matrix *mat* to the companion matrix of the polynomial 663 *poly* which must have degree *n*. 664 The underscore method reads `n + 1` input coefficients. 665 666Special functions 667------------------------------------------------------------------------------- 668 669.. function:: void arb_mat_exp_taylor_sum(arb_mat_t S, const arb_mat_t A, slong N, slong prec) 670 671 Sets *S* to the truncated exponential Taylor series `S = \sum_{k=0}^{N-1} A^k / k!`. 672 Uses rectangular splitting to compute the sum using `O(\sqrt{N})` 673 matrix multiplications. The recurrence relation for factorials 674 is used to get scalars that are small integers instead of full 675 factorials. As in [Joh2014b]_, all divisions are postponed to 676 the end by computing partial factorials of length `O(\sqrt{N})`. 677 The scalars could be reduced by doing more divisions, but this 678 appears to be slower in most cases. 679 680.. function:: void arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec) 681 682 Sets *B* to the exponential of the matrix *A*, defined by the Taylor series 683 684 .. math :: 685 686 \exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}. 687 688 The function is evaluated as `\exp(A/2^r)^{2^r}`, where `r` is chosen 689 to give rapid convergence. 690 691 The elementwise error when truncating the Taylor series after *N* 692 terms is bounded by the error in the infinity norm, for which we have 693 694 .. math :: 695 \left\|\exp(2^{-r}A) - \sum_{k=0}^{N-1} 696 \frac{\left(2^{-r} A\right)^k}{k!} \right\|_{\infty} = 697 \left\|\sum_{k=N}^{\infty} \frac{\left(2^{-r} A\right)^k}{k!}\right\|_{\infty} \le 698 \sum_{k=N}^{\infty} \frac{(2^{-r} \|A\|_{\infty})^k}{k!}. 699 700 We bound the sum on the right using :func:`mag_exp_tail`. 701 Truncation error is not added to entries whose values are determined 702 by the sparsity structure of `A`. 703 704.. function:: void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec) 705 706 Sets *trace* to the trace of the matrix, i.e. the sum of entries on the 707 main diagonal of *mat*. The matrix is required to be square. 708 709.. function:: void _arb_mat_diag_prod(arb_t res, const arb_mat_t mat, slong a, slong b, slong prec) 710 711.. function:: void arb_mat_diag_prod(arb_t res, const arb_mat_t mat, slong prec) 712 713 Sets *res* to the product of the entries on the main diagonal of *mat*. 714 The underscore method computes the product of the entries between 715 index *a* inclusive and *b* exclusive (the indices must be in range). 716 717Sparsity structure 718------------------------------------------------------------------------------- 719 720.. function:: void arb_mat_entrywise_is_zero(fmpz_mat_t dest, const arb_mat_t src) 721 722 Sets each entry of *dest* to indicate whether the corresponding 723 entry of *src* is certainly zero. 724 If the entry of *src* at row `i` and column `j` is zero according to 725 :func:`arb_is_zero` then the entry of *dest* at that row and column 726 is set to one, otherwise that entry of *dest* is set to zero. 727 728.. function:: void arb_mat_entrywise_not_is_zero(fmpz_mat_t dest, const arb_mat_t src) 729 730 Sets each entry of *dest* to indicate whether the corresponding 731 entry of *src* is not certainly zero. 732 This the complement of :func:`arb_mat_entrywise_is_zero`. 733 734.. function:: slong arb_mat_count_is_zero(const arb_mat_t mat) 735 736 Returns the number of entries of *mat* that are certainly zero 737 according to :func:`arb_is_zero`. 738 739.. function:: slong arb_mat_count_not_is_zero(const arb_mat_t mat) 740 741 Returns the number of entries of *mat* that are not certainly zero. 742 743Component and error operations 744------------------------------------------------------------------------------- 745 746.. function:: void arb_mat_get_mid(arb_mat_t B, const arb_mat_t A) 747 748 Sets the entries of *B* to the exact midpoints of the entries of *A*. 749 750.. function:: void arb_mat_add_error_mag(arb_mat_t mat, const mag_t err) 751 752 Adds *err* in-place to the radii of the entries of *mat*. 753 754Eigenvalues and eigenvectors 755------------------------------------------------------------------------------- 756 757To compute eigenvalues and eigenvectors, one can convert to an 758:type:`acb_mat_t` and use the functions in :ref:`acb_mat.h: Eigenvalues and eigenvectors<acb-mat-eigenvalues>`. 759In the future dedicated methods for real matrices will be added here. 760