Commit 4ec5f73fa4d5cb73caf7f101750cf0a5ce3a4255 - flint-arb

New upstream version 2.16.0 Julien Puydt 5 years ago

59 changed file(s) with 4757 addition(s) and 492 deletion(s). Raw diff Collapse all Expand all

+11

-6

README.md less more

71	71
72	72	For more example programs, see: http://arblib.org/examples.html
73	73
74		## General features
	74	## Features
75	75
76	76	Besides basic arithmetic, Arb allows working with univariate
77	77	polynomials, truncated power series, and matrices
78	78	over both real and complex numbers.
79	79
80	80	Basic linear algebra is supported, including matrix multiplication,
81		determinant, inverse, nonsingular solving and matrix exponential.
	81	determinant, inverse, nonsingular solving, matrix exponential,
	82	and computation of eigenvalues and eigenvectors.
82	83
83		Support for polynomial and power series is quite extensive,
	84	Support for polynomials and power series is quite extensive,
84	85	including methods for composition, reversion, product trees,
85	86	multipoint evaluation and interpolation, complex root isolation,
86	87	and transcendental functions of power series.

109	110	Various low-level optimizations have also been implemented
110	111	to reduce overhead at precisions of just a few machine
111	112	words. Most operations on polynomials and power series
112		use asymptotically fast FFT multiplication.
	113	use asymptotically fast FFT multiplication based on FLINT.
	114	Similarly, most operations on large matrices take advantage
	115	of the fast integer matrix multiplication in FLINT.
113	116
114	117	For basic arithmetic, Arb should generally be around as fast
115	118	as MPFR (http://mpfr.org), though it can be a bit slower

153	156	Julia interface to Arb. The Nemo installation script will
154	157	create a local installation of Arb along with other dependencies.
155	158
156		An experimental standalone Python interface to FLINT and Arb is also available
	159	A standalone Python interface to FLINT and Arb is also available
157	160	(<https://github.com/fredrik-johansson/python-flint>).
158	161
159	162	A separate wrapper of transcendental functions for use with the

161	164	(<https://github.com/fredrik-johansson/arbcmath>).
162	165
163	166	Other third-party wrappers include:
	167	* A Julia interface: https://github.com/JeffreySarnoff/ArbNumerics.jl
	168	* Another Julia interface: https://github.com/JuliaArbTypes/ArbFloats.jl
164	169	* Java wrapper using JNA: https://github.com/crowlogic/arb/
165		* Another Julia interface: https://github.com/JuliaArbTypes/ArbFloats.jl
	170

-0

acb.h less more

65	65	acb_get_imag(arb_t im, const acb_t z)
66	66	{
67	67	arb_set(im, acb_imagref(z));
	68	}
	69
	70	ACB_INLINE void
	71	acb_get_mid(acb_t res, const acb_t x)
	72	{
	73	arb_get_mid_arb(acb_realref(res), acb_realref(x));
	74	arb_get_mid_arb(acb_imagref(res), acb_imagref(x));
68	75	}
69	76
70	77	ACB_INLINE int

+1021

-0

acb_mat/approx_eig_qr.c less more

	0	/*
	1	Copyright 2013 Timo Hartmann
	2	Copyright 2018 Fredrik Johansson
	3
	4	This file is part of Arb.
	5
	6	Arb is free software: you can redistribute it and/or modify it under
	7	the terms of the GNU Lesser General Public License (LGPL) as published
	8	by the Free Software Foundation; either version 2.1 of the License, or
	9	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	10	*/
	11
	12	/*
	13	Adapted from eigen.py in mpmath (written by Timo Hartmann)
	14
	15	Todo items present in the original code:
	16
	17	- Implement balancing
	18	- Aggressive early deflation
	19	*/
	20
	21	#include "acb_mat.h"
	22
	23	static void
	24	acb_approx_mag(mag_t res, const acb_t x)
	25	{
	26	mag_t t;
	27	mag_init(t);
	28	arf_get_mag(res, arb_midref(acb_realref(x)));
	29	arf_get_mag(t, arb_midref(acb_imagref(x)));
	30	mag_hypot(res, res, t);
	31	mag_clear(t);
	32	}
	33
	34	static void
	35	acb_approx_mul(acb_t res, const acb_t x, const acb_t y, slong prec)
	36	{
	37	arf_complex_mul(arb_midref(acb_realref(res)), arb_midref(acb_imagref(res)),
	38	arb_midref(acb_realref(x)), arb_midref(acb_imagref(x)),
	39	arb_midref(acb_realref(y)), arb_midref(acb_imagref(y)), prec, ARF_RND_DOWN);
	40	}
	41
	42	static void
	43	acb_approx_add(acb_t res, const acb_t x, const acb_t y, slong prec)
	44	{
	45	arf_add(arb_midref(acb_realref(res)), arb_midref(acb_realref(x)), arb_midref(acb_realref(y)), prec, ARF_RND_DOWN);
	46	arf_add(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)), arb_midref(acb_imagref(y)), prec, ARF_RND_DOWN);
	47	}
	48
	49	static void
	50	acb_approx_sub(acb_t res, const acb_t x, const acb_t y, slong prec)
	51	{
	52	arf_sub(arb_midref(acb_realref(res)), arb_midref(acb_realref(x)), arb_midref(acb_realref(y)), prec, ARF_RND_DOWN);
	53	arf_sub(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)), arb_midref(acb_imagref(y)), prec, ARF_RND_DOWN);
	54	}
	55
	56	static void
	57	acb_approx_set(acb_t res, const acb_t x)
	58	{
	59	arf_set(arb_midref(acb_realref(res)), arb_midref(acb_realref(x)));
	60	arf_set(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)));
	61	}
	62
	63	static void
	64	acb_approx_div_arb(acb_t res, const acb_t x, const arb_t y, slong prec)
	65	{
	66	arf_div(arb_midref(acb_realref(res)), arb_midref(acb_realref(x)), arb_midref(y), prec, ARF_RND_DOWN);
	67	arf_div(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)), arb_midref(y), prec, ARF_RND_DOWN);
	68	}
	69
	70	static void
	71	acb_approx_inv(acb_t z, const acb_t x, slong prec)
	72	{
	73	arf_set(arb_midref(acb_realref(z)), arb_midref(acb_realref(x)));
	74	arf_set(arb_midref(acb_imagref(z)), arb_midref(acb_imagref(x)));
	75
	76	mag_zero(arb_radref(acb_realref(z)));
	77	mag_zero(arb_radref(acb_imagref(z)));
	78
	79	acb_inv(z, z, prec);
	80
	81	mag_zero(arb_radref(acb_realref(z)));
	82	mag_zero(arb_radref(acb_imagref(z)));
	83	}
	84
	85	static void
	86	acb_approx_div(acb_t z, const acb_t x, const acb_t y, slong prec)
	87	{
	88	acb_t t;
	89	acb_init(t);
	90	acb_approx_inv(t, y, prec);
	91	acb_approx_mul(z, x, t, prec);
	92	acb_clear(t);
	93	}
	94
	95	void
	96	acb_mat_approx_qr_step(acb_mat_t A, acb_mat_t Q, slong n0, slong n1, const acb_t shift, slong prec)
	97	{
	98	slong j, k, n;
	99	acb_t c, s, negs, cc, cs, negcs, t;
	100	acb_struct v1[2];
	101	acb_struct v1neg[2];
	102	acb_struct v2[2];
	103	acb_struct v2neg[2];
	104	acb_struct v3[2];
	105	arb_t v, u;
	106
	107	n = acb_mat_nrows(A);
	108
	109	acb_init(c);
	110	acb_init(s);
	111	acb_init(negs);
	112	acb_init(cc);
	113	acb_init(cs);
	114	acb_init(negcs);
	115	acb_init(t);
	116	arb_init(v);
	117	arb_init(u);
	118
	119	/* Calculate Givens rotation */
	120	acb_approx_sub(c, acb_mat_entry(A, n0, n0), shift, prec);
	121	acb_approx_set(s, acb_mat_entry(A, n0 + 1, n0));
	122
	123	arf_sosq(arb_midref(v), arb_midref(acb_realref(c)), arb_midref(acb_imagref(c)), prec, ARF_RND_DOWN);
	124	arf_sosq(arb_midref(u), arb_midref(acb_realref(s)), arb_midref(acb_imagref(s)), prec, ARF_RND_DOWN);
	125	arf_add(arb_midref(v), arb_midref(v), arb_midref(u), prec, ARF_RND_DOWN);
	126	arf_sqrt(arb_midref(v), arb_midref(v), prec, ARF_RND_DOWN);
	127
	128	if (arb_is_zero(v))
	129	{
	130	arb_one(v);
	131	acb_one(c);
	132	acb_zero(s);
	133	}
	134	else
	135	{
	136	acb_approx_div_arb(c, c, v, prec);
	137	acb_approx_div_arb(s, s, v, prec);
	138	}
	139
	140	acb_conj(cc, c);
	141	acb_conj(cs, s);
	142	acb_neg(negs, s);
	143	acb_neg(negcs, cs);
	144
	145	v1[0] = *c;
	146	v1[1] = *s;
	147	v1neg[0] = *c;
	148	v1neg[1] = *negs;
	149	v2[0] = *cc;
	150	v2[1] = *cs;
	151	v2neg[0] = *cc;
	152	v2neg[1] = *negcs;
	153
	154	/* Apply Givens rotation from the left */
	155	for (k = n0; k < n; k++)
	156	{
	157	v3[0] = *acb_mat_entry(A, n0, k);
	158	v3[1] = *acb_mat_entry(A, n0 + 1, k);
	159
	160	/* x = A[n0 ,k] */
	161	/* y = A[n0+1,k] */
	162	/* A[n0, k] = cc * x + cs * y */
	163	/* A[n0 + 1, k] = c * y - s * x */
	164
	165	acb_approx_dot(t, NULL, 0, v2, 1, v3, 1, 2, prec);
	166	acb_approx_dot(acb_mat_entry(A, n0 + 1, k), NULL, 0, v1neg, 1, v3 + 1, -1, 2, prec);
	167	acb_swap(t, acb_mat_entry(A, n0, k));
	168	}
	169
	170	/* Apply Givens rotation from the right */
	171	for (k = 0; k < FLINT_MIN(n1, n0 + 3); k++)
	172	{
	173	/* x = A[k,n0 ] */
	174	/* y = A[k,n0+1] */
	175	/* A[k,n0 ] = c * x + s * y */
	176	/* A[k,n0+1] = cc * y - cs * x */
	177
	178	v3[0] = *acb_mat_entry(A, k, n0);
	179	v3[1] = *acb_mat_entry(A, k, n0 + 1);
	180
	181	acb_approx_dot(t, NULL, 0, v1, 1, v3, 1, 2, prec);
	182	acb_approx_dot(acb_mat_entry(A, k, n0 + 1), NULL, 0, v2neg, 1, v3 + 1, -1, 2, prec);
	183	acb_swap(t, acb_mat_entry(A, k, n0));
	184	}
	185
	186	if (Q != NULL)
	187	{
	188	for (k = 0; k < n; k++)
	189	{
	190	/* x = Q[k,n0 ] */
	191	/* y = Q[k,n0+1] */
	192	/* Q[k,n0 ] = c * x + s * y */
	193	/* Q[k,n0+1] = cc * y - cs * x */
	194
	195	v3[0] = *acb_mat_entry(Q, k, n0);
	196	v3[1] = *acb_mat_entry(Q, k, n0 + 1);
	197
	198	acb_approx_dot(t, NULL, 0, v1, 1, v3, 1, 2, prec);
	199	acb_approx_dot(acb_mat_entry(Q, k, n0 + 1), NULL, 0, v2neg, 1, v3 + 1, -1, 2, prec);
	200	acb_swap(t, acb_mat_entry(Q, k, n0));
	201	}
	202	}
	203
	204	for (j = n0; j < n1 - 2; j++)
	205	{
	206	/* Calculate Givens rotation */
	207	acb_set(c, acb_mat_entry(A, j + 1, j));
	208	acb_set(s, acb_mat_entry(A, j + 2, j));
	209
	210	arf_sosq(arb_midref(v), arb_midref(acb_realref(c)), arb_midref(acb_imagref(c)), prec, ARF_RND_DOWN);
	211	arf_sosq(arb_midref(u), arb_midref(acb_realref(s)), arb_midref(acb_imagref(s)), prec, ARF_RND_DOWN);
	212	arf_add(arb_midref(v), arb_midref(v), arb_midref(u), prec, ARF_RND_DOWN);
	213	arf_sqrt(arb_midref(v), arb_midref(v), prec, ARF_RND_DOWN);
	214
	215	if (arb_is_zero(v))
	216	{
	217	acb_zero(acb_mat_entry(A, j + 1, j));
	218	arb_one(v);
	219	acb_one(c);
	220	acb_zero(s);
	221	}
	222	else
	223	{
	224	acb_set_arb(acb_mat_entry(A, j + 1, j), v);
	225	acb_approx_div_arb(c, c, v, prec);
	226	acb_approx_div_arb(s, s, v, prec);
	227	}
	228
	229	acb_zero(acb_mat_entry(A, j + 2, j));
	230
	231	acb_conj(cc, c);
	232	acb_conj(cs, s);
	233	acb_neg(negs, s);
	234	acb_neg(negcs, cs);
	235
	236	v1[0] = *c;
	237	v1[1] = *s;
	238	v1neg[0] = *c;
	239	v1neg[1] = *negs;
	240	v2[0] = *cc;
	241	v2[1] = *cs;
	242	v2neg[0] = *cc;
	243	v2neg[1] = *negcs;
	244
	245	/* Apply Givens rotation from the left */
	246	for (k = j + 1; k < n; k++)
	247	{
	248	v3[0] = *acb_mat_entry(A, j + 1, k);
	249	v3[1] = *acb_mat_entry(A, j + 2, k);
	250
	251	/* x = A[j+1, k] */
	252	/* y = A[j+2, k] */
	253	/* A[j + 1, k] = cc * x + cs * y */
	254	/* A[j + 2, k] = c * y - s * x */
	255
	256	acb_approx_dot(t, NULL, 0, v2, 1, v3, 1, 2, prec);
	257	acb_approx_dot(acb_mat_entry(A, j + 2, k), NULL, 0, v1neg, 1, v3 + 1, -1, 2, prec);
	258	acb_swap(t, acb_mat_entry(A, j + 1, k));
	259	}
	260
	261	/* Apply Givens rotation from the right */
	262	for (k = 0; k < FLINT_MIN(n1, j + 4); k++)
	263	{
	264	/* x = A[k,j+1] */
	265	/* y = A[k,j+2] */
	266	/* A[k,j+1] = c * x + s * y */
	267	/* A[k,j+2] = cc * y - cs * x */
	268
	269	v3[0] = *acb_mat_entry(A, k, j + 1);
	270	v3[1] = *acb_mat_entry(A, k, j + 2);
	271
	272	acb_approx_dot(t, NULL, 0, v1, 1, v3, 1, 2, prec);
	273	acb_approx_dot(acb_mat_entry(A, k, j + 2), NULL, 0, v2neg, 1, v3 + 1, -1, 2, prec);
	274	acb_swap(t, acb_mat_entry(A, k, j + 1));
	275	}
	276
	277	if (Q != NULL)
	278	{
	279	for (k = 0; k < n; k++)
	280	{
	281	/* x = Q[k,j+1] */
	282	/* y = Q[k,j+2] */
	283	/* Q[k,j+1] = c * x + s * y */
	284	/* Q[k,j+2] = cc * y - cs * x */
	285
	286	v3[0] = *acb_mat_entry(Q, k, j + 1);
	287	v3[1] = *acb_mat_entry(Q, k, j + 2);
	288
	289	acb_approx_dot(t, NULL, 0, v1, 1, v3, 1, 2, prec);
	290	acb_approx_dot(acb_mat_entry(Q, k, j + 2), NULL, 0, v2neg, 1, v3 + 1, -1, 2, prec);
	291	acb_swap(t, acb_mat_entry(Q, k, j + 1));
	292	}
	293	}
	294	}
	295
	296	acb_clear(c);
	297	acb_clear(s);
	298	acb_clear(negs);
	299	acb_clear(cc);
	300	acb_clear(cs);
	301	acb_clear(negcs);
	302	acb_clear(t);
	303	arb_clear(v);
	304	arb_clear(u);
	305	}
	306
	307	void
	308	acb_mat_approx_hessenberg_reduce_0(acb_mat_t A, acb_ptr T, slong prec)
	309	{
	310	slong i, j, k, n;
	311	arf_t scale, scale_inv, tt, H, G, f;
	312	acb_ptr V1, V2;
	313	acb_t ff, GG, TT;
	314	acb_t F;
	315
	316	n = acb_mat_nrows(A);
	317	if (n <= 2)
	318	return;
	319
	320	arf_init(scale);
	321	arf_init(scale_inv);
	322	arf_init(tt);
	323	arf_init(H);
	324	arf_init(G);
	325	arf_init(f);
	326	acb_init(F);
	327	V1 = _acb_vec_init(n + 1);
	328	V2 = _acb_vec_init(n + 1);
	329	acb_init(ff);
	330	acb_init(GG);
	331	acb_init(TT);
	332
	333	for (i = n - 1; i >= 2; i--)
	334	{
	335	/* Scale the vector (todo: is this needed?) */
	336	arf_zero(scale);
	337	for (k = 0; k < i; k++)
	338	{
	339	arf_abs(tt, arb_midref(acb_realref(acb_mat_entry(A, i, k))));
	340	arf_add(scale, scale, tt, prec, ARF_RND_DOWN);
	341	arf_abs(tt, arb_midref(acb_imagref(acb_mat_entry(A, i, k))));
	342	arf_add(scale, scale, tt, prec, ARF_RND_DOWN);
	343	}
	344
	345	arf_ui_div(scale_inv, 1, scale, prec, ARF_RND_DOWN);
	346
	347	if (arf_is_zero(scale))
	348	{
	349	acb_zero(T + i);
	350	acb_zero(acb_mat_entry(A, i, i - 1));
	351	continue;
	352	}
	353
	354	/* Calculate parameters for Householder transformation. */
	355	arf_zero(H);
	356	for (k = 0; k < i; k++)
	357	{
	358	arf_ptr Aikr, Aiki;
	359
	360	Aikr = arb_midref(acb_realref(acb_mat_entry(A, i, k)));
	361	Aiki = arb_midref(acb_imagref(acb_mat_entry(A, i, k)));
	362
	363	arf_mul(Aikr, Aikr, scale_inv, prec, ARF_RND_DOWN);
	364	arf_mul(Aiki, Aiki, scale_inv, prec, ARF_RND_DOWN);
	365
	366	arf_addmul(H, Aikr, Aikr, prec, ARF_RND_DOWN);
	367	arf_addmul(H, Aiki, Aiki, prec, ARF_RND_DOWN);
	368	}
	369
	370	acb_set(F, acb_mat_entry(A, i, i - 1));
	371
	372	/* f = abs(F) */
	373	arf_mul(f, arb_midref(acb_realref(F)), arb_midref(acb_realref(F)), prec, ARF_RND_DOWN);
	374	arf_addmul(f, arb_midref(acb_imagref(F)), arb_midref(acb_imagref(F)), prec, ARF_RND_DOWN);
	375	arf_sqrt(f, f, prec, ARF_RND_DOWN);
	376
	377	arf_sqrt(G, H, prec, ARF_RND_DOWN);
	378
	379	/* A[i,i-1] = -G scale */
	380	arf_mul(arb_midref(acb_realref(acb_mat_entry(A, i, i - 1))), G, scale, prec, ARF_RND_DOWN);
	381	arf_neg(arb_midref(acb_realref(acb_mat_entry(A, i, i - 1))), arb_midref(acb_realref(acb_mat_entry(A, i, i - 1))));
	382	arf_zero(arb_midref(acb_imagref(acb_mat_entry(A, i, i - 1))));
	383
	384	if (arf_is_zero(f))
	385	{
	386	arb_set_arf(acb_realref(T + i), G);
	387	arb_zero(acb_imagref(T + i));
	388	}
	389	else
	390	{
	391	/* ff = F / f */
	392	arf_div(arb_midref(acb_realref(ff)), arb_midref(acb_realref(F)), f, prec, ARF_RND_DOWN);
	393	arf_div(arb_midref(acb_imagref(ff)), arb_midref(acb_imagref(F)), f, prec, ARF_RND_DOWN);
	394
	395	/* T[i] = F + G ff */
	396	acb_set(T + i, F);
	397	arf_addmul(arb_midref(acb_realref(T + i)), arb_midref(acb_realref(ff)), G, prec, ARF_RND_DOWN);
	398	arf_addmul(arb_midref(acb_imagref(T + i)), arb_midref(acb_imagref(ff)), G, prec, ARF_RND_DOWN);
	399
	400	/* A[i,i-1] = ff /
	401	acb_approx_mul(acb_mat_entry(A, i, i - 1), acb_mat_entry(A, i, i - 1), ff, prec);
	402	}
	403
	404	arf_addmul(H, G, f, prec, ARF_RND_DOWN);
	405	arf_rsqrt(H, H, prec, ARF_RND_DOWN);
	406
	407	arf_mul(arb_midref(acb_realref(T + i)), arb_midref(acb_realref(T + i)), H, prec, ARF_RND_DOWN);
	408	arf_mul(arb_midref(acb_imagref(T + i)), arb_midref(acb_imagref(T + i)), H, prec, ARF_RND_DOWN);
	409
	410	for (k = 0; k < i - 1; k++)
	411	{
	412	arf_mul(arb_midref(acb_realref(acb_mat_entry(A, i, k))),
	413	arb_midref(acb_realref(acb_mat_entry(A, i, k))),
	414	H, prec, ARF_RND_DOWN);
	415	arf_mul(arb_midref(acb_imagref(acb_mat_entry(A, i, k))),
	416	arb_midref(acb_imagref(acb_mat_entry(A, i, k))),
	417	H, prec, ARF_RND_DOWN);
	418	}
	419
	420	/* todo: optimize copies below */
	421	/* todo: conj mid etc... */
	422
	423	/* Apply Householder transformation (from the right). */
	424	for (j = 0; j < i; j++)
	425	{
	426	acb_conj(V1, T + i);
	427	acb_set(V2, acb_mat_entry(A, j, i - 1));
	428	for (k = 0; k < i - 1; k++)
	429	{
	430	acb_conj(V1 + k + 1, acb_mat_entry(A, i, k));
	431	acb_set(V2 + k + 1, acb_mat_entry(A, j, k));
	432	}
	433	acb_approx_dot(GG, NULL, 0, V1, 1, V2, 1, i, prec);
	434
	435	acb_approx_mul(TT, GG, T + i, prec);
	436	acb_approx_sub(acb_mat_entry(A, j, i - 1),
	437	acb_mat_entry(A, j, i - 1), TT, prec);
	438	for (k = 0; k < i - 1; k++)
	439	{
	440	acb_approx_mul(TT, GG, acb_mat_entry(A, i, k), prec);
	441	acb_approx_sub(acb_mat_entry(A, j, k),
	442	acb_mat_entry(A, j, k), TT, prec);
	443	}
	444	}
	445
	446	for (j = 0; j < n; j++)
	447	{
	448	acb_set(V1, T + i);
	449	acb_set(V2, acb_mat_entry(A, i - 1, j));
	450	for (k = 0; k < i - 1; k++)
	451	{
	452	acb_set(V1 + k + 1, acb_mat_entry(A, i, k));
	453	acb_set(V2 + k + 1, acb_mat_entry(A, k, j));
	454	}
	455	acb_approx_dot(GG, NULL, 0, V1, 1, V2, 1, i, prec);
	456
	457	acb_conj(TT, T + i);
	458	acb_approx_mul(TT, GG, TT, prec);
	459	acb_approx_sub(acb_mat_entry(A, i - 1, j),
	460	acb_mat_entry(A, i - 1, j), TT, prec);
	461	for (k = 0; k < i - 1; k++)
	462	{
	463	acb_conj(TT, acb_mat_entry(A, i, k));
	464	acb_approx_mul(TT, GG, TT, prec);
	465	acb_approx_sub(acb_mat_entry(A, k, j),
	466	acb_mat_entry(A, k, j), TT, prec);
	467	}
	468	}
	469	}
	470
	471	arf_clear(scale);
	472	arf_clear(scale_inv);
	473	arf_clear(tt);
	474	arf_clear(H);
	475	arf_clear(G);
	476	arf_clear(f);
	477	acb_clear(F);
	478	_acb_vec_clear(V1, n + 1);
	479	_acb_vec_clear(V2, n + 1);
	480	acb_clear(ff);
	481	acb_clear(GG);
	482	acb_clear(TT);
	483	}
	484
	485	void
	486	acb_mat_approx_hessenberg_reduce_1(acb_mat_t A, acb_srcptr T, slong prec)
	487	{
	488	slong i, j, k, n;
	489	acb_t G, t;
	490
	491	n = acb_mat_nrows(A);
	492
	493	if (n <= 1)
	494	{
	495	if (n == 1)
	496	acb_one(acb_mat_entry(A, 0, 0));
	497	return;
	498	}
	499
	500	acb_one(acb_mat_entry(A, 0, 0));
	501	acb_one(acb_mat_entry(A, 1, 1));
	502	acb_zero(acb_mat_entry(A, 0, 1));
	503	acb_zero(acb_mat_entry(A, 1, 0));
	504
	505	acb_init(G);
	506	acb_init(t);
	507
	508	for (i = 2; i < n; i++)
	509	{
	510	if (!acb_is_zero(T + i))
	511	{
	512	/* todo: rewrite using approx_dot */
	513	for (j = 0; j < i; j++)
	514	{
	515	acb_approx_mul(G, T + i, acb_mat_entry(A, i - 1, j), prec);
	516	for (k = 0; k < i - 1; k++)
	517	{
	518	acb_approx_mul(t, acb_mat_entry(A, i, k), acb_mat_entry(A, k, j), prec);
	519	acb_approx_add(G, G, t, prec);
	520	}
	521
	522	acb_conj(t, T + i);
	523	acb_approx_mul(t, G, t, prec);
	524	acb_approx_sub(acb_mat_entry(A, i - 1, j), acb_mat_entry(A, i - 1, j), t, prec);
	525	for (k = 0; k < i - 1; k++)
	526	{
	527	acb_conj(t, acb_mat_entry(A, i, k));
	528	acb_approx_mul(t, G, t, prec);
	529	acb_approx_sub(acb_mat_entry(A, k, j), acb_mat_entry(A, k, j), t, prec);
	530	}
	531	}
	532	}
	533
	534	acb_one(acb_mat_entry(A, i, i));
	535	for (j = 0; j < i; j++)
	536	{
	537	acb_zero(acb_mat_entry(A, j, i));
	538	acb_zero(acb_mat_entry(A, i, j));
	539	}
	540	}
	541
	542	acb_clear(G);
	543	acb_clear(t);
	544	}
	545
	546	/* Right eigenvectors of a triu matrix. No aliasing. */
	547	void
	548	acb_mat_approx_eig_triu_r(acb_mat_t ER, const acb_mat_t A, slong prec)
	549	{
	550	slong i, j, k, n;
	551	mag_t tm, smin, unfl, simin, smlnum, rmax;
	552	acb_t r, s, t;
	553
	554	n = acb_mat_nrows(A);
	555
	556	acb_mat_one(ER);
	557
	558	acb_init(r);
	559	acb_init(s);
	560	acb_init(t);
	561	mag_init(tm);
	562	mag_init(smin);
	563	mag_init(smlnum);
	564	mag_init(unfl);
	565	mag_init(simin);
	566	mag_init(rmax);
	567
	568	mag_set_ui_2exp_si(unfl, 1, -30 * prec);
	569	mag_mul_ui(smlnum, unfl, n);
	570	mag_mul_2exp_si(smlnum, smlnum, prec);
	571	mag_set_ui_2exp_si(simin, 1, prec / 2);
	572	mag_one(rmax);
	573
	574	for (i = 1; i < n; i++)
	575	{
	576	acb_set(s, acb_mat_entry(A, i, i));
	577
	578	/* smin = max(eps * abs(s), smlnum) */
	579	acb_approx_mag(smin, s);
	580	mag_mul_2exp_si(smin, smin, -prec);
	581	mag_max(smin, smin, smlnum);
	582
	583	for (j = i - 1; j >= 0; j--)
	584	{
	585	acb_approx_dot(r, NULL, 0, A->rows[j] + j + 1, 1, ER->rows[i] + j + 1, 1, i - j, prec);
	586	acb_approx_sub(t, acb_mat_entry(A, j, j), s, prec);
	587
	588	/* if abs(t) < smin: t = smin */
	589	acb_approx_mag(tm, t);
	590	if (mag_cmp(tm, smin) < 0)
	591	{
	592	acb_zero(t);
	593	arf_set_mag(arb_midref(acb_realref(t)), smin);
	594	}
	595
	596	acb_approx_div(acb_mat_entry(ER, i, j), r, t, prec);
	597	acb_neg(acb_mat_entry(ER, i, j), acb_mat_entry(ER, i, j));
	598
	599	acb_approx_mag(tm, r);
	600	mag_max(rmax, rmax, tm);
	601	if (mag_cmp(rmax, simin) > 0)
	602	{
	603	arb_t b;
	604	arb_init(b);
	605	arf_set_mag(arb_midref(b), rmax);
	606
	607	for (k = j; k < i + 1; k++)
	608	{
	609	acb_approx_div_arb(acb_mat_entry(ER, i, k),
	610	acb_mat_entry(ER, i, k), b, prec);
	611	}
	612
	613	mag_one(rmax);
	614	arb_clear(b);
	615	}
	616	}
	617
	618	if (mag_cmp_2exp_si(rmax, 0) != 0)
	619	{
	620	arb_t b;
	621	arb_init(b);
	622	arf_set_mag(arb_midref(b), rmax);
	623
	624	for (k = 0; k < i + 1; k++)
	625	{
	626	acb_approx_div_arb(acb_mat_entry(ER, i, k),
	627	acb_mat_entry(ER, i, k), b, prec);
	628	}
	629
	630	arb_clear(b);
	631	}
	632	}
	633
	634	acb_mat_transpose(ER, ER);
	635
	636	acb_clear(r);
	637	acb_clear(s);
	638	acb_clear(t);
	639	mag_clear(tm);
	640	mag_clear(smin);
	641	mag_clear(smlnum);
	642	mag_clear(unfl);
	643	mag_clear(simin);
	644	mag_clear(rmax);
	645	}
	646
	647	/* Left eigenvectors of a triu matrix. No aliasing. */
	648	void
	649	acb_mat_approx_eig_triu_l(acb_mat_t EL, const acb_mat_t A, slong prec)
	650	{
	651	slong i, j, k, n;
	652	mag_t tm, smin, unfl, simin, smlnum, rmax;
	653	acb_t r, s, t;
	654	acb_mat_t AT;
	655
	656	n = acb_mat_nrows(A);
	657	acb_mat_init(AT, n, n);
	658
	659	acb_mat_one(EL);
	660	acb_mat_transpose(AT, A);
	661
	662	acb_init(r);
	663	acb_init(s);
	664	acb_init(t);
	665	mag_init(tm);
	666	mag_init(smin);
	667	mag_init(smlnum);
	668	mag_init(unfl);
	669	mag_init(simin);
	670	mag_init(rmax);
	671
	672	mag_set_ui_2exp_si(unfl, 1, -30 * prec);
	673	mag_mul_ui(smlnum, unfl, n);
	674	mag_mul_2exp_si(smlnum, smlnum, prec);
	675	mag_set_ui_2exp_si(simin, 1, prec / 2);
	676	mag_one(rmax);
	677
	678	for (i = 0; i < n - 1; i++)
	679	{
	680	acb_set(s, acb_mat_entry(AT, i, i));
	681
	682	/* smin = max(eps * abs(s), smlnum) */
	683	acb_approx_mag(smin, s);
	684	mag_mul_2exp_si(smin, smin, -prec);
	685	mag_max(smin, smin, smlnum);
	686
	687	for (j = i + 1; j < n; j++)
	688	{
	689	acb_approx_dot(r, NULL, 0, EL->rows[i] + i, 1, AT->rows[j] + i, 1, j - i, prec);
	690	acb_approx_sub(t, acb_mat_entry(AT, j, j), s, prec);
	691
	692	/* if abs(t) < smin: t = smin */
	693	acb_approx_mag(tm, t);
	694	if (mag_cmp(tm, smin) < 0)
	695	{
	696	acb_zero(t);
	697	arf_set_mag(arb_midref(acb_realref(t)), smin);
	698	}
	699
	700	acb_approx_div(acb_mat_entry(EL, i, j), r, t, prec);
	701	acb_neg(acb_mat_entry(EL, i, j), acb_mat_entry(EL, i, j));
	702
	703	acb_approx_mag(tm, r);
	704	mag_max(rmax, rmax, tm);
	705	if (mag_cmp(rmax, simin) > 0)
	706	{
	707	arb_t b;
	708	arb_init(b);
	709	arf_set_mag(arb_midref(b), rmax);
	710
	711	for (k = i; k < j + 1; k++)
	712	{
	713	acb_approx_div_arb(acb_mat_entry(EL, i, k),
	714	acb_mat_entry(EL, i, k), b, prec);
	715	}
	716
	717	mag_one(rmax);
	718	arb_clear(b);
	719	}
	720	}
	721
	722	if (mag_cmp_2exp_si(rmax, 0) != 0)
	723	{
	724	arb_t b;
	725	arb_init(b);
	726	arf_set_mag(arb_midref(b), rmax);
	727
	728	for (k = i; k < n; k++)
	729	{
	730	acb_approx_div_arb(acb_mat_entry(EL, i, k),
	731	acb_mat_entry(EL, i, k), b, prec);
	732	}
	733
	734	arb_clear(b);
	735	}
	736	}
	737
	738	acb_clear(r);
	739	acb_clear(s);
	740	acb_clear(t);
	741	mag_clear(tm);
	742	mag_clear(smin);
	743	mag_clear(smlnum);
	744	mag_clear(unfl);
	745	mag_clear(simin);
	746	mag_clear(rmax);
	747	}
	748
	749	int
	750	acb_mat_approx_hessenberg_qr(acb_mat_t A, acb_mat_t Q, const mag_t tol, slong maxiter, slong prec)
	751	{
	752	slong n, i, j, k, n0, n1, iter, total_iter;
	753	mag_t norm, tm, um, eps, ts;
	754	acb_t shift, s, t, a, b;
	755	int result;
	756
	757	n = acb_mat_nrows(A);
	758
	759	if (n <= 1)
	760	return 1;
	761
	762	mag_init(norm);
	763	mag_init(tm);
	764	mag_init(um);
	765	mag_init(eps);
	766	mag_init(ts);
	767	acb_init(shift);
	768	acb_init(s);
	769	acb_init(t);
	770	acb_init(a);
	771	acb_init(b);
	772
	773	for (i = 0; i < n; i++)
	774	{
	775	for (j = 0; j < FLINT_MIN(i + 2, n); j++)
	776	{
	777	arf_get_mag(tm, arb_midref(acb_realref(acb_mat_entry(A, j, i))));
	778	arf_get_mag(um, arb_midref(acb_imagref(acb_mat_entry(A, j, i))));
	779	mag_addmul(norm, tm, tm);
	780	mag_addmul(norm, um, um);
	781	}
	782	}
	783
	784	mag_sqrt(norm, norm);
	785	mag_div_ui(norm, norm, n);
	786
	787	if (mag_is_zero(norm))
	788	return 1;
	789
	790	if (mag_is_inf(norm))
	791	return 0;
	792
	793	if (tol == NULL)
	794	{
	795	mag_one(eps);
	796	mag_mul_2exp_si(eps, eps, -prec);
	797	mag_div_ui(eps, eps, 100 * n);
	798	}
	799	else
	800	{
	801	mag_set(eps, tol);
	802	}
	803
	804	if (maxiter <= 0)
	805	{
	806	maxiter = 14 * n;
	807	maxiter += prec / 10;
	808	}
	809
	810	/* The active submatrix is A[n0:n1,n0:n1]. */
	811	n0 = 0;
	812	n1 = n;
	813
	814	iter = total_iter = 0;
	815	result = 0;
	816
	817	while (1)
	818	{
	819	k = n0;
	820
	821	/* flint_printf("total_iter %wd active %wd\n", total_iter, n1 - n0); */
	822
	823	while (k + 1 < n1)
	824	{
	825	mag_zero(ts);
	826	arf_get_mag(tm, arb_midref(acb_realref(acb_mat_entry(A, k, k))));
	827	mag_add(ts, ts, tm);
	828	arf_get_mag(tm, arb_midref(acb_imagref(acb_mat_entry(A, k, k))));
	829	mag_add(ts, ts, tm);
	830	arf_get_mag(tm, arb_midref(acb_realref(acb_mat_entry(A, k + 1, k + 1))));
	831	mag_add(ts, ts, tm);
	832	arf_get_mag(tm, arb_midref(acb_imagref(acb_mat_entry(A, k + 1, k + 1))));
	833	mag_add(ts, ts, tm);
	834
	835	/* if s < eps * norm, s = norm */
	836	mag_mul(tm, eps, norm);
	837	if (mag_cmp(ts, tm) < 0)
	838	mag_set(ts, norm);
	839
	840	/* if abs(A[k+1,k]) < eps * s, break */
	841	arf_get_mag(tm, arb_midref(acb_realref(acb_mat_entry(A, k + 1, k))));
	842	arf_get_mag(um, arb_midref(acb_imagref(acb_mat_entry(A, k + 1, k))));
	843	mag_hypot(tm, tm, um);
	844	mag_mul(um, eps, ts);
	845	if (mag_cmp(tm, um) < 0)
	846	break;
	847
	848	k++;
	849	}
	850
	851	/* Deflation found at position (k+1, k). */
	852	if (k + 1 < n1)
	853	{
	854	acb_zero(acb_mat_entry(A, k + 1, k));
	855	n0 = k + 1;
	856	iter = 0;
	857
	858	if (n0 + 1 >= n1)
	859	{
	860	/* Block of size at most two has converged. */
	861	n0 = 0;
	862	n1 = k + 1;
	863	if (n1 < 2)
	864	{
	865	/* QR algorithm has converged. */
	866	result = 1;
	867	goto cleanup;
	868	}
	869	}
	870	}
	871	else
	872	{
	873	if (iter % 30 == 10)
	874	{
	875	/* Exceptional shift */
	876	acb_set(shift, acb_mat_entry(A, n1 - 1, n1 - 2));
	877	}
	878	else if (iter % 30 == 20)
	879	{
	880	/* Exceptional shift */
	881	acb_abs(acb_realref(shift), acb_mat_entry(A, n1 - 1, n1 - 2), prec);
	882	arb_zero(acb_imagref(shift));
	883	}
	884	else if (iter % 30 == 29)
	885	{
	886	/* Exceptional shift */
	887	acb_zero(shift);
	888	arf_set_mag(arb_midref(acb_realref(shift)), norm);
	889	}
	890	else
	891	{
	892	acb_approx_add(t, acb_mat_entry(A, n1 - 2, n1 - 2), acb_mat_entry(A, n1 - 1, n1 - 1), prec);
	893	acb_approx_sub(a, acb_mat_entry(A, n1 - 1, n1 - 1), acb_mat_entry(A, n1 - 2, n1 - 2), prec);
	894	acb_approx_mul(a, a, a, prec);
	895	acb_approx_mul(b, acb_mat_entry(A, n1 - 1, n1 - 2), acb_mat_entry(A, n1 - 2, n1 - 1), prec);
	896	acb_mul_2exp_si(b, b, 2);
	897	acb_approx_add(s, a, b, prec);
	898
	899	if (arb_is_positive(acb_realref(s)))
	900	{
	901	acb_sqrt(s, s, prec);
	902	acb_get_mid(s, s);
	903	}
	904	else
	905	{
	906	acb_neg(s, s);
	907	acb_sqrt(s, s, prec);
	908	acb_get_mid(s, s);
	909	acb_mul_onei(s, s);
	910	}
	911
	912	acb_approx_add(a, t, s, prec);
	913	acb_approx_sub(b, t, s, prec);
	914	acb_mul_2exp_si(a, a, -1);
	915	acb_mul_2exp_si(b, b, -1);
	916
	917	acb_approx_sub(s, acb_mat_entry(A, n1 - 1, n1 - 1), a, prec);
	918	acb_approx_sub(t, acb_mat_entry(A, n1 - 1, n1 - 1), b, prec);
	919	acb_get_mag(tm, s);
	920	acb_get_mag(um, t);
	921
	922	if (mag_cmp(tm, um) > 0)
	923	acb_set(shift, b);
	924	else
	925	acb_set(shift, a);
	926	}
	927
	928	mag_zero(arb_radref(acb_realref(shift)));
	929	mag_zero(arb_radref(acb_imagref(shift)));
	930
	931	iter++;
	932	total_iter++;
	933
	934	acb_mat_approx_qr_step(A, Q, n0, n1, shift, prec);
	935
	936	if (iter > maxiter)
	937	{
	938	result = 0;
	939	goto cleanup;
	940	}
	941	}
	942	}
	943
	944	cleanup:
	945	mag_clear(norm);
	946	mag_clear(tm);
	947	mag_clear(um);
	948	mag_clear(ts);
	949	mag_clear(eps);
	950	acb_clear(shift);
	951	acb_clear(s);
	952	acb_clear(t);
	953	acb_clear(a);
	954	acb_clear(b);
	955
	956	return result;
	957	}
	958
	959	int
	960	acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)
	961	{
	962	slong n, i, j;
	963	acb_ptr T;
	964	acb_mat_t Acopy, Q;
	965	int result;
	966
	967	n = acb_mat_nrows(A);
	968
	969	T = _acb_vec_init(n);
	970	acb_mat_init(Acopy, n, n);
	971	acb_mat_get_mid(Acopy, A);
	972
	973	acb_mat_approx_hessenberg_reduce_0(Acopy, T, prec);
	974
	975	if (L != NULL \|\| R != NULL)
	976	{
	977	acb_mat_init(Q, n, n);
	978	acb_mat_set(Q, Acopy);
	979	acb_mat_approx_hessenberg_reduce_1(Q, T, prec);
	980	}
	981
	982	for (i = 0; i < n; i++)
	983	for (j = i + 2; j < n; j++)
	984	acb_zero(acb_mat_entry(Acopy, j, i));
	985
	986	result = acb_mat_approx_hessenberg_qr(Acopy,
	987	(L != NULL \|\| R != NULL) ? Q : NULL, tol, maxiter, prec);
	988
	989	for (i = 0; i < n; i++)
	990	acb_get_mid(E + i, acb_mat_entry(Acopy, i, i));
	991
	992	if (R != NULL)
	993	{
	994	acb_mat_t ER;
	995	acb_mat_init(ER, n, n);
	996	acb_mat_approx_eig_triu_r(ER, Acopy, prec);
	997	acb_mat_approx_mul(R, Q, ER, prec);
	998	acb_mat_get_mid(R, R); /* zero radii */
	999	acb_mat_clear(ER);
	1000	}
	1001
	1002	if (L != NULL)
	1003	{
	1004	acb_mat_t EL;
	1005	acb_mat_init(EL, n, n);
	1006	acb_mat_approx_eig_triu_l(EL, Acopy, prec);
	1007	acb_mat_conjugate_transpose(Q, Q);
	1008	acb_mat_approx_mul(L, EL, Q, prec);
	1009	acb_mat_get_mid(L, L); /* zero radii */
	1010	acb_mat_clear(EL);
	1011	}
	1012
	1013	if (L != NULL \|\| R != NULL)
	1014	acb_mat_clear(Q);
	1015
	1016	_acb_vec_clear(T, n);
	1017	acb_mat_clear(Acopy);
	1018
	1019	return result;
	1020	}

+30

-0

acb_mat/approx_inv.c less more

	0	/*
	1	Copyright (C) 2012 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_approx_inv(acb_mat_t X, const acb_mat_t A, slong prec)
	15	{
	16	if (X == A)
	17	{
	18	int r;
	19	acb_mat_t T;
	20	acb_mat_init(T, acb_mat_nrows(A), acb_mat_ncols(A));
	21	r = acb_mat_approx_inv(T, A, prec);
	22	acb_mat_swap(T, X);
	23	acb_mat_clear(T);
	24	return r;
	25	}
	26
	27	acb_mat_one(X);
	28	return acb_mat_approx_solve(X, A, X, prec);
	29	}

-14

acb_mat/approx_mul.c less more

83	83	}
84	84	}
85	85
86		int
87		acb_mat_is_exact(const acb_mat_t A)
88		{
89		slong i, j;
90
91		for (i = 0; i < acb_mat_nrows(A); i++)
92		for (j = 0; j < acb_mat_ncols(A); j++)
93		if (!mag_is_zero(arb_radref(acb_realref(acb_mat_entry(A, i, j)))) \|\|
94		!mag_is_zero(arb_radref(acb_imagref(acb_mat_entry(A, i, j)))))
95		return 0;
96
97		return 1;
98		}
99
100	86	void
101	87	acb_mat_approx_mul(acb_mat_t C, const acb_mat_t A, const acb_mat_t B, slong prec)
102	88	{

+49

-0

acb_mat/companion.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	void
	14	_acb_mat_companion(acb_mat_t A, acb_srcptr poly, slong prec)
	15	{
	16	slong i, j, n;
	17	acb_t c;
	18
	19	n = acb_mat_nrows(A);
	20
	21	if (n == 0)
	22	return;
	23
	24	for (i = 0; i < n - 1; i++)
	25	for (j = 0; j < n; j++)
	26	acb_set_ui(acb_mat_entry(A, i, j), (i + 1) == j);
	27
	28	acb_init(c);
	29	acb_inv(c, poly + n, prec);
	30	acb_neg(c, c);
	31	for (j = 0; j < n; j++)
	32	acb_mul(acb_mat_entry(A, n - 1, j), poly + j, c, prec);
	33	acb_clear(c);
	34	}
	35
	36	void
	37	acb_mat_companion(acb_mat_t A, const acb_poly_t poly, slong prec)
	38	{
	39	slong n = acb_mat_nrows(A);
	40
	41	if (n != acb_poly_degree(poly) \|\| n != acb_mat_ncols(A))
	42	{
	43	flint_printf("acb_mat_companion: incompatible dimensions!\n");
	44	flint_abort();
	45	}
	46
	47	_acb_mat_companion(A, poly->coeffs, prec);
	48	}

-75

acb_mat/det.c less more

38	38	acb_clear(a);
39	39	}
40	40
41		int
42		acb_mat_is_finite(const acb_mat_t A)
43		{
44		slong i, j, n, m;
45
46		n = acb_mat_nrows(A);
47		m = acb_mat_ncols(A);
48
49		for (i = 0; i < n; i++)
50		for (j = 0; j < m; j++)
51		if (!acb_is_finite(acb_mat_entry(A, i, j)))
52		return 0;
53
54		return 1;
55		}
56
57		int
58		acb_mat_is_triu(const acb_mat_t A)
59		{
60		slong i, j, n, m;
61
62		n = acb_mat_nrows(A);
63		m = acb_mat_ncols(A);
64
65		for (i = 1; i < n; i++)
66		for (j = 0; j < FLINT_MIN(i, m); j++)
67		if (!acb_is_zero(acb_mat_entry(A, i, j)))
68		return 0;
69
70		return 1;
71		}
72
73		int
74		acb_mat_is_tril(const acb_mat_t A)
75		{
76		slong i, j, n, m;
77
78		n = acb_mat_nrows(A);
79		m = acb_mat_ncols(A);
80
81		for (i = 0; i < n; i++)
82		for (j = i + 1; j < m; j++)
83		if (!acb_is_zero(acb_mat_entry(A, i, j)))
84		return 0;
85
86		return 1;
87		}
88
89		void
90		acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong a, slong b, slong prec)
91		{
92		if (b - a == 0)
93		acb_one(res);
94		else if (b - a == 1)
95		acb_set_round(res, acb_mat_entry(A, a, a), prec);
96		else if (b - a == 2)
97		acb_mul(res, acb_mat_entry(A, a, a), acb_mat_entry(A, a + 1, a + 1), prec);
98		else if (b - a == 3)
99		{
100		acb_mul(res, acb_mat_entry(A, a, a), acb_mat_entry(A, a + 1, a + 1), prec);
101		acb_mul(res, res, acb_mat_entry(A, a + 2, a + 2), prec);
102		}
103		else
104		{
105		acb_t t;
106		acb_init(t);
107		acb_mat_diag_prod(t, A, a, a + (b - a) / 2, prec);
108		acb_mat_diag_prod(res, A, a + (b - a) / 2, b, prec);
109		acb_mul(res, res, t, prec);
110		acb_clear(t);
111		}
112		}
113
114	41	void
115	42	acb_mat_det(acb_t det, const acb_mat_t A, slong prec)
116	43	{

142	69	}
143	70	else if (acb_mat_is_tril(A) \|\| acb_mat_is_triu(A))
144	71	{
145		acb_mat_diag_prod(det, A, 0, n, prec);
	72	acb_mat_diag_prod(det, A, prec);
146	73	}
147	74	else if (n == 3)
148	75	{

157	84	acb_mat_det_precond(det, A, prec);
158	85	}
159	86	}
160

-4

acb_mat/det_lu.c less more

78	78	}
79	79
80	80	void
81		acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong a, slong b, slong prec);
82
83		void
84	81	acb_mat_det_lu_inplace(acb_t det, acb_mat_t A, slong prec)
85	82	{
86	83	slong i, n, sign, rank;

91	88	sign = (rank < 0) ? -1 : 1;
92	89	rank = FLINT_ABS(rank);
93	90
94		acb_mat_diag_prod(det, A, 0, rank, prec);
	91	_acb_mat_diag_prod(det, A, 0, rank, prec);
95	92	acb_mul_si(det, det, sign, prec);
96	93
97	94	/* bound unreduced part using Hadamard's inequality */

-5

acb_mat/det_precond.c less more

74	74	}
75	75
76	76	void
77		acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong a, slong b, slong prec);
78
79		void
80	77	acb_mat_det_precond(acb_t det, const acb_mat_t A, slong prec)
81	78	{
82	79	acb_mat_t LU, Linv, Uinv;

112	109	acb_mat_one(Uinv);
113	110	acb_mat_approx_solve_triu(Uinv, LU, Uinv, 0, prec);
114	111
115		acb_mat_diag_prod(detU, Uinv, 0, n, prec);
	112	acb_mat_diag_prod(detU, Uinv, prec);
116	113
117	114	acb_mat_mul(LU, A, Uinv, prec);
118	115	_apply_permutation(LU, P, n);

157	154	_perm_clear(P);
158	155	acb_mat_clear(LU);
159	156	}
160

+48

-0

acb_mat/diag_prod.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	void
	14	_acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong a, slong b, slong prec)
	15	{
	16	if (b - a == 0)
	17	acb_one(res);
	18	else if (b - a == 1)
	19	acb_set_round(res, acb_mat_entry(A, a, a), prec);
	20	else if (b - a == 2)
	21	acb_mul(res, acb_mat_entry(A, a, a), acb_mat_entry(A, a + 1, a + 1), prec);
	22	else if (b - a == 3)
	23	{
	24	acb_mul(res, acb_mat_entry(A, a, a), acb_mat_entry(A, a + 1, a + 1), prec);
	25	acb_mul(res, res, acb_mat_entry(A, a + 2, a + 2), prec);
	26	}
	27	else
	28	{
	29	acb_t t;
	30	acb_init(t);
	31	_acb_mat_diag_prod(t, A, a, a + (b - a) / 2, prec);
	32	_acb_mat_diag_prod(res, A, a + (b - a) / 2, b, prec);
	33	acb_mul(res, res, t, prec);
	34	acb_clear(t);
	35	}
	36	}
	37
	38	void
	39	acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong prec)
	40	{
	41	slong m, n;
	42
	43	m = acb_mat_nrows(A);
	44	n = acb_mat_nrows(A);
	45
	46	_acb_mat_diag_prod(res, A, 0, FLINT_MIN(m, n), prec);
	47	}

+411

-0

acb_mat/eig_enclosure_rump.c less more

	0	/*
	1	Copyright 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	/*
	14	Follows section 13.4 of Siegfried M. Rump, "Verication methods:
	15	Rigorous results using floating-point arithmetic",
	16	Acta Numerica 19 (2010), pp. 287 - 449, implemented as
	17	verifyeig() in INTLAB.
	18
	19	Cheat sheet for the formulas in Rump's paper.
	20	Assuming U = first n-k indices, V = last k indices.
	21
	22	U U^T M = selects first n-k rows from M, zeroing rest
	23	V V^T M = selects last k rows from M, zeroing rest
	24	M U U^T = selects first n-k columns from M, zeroing rest
	25	M V V^T = selects last k columns from M, zeroing rest
	26
	27	U^T M = selects first n-k rows from M, truncating matrix
	28	V^T M = selects last k rows from M, truncating matrix
	29	M U = selects first n-k columns from M, truncating matrix
	30	M V = selects last k columns from M, truncating matrix
	31
	32	X V^T = extends X to n x n matrix (placing X on right)
	33	*/
	34
	35	static void
	36	acb_approx_neg(acb_t res, const acb_t x)
	37	{
	38	arf_neg(arb_midref(acb_realref(res)), arb_midref(acb_realref(x)));
	39	arf_neg(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)));
	40	}
	41
	42	static void
	43	acb_approx_sub(acb_t res, const acb_t x, const acb_t y, slong prec)
	44	{
	45	arf_sub(arb_midref(acb_realref(res)), arb_midref(acb_realref(x)), arb_midref(acb_realref(y)), prec, ARF_RND_DOWN);
	46	arf_sub(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)), arb_midref(acb_imagref(y)), prec, ARF_RND_DOWN);
	47	}
	48
	49	/* todo: separate out */
	50	void
	51	acb_mat_bound_max_norm(mag_t res, const acb_mat_t A)
	52	{
	53	mag_t t;
	54	slong i, j;
	55
	56	mag_init(t);
	57	mag_zero(res);
	58
	59	for (i = 0; i < acb_mat_nrows(A); i++)
	60	{
	61	for (j = 0; j < acb_mat_ncols(A); j++)
	62	{
	63	acb_get_mag(t, acb_mat_entry(A, i, j));
	64	mag_max(res, res, t);
	65	}
	66	}
	67
	68	mag_clear(t);
	69	}
	70
	71	static void
	72	arb_mat_nonnegative_eig_bound(mag_t eps, const arb_mat_t M, slong prec)
	73	{
	74	/* Cheap, but poor for defective eigenvalues */
	75	arb_mat_bound_frobenius_norm(eps, M);
	76
	77	/* Use Perron root bound. TODO: do something direct for k = 2. */
	78	if (1)
	79	{
	80	acb_mat_t A, R, E;
	81	arb_mat_t V, MV;
	82	mag_t tm, um, vbound;
	83	slong i, j, k;
	84
	85	k = arb_mat_nrows(M);
	86
	87	acb_mat_init(A, k, k);
	88	acb_mat_init(R, k, k);
	89	acb_mat_init(E, 1, k);
	90	arb_mat_init(V, k, k);
	91	arb_mat_init(MV, k, k);
	92	mag_init(tm);
	93	mag_init(um);
	94	mag_init(vbound);
	95
	96	acb_mat_set_arb_mat(A, M);
	97	/* TODO: could probably lower precision if precision is very high? */
	98	acb_mat_approx_eig_qr(acb_mat_entry(E, 0, 0), NULL, R, A, NULL, 0, prec);
	99
	100	for (i = 0; i < k; i++)
	101	{
	102	for (j = 0; j < k; j++)
	103	{
	104	acb_get_mag(tm, acb_mat_entry(R, i, j));
	105	arf_set_mag(arb_midref(arb_mat_entry(V, i, j)), tm);
	106	}
	107	}
	108
	109	arb_mat_mul(MV, M, V, MAG_BITS);
	110
	111	for (j = 0; j < k; j++)
	112	{
	113	mag_zero(vbound);
	114
	115	for (i = 0; i < k; i++)
	116	{
	117	arb_get_mag(tm, arb_mat_entry(MV, i, j));
	118	arb_get_mag_lower(um, arb_mat_entry(V, i, j));
	119	mag_div(tm, tm, um);
	120	mag_max(vbound, vbound, tm);
	121	}
	122
	123	mag_min(eps, eps, vbound);
	124	}
	125
	126	acb_mat_clear(A);
	127	acb_mat_clear(R);
	128	acb_mat_clear(E);
	129	arb_mat_clear(V);
	130	arb_mat_clear(MV);
	131	mag_clear(tm);
	132	mag_clear(um);
	133	mag_clear(vbound);
	134	}
	135	}
	136
	137	static void
	138	acb_approx_mag(mag_t res, const acb_t x)
	139	{
	140	mag_t t;
	141	mag_init(t);
	142	arf_get_mag(res, arb_midref(acb_realref(x)));
	143	arf_get_mag(t, arb_midref(acb_imagref(x)));
	144	mag_hypot(res, res, t);
	145	mag_clear(t);
	146	}
	147
	148	/* Extract k largest rows to freeze */
	149	static void
	150	partition_X_sorted(slong * u, slong * v, const acb_mat_t X, slong prec)
	151	{
	152	slong i, j, n, k, c;
	153	slong * row_idx;
	154	mag_ptr row_mag;
	155	mag_t t;
	156
	157	n = acb_mat_nrows(X);
	158	k = acb_mat_ncols(X);
	159
	160	row_mag = _mag_vec_init(n);
	161	row_idx = flint_malloc(sizeof(slong) * n);
	162	mag_init(t);
	163
	164	for (i = 0; i < n; i++)
	165	{
	166	row_idx[i] = i;
	167
	168	for (j = 0; j < k; j++)
	169	{
	170	acb_approx_mag(t, acb_mat_entry(X, i, j));
	171	mag_add(row_mag + i, row_mag + i, t);
	172	}
	173	}
	174
	175	/* Bubble sort... */
	176	for (i = 0; i < n - 1; i++)
	177	{
	178	for (j = 0; j < n - i - 1; j++)
	179	{
	180	if (mag_cmp(row_mag + j, row_mag + j + 1) > 0)
	181	{
	182	mag_swap(row_mag + j, row_mag + j + 1);
	183	c = row_idx[j];
	184	row_idx[j] = row_idx[j + 1];
	185	row_idx[j + 1] = c;
	186	}
	187	}
	188	}
	189
	190	/* Not frozen rows of the approximation. */
	191	for (i = 0; i < n - k; i++)
	192	u[i] = row_idx[i];
	193
	194	/* Frozen rows of the approximation. */
	195	for (i = 0; i < k; i++)
	196	v[i] = row_idx[n - k + i];
	197
	198	_mag_vec_clear(row_mag, n);
	199	flint_free(row_idx);
	200	mag_clear(t);
	201	}
	202
	203	static void
	204	partition_X_trivial(slong * u, slong * v, const acb_mat_t X, slong prec)
	205	{
	206	slong n, k, i;
	207
	208	n = acb_mat_nrows(X);
	209	k = acb_mat_ncols(X);
	210
	211	/* Not frozen rows of the approximation. */
	212	for (i = 0; i < n - k; i++)
	213	u[i] = i;
	214	/* Frozen rows of the approximation. */
	215	for (i = 0; i < k; i++)
	216	v[i] = n - k + i;
	217	}
	218
	219	void
	220	acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t X, const acb_mat_t A,
	221	const acb_t lambda_approx, const acb_mat_t X_approx, slong prec)
	222	{
	223	slong n, k, i, j, iter, maxiter;
	224	slong u, v;
	225	acb_mat_t R, I, T, Y, Y0, UY, VY, Yeps;
	226	mag_t eps;
	227
	228	n = acb_mat_nrows(A);
	229	k = acb_mat_ncols(X_approx);
	230
	231	if (k < 1 \|\| k > n \|\| n != acb_mat_nrows(X_approx) \|\| n != acb_mat_ncols(A))
	232	{
	233	flint_printf("bad matrix dimensions in acb_mat_eig_enclosure_rump\n");
	234	flint_abort();
	235	}
	236
	237	/* Not frozen rows of the approximation. */
	238	u = flint_malloc(sizeof(slong) * (n - k));
	239	/* Frozen rows of the approximation. */
	240	v = flint_malloc(sizeof(slong) * k);
	241
	242	if (k == 1)
	243	partition_X_sorted(u, v, X_approx, prec);
	244	else
	245	partition_X_trivial(u, v, X_approx, prec);
	246
	247	mag_init(eps);
	248	acb_mat_init(R, n, n);
	249	acb_mat_init(UY, n, k);
	250	acb_mat_init(VY, k, k);
	251	acb_mat_init(T, n, n);
	252	acb_mat_init(Y, n, k);
	253	acb_mat_init(Y0, n, k);
	254	acb_mat_init(Yeps, n, k);
	255
	256	/* Preconditioner:
	257	R ~= ((A - lambda_approx I) U U^T - X_approx V^T)^(-1) */
	258	acb_mat_get_mid(R, A);
	259	for (i = 0; i < n; i++)
	260	acb_approx_sub(acb_mat_entry(R, i, i),
	261	acb_mat_entry(R, i, i), lambda_approx, prec);
	262
	263	for (i = 0; i < n; i++)
	264	for (j = 0; j < k; j++)
	265	acb_approx_neg(acb_mat_entry(R, i, v[j]),
	266	acb_mat_entry(X_approx, i, j));
	267
	268	acb_mat_init(I, n, n);
	269	acb_mat_one(I);
	270	acb_mat_approx_solve(R, R, I, prec);
	271	acb_mat_clear(I);
	272
	273	/* T = I - R * ((A - lambda_approx I) U U^T - X_approx V^T) */
	274	/* Y = Y_0 = -R * ((A - lambda_approx I) X_approx) */
	275	acb_mat_set(T, A);
	276	for (i = 0; i < n; i++)
	277	acb_sub(acb_mat_entry(T, i, i), acb_mat_entry(T, i, i), lambda_approx, prec);
	278
	279	acb_mat_mul(Y0, T, X_approx, prec);
	280	acb_mat_mul(Y0, R, Y0, prec);
	281	acb_mat_neg(Y0, Y0);
	282	acb_mat_set(Y, Y0);
	283
	284	for (i = 0; i < n; i++)
	285	for (j = 0; j < k; j++)
	286	acb_neg(acb_mat_entry(T, i, v[j]), acb_mat_entry(X_approx, i, j));
	287
	288	acb_mat_mul(T, R, T, prec);
	289	acb_mat_neg(T, T);
	290	for (i = 0; i < n; i++)
	291	acb_add_ui(acb_mat_entry(T, i, i), acb_mat_entry(T, i, i), 1, prec);
	292
	293	/* Iteration with epsilon-inflation */
	294	/* Y represents the error with respect to lambda_approx and X_approx */
	295	/* TODO: what number of iterations is actually reasonable? */
	296	/* TODO: what size of epsilon is actually reasonable? */
	297	maxiter = 5 + FLINT_BIT_COUNT(prec);
	298	for (iter = 0; iter < maxiter; iter++)
	299	{
	300	/* Inflate Y. TODO: make it elementwise? */
	301	acb_mat_bound_max_norm(eps, Y);
	302	if (mag_is_zero(eps))
	303	mag_set_ui_2exp_si(eps, 1, -20 * prec);
	304	mag_mul_2exp_si(eps, eps, -3 + 2 * iter);
	305	/* if (iter > 3)
	306	mag_mul_2exp_si(eps, eps, (prec / 2) * (iter - 3) / (maxiter - 3)); */
	307
	308	acb_mat_add_error_mag(Y, eps);
	309	acb_mat_set(Yeps, Y);
	310
	311	/* Y = Y0 + T Y + R ((U U^T Y) V^T Y) */
	312	acb_mat_zero(UY);
	313	acb_mat_zero(VY);
	314
	315	/* U U^T Y -- zero the rows at indices v. */
	316	acb_mat_set(UY, Y);
	317	for (i = 0; i < k; i++)
	318	for (j = 0; j < k; j++)
	319	acb_zero(acb_mat_entry(UY, v[i], j));
	320
	321	/* V^T Y -- extract rows at indices v */
	322	for (i = 0; i < k; i++)
	323	for (j = 0; j < k; j++)
	324	acb_set(acb_mat_entry(VY, i, j), acb_mat_entry(Y, v[i], j));
	325
	326	acb_mat_mul(UY, UY, VY, prec);
	327	acb_mat_mul(UY, R, UY, prec);
	328
	329	acb_mat_mul(Y, T, Y, prec);
	330	acb_mat_add(Y, Y, UY, prec);
	331	acb_mat_add(Y, Y, Y0, prec);
	332
	333	if (acb_mat_contains(Yeps, Y))
	334	{
	335	acb_get_mid(lambda, lambda_approx);
	336
	337	if (J != NULL)
	338	{
	339	/* J = lambda_approx I_k + V^T Y */
	340	for (i = 0; i < k; i++)
	341	for (j = 0; j < k; j++)
	342	acb_set(acb_mat_entry(J, i, j), acb_mat_entry(Y, v[i], j));
	343
	344	for (i = 0; i < k; i++)
	345	acb_add(acb_mat_entry(J, i, i), acb_mat_entry(J, i, i), lambda, prec);
	346	}
	347
	348	/* The correction for the frozen rows corresponds
	349	to the eigenvalue. */
	350	if (k == 1)
	351	{
	352	/* Just one eigenvalue. */
	353	acb_get_mag(eps, acb_mat_entry(Y, v[0], 0));
	354	}
	355	else
	356	{
	357	/* Inclusion of eigenvalues of lambda_approx I_k + V^T Y. */
	358	arb_mat_t M;
	359	arb_mat_init(M, k, k);
	360
	361	/* Extract rows of Y corresponding to the eigenvalue correction. */
	362	for (i = 0; i < k; i++)
	363	{
	364	for (j = 0; j < k; j++)
	365	{
	366	acb_get_mag(eps, acb_mat_entry(Y, v[i], j));
	367	arf_set_mag(arb_midref(arb_mat_entry(M, i, j)), eps);
	368	}
	369	}
	370
	371	arb_mat_nonnegative_eig_bound(eps, M, prec);
	372	arb_mat_clear(M);
	373	}
	374
	375	/* Error bound for eigenvalues. */
	376	acb_add_error_mag(lambda, eps);
	377
	378	acb_mat_get_mid(X, X_approx);
	379	/* Error bounds for eigenvectors. */
	380	/* Update the not frozen rows of the eigenvectors. */
	381	for (i = 0; i < n - k; i++)
	382	{
	383	for (j = 0; j < k; j++)
	384	acb_add(acb_mat_entry(X, u[i], j),
	385	acb_mat_entry(X, u[i], j),
	386	acb_mat_entry(Y, u[i], j), prec);
	387	}
	388
	389	goto cleanup;
	390	}
	391	}
	392
	393	/* We failed to find an enclosure. */
	394	acb_indeterminate(lambda);
	395	acb_mat_indeterminate(X);
	396	if (J != NULL)
	397	acb_mat_indeterminate(J);
	398
	399	cleanup:
	400	acb_mat_clear(R);
	401	acb_mat_clear(T);
	402	acb_mat_clear(Y);
	403	acb_mat_clear(Y0);
	404	acb_mat_clear(Yeps);
	405	acb_mat_clear(UY);
	406	acb_mat_clear(VY);
	407	mag_clear(eps);
	408	flint_free(u);
	409	flint_free(v);
	410	}

+63

-0

acb_mat/eig_global_enclosure.c less more

	0	/*
	1	Copyright 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	void
	14	acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)
	15	{
	16	acb_mat_t Y, R1, R2;
	17	slong i, j, n;
	18	mag_t r1, r2;
	19
	20	n = acb_mat_nrows(A);
	21
	22	acb_mat_init(Y, n, n);
	23	acb_mat_init(R1, n, n);
	24	acb_mat_init(R2, n, n);
	25	mag_init(r1);
	26	mag_init(r2);
	27
	28	/* Y ~= inv(R) */
	29	acb_mat_one(R1);
	30	acb_mat_approx_solve(Y, R, R1, prec);
	31
	32	/* R2 = YR - I /
	33	acb_mat_mul(R2, Y, R, prec);
	34	for (i = 0; i < n; i++)
	35	acb_sub_ui(acb_mat_entry(R2, i, i), acb_mat_entry(R2, i, i), 1, prec);
	36
	37	acb_mat_bound_inf_norm(r2, R2);
	38
	39	if (mag_cmp_2exp_si(r2, 0) < 0)
	40	{
	41	/* R1 = Y(AR - RD) /
	42	acb_mat_mul(R2, A, R, prec);
	43	for (i = 0; i < n; i++)
	44	for (j = 0; j < n; j++)
	45	acb_submul(acb_mat_entry(R2, i, j), acb_mat_entry(R, i, j), E + j, prec);
	46	acb_mat_mul(R1, Y, R2, prec);
	47
	48	acb_mat_bound_inf_norm(r1, R1);
	49	mag_geom_series(r2, r2, 0);
	50	mag_mul(eps, r1, r2);
	51	}
	52	else
	53	{
	54	mag_inf(eps);
	55	}
	56
	57	acb_mat_clear(R1);
	58	acb_mat_clear(R2);
	59	acb_mat_clear(Y);
	60	mag_clear(r1);
	61	mag_clear(r2);
	62	}

+33

-0

acb_mat/eig_multiple.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_eig_multiple(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	15	{
	16	slong n;
	17	acb_ptr F;
	18	int success;
	19
	20	n = arb_mat_nrows(A);
	21	F = _acb_vec_init(n);
	22
	23	success = acb_mat_eig_simple_vdhoeven_mourrain(F, NULL, NULL, A, E_approx, R_approx, prec);
	24
	25	if (!success)
	26	success = acb_mat_eig_multiple_rump(F, A, E_approx, R_approx, prec);
	27
	28	_acb_vec_set(E, F, n);
	29	_acb_vec_clear(F, n);
	30
	31	return success;
	32	}

+181

-0

acb_mat/eig_multiple_rump.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	static void
	14	acb_approx_mag(mag_t res, const acb_t x)
	15	{
	16	mag_t t;
	17	mag_init(t);
	18	arf_get_mag(res, arb_midref(acb_realref(x)));
	19	arf_get_mag(t, arb_midref(acb_imagref(x)));
	20	mag_hypot(res, res, t);
	21	mag_clear(t);
	22	}
	23
	24	static int
	25	close(const acb_t x, const acb_t y, const mag_t eps)
	26	{
	27	arf_t t;
	28	mag_t a, b;
	29	int result;
	30
	31	mag_init(a);
	32	mag_init(b);
	33	arf_init(t);
	34	arf_sub(t, arb_midref(acb_realref(x)), arb_midref(acb_realref(y)), MAG_BITS, ARF_RND_UP);
	35	arf_get_mag(a, t);
	36	arf_sub(t, arb_midref(acb_imagref(x)), arb_midref(acb_imagref(y)), MAG_BITS, ARF_RND_UP);
	37	arf_get_mag(b, t);
	38
	39	mag_hypot(a, a, b);
	40	result = (mag_cmp(a, eps) <= 0);
	41
	42	mag_clear(a);
	43	mag_clear(b);
	44	arf_clear(t);
	45
	46	return result;
	47	}
	48
	49	int
	50	acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	51	{
	52	slong c, i, j, k, n;
	53	acb_mat_t X;
	54	acb_ptr F;
	55	int result;
	56	slong iter;
	57	mag_t escale, eps, tm, um;
	58	slong ** cluster;
	59	slong * cluster_size;
	60	slong num_clusters;
	61
	62	n = acb_mat_nrows(A);
	63
	64	if (n == 0)
	65	return 1;
	66
	67	cluster = flint_malloc(sizeof(slong ) n);
	68	for (i = 0; i < n; i++)
	69	cluster[i] = flint_malloc(sizeof(slong) * n);
	70	cluster_size = flint_malloc(sizeof(slong) * n);
	71
	72	mag_init(eps);
	73	mag_init(escale);
	74	mag_init(tm);
	75	mag_init(um);
	76
	77	mag_zero(escale);
	78	for (i = 0; i < n; i++)
	79	{
	80	acb_approx_mag(tm, E_approx + i);
	81	mag_max(escale, escale, tm);
	82	}
	83
	84	/* todo: when num_clusters = 1, could use fallback global enclosure */
	85
	86	/* todo: initial clustering could be allowed to be zero */
	87	/* M 2^(-0.75p), M 2^(-0.5p), M 2^(-0.25p), ... */
	88	for (iter = 0; iter < 2; iter++)
	89	{
	90	mag_mul_2exp_si(eps, escale, -prec + (iter + 1) * prec/4);
	91
	92	/* Group the eigenvalue approximations. */
	93	num_clusters = 0;
	94	for (i = 0; i < n; i++)
	95	{
	96	int new_cluster = 1;
	97
	98	for (j = 0; j < num_clusters && new_cluster; j++)
	99	{
	100	if (close(E_approx + i, E_approx + cluster[j][0], eps))
	101	{
	102	cluster[j][cluster_size[j]] = i;
	103	cluster_size[j]++;
	104	new_cluster = 0;
	105	}
	106	}
	107
	108	if (new_cluster)
	109	{
	110	cluster[num_clusters][0] = i;
	111	cluster_size[num_clusters] = 1;
	112	num_clusters++;
	113	}
	114	}
	115
	116	result = 1;
	117
	118	F = _acb_vec_init(num_clusters);
	119
	120	for (c = 0; c < num_clusters && result; c++)
	121	{
	122	k = cluster_size[c];
	123
	124	acb_mat_init(X, n, k);
	125
	126	for (i = 0; i < n; i++)
	127	for (j = 0; j < k; j++)
	128	acb_set(acb_mat_entry(X, i, j), acb_mat_entry(R_approx, i, cluster[c][j]));
	129
	130	acb_mat_eig_enclosure_rump(F + c, NULL, X, A, E_approx + cluster[c][0], X, prec);
	131
	132	if (!acb_is_finite(F + c))
	133	result = 0;
	134
	135	acb_mat_clear(X);
	136	}
	137
	138	for (i = 0; i < num_clusters; i++)
	139	{
	140	for (j = i + 1; j < num_clusters; j++)
	141	{
	142	if (acb_overlaps(F + i, F + j))
	143	result = 0;
	144	}
	145	}
	146
	147	if (result)
	148	{
	149	i = 0;
	150	for (c = 0; c < num_clusters; c++)
	151	{
	152	for (j = 0; j < cluster_size[c]; j++)
	153	{
	154	acb_set(E + i, F + c);
	155	i++;
	156	}
	157	}
	158	}
	159
	160	_acb_vec_clear(F, num_clusters);
	161
	162	if (result)
	163	break;
	164	}
	165
	166	if (!result)
	167	_acb_vec_indeterminate(E, n);
	168
	169	for (i = 0; i < n; i++)
	170	flint_free(cluster[i]);
	171	flint_free(cluster);
	172	flint_free(cluster_size);
	173
	174	mag_clear(eps);
	175	mag_clear(escale);
	176	mag_clear(tm);
	177	mag_clear(um);
	178
	179	return result;
	180	}

+19

-0

acb_mat/eig_simple.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_eig_simple(acb_ptr E, acb_mat_t L, acb_mat_t R,
	15	const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	16	{
	17	return acb_mat_eig_simple_vdhoeven_mourrain(E, L, R, A, E_approx, R_approx, prec);
	18	}

+81

-0

acb_mat/eig_simple_rump.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_eig_simple_rump(acb_ptr E, acb_mat_t L, acb_mat_t R,
	15	const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	16	{
	17	slong i, j, n;
	18	acb_mat_t X, R2;
	19	int result;
	20
	21	n = acb_mat_nrows(A);
	22
	23	if (n == 0)
	24	return 1;
	25
	26	if (n == 1)
	27	{
	28	acb_set_round(E, acb_mat_entry(A, 0, 0), prec);
	29	if (L != NULL)
	30	acb_one(acb_mat_entry(L, 0, 0));
	31	if (R != NULL)
	32	acb_one(acb_mat_entry(R, 0, 0));
	33	return 1;
	34	}
	35
	36	acb_mat_init(X, n, 1);
	37	acb_mat_init(R2, n, n);
	38
	39	result = 1;
	40
	41	for (i = 0; i < n && result; i++)
	42	{
	43	for (j = 0; j < n; j++)
	44	acb_set(acb_mat_entry(X, j, 0), acb_mat_entry(R_approx, j, i));
	45
	46	acb_mat_eig_enclosure_rump(E + i, NULL, X, A, E_approx + i, X, prec);
	47
	48	if (!acb_is_finite(E + i))
	49	result = 0;
	50
	51	for (j = 0; j < i; j++)
	52	if (acb_overlaps(E + i, E + j))
	53	result = 0;
	54
	55	for (j = 0; j < n; j++)
	56	acb_set(acb_mat_entry(R2, j, i), acb_mat_entry(X, j, 0));
	57	}
	58
	59	if (R != NULL)
	60	{
	61	if (result)
	62	acb_mat_set(R, R2);
	63	else
	64	acb_mat_indeterminate(R);
	65	}
	66
	67	if (L != NULL)
	68	{
	69	if (!result \|\| !acb_mat_inv(L, R, prec))
	70	acb_mat_indeterminate(L);
	71	}
	72
	73	if (!result)
	74	_acb_vec_indeterminate(E, n);
	75
	76	acb_mat_clear(X);
	77	acb_mat_clear(R2);
	78
	79	return result;
	80	}

+252

-0

acb_mat/eig_simple_vdhoeven_mourrain.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	/* todo: move out */
	14	void
	15	acb_mat_inf_norm(arb_t res, const acb_mat_t A, slong prec)
	16	{
	17	slong i, j, m, n;
	18	arb_t s, t;
	19
	20	m = acb_mat_nrows(A);
	21	n = acb_mat_nrows(A);
	22
	23	if (m == 0 \|\| n == 0)
	24	{
	25	arb_zero(res);
	26	return;
	27	}
	28
	29	arb_init(s);
	30	arb_init(t);
	31
	32	arb_zero(res);
	33
	34	for (i = 0; i < m; i++)
	35	{
	36	acb_abs(s, acb_mat_entry(A, i, 0), prec);
	37
	38	for (j = 1; j < n; j++)
	39	{
	40	acb_abs(t, acb_mat_entry(A, i, j), prec);
	41	arb_add(s, s, t, prec);
	42	}
	43
	44	arb_max(res, res, s, prec);
	45	}
	46
	47	arb_clear(s);
	48	arb_clear(t);
	49	}
	50
	51	static void
	52	diagonal_certify(arb_t epsilon, arb_t eta1, arb_t eta2, const acb_mat_t D, const acb_mat_t H, slong prec)
	53	{
	54	arb_t mu, sigma, alpha, t, u, v;
	55	acb_t d;
	56	slong i, j, n;
	57
	58	arb_init(mu);
	59	arb_init(sigma);
	60	arb_init(alpha);
	61	arb_init(t);
	62	arb_init(u);
	63	arb_init(v);
	64	acb_init(d);
	65
	66	n = acb_mat_nrows(D);
	67
	68	/* D = diagonal matrix; H = off diagonal matrix */
	69
	70	/* mu = \|\|D\|\| */
	71	acb_mat_inf_norm(mu, D, prec);
	72
	73	/* sigma = sigma(D) = separation number */
	74	arb_pos_inf(sigma);
	75	for (i = 0; i < n; i++)
	76	{
	77	for (j = i + 1; j < n; j++)
	78	{
	79	acb_sub(d, acb_mat_entry(D, i, i), acb_mat_entry(D, j, j), prec);
	80	acb_abs(t, d, prec);
	81	arb_min(sigma, sigma, t, prec);
	82	}
	83	}
	84
	85	/* eta1 = \|\|Delta(H)\|\| Delta = diagonal projection */
	86	arb_zero(eta1);
	87
	88	/* eta2 = \|\|Omega(H)\|\| Omega = off-diagonal projection */
	89	acb_mat_inf_norm(eta2, H, prec);
	90
	91	/* alpha = min(sigma / (6 mu), 1/4) */
	92	arb_div(t, sigma, mu, prec);
	93	arb_div_ui(t, t, 6, prec);
	94	arb_set_d(u, 0.25);
	95	arb_min(alpha, t, u, prec);
	96
	97	arb_add(t, eta1, eta2, prec);
	98	arb_mul(u, alpha, mu, prec);
	99	arb_mul_2exp_si(u, u, -3);
	100
	101	arb_mul(v, alpha, sigma, prec);
	102	arb_mul_2exp_si(v, v, -3);
	103
	104	if (arb_le(t, u) && arb_le(eta2, v))
	105	{
	106	arb_div(epsilon, eta2, sigma, prec);
	107	arb_mul_ui(epsilon, epsilon, 3, prec);
	108	}
	109	else
	110	{
	111	arb_indeterminate(epsilon);
	112	}
	113
	114	arb_clear(mu);
	115	arb_clear(sigma);
	116	arb_clear(alpha);
	117	arb_clear(t);
	118	arb_clear(u);
	119	arb_clear(v);
	120	acb_clear(d);
	121	}
	122
	123	int
	124	acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E,
	125	acb_mat_t L, acb_mat_t R, const acb_mat_t A,
	126	acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	127	{
	128	acb_mat_t D, T, AT;
	129	int result;
	130	slong i, j, n;
	131
	132	result = 0;
	133	n = acb_mat_nrows(A);
	134
	135	if (n == 0)
	136	return 1;
	137
	138	if (n == 1)
	139	{
	140	acb_set_round(E, acb_mat_entry(A, 0, 0), prec);
	141	if (L != NULL)
	142	acb_one(acb_mat_entry(L, 0, 0));
	143	if (R != NULL)
	144	acb_one(acb_mat_entry(R, 0, 0));
	145	return 1;
	146	}
	147
	148	acb_mat_init(D, n, n);
	149	acb_mat_init(T, n, n);
	150	acb_mat_init(AT, n, n);
	151
	152	/* T D = A T */
	153	acb_mat_get_mid(T, R_approx);
	154	acb_mat_mul(AT, A, T, prec);
	155
	156	if (acb_mat_solve(D, T, AT, prec))
	157	{
	158	acb_mat_t DD, DH;
	159	arb_t epsilon, eta1, eta2;
	160
	161	arb_init(epsilon);
	162	arb_init(eta1);
	163	arb_init(eta2);
	164	acb_mat_init(DD, n, n);
	165	acb_mat_init(DH, n, n);
	166
	167	for (i = 0; i < n; i++)
	168	acb_set(acb_mat_entry(DD, i, i), acb_mat_entry(D, i, i));
	169
	170	for (i = 0; i < n; i++)
	171	for (j = 0; j < n; j++)
	172	if (i != j)
	173	acb_set(acb_mat_entry(DH, i, j), acb_mat_entry(D, i, j));
	174
	175	diagonal_certify(epsilon, eta1, eta2, DD, DH, 2 * MAG_BITS);
	176
	177	if (arb_is_finite(epsilon))
	178	{
	179	for (i = 0; i < n; i++)
	180	{
	181	/* note: paper actually uses D_c which would be better? */
	182	acb_set(E + i, acb_mat_entry(D, i, i));
	183	arb_add_error(acb_realref(E + i), eta2);
	184	arb_add_error(acb_imagref(E + i), eta2);
	185	}
	186
	187	result = 1;
	188
	189	/* unlikely */
	190	for (i = 0; i < n; i++)
	191	for (j = i + 1; j < n; j++)
	192	if (acb_overlaps(E + i, E + j))
	193	result = 0;
	194
	195	if (result && (R != NULL \|\| L != NULL))
	196	{
	197	mag_t ep, em;
	198	mag_init(ep);
	199	mag_init(em);
	200
	201	arb_get_mag(ep, epsilon);
	202	acb_mat_zero(D);
	203	acb_mat_add_error_mag(D, ep);
	204	acb_mat_mul(D, T, D, MAG_BITS);
	205
	206	for (i = 0; i < n; i++)
	207	{
	208	for (j = 0; j < n; j++)
	209	{
	210	acb_get_mag(ep, acb_mat_entry(D, i, j));
	211	acb_add_error_mag(acb_mat_entry(T, i, j), ep);
	212	}
	213	}
	214
	215	if (R != NULL)
	216	acb_mat_set(R, T);
	217
	218	if (L != NULL)
	219	{
	220	if (!acb_mat_inv(L, T, prec))
	221	acb_mat_indeterminate(L);
	222	}
	223
	224	mag_clear(ep);
	225	mag_clear(em);
	226	}
	227	}
	228
	229	acb_mat_clear(DD);
	230	acb_mat_clear(DH);
	231	arb_clear(epsilon);
	232	arb_clear(eta1);
	233	arb_clear(eta2);
	234	}
	235
	236	if (!result)
	237	{
	238	for (i = 0; i < n; i++)
	239	acb_indeterminate(E + i);
	240	if (L != NULL)
	241	acb_mat_indeterminate(L);
	242	if (R != NULL)
	243	acb_mat_indeterminate(R);
	244	}
	245
	246	acb_mat_clear(D);
	247	acb_mat_clear(T);
	248	acb_mat_clear(AT);
	249
	250	return result;
	251	}

+22

-0

acb_mat/indeterminate.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	void
	14	acb_mat_indeterminate(acb_mat_t A)
	15	{
	16	slong i, j;
	17
	18	for (i = 0; i < acb_mat_nrows(A); i++)
	19	for (j = 0; j < acb_mat_ncols(A); j++)
	20	acb_indeterminate(acb_mat_entry(A, i, j));
	21	}

+26

-0

acb_mat/is_exact.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_is_exact(const acb_mat_t A)
	15	{
	16	slong i, j;
	17
	18	for (i = 0; i < acb_mat_nrows(A); i++)
	19	for (j = 0; j < acb_mat_ncols(A); j++)
	20	if (!mag_is_zero(arb_radref(acb_realref(acb_mat_entry(A, i, j)))) \|\|
	21	!mag_is_zero(arb_radref(acb_imagref(acb_mat_entry(A, i, j)))))
	22	return 0;
	23
	24	return 1;
	25	}

+28

-0

acb_mat/is_finite.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_is_finite(const acb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = acb_mat_nrows(A);
	19	m = acb_mat_ncols(A);
	20
	21	for (i = 0; i < n; i++)
	22	for (j = 0; j < m; j++)
	23	if (!acb_is_finite(acb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+28

-0

acb_mat/is_tril.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_is_tril(const acb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = acb_mat_nrows(A);
	19	m = acb_mat_ncols(A);
	20
	21	for (i = 0; i < n; i++)
	22	for (j = i + 1; j < m; j++)
	23	if (!acb_is_zero(acb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+28

-0

acb_mat/is_triu.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_is_triu(const acb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = acb_mat_nrows(A);
	19	m = acb_mat_ncols(A);
	20
	21	for (i = 1; i < n; i++)
	22	for (j = 0; j < FLINT_MIN(i, m); j++)
	23	if (!acb_is_zero(acb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+28

-0

acb_mat/is_zero.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	acb_mat_is_zero(const acb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = acb_mat_nrows(A);
	19	m = acb_mat_ncols(A);
	20
	21	for (i = 0; i < n; i++)
	22	for (j = 0; j < m; j++)
	23	if (!acb_is_zero(acb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+58

-0

acb_mat/randtest_eig.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	void
	14	acb_mat_randtest_eig(acb_mat_t A, flint_rand_t state, acb_srcptr E, slong prec)
	15	{
	16	slong n, i, j, ebits;
	17	acb_mat_t U, Q;
	18
	19	n = acb_mat_nrows(A);
	20	ebits = 1 + n_randint(state, 5);
	21
	22	acb_mat_init(U, n, n);
	23	acb_mat_init(Q, n, n);
	24
	25	/* Skew-Hermitian matrix */
	26	acb_mat_randtest(Q, state, prec, 1);
	27	if (n_randint(state, 2))
	28	acb_mat_get_mid(Q, Q);
	29	for (i = 0; i < n; i++)
	30	{
	31	for (j = i + 1; j < n; j++)
	32	{
	33	acb_neg(acb_mat_entry(Q, i, j), acb_mat_entry(Q, j, i));
	34	acb_conj(acb_mat_entry(Q, i, j), acb_mat_entry(Q, i, j));
	35	}
	36	arb_zero(acb_realref(acb_mat_entry(Q, i, i)));
	37	}
	38
	39	acb_mat_exp(Q, Q, prec);
	40
	41	acb_mat_randtest(U, state, prec, ebits);
	42	if (n_randint(state, 2))
	43	acb_mat_get_mid(U, U);
	44	for (i = 0; i < n; i++)
	45	for (j = 0; j < i; j++)
	46	acb_zero(acb_mat_entry(U, i, j));
	47
	48	for (i = 0; i < n; i++)
	49	acb_set(acb_mat_entry(U, i, i), E + i);
	50
	51	acb_mat_mul(U, Q, U, prec);
	52	acb_mat_conjugate_transpose(Q, Q);
	53	acb_mat_mul(A, U, Q, prec);
	54
	55	acb_mat_clear(U);
	56	acb_mat_clear(Q);
	57	}

+194

-0

acb_mat/test/t-approx_eig_qr.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("approx_eig_qr....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	/* Test random & DFT matrices */
	24	for (iter = 0; iter < 200 * arb_test_multiplier(); iter++)
	25	{
	26	acb_mat_t A, L, R;
	27	acb_ptr E;
	28	acb_t t;
	29	mag_t b;
	30	slong i, j, n, prec, goal, c0, c1, c2, c3;
	31	int wantL, wantR, result, dft;
	32
	33	dft = n_randint(state, 2);
	34	if (dft)
	35	n = n_randint(state, 30);
	36	else
	37	n = n_randint(state, 15);
	38	goal = 2 + n_randint(state, 100);
	39	wantL = n_randint(state, 2);
	40	wantR = n_randint(state, 2);
	41
	42	acb_mat_init(A, n, n);
	43	acb_mat_init(L, n, n);
	44	acb_mat_init(R, n, n);
	45	acb_init(t);
	46	mag_init(b);
	47	E = _acb_vec_init(n);
	48
	49	for (prec = 32; ; prec *= 2)
	50	{
	51	if (dft)
	52	{
	53	acb_mat_dft(A, 0, prec);
	54	}
	55	else
	56	{
	57	acb_mat_randtest(A, state, 2 + n_randint(state, 200), 5);
	58	acb_mat_get_mid(A, A);
	59	}
	60
	61	acb_mat_approx_eig_qr(E, wantL ? L : NULL, wantR ? R : NULL, A, NULL, 0, prec);
	62
	63	if (dft)
	64	{
	65	/* Verify the known eigenvalues + multiplicities */
	66	c0 = c1 = c2 = c3 = 0;
	67
	68	for (i = 0; i < n; i++)
	69	{
	70	acb_set_d_d(t, 1.0, 0.0);
	71	acb_sub(t, t, E + i, prec);
	72	acb_get_mag(b, t);
	73	c0 += (mag_cmp_2exp_si(b, -goal) < 0);
	74
	75	acb_set_d_d(t, -1.0, 0.0);
	76	acb_sub(t, t, E + i, prec);
	77	acb_get_mag(b, t);
	78	c1 += (mag_cmp_2exp_si(b, -goal) < 0);
	79
	80	acb_set_d_d(t, 0.0, 1.0);
	81	acb_sub(t, t, E + i, prec);
	82	acb_get_mag(b, t);
	83	c2 += (mag_cmp_2exp_si(b, -goal) < 0);
	84
	85	acb_set_d_d(t, 0.0, -1.0);
	86	acb_sub(t, t, E + i, prec);
	87	acb_get_mag(b, t);
	88	c3 += (mag_cmp_2exp_si(b, -goal) < 0);
	89	}
	90
	91	result = (n == 0 \|\| (c0 == (n+4)/4 && c1 == (n+2)/4 && c2 == (n-1)/4 && c3 == (n+1)/4));
	92	}
	93	else
	94	{
	95	result = 1;
	96	}
	97
	98	if (result && wantL)
	99	{
	100	acb_mat_t LA, D;
	101	acb_mat_init(LA, n, n);
	102	acb_mat_init(D, n, n);
	103
	104	/* Check LA - lambda L = 0 */
	105	acb_mat_approx_mul(LA, L, A, prec);
	106	for (i = 0; i < n; i++)
	107	acb_set(acb_mat_entry(D, i, i), E + i);
	108	acb_mat_approx_mul(D, D, L, prec);
	109	acb_mat_sub(LA, LA, D, prec);
	110
	111	for (i = 0; i < n; i++)
	112	{
	113	for (j = 0; j < n; j++)
	114	{
	115	acb_get_mag(b, acb_mat_entry(LA, i, j));
	116	result = result && (mag_cmp_2exp_si(b, -goal) < 0);
	117	}
	118	}
	119
	120	acb_mat_clear(LA);
	121	acb_mat_clear(D);
	122	}
	123
	124	if (result && wantR)
	125	{
	126	acb_mat_t AR, D;
	127	acb_mat_init(AR, n, n);
	128	acb_mat_init(D, n, n);
	129
	130	/* Check AR - R lambda = 0 */
	131	acb_mat_approx_mul(AR, A, R, prec);
	132	for (i = 0; i < n; i++)
	133	acb_set(acb_mat_entry(D, i, i), E + i);
	134	acb_mat_approx_mul(D, R, D, prec);
	135	acb_mat_sub(AR, AR, D, prec);
	136
	137	for (i = 0; i < n; i++)
	138	{
	139	for (j = 0; j < n; j++)
	140	{
	141	acb_get_mag(b, acb_mat_entry(AR, i, j));
	142	result = result && (mag_cmp_2exp_si(b, -goal) < 0);
	143	}
	144	}
	145
	146	acb_mat_clear(AR);
	147	acb_mat_clear(D);
	148	}
	149
	150	if (result)
	151	break;
	152
	153	if (prec > 2000)
	154	{
	155	flint_printf("FAIL (convergence, dft = %d)\n\n", dft);
	156	flint_printf("n = %wd\n\n", n);
	157	acb_mat_printd(A, 10);
	158	flint_printf("\n\n");
	159	for (i = 0; i < n; i++)
	160	{
	161	acb_printn(E + i, 50, 0);
	162	flint_printf("\n");
	163	}
	164	flint_printf("\n");
	165	if (wantL)
	166	{
	167	flint_printf("L = \n");
	168	acb_mat_printd(L, 10);
	169	flint_printf("\n\n");
	170	}
	171	if (wantR)
	172	{
	173	flint_printf("R = \n");
	174	acb_mat_printd(R, 10);
	175	flint_printf("\n\n");
	176	}
	177	flint_abort();
	178	}
	179	}
	180
	181	acb_mat_clear(A);
	182	acb_mat_clear(L);
	183	acb_mat_clear(R);
	184	_acb_vec_clear(E, n);
	185	acb_clear(t);
	186	mag_clear(b);
	187	}
	188
	189	flint_randclear(state);
	190	flint_cleanup();
	191	flint_printf("PASS\n");
	192	return EXIT_SUCCESS;
	193	}

+65

-0

acb_mat/test/t-companion.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int
	14	main(void)
	15	{
	16	slong iter;
	17	flint_rand_t state;
	18
	19	flint_printf("companion....");
	20	fflush(stdout);
	21
	22	flint_randinit(state);
	23
	24	for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
	25	{
	26	acb_mat_t A;
	27	acb_poly_t f, g;
	28	slong n, prec;
	29
	30	acb_poly_init(f);
	31	acb_poly_init(g);
	32
	33	do {
	34	acb_poly_randtest(f, state, 1 + n_randint(state, 8), 1 + n_randint(state, 1000), 10);
	35	} while (acb_poly_degree(f) < 0);
	36
	37	n = acb_poly_degree(f);
	38	prec = 2 + n_randint(state, 200);
	39
	40	acb_mat_init(A, n, n);
	41	acb_mat_randtest(A, state, 1 + n_randint(state, 1000), 10);
	42	acb_mat_companion(A, f, prec);
	43	acb_mat_charpoly(g, A, prec);
	44	acb_poly_scalar_mul(g, g, acb_poly_get_coeff_ptr(f, n), prec);
	45
	46	if (!acb_poly_contains(g, f))
	47	{
	48	flint_printf("FAIL\n");
	49	flint_printf("A:\n"), acb_mat_printd(A, 15), flint_printf("\n");
	50	flint_printf("f:\n"), acb_poly_printd(f, 15), flint_printf("\n");
	51	flint_printf("g:\n"), acb_poly_printd(g, 15), flint_printf("\n");
	52	flint_abort();
	53	}
	54
	55	acb_mat_clear(A);
	56	acb_poly_clear(f);
	57	acb_poly_clear(g);
	58	}
	59
	60	flint_randclear(state);
	61	flint_cleanup();
	62	flint_printf("PASS\n");
	63	return 0;
	64	}

+247

-0

acb_mat/test/t-eig_enclosure_rump.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("eig_enclosure_rump....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	/* Test random matrices */
	24	for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
	25	{
	26	acb_mat_t A, X, R, AR, J, RJ, Z, Z0;
	27	acb_ptr E, F;
	28	acb_t b, lambda;
	29	slong i, j, h, k, n, iter2, prec, found_eigenvalue;
	30
	31	n = 1 + n_randint(state, 7);
	32	prec = 2 + n_randint(state, 200);
	33
	34	acb_mat_init(A, n, n);
	35	acb_mat_init(X, n, n);
	36	acb_init(lambda);
	37	acb_init(b);
	38	E = _acb_vec_init(n);
	39	F = _acb_vec_init(n);
	40
	41	if (n_randint(state, 2))
	42	{
	43	for (i = 0; i < n; i++)
	44	acb_randtest(E + i, state, prec, 3);
	45	}
	46	else
	47	{
	48	/* Randomly repeat eigenvalues. */
	49	for (i = 0; i < n; i++)
	50	{
	51	if (i == 0 \|\| n_randint(state, 2))
	52	acb_randtest(E + i, state, prec, 3);
	53	else
	54	acb_set(E + i, E + n_randint(state, i));
	55	}
	56	}
	57
	58	if (n_randint(state, 2))
	59	{
	60	for (i = 0; i < n; i++)
	61	acb_get_mid(E + i, E + i);
	62	}
	63
	64	acb_mat_randtest_eig(A, state, E, prec);
	65	acb_mat_approx_eig_qr(F, NULL, X, A, NULL, 0, prec);
	66
	67	/* Perturb F further. */
	68	if (n_randint(state, 4) == 0)
	69	{
	70	for (i = 0; i < n; i++)
	71	{
	72	acb_randtest(b, state, prec, 1);
	73	acb_mul_2exp_si(b, b, -n_randint(state, prec));
	74	acb_add(F + i, F + i, b, prec);
	75	}
	76	}
	77
	78	/* Perturb X further. */
	79	if (n_randint(state, 10) == 0)
	80	{
	81	j = n_randint(state, n);
	82
	83	for (i = 0; i < n; i++)
	84	{
	85	acb_randtest(b, state, prec, 1);
	86	acb_mul_2exp_si(b, b, -10 - n_randint(state, prec));
	87	acb_add(acb_mat_entry(X, i, j), acb_mat_entry(X, i, j), b, prec);
	88	}
	89	}
	90
	91	/* Test k = 1 */
	92	if (1)
	93	{
	94	acb_mat_init(R, n, 1);
	95	acb_mat_init(AR, n, 1);
	96	acb_mat_init(Z, n, 1);
	97	acb_mat_init(Z0, n, 1);
	98
	99	for (j = 0; j < n; j++)
	100	{
	101	acb_set(lambda, F + j);
	102	for (i = 0; i < n; i++)
	103	acb_set(acb_mat_entry(R, i, 0), acb_mat_entry(X, i, j));
	104	acb_mat_eig_enclosure_rump(lambda, NULL, R, A, lambda, R, prec);
	105
	106	acb_mat_mul(AR, A, R, prec);
	107	acb_mat_neg(Z, AR);
	108	acb_mat_scalar_addmul_acb(Z, R, lambda, prec);
	109
	110	if (!acb_mat_contains(Z, Z0))
	111	{
	112	flint_printf("FAIL: not containing zero!\n\n");
	113	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	114	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	115	flint_printf("lambda = \n"); acb_printd(lambda, 20); flint_printf("\n\n");
	116	flint_printf("Z = \n"); acb_mat_printd(Z, 20); flint_printf("\n\n");
	117	flint_printf("E = \n");
	118	for (j = 0; j < n; j++)
	119	{
	120	acb_printd(E + j, 20);
	121	flint_printf("\n");
	122	}
	123	flint_abort();
	124	}
	125
	126	found_eigenvalue = 0;
	127	for (j = 0; j < n; j++)
	128	{
	129	if (acb_contains(lambda, E + j))
	130	found_eigenvalue++;
	131	}
	132
	133	if (found_eigenvalue == 0)
	134	{
	135	flint_printf("FAIL: eigenvalue not found\n\n");
	136	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	137	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	138	flint_printf("lambda = \n"); acb_printd(lambda, 20); flint_printf("\n\n");
	139	flint_printf("Z = \n"); acb_mat_printd(Z, 20); flint_printf("\n\n");
	140	flint_printf("E = \n");
	141	for (j = 0; j < n; j++)
	142	{
	143	acb_printd(E + j, 20);
	144	flint_printf("\n");
	145	}
	146	flint_abort();
	147	}
	148	}
	149
	150	acb_mat_clear(R);
	151	acb_mat_clear(AR);
	152	acb_mat_clear(Z);
	153	acb_mat_clear(Z0);
	154	}
	155
	156	/* Test k > 1 */
	157	for (iter2 = 1; iter2 < n; iter2++)
	158	{
	159	k = n_randint(state, n + 1);
	160	k = FLINT_MAX(k, 2);
	161
	162	acb_mat_init(R, n, k);
	163	acb_mat_init(AR, n, k);
	164	acb_mat_init(Z, n, k);
	165	acb_mat_init(Z0, n, k);
	166	acb_mat_init(J, k, k);
	167	acb_mat_init(RJ, n, k);
	168
	169	/* Random selection */
	170	for (h = 0; h < k; h++)
	171	{
	172	j = n_randint(state, n);
	173	if (h == 0 \|\| n_randint(state, 2))
	174	acb_set(lambda, F + j);
	175	for (i = 0; i < n; i++)
	176	acb_set(acb_mat_entry(R, i, h), acb_mat_entry(X, i, j));
	177	}
	178
	179	acb_mat_eig_enclosure_rump(lambda, J, R, A, lambda, R, prec);
	180
	181	/* AY = YJ */
	182	acb_mat_mul(AR, A, R, prec);
	183	acb_mat_mul(RJ, R, J, prec);
	184	acb_mat_sub(Z, AR, RJ, prec);
	185
	186	if (!acb_mat_contains(Z, Z0))
	187	{
	188	flint_printf("FAIL: not containing zero! (k = %wd, prec = %wd)\n\n", k, prec);
	189	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	190	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	191	flint_printf("lambda = \n"); acb_printd(lambda, 20); flint_printf("\n\n");
	192	flint_printf("J = \n"); acb_mat_printd(J, 20); flint_printf("\n\n");
	193	flint_printf("Z = \n"); acb_mat_printd(Z, 20); flint_printf("\n\n");
	194	flint_printf("E = \n");
	195	for (j = 0; j < n; j++)
	196	{
	197	acb_printd(E + j, 20);
	198	flint_printf("\n");
	199	}
	200	flint_abort();
	201	}
	202
	203	found_eigenvalue = 0;
	204	for (j = 0; j < n; j++)
	205	{
	206	if (acb_contains(lambda, E + j))
	207	found_eigenvalue++;
	208	}
	209
	210	if (found_eigenvalue < k)
	211	{
	212	flint_printf("FAIL: eigenvalue not found (k = %wd, found = %wd)\n\n", k, found_eigenvalue);
	213	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	214	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	215	flint_printf("lambda = \n"); acb_printd(lambda, 20); flint_printf("\n\n");
	216	flint_printf("Z = \n"); acb_mat_printd(Z, 20); flint_printf("\n\n");
	217	flint_printf("E = \n");
	218	for (j = 0; j < n; j++)
	219	{
	220	acb_printd(E + j, 20);
	221	flint_printf("\n");
	222	}
	223	flint_abort();
	224	}
	225
	226	acb_mat_clear(R);
	227	acb_mat_clear(AR);
	228	acb_mat_clear(Z);
	229	acb_mat_clear(Z0);
	230	acb_mat_clear(J);
	231	acb_mat_clear(RJ);
	232	}
	233
	234	acb_mat_clear(A);
	235	acb_mat_clear(X);
	236	acb_clear(lambda);
	237	acb_clear(b);
	238	_acb_vec_clear(E, n);
	239	_acb_vec_clear(F, n);
	240	}
	241
	242	flint_randclear(state);
	243	flint_cleanup();
	244	flint_printf("PASS\n");
	245	return EXIT_SUCCESS;
	246	}

+110

-0

acb_mat/test/t-eig_global_enclosure.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("eig_global_enclosure....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 250 * arb_test_multiplier(); iter++)
	24	{
	25	acb_mat_t A, X;
	26	acb_ptr E, F, B;
	27	mag_t eps;
	28	slong i, j, n, prec;
	29	int success;
	30
	31	n = n_randint(state, 8);
	32	prec = 2 + n_randint(state, 100);
	33
	34	acb_mat_init(A, n, n);
	35	acb_mat_init(X, n, n);
	36	E = _acb_vec_init(n);
	37	F = _acb_vec_init(n);
	38	B = _acb_vec_init(n);
	39	mag_init(eps);
	40
	41	for (i = 0; i < n; i++)
	42	acb_randtest(E + i, state, prec, 3);
	43
	44	acb_mat_randtest_eig(A, state, E, prec);
	45	acb_mat_approx_eig_qr(F, NULL, X, A, NULL, 0, prec);
	46
	47	/* perturb F further */
	48	if (n_randint(state, 2))
	49	{
	50	for (i = 0; i < n; i++)
	51	{
	52	acb_randtest(B, state, prec, 1);
	53	acb_mul_2exp_si(B, B, -n_randint(state, prec));
	54	acb_add(F + i, F + i, B, prec);
	55	}
	56	}
	57
	58	acb_mat_eig_global_enclosure(eps, A, F, X, prec);
	59
	60	_acb_vec_set(B, F, n);
	61	for (i = 0; i < n; i++)
	62	acb_add_error_mag(B + i, eps);
	63
	64	for (i = 0; i < n; i++)
	65	{
	66	success = 0;
	67
	68	for (j = 0; j < n; j++)
	69	{
	70	if (acb_contains(B + j, E + i))
	71	{
	72	success = 1;
	73	break;
	74	}
	75	}
	76
	77	if (!success)
	78	{
	79	flint_printf("FAIL\n\n");
	80	flint_printf("A =\n");
	81	acb_mat_printd(A, 20);
	82	flint_printf("\n\ni = %wd\n\n", i);
	83	flint_printf("\nE =\n");
	84	for (j = 0; j < n; j++)
	85	{
	86	acb_printn(E + j, 20, 0); flint_printf("\n");
	87	}
	88	flint_printf("\nB =\n");
	89	for (j = 0; j < n; j++)
	90	{
	91	acb_printn(B + j, 20, 0); flint_printf("\n");
	92	}
	93	flint_abort();
	94	}
	95	}
	96
	97	acb_mat_clear(A);
	98	acb_mat_clear(X);
	99	_acb_vec_clear(E, n);
	100	_acb_vec_clear(F, n);
	101	_acb_vec_clear(B, n);
	102	mag_clear(eps);
	103	}
	104
	105	flint_randclear(state);
	106	flint_cleanup();
	107	flint_printf("PASS\n");
	108	return EXIT_SUCCESS;
	109	}

+238

-0

acb_mat/test/t-eig_multiple.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("eig_multiple....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_mat_t A, R;
	26	acb_ptr E, F;
	27	acb_t b;
	28	slong i, j, n, prec, count, count2;
	29	int result;
	30
	31	n = n_randint(state, 8);
	32	prec = 2 + n_randint(state, 200);
	33
	34	acb_init(b);
	35	acb_mat_init(A, n, n);
	36	acb_mat_init(R, n, n);
	37	E = _acb_vec_init(n);
	38	F = _acb_vec_init(n);
	39
	40	if (n_randint(state, 10) != 0)
	41	{
	42	for (i = 0; i < n; i++)
	43	acb_randtest(E + i, state, prec, 2);
	44	}
	45	else
	46	{
	47	/* Randomly repeat eigenvalues. */
	48	for (i = 0; i < n; i++)
	49	{
	50	if (i == 0 \|\| n_randint(state, 2))
	51	acb_randtest(E + i, state, prec, 2);
	52	else
	53	acb_set(E + i, E + n_randint(state, i));
	54	}
	55	}
	56
	57	if (n_randint(state, 2))
	58	{
	59	for (i = 0; i < n; i++)
	60	acb_get_mid(E + i, E + i);
	61	}
	62
	63	acb_mat_randtest_eig(A, state, E, prec);
	64	acb_mat_approx_eig_qr(F, NULL, R, A, NULL, 0, prec);
	65
	66	/* Perturb F further. */
	67	if (n_randint(state, 10) == 0)
	68	{
	69	for (i = 0; i < n; i++)
	70	{
	71	acb_randtest(b, state, prec, 1);
	72	acb_mul_2exp_si(b, b, -n_randint(state, prec));
	73	acb_add(F + i, F + i, b, prec);
	74	}
	75	}
	76
	77	/* Perturb R further. */
	78	if (n_randint(state, 10) == 0)
	79	{
	80	j = n_randint(state, n);
	81
	82	for (i = 0; i < n; i++)
	83	{
	84	acb_randtest(b, state, prec, 1);
	85	acb_mul_2exp_si(b, b, -10 - n_randint(state, prec));
	86	acb_add(acb_mat_entry(R, i, j), acb_mat_entry(R, i, j), b, prec);
	87	}
	88	}
	89
	90	if (n_randint(state, 2))
	91	result = acb_mat_eig_multiple_rump(F, A, E, R, prec);
	92	else
	93	result = acb_mat_eig_multiple(F, A, E, R, prec);
	94
	95	if (result)
	96	{
	97	count = 0;
	98	count2 = 0;
	99
	100	for (i = 0; i < n; i++)
	101	{
	102	for (j = 0; j < n; j++)
	103	{
	104	if (j == 0 \|\| !acb_equal(F + j, F + j - 1))
	105	count += acb_contains(F + j, E + i);
	106	}
	107
	108	for (j = 0; j < n; j++)
	109	{
	110	if (j == 0 \|\| !acb_equal(F + j, F + j - 1))
	111	count2 += acb_overlaps(F + j, E + i);
	112	}
	113	}
	114
	115	if (count != n \|\| count2 != n)
	116	{
	117	flint_printf("FAIL: count\n\n");
	118	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	119	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	120	flint_printf("count = %wd, count2 = %wd\n\n", count, count2);
	121	flint_printf("E = \n");
	122	for (j = 0; j < n; j++)
	123	{
	124	acb_printd(E + j, 20);
	125	flint_printf("\n");
	126	}
	127	flint_printf("F = \n");
	128	for (j = 0; j < n; j++)
	129	{
	130	acb_printd(F + j, 20);
	131	flint_printf("\n");
	132	}
	133	flint_abort();
	134	}
	135	}
	136
	137	acb_mat_clear(A);
	138	acb_mat_clear(R);
	139	_acb_vec_clear(E, n);
	140	_acb_vec_clear(F, n);
	141	acb_clear(b);
	142	}
	143
	144	/* Test convergence for DFT matrices */
	145	for (iter = 0; iter < 50 * arb_test_multiplier(); iter++)
	146	{
	147	acb_mat_t A, R, QC;
	148	acb_ptr E;
	149	acb_t t;
	150	fmpq_mat_t Q, Qinv;
	151	slong i, n, c0, c1, c2, c3;
	152	slong prec;
	153	int algorithm, result;
	154
	155	n = n_randint(state, 30);
	156	algorithm = n_randint(state, 2);
	157	acb_mat_init(A, n, n);
	158	acb_mat_init(R, n, n);
	159	E = _acb_vec_init(n);
	160	acb_init(t);
	161	acb_mat_init(QC, n, n);
	162	fmpq_mat_init(Q, n, n);
	163	fmpq_mat_init(Qinv, n, n);
	164
	165	/* The current algorithm is not robust enough. */
	166	#if 0
	167	do {
	168	fmpq_mat_randtest(Q, state, 2 + n_randint(state, 10));
	169	} while (!fmpq_mat_inv(Qinv, Q));
	170	#else
	171	fmpq_mat_one(Q);
	172	fmpq_mat_one(Qinv);
	173	#endif
	174
	175	for (prec = 32; ; prec *= 2)
	176	{
	177	if (prec > 10000)
	178	{
	179	flint_printf("FAIL: unsuccessful, prec > 10000\n\n");
	180	flint_printf("algorithm = %d, iter %wd\n\n", algorithm, iter);
	181	flint_abort();
	182	}
	183
	184	acb_mat_dft(A, 0, prec);
	185
	186	#if 0
	187	acb_mat_set_fmpq_mat(QC, Q, prec);
	188	acb_mat_mul(A, A, QC, prec);
	189	acb_mat_set_fmpq_mat(QC, Qinv, prec);
	190	acb_mat_mul(A, QC, A, prec);
	191	#endif
	192
	193	acb_mat_approx_eig_qr(E, NULL, R, A, NULL, 0, prec);
	194
	195	if (algorithm == 0)
	196	result = acb_mat_eig_multiple_rump(E, A, E, R, prec);
	197	else
	198	result = acb_mat_eig_multiple(E, A, E, R, prec);
	199
	200	/* Verify the known eigenvalues + multiplicities */
	201	if (result)
	202	{
	203	c0 = c1 = c2 = c3 = 0;
	204
	205	for (i = 0; i < n; i++)
	206	{
	207	acb_set_d_d(t, 1.0, 0.0);
	208	c0 += acb_contains(E + i, t);
	209	acb_set_d_d(t, -1.0, 0.0);
	210	c1 += acb_contains(E + i, t);
	211	acb_set_d_d(t, 0.0, 1.0);
	212	c2 += acb_contains(E + i, t);
	213	acb_set_d_d(t, 0.0, -1.0);
	214	c3 += acb_contains(E + i, t);
	215	}
	216
	217	result = (n == 0 \|\| (c0 == (n+4)/4 && c1 == (n+2)/4 && c2 == (n-1)/4 && c3 == (n+1)/4));
	218	}
	219
	220	if (result)
	221	break;
	222	}
	223
	224	acb_mat_clear(A);
	225	acb_mat_clear(R);
	226	acb_mat_clear(QC);
	227	_acb_vec_clear(E, n);
	228	acb_clear(t);
	229	fmpq_mat_clear(Q);
	230	fmpq_mat_clear(Qinv);
	231	}
	232
	233	flint_randclear(state);
	234	flint_cleanup();
	235	flint_printf("PASS\n");
	236	return EXIT_SUCCESS;
	237	}

+280

-0

acb_mat/test/t-eig_simple.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_mat.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("eig_simple....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 5000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_mat_t A, L, R, LAR, D;
	26	acb_ptr E, F;
	27	acb_t b;
	28	slong i, j, n, prec, count, count2;
	29	int result, algorithm;
	30
	31	n = n_randint(state, 8);
	32	prec = 2 + n_randint(state, 200);
	33	algorithm = n_randint(state, 3);
	34
	35	acb_mat_init(A, n, n);
	36	acb_mat_init(L, n, n);
	37	acb_mat_init(R, n, n);
	38	acb_mat_init(LAR, n, n);
	39	acb_mat_init(D, n, n);
	40	acb_init(b);
	41	E = _acb_vec_init(n);
	42	F = _acb_vec_init(n);
	43
	44	if (n_randint(state, 10) != 0)
	45	{
	46	for (i = 0; i < n; i++)
	47	acb_randtest(E + i, state, prec, 2);
	48	}
	49	else
	50	{
	51	/* Randomly repeat eigenvalues. */
	52	for (i = 0; i < n; i++)
	53	{
	54	if (i == 0 \|\| n_randint(state, 2))
	55	acb_randtest(E + i, state, prec, 2);
	56	else
	57	acb_set(E + i, E + n_randint(state, i));
	58	}
	59	}
	60
	61	if (n_randint(state, 2))
	62	{
	63	for (i = 0; i < n; i++)
	64	acb_get_mid(E + i, E + i);
	65	}
	66
	67	acb_mat_randtest_eig(A, state, E, prec);
	68	acb_mat_approx_eig_qr(F, NULL, R, A, NULL, 0, prec);
	69
	70	/* Perturb F further. */
	71	if (n_randint(state, 10) == 0)
	72	{
	73	for (i = 0; i < n; i++)
	74	{
	75	acb_randtest(b, state, prec, 1);
	76	acb_mul_2exp_si(b, b, -n_randint(state, prec));
	77	acb_add(F + i, F + i, b, prec);
	78	}
	79	}
	80
	81	/* Perturb R further. */
	82	if (n_randint(state, 10) == 0)
	83	{
	84	j = n_randint(state, n);
	85
	86	for (i = 0; i < n; i++)
	87	{
	88	acb_randtest(b, state, prec, 1);
	89	acb_mul_2exp_si(b, b, -10 - n_randint(state, prec));
	90	acb_add(acb_mat_entry(R, i, j), acb_mat_entry(R, i, j), b, prec);
	91	}
	92	}
	93
	94	if (algorithm == 0)
	95	result = acb_mat_eig_simple(F, L, R, A, E, R, prec);
	96	else if (algorithm == 1)
	97	result = acb_mat_eig_simple_rump(F, L, R, A, E, R, prec);
	98	else
	99	result = acb_mat_eig_simple_vdhoeven_mourrain(F, L, R, A, E, R, prec);
	100
	101	acb_mat_mul(LAR, L, A, prec);
	102	acb_mat_mul(LAR, LAR, R, prec);
	103	for (i = 0; i < n; i++)
	104	acb_set(acb_mat_entry(D, i, i), F + i);
	105
	106	if (!acb_mat_overlaps(LAR, D))
	107	{
	108	flint_printf("FAIL: overlap\n\n");
	109	flint_printf("algorithm = %d\n\n", algorithm);
	110	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	111	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	112	flint_printf("L = \n"); acb_mat_printd(L, 20); flint_printf("\n\n");
	113	flint_printf("D = \n"); acb_mat_printd(D, 20); flint_printf("\n\n");
	114	flint_printf("LAR = \n"); acb_mat_printd(LAR, 20); flint_printf("\n\n");
	115	flint_abort();
	116	}
	117
	118	if (result)
	119	{
	120	for (i = 0; i < n; i++)
	121	{
	122	count = 0;
	123	for (j = 0; j < n; j++)
	124	count += acb_contains(F + i, E + j);
	125
	126	count2 = 0;
	127	for (j = 0; j < n; j++)
	128	count2 += acb_overlaps(F + i, E + j);
	129
	130	if (count != 1 \|\| count2 != 1)
	131	{
	132	flint_printf("FAIL: count\n\n");
	133	flint_printf("algorithm = %d\n\n", algorithm);
	134	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	135	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	136	flint_printf("L = \n"); acb_mat_printd(L, 20); flint_printf("\n\n");
	137	flint_printf("D = \n"); acb_mat_printd(D, 20); flint_printf("\n\n");
	138	flint_printf("LAR = \n"); acb_mat_printd(LAR, 20); flint_printf("\n\n");
	139	flint_printf("i = %wd, count = %wd, count2 = %wd\n\n", i, count, count2);
	140	flint_abort();
	141	}
	142	}
	143	}
	144
	145	acb_mat_clear(A);
	146	acb_mat_clear(L);
	147	acb_mat_clear(R);
	148	acb_mat_clear(LAR);
	149	acb_mat_clear(D);
	150	acb_clear(b);
	151	_acb_vec_clear(E, n);
	152	_acb_vec_clear(F, n);
	153	}
	154
	155	/* Test convergence, given companion matrices */
	156	for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
	157	{
	158	acb_mat_t A, R, QC;
	159	acb_ptr E;
	160	acb_ptr roots;
	161	fmpq_mat_t Q, Qinv;
	162	acb_poly_t f;
	163	slong i, j, n, prec, count, count2;
	164	int algorithm, success;
	165
	166	algorithm = n_randint(state, 3);
	167	n = n_randint(state, 10);
	168	roots = _acb_vec_init(n);
	169	E = _acb_vec_init(n);
	170	acb_poly_init(f);
	171	acb_mat_init(A, n, n);
	172	acb_mat_init(R, n, n);
	173	fmpq_mat_init(Q, n, n);
	174	fmpq_mat_init(Qinv, n, n);
	175	acb_mat_init(QC, n, n);
	176
	177	for (i = 0; i < n; i++)
	178	{
	179	new_root:
	180	acb_randtest(roots + i, state, 2 + n_randint(state, 100), 4);
	181	acb_get_mid(roots + i, roots + i);
	182
	183	for (j = 0; j < i; j++)
	184	if (acb_equal(roots + i, roots + j))
	185	goto new_root;
	186	}
	187
	188	do {
	189	fmpq_mat_randtest(Q, state, 2 + n_randint(state, 100));
	190	} while (!fmpq_mat_inv(Qinv, Q));
	191
	192	success = 0;
	193
	194	for (prec = 32; !success; prec *= 2)
	195	{
	196	if (prec > 10000)
	197	{
	198	flint_printf("FAIL: unsuccessful, prec > 10000\n\n");
	199	flint_printf("algorithm = %d, iter %wd\n\n", algorithm, iter);
	200	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	201	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	202	flint_printf("roots = \n");
	203	for (j = 0; j < n; j++)
	204	{
	205	acb_printd(roots + j, 20);
	206	flint_printf("\n");
	207	}
	208	flint_abort();
	209	}
	210
	211	acb_poly_product_roots(f, roots, n, prec);
	212
	213	acb_mat_companion(A, f, prec);
	214	acb_mat_set_fmpq_mat(QC, Q, prec);
	215	acb_mat_mul(A, A, QC, prec);
	216	acb_mat_set_fmpq_mat(QC, Qinv, prec);
	217	acb_mat_mul(A, QC, A, prec);
	218
	219	acb_mat_approx_eig_qr(E, NULL, R, A, NULL, 0, prec);
	220
	221	if (algorithm == 0)
	222	success = acb_mat_eig_simple(E, NULL, NULL, A, E, R, prec);
	223	else if (algorithm == 1)
	224	success = acb_mat_eig_simple_rump(E, NULL, NULL, A, E, R, prec);
	225	else
	226	success = acb_mat_eig_simple_vdhoeven_mourrain(E, NULL, NULL, A, E, R, prec);
	227
	228	if (success)
	229	{
	230	for (i = 0; i < n; i++)
	231	{
	232	count = 0;
	233	for (j = 0; j < n; j++)
	234	count += acb_contains(E + i, roots + j);
	235
	236	count2 = 0;
	237	for (j = 0; j < n; j++)
	238	count2 += acb_overlaps(E + i, roots + j);
	239
	240	if (count != 1 \|\| count2 != 1)
	241	{
	242	flint_printf("FAIL: count\n\n");
	243	flint_printf("algorithm = %d\n\n", algorithm);
	244	flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
	245	flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
	246	flint_printf("i = %wd, count = %wd, count2 = %wd\n\n", i, count, count2);
	247	flint_printf("roots = \n");
	248	for (j = 0; j < n; j++)
	249	{
	250	acb_printd(roots + j, 20);
	251	flint_printf("\n");
	252	}
	253	flint_printf("E = \n");
	254	for (j = 0; j < n; j++)
	255	{
	256	acb_printd(E + j, 20);
	257	flint_printf("\n");
	258	}
	259	flint_abort();
	260	}
	261	}
	262	}
	263	}
	264
	265	fmpq_mat_clear(Q);
	266	fmpq_mat_clear(Qinv);
	267	acb_mat_clear(QC);
	268	acb_mat_clear(A);
	269	acb_mat_clear(R);
	270	acb_poly_clear(f);
	271	_acb_vec_clear(roots, n);
	272	_acb_vec_clear(E, n);
	273	}
	274
	275	flint_randclear(state);
	276	flint_cleanup();
	277	flint_printf("PASS\n");
	278	return EXIT_SUCCESS;
	279	}

+44

-3

acb_mat.h less more

93	93
94	94	void acb_mat_randtest(acb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits);
95	95
	96	void acb_mat_randtest_eig(acb_mat_t A, flint_rand_t state, acb_srcptr E, slong prec);
	97
96	98	/* I/O */
97	99
98	100	void acb_mat_fprintd(FILE * file, const acb_mat_t mat, slong digits);

131	133	acb_mat_is_square(const acb_mat_t mat)
132	134	{
133	135	return (mat->r == mat->c);
	136	}
	137
	138	int acb_mat_is_exact(const acb_mat_t mat);
	139
	140	int acb_mat_is_zero(const acb_mat_t mat);
	141	int acb_mat_is_finite(const acb_mat_t mat);
	142	int acb_mat_is_triu(const acb_mat_t mat);
	143	int acb_mat_is_tril(const acb_mat_t mat);
	144
	145	ACB_MAT_INLINE int
	146	acb_mat_is_diag(const acb_mat_t mat)
	147	{
	148	return acb_mat_is_tril(mat) && acb_mat_is_triu(mat);
134	149	}
135	150
136	151	/* Radius and interval operations */

167	182	void acb_mat_one(acb_mat_t mat);
168	183
169	184	void acb_mat_ones(acb_mat_t mat);
	185
	186	void acb_mat_indeterminate(acb_mat_t mat);
170	187
171	188	void acb_mat_dft(acb_mat_t res, int kind, slong prec);
172	189

394	411	int acb_mat_approx_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec);
395	412	void acb_mat_approx_solve_lu_precomp(acb_mat_t X, const slong * perm, const acb_mat_t A, const acb_mat_t B, slong prec);
396	413	int acb_mat_approx_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec);
	414	int acb_mat_approx_inv(acb_mat_t X, const acb_mat_t A, slong prec);
397	415
398	416	int acb_mat_inv(acb_mat_t X, const acb_mat_t A, slong prec);
399	417

401	419	void acb_mat_det_precond(acb_t det, const acb_mat_t A, slong prec);
402	420	void acb_mat_det(acb_t det, const acb_mat_t A, slong prec);
403	421
	422	/* Eigenvalues and eigenvectors */
	423
	424	int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec);
	425
	426	void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec);
	427
	428	void acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t X, const acb_mat_t A,
	429	const acb_t lambda_approx, const acb_mat_t X_approx, slong prec);
	430
	431	int acb_mat_eig_simple_rump(acb_ptr E, acb_mat_t L, acb_mat_t R,
	432	const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
	433	int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, acb_mat_t L, acb_mat_t R,
	434	const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
	435	int acb_mat_eig_simple(acb_ptr E, acb_mat_t L, acb_mat_t R,
	436	const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
	437
	438	int acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
	439	int acb_mat_eig_multiple(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
	440
404	441	/* Special functions */
405	442
406	443	void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec);
407	444
408	445	void acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec);
409	446
410		void _acb_mat_charpoly(acb_ptr cp, const acb_mat_t mat, slong prec);
411
412		void acb_mat_charpoly(acb_poly_t cp, const acb_mat_t mat, slong prec);
	447	void _acb_mat_charpoly(acb_ptr poly, const acb_mat_t mat, slong prec);
	448	void acb_mat_charpoly(acb_poly_t poly, const acb_mat_t mat, slong prec);
	449	void _acb_mat_companion(acb_mat_t mat, acb_srcptr poly, slong prec);
	450	void acb_mat_companion(acb_mat_t mat, const acb_poly_t poly, slong prec);
413	451
414	452	void acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec);
	453
	454	void _acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong a, slong b, slong prec);
	455	void acb_mat_diag_prod(acb_t res, const acb_mat_t A, slong prec);
415	456
416	457	ACB_MAT_INLINE slong
417	458	acb_mat_allocated_bytes(const acb_mat_t x)

-1

arb/version.c less more

10	10
11	11	#include "arb.h"
12	12
13		const char * arb_version = "2.15.1";
	13	const char * arb_version = "2.16.0";
14	14

+25

-3

arb_mat/addmul_rad_mag_fast.c less more

20	20	slack for accumulated sums. */
21	21	#define DOUBLE_MAX_OFFSET 900
22	22
	23	static __inline__ double dot8(const double * A, const double * B)
	24	{
	25	return ((A[0] * B[0] + A[1] * B[1]) + (A[2] * B[2] + A[3] * B[3])) +
	26	((A[4] * B[4] + A[5] * B[5]) + (A[6] * B[6] + A[7] * B[7]));
	27	}
	28
23	29	/* Upper bound of matrix product, assuming nonnegative entries and
24	30	no overflow/underflow. B is pre-transposed. Straightforward blocked
25	31	implementation; could use BLAS, but this matrix product is rarely going
26	32	to be the bottleneck. */
	33
27	34	static void
28	35	_d_mat_addmul(double * C, const double * A, const double * B, slong ar, slong ac, slong bc)
29	36	{

42	49	{
43	50	for (j = jj; j < FLINT_MIN(jj + BLOCK_SIZE, bc); j++)
44	51	{
45		t = 0.0;
46		for (k = kk; k < FLINT_MIN(kk + BLOCK_SIZE, ac); k++)
47		t += A[i * ac + k] * B[j * ac + k];
	52	if (BLOCK_SIZE == 32 && kk + BLOCK_SIZE <= ac)
	53	{
	54	double t0, t1, t2, t3;
	55
	56	t0 = dot8(A + i * ac + kk + 0, B + j * ac + kk + 0);
	57	t1 = dot8(A + i * ac + kk + 8, B + j * ac + kk + 8);
	58	t2 = dot8(A + i * ac + kk + 16, B + j * ac + kk + 16);
	59	t3 = dot8(A + i * ac + kk + 24, B + j * ac + kk + 24);
	60
	61	t = (t0 + t1) + (t2 + t3);
	62	}
	63	else
	64	{
	65	t = 0.0;
	66
	67	for (k = kk; k < FLINT_MIN(kk + BLOCK_SIZE, ac); k++)
	68	t += A[i * ac + k] * B[j * ac + k];
	69	}
48	70
49	71	C[i * bc + j] += t;
50	72	}

+30

-0

arb_mat/approx_inv.c less more

	0	/*
	1	Copyright (C) 2012 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec)
	15	{
	16	if (X == A)
	17	{
	18	int r;
	19	arb_mat_t T;
	20	arb_mat_init(T, arb_mat_nrows(A), arb_mat_ncols(A));
	21	r = arb_mat_approx_inv(T, A, prec);
	22	arb_mat_swap(T, X);
	23	arb_mat_clear(T);
	24	return r;
	25	}
	26
	27	arb_mat_one(X);
	28	return arb_mat_approx_solve(X, A, X, prec);
	29	}

-13

arb_mat/approx_mul.c less more

79	79	}
80	80	}
81	81
82		int
83		arb_mat_is_exact(const arb_mat_t A)
84		{
85		slong i, j;
86
87		for (i = 0; i < arb_mat_nrows(A); i++)
88		for (j = 0; j < arb_mat_ncols(A); j++)
89		if (!mag_is_zero(arb_radref(arb_mat_entry(A, i, j))))
90		return 0;
91
92		return 1;
93		}
94
95	82	void
96	83	arb_mat_approx_mul(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)
97	84	{

+49

-0

arb_mat/companion.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	void
	14	_arb_mat_companion(arb_mat_t A, arb_srcptr poly, slong prec)
	15	{
	16	slong i, j, n;
	17	arb_t c;
	18
	19	n = arb_mat_nrows(A);
	20
	21	if (n == 0)
	22	return;
	23
	24	for (i = 0; i < n - 1; i++)
	25	for (j = 0; j < n; j++)
	26	arb_set_ui(arb_mat_entry(A, i, j), (i + 1) == j);
	27
	28	arb_init(c);
	29	arb_inv(c, poly + n, prec);
	30	arb_neg(c, c);
	31	for (j = 0; j < n; j++)
	32	arb_mul(arb_mat_entry(A, n - 1, j), poly + j, c, prec);
	33	arb_clear(c);
	34	}
	35
	36	void
	37	arb_mat_companion(arb_mat_t A, const arb_poly_t poly, slong prec)
	38	{
	39	slong n = arb_mat_nrows(A);
	40
	41	if (n != arb_poly_degree(poly) \|\| n != arb_mat_ncols(A))
	42	{
	43	flint_printf("arb_mat_companion: incompatible dimensions!\n");
	44	flint_abort();
	45	}
	46
	47	_arb_mat_companion(A, poly->coeffs, prec);
	48	}

-52

arb_mat/det.c less more

38	38	arb_clear(a);
39	39	}
40	40
41		int
42		arb_mat_is_finite(const arb_mat_t A)
43		{
44		slong i, j, n, m;
45
46		n = arb_mat_nrows(A);
47		m = arb_mat_ncols(A);
48
49		for (i = 0; i < n; i++)
50		for (j = 0; j < m; j++)
51		if (!arb_is_finite(arb_mat_entry(A, i, j)))
52		return 0;
53
54		return 1;
55		}
56
57		int
58		arb_mat_is_triu(const arb_mat_t A)
59		{
60		slong i, j, n, m;
61
62		n = arb_mat_nrows(A);
63		m = arb_mat_ncols(A);
64
65		for (i = 1; i < n; i++)
66		for (j = 0; j < FLINT_MIN(i, m); j++)
67		if (!arb_is_zero(arb_mat_entry(A, i, j)))
68		return 0;
69
70		return 1;
71		}
72
73		int
74		arb_mat_is_tril(const arb_mat_t A)
75		{
76		slong i, j, n, m;
77
78		n = arb_mat_nrows(A);
79		m = arb_mat_ncols(A);
80
81		for (i = 0; i < n; i++)
82		for (j = i + 1; j < m; j++)
83		if (!arb_is_zero(arb_mat_entry(A, i, j)))
84		return 0;
85
86		return 1;
87		}
88
89	41	void
90	42	arb_mat_det(arb_t det, const arb_mat_t A, slong prec)
91	43	{
92		slong k, n;
	44	slong n;
93	45
94	46	if (!arb_mat_is_square(A))
95	47	{

117	69	}
118	70	else if (arb_mat_is_tril(A) \|\| arb_mat_is_triu(A))
119	71	{
120		arb_set(det, arb_mat_entry(A, 0, 0));
121		for (k = 1; k < n; k++)
122		arb_mul(det, det, arb_mat_entry(A, k, k), prec);
	72	arb_mat_diag_prod(det, A, prec);
123	73	}
124	74	else if (n == 3)
125	75	{

-1

arb_mat/det_lu.c less more

139	139	arb_mat_clear(T);
140	140	}
141	141	}
142

-5

arb_mat/det_precond.c less more

77	77	{
78	78	arb_mat_t LU, Linv, Uinv;
79	79	arb_t detU;
80		slong i, n;
	80	slong n;
81	81	slong *P;
82	82
83	83	n = arb_mat_nrows(A);

108	108	arb_mat_one(Uinv);
109	109	arb_mat_approx_solve_triu(Uinv, LU, Uinv, 0, prec);
110	110
111		arb_set(detU, arb_mat_entry(Uinv, 0, 0));
112		for (i = 1; i < n; i++)
113		arb_mul(detU, detU, arb_mat_entry(Uinv, i, i), prec);
	111	arb_mat_diag_prod(detU, Uinv, prec);
114	112
115	113	arb_mat_mul(LU, A, Uinv, prec);
116	114	_apply_permutation(LU, P, n);

140	138	_perm_clear(P);
141	139	arb_mat_clear(LU);
142	140	}
143

+55

-0

arb_mat/diag_prod.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	void
	14	_arb_mat_diag_prod(arb_t res, const arb_mat_t A, slong a, slong b, slong prec)
	15	{
	16	if (b - a == 0)
	17	{
	18	arb_one(res);
	19	}
	20	else if (b - a == 1)
	21	{
	22	arb_set_round(res, arb_mat_entry(A, a, a), prec);
	23	}
	24	else
	25	{
	26	slong i;
	27
	28	arb_mul(res, arb_mat_entry(A, a, a), arb_mat_entry(A, a + 1, a + 1), prec);
	29	for (i = a + 2; i < b; i++)
	30	arb_mul(res, res, arb_mat_entry(A, i, i), prec);
	31
	32	/* no advantage? */
	33
	34	/*
	35	arb_t t;
	36	arb_init(t);
	37	_arb_mat_diag_prod(t, A, a, a + (b - a) / 2, prec);
	38	_arb_mat_diag_prod(res, A, a + (b - a) / 2, b, prec);
	39	arb_mul(res, res, t, prec);
	40	arb_clear(t);
	41	*/
	42	}
	43	}
	44
	45	void
	46	arb_mat_diag_prod(arb_t res, const arb_mat_t A, slong prec)
	47	{
	48	slong m, n;
	49
	50	m = arb_mat_nrows(A);
	51	n = arb_mat_nrows(A);
	52
	53	_arb_mat_diag_prod(res, A, 0, FLINT_MIN(m, n), prec);
	54	}

+22

-0

arb_mat/indeterminate.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	void
	14	arb_mat_indeterminate(arb_mat_t A)
	15	{
	16	slong i, j;
	17
	18	for (i = 0; i < arb_mat_nrows(A); i++)
	19	for (j = 0; j < arb_mat_ncols(A); j++)
	20	arb_indeterminate(arb_mat_entry(A, i, j));
	21	}

+25

-0

arb_mat/is_exact.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	arb_mat_is_exact(const arb_mat_t A)
	15	{
	16	slong i, j;
	17
	18	for (i = 0; i < arb_mat_nrows(A); i++)
	19	for (j = 0; j < arb_mat_ncols(A); j++)
	20	if (!mag_is_zero(arb_radref(arb_mat_entry(A, i, j))))
	21	return 0;
	22
	23	return 1;
	24	}

+28

-0

arb_mat/is_finite.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	arb_mat_is_finite(const arb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = arb_mat_nrows(A);
	19	m = arb_mat_ncols(A);
	20
	21	for (i = 0; i < n; i++)
	22	for (j = 0; j < m; j++)
	23	if (!arb_is_finite(arb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+28

-0

arb_mat/is_tril.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	arb_mat_is_tril(const arb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = arb_mat_nrows(A);
	19	m = arb_mat_ncols(A);
	20
	21	for (i = 0; i < n; i++)
	22	for (j = i + 1; j < m; j++)
	23	if (!arb_is_zero(arb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+28

-0

arb_mat/is_triu.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	arb_mat_is_triu(const arb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = arb_mat_nrows(A);
	19	m = arb_mat_ncols(A);
	20
	21	for (i = 1; i < n; i++)
	22	for (j = 0; j < FLINT_MIN(i, m); j++)
	23	if (!arb_is_zero(arb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+28

-0

arb_mat/is_zero.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	arb_mat_is_zero(const arb_mat_t A)
	15	{
	16	slong i, j, n, m;
	17
	18	n = arb_mat_nrows(A);
	19	m = arb_mat_ncols(A);
	20
	21	for (i = 0; i < n; i++)
	22	for (j = 0; j < m; j++)
	23	if (!arb_is_zero(arb_mat_entry(A, i, j)))
	24	return 0;
	25
	26	return 1;
	27	}

+152

-146

arb_mat/mul_block.c less more

29	29	return 1;
30	30	}
31	31
32		int arb_mat_is_zero(const arb_mat_t A)
	32	/* allow changing this from the test code */
	33	ARB_DLL slong arb_mat_mul_block_min_block_size = 0;
	34
	35	void
	36	arb_mat_mid_addmul_block_fallback(arb_mat_t C,
	37	const arb_mat_t A, const arb_mat_t B,
	38	slong block_start,
	39	slong block_end,
	40	slong prec)
33	41	{
34		slong i, j, M, N;
	42	slong M, P, n;
	43	slong i, j;
	44	arb_ptr tmpA, tmpB;
35	45
36	46	M = arb_mat_nrows(A);
37		N = arb_mat_ncols(A);
	47	P = arb_mat_ncols(B);
	48
	49	n = block_end - block_start;
	50
	51	tmpA = flint_malloc(sizeof(arb_struct) * (M * n + P * n));
	52	tmpB = tmpA + M * n;
38	53
39	54	for (i = 0; i < M; i++)
40	55	{
41		for (j = 0; j < N; j++)
42		{
43		if (!arb_is_zero(arb_mat_entry(A, i, j)))
44		return 0;
45		}
46		}
47
48		return 1;
	56	for (j = 0; j < n; j++)
	57	{
	58	arb_midref(tmpA + i n + j) = *arb_midref(arb_mat_entry(A, i, block_start + j));
	59	mag_init(arb_radref(tmpA + i * n + j));
	60	}
	61	}
	62
	63	for (i = 0; i < P; i++)
	64	{
	65	for (j = 0; j < n; j++)
	66	{
	67	arb_midref(tmpB + i n + j) = *arb_midref(arb_mat_entry(B, block_start + j, i));
	68	mag_init(arb_radref(tmpB + i * n + j));
	69	}
	70	}
	71
	72	for (i = 0; i < M; i++)
	73	{
	74	for (j = 0; j < P; j++)
	75	{
	76	arb_dot(arb_mat_entry(C, i, j),
	77	(block_start == 0) ? NULL : arb_mat_entry(C, i, j), 0,
	78	tmpA + i * n, 1, tmpB + j * n, 1, n, prec);
	79	}
	80	}
	81
	82	flint_free(tmpA);
49	83	}
50
51		#define MIN_BLOCK_SIZE 8
52	84
53	85	void
54	86	arb_mat_mid_addmul_block_prescaled(arb_mat_t C,

60	92	slong prec)
61	93	{
62	94	slong M, P, n;
63		slong i, j, k;
	95	slong i, j;
	96	slong M0, M1, P0, P1, Mstep, Pstep;
64	97	int inexact;
65	98
66	99	/* flint_printf("block mul from %wd to %wd\n", block_start, block_end); */

70	103
71	104	n = block_end - block_start;
72	105
73		if (n < MIN_BLOCK_SIZE)
74		{
75		/* Fall back to naive multiplication of the midpoints. */
76		for (i = 0; i < M; i++)
77		{
78		for (j = 0; j < P; j++)
79		{
80		for (k = 0; k < n; k++)
81		{
82		arb_ptr Cij = arb_mat_entry(C, i, j);
83
84		/* are we writing this Cij for the first time? */
85		if (block_start + k == 0)
	106	/* Create sub-blocks to keep matrices nearly square. Necessary? */
	107	#if 1
	108	Mstep = (M < 2 * n) ? M : n;
	109	Pstep = (P < 2 * n) ? P : n;
	110	#else
	111	Mstep = M;
	112	Pstep = P;
	113	#endif
	114
	115	for (M0 = 0; M0 < M; M0 += Mstep)
	116	{
	117	for (P0 = 0; P0 < P; P0 += Pstep)
	118	{
	119	fmpz_mat_t AA, BB, CC;
	120	arb_t t;
	121	fmpz_t e;
	122
	123	M1 = FLINT_MIN(M0 + Mstep, M);
	124	P1 = FLINT_MIN(P0 + Pstep, P);
	125
	126	fmpz_mat_init(AA, M1 - M0, n);
	127	fmpz_mat_init(BB, n, P1 - P0);
	128	fmpz_mat_init(CC, M1 - M0, P1 - P0);
	129
	130	/* Convert to fixed-point matrices. */
	131	for (i = M0; i < M1; i++)
	132	{
	133	if (A_min[i] == WORD_MIN) /* only zeros in this row */
	134	continue;
	135
	136	for (j = 0; j < n; j++)
	137	{
	138	inexact = arf_get_fmpz_fixed_si(fmpz_mat_entry(AA, i - M0, j),
	139	arb_midref(arb_mat_entry(A, i, block_start + j)), A_min[i]);
	140
	141	if (inexact)
86	142	{
87		inexact = arf_mul(arb_midref(Cij),
88		arb_midref(arb_mat_entry(A, i, block_start + k)),
89		arb_midref(arb_mat_entry(B, block_start + k, j)), prec, ARB_RND);
90
91		/* important: must not use unsafe fast mag method here because C might not be demoted */
92		if (inexact)
93		arf_mag_set_ulp(arb_radref(Cij), arb_midref(Cij), prec);
94		else
95		mag_zero(arb_radref(Cij));
	143	flint_printf("matrix multiplication: bad exponent!\n");
	144	flint_abort();
	145	}
	146	}
	147	}
	148
	149	for (i = P0; i < P1; i++)
	150	{
	151	if (B_min[i] == WORD_MIN) /* only zeros in this column */
	152	continue;
	153
	154	for (j = 0; j < n; j++)
	155	{
	156	inexact = arf_get_fmpz_fixed_si(fmpz_mat_entry(BB, j, i - P0),
	157	arb_midref(arb_mat_entry(B, block_start + j, i)), B_min[i]);
	158
	159	if (inexact)
	160	{
	161	flint_printf("matrix multiplication: bad exponent!\n");
	162	flint_abort();
	163	}
	164	}
	165	}
	166
	167	/* The main multiplication */
	168	fmpz_mat_mul(CC, AA, BB);
	169	/* flint_printf("bits %wd %wd %wd\n", fmpz_mat_max_bits(CC),
	170	fmpz_mat_max_bits(AA), fmpz_mat_max_bits(BB)); */
	171
	172	fmpz_mat_clear(AA);
	173	fmpz_mat_clear(BB);
	174
	175	arb_init(t);
	176
	177	/* Add to the result matrix */
	178	for (i = M0; i < M1; i++)
	179	{
	180	for (j = P0; j < P1; j++)
	181	{
	182	*e = A_min[i] + B_min[j];
	183
	184	/* The first time we write this Cij */
	185	if (block_start == 0)
	186	{
	187	arb_set_round_fmpz_2exp(arb_mat_entry(C, i, j),
	188	fmpz_mat_entry(CC, i - M0, j - P0), e, prec);
96	189	}
97	190	else
98	191	{
99		inexact = arf_addmul(arb_midref(Cij),
100		arb_midref(arb_mat_entry(A, i, block_start + k)),
101		arb_midref(arb_mat_entry(B, block_start + k, j)), prec, ARB_RND);
102
103		/* here the unsafe method is ok */
104		if (inexact)
105		arf_mag_fast_add_ulp(arb_radref(Cij), arb_radref(Cij), arb_midref(Cij), prec);
	192	arb_set_round_fmpz_2exp(t, fmpz_mat_entry(CC, i - M0, j - P0), e, prec);
	193	arb_add(arb_mat_entry(C, i, j), arb_mat_entry(C, i, j), t, prec);
106	194	}
107	195	}
108	196	}
109		}
110		}
111		else /* Multiply using fmpz_mat_mul */
112		{
113		slong M0, M1, P0, P1, Mstep, Pstep;
114
115		/* Create sub-blocks to keep matrices nearly square. Necessary? */
116		#if 1
117		Mstep = (M < 2 * n) ? M : n;
118		Pstep = (P < 2 * n) ? P : n;
119		#else
120		Mstep = M;
121		Pstep = P;
122		#endif
123
124		for (M0 = 0; M0 < M; M0 += Mstep)
125		{
126		for (P0 = 0; P0 < P; P0 += Pstep)
127		{
128		fmpz_mat_t AA, BB, CC;
129		arb_t t;
130		fmpz_t e;
131
132		M1 = FLINT_MIN(M0 + Mstep, M);
133		P1 = FLINT_MIN(P0 + Pstep, P);
134
135		fmpz_mat_init(AA, M1 - M0, n);
136		fmpz_mat_init(BB, n, P1 - P0);
137		fmpz_mat_init(CC, M1 - M0, P1 - P0);
138
139		/* Convert to fixed-point matrices. */
140		for (i = M0; i < M1; i++)
141		{
142		if (A_min[i] == WORD_MIN) /* only zeros in this row */
143		continue;
144
145		for (j = 0; j < n; j++)
146		{
147		inexact = arf_get_fmpz_fixed_si(fmpz_mat_entry(AA, i - M0, j),
148		arb_midref(arb_mat_entry(A, i, block_start + j)), A_min[i]);
149
150		if (inexact)
151		{
152		flint_printf("matrix multiplication: bad exponent!\n");
153		flint_abort();
154		}
155		}
156		}
157
158		for (i = P0; i < P1; i++)
159		{
160		if (B_min[i] == WORD_MIN) /* only zeros in this column */
161		continue;
162
163		for (j = 0; j < n; j++)
164		{
165		inexact = arf_get_fmpz_fixed_si(fmpz_mat_entry(BB, j, i - P0),
166		arb_midref(arb_mat_entry(B, block_start + j, i)), B_min[i]);
167
168		if (inexact)
169		{
170		flint_printf("matrix multiplication: bad exponent!\n");
171		flint_abort();
172		}
173		}
174		}
175
176		/* The main multiplication */
177		fmpz_mat_mul(CC, AA, BB);
178		/* flint_printf("bits %wd %wd %wd\n", fmpz_mat_max_bits(CC),
179		fmpz_mat_max_bits(AA), fmpz_mat_max_bits(BB)); */
180
181		fmpz_mat_clear(AA);
182		fmpz_mat_clear(BB);
183
184		arb_init(t);
185
186		/* Add to the result matrix */
187		for (i = M0; i < M1; i++)
188		{
189		for (j = P0; j < P1; j++)
190		{
191		*e = A_min[i] + B_min[j];
192
193		/* The first time we write this Cij */
194		if (block_start == 0)
195		{
196		arb_set_round_fmpz_2exp(arb_mat_entry(C, i, j),
197		fmpz_mat_entry(CC, i - M0, j - P0), e, prec);
198		}
199		else
200		{
201		arb_set_round_fmpz_2exp(t, fmpz_mat_entry(CC, i - M0, j - P0), e, prec);
202		arb_add(arb_mat_entry(C, i, j), arb_mat_entry(C, i, j), t, prec);
203		}
204		}
205		}
206		arb_clear(t);
207
208		fmpz_mat_clear(CC);
209		}
	197	arb_clear(t);
	198
	199	fmpz_mat_clear(CC);
210	200	}
211	201	}
212	202	}

221	211	slong A_bot, B_bot;
222	212	slong block_start, block_end, i, j, bot, top, max_height;
223	213	slong b, A_max_bits, B_max_bits;
	214	slong min_block_size;
224	215	arb_srcptr t;
225	216	int A_exact, B_exact;
226	217	double A_density, B_density;

348	339	return;
349	340	}
350	341
	342	if (arb_mat_mul_block_min_block_size != 0)
	343	min_block_size = arb_mat_mul_block_min_block_size;
	344	else
	345	min_block_size = 30;
	346
351	347	block_start = 0;
352	348	while (block_start < N)
353	349	{

444	440	}
445	441
446	442	blocks_built:
447		arb_mat_mid_addmul_block_prescaled(C, A, B,
448		block_start, block_end, A_min, B_min, prec);
	443	if (block_end - block_start < min_block_size)
	444	{
	445	block_end = FLINT_MIN(N, block_start + min_block_size);
	446
	447	arb_mat_mid_addmul_block_fallback(C, A, B,
	448	block_start, block_end, prec);
	449	}
	450	else
	451	{
	452	arb_mat_mid_addmul_block_prescaled(C, A, B,
	453	block_start, block_end, A_min, B_min, prec);
	454	}
449	455
450	456	block_start = block_end;
451	457	}

+65

-0

arb_mat/test/t-companion.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_mat.h"
	12
	13	int
	14	main(void)
	15	{
	16	slong iter;
	17	flint_rand_t state;
	18
	19	flint_printf("companion....");
	20	fflush(stdout);
	21
	22	flint_randinit(state);
	23
	24	for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
	25	{
	26	arb_mat_t A;
	27	arb_poly_t f, g;
	28	slong n, prec;
	29
	30	arb_poly_init(f);
	31	arb_poly_init(g);
	32
	33	do {
	34	arb_poly_randtest(f, state, 1 + n_randint(state, 8), 1 + n_randint(state, 1000), 10);
	35	} while (arb_poly_degree(f) < 0);
	36
	37	n = arb_poly_degree(f);
	38	prec = 2 + n_randint(state, 200);
	39
	40	arb_mat_init(A, n, n);
	41	arb_mat_randtest(A, state, 1 + n_randint(state, 1000), 10);
	42	arb_mat_companion(A, f, prec);
	43	arb_mat_charpoly(g, A, prec);
	44	arb_poly_scalar_mul(g, g, arb_poly_get_coeff_ptr(f, n), prec);
	45
	46	if (!arb_poly_contains(g, f))
	47	{
	48	flint_printf("FAIL\n");
	49	flint_printf("A:\n"), arb_mat_printd(A, 15), flint_printf("\n");
	50	flint_printf("f:\n"), arb_poly_printd(f, 15), flint_printf("\n");
	51	flint_printf("g:\n"), arb_poly_printd(g, 15), flint_printf("\n");
	52	flint_abort();
	53	}
	54
	55	arb_mat_clear(A);
	56	arb_poly_clear(f);
	57	arb_poly_clear(g);
	58	}
	59
	60	flint_randclear(state);
	61	flint_cleanup();
	62	flint_printf("PASS\n");
	63	return 0;
	64	}

-0

arb_mat/test/t-mul_block.c less more

9	9	*/
10	10
11	11	#include "arb_mat.h"
	12
	13	ARB_DLL extern slong arb_mat_mul_block_min_block_size;
12	14
13	15	int main()
14	16	{

31	33	rbits1 = 2 + n_randint(state, 200);
32	34	rbits2 = 2 + n_randint(state, 200);
33	35	rbits3 = 2 + n_randint(state, 200);
	36
	37	arb_mat_mul_block_min_block_size = n_randint(state, 10);
34	38
35	39	m = n_randint(state, 10);
36	40	n = n_randint(state, 10);

114	118	n = n_randint(state, 40);
115	119	p = n_randint(state, 40);
116	120
	121	arb_mat_mul_block_min_block_size = n_randint(state, 10);
	122
117	123	if (n_randint(state, 4) == 0)
118	124	{
119	125	exp1 = 4 + n_randint(state, FLINT_BITS);

+23

-3

arb_mat.h less more

125	125	return (mat->r == mat->c);
126	126	}
127	127
	128	int arb_mat_is_exact(const arb_mat_t A);
	129
	130	int arb_mat_is_zero(const arb_mat_t mat);
	131	int arb_mat_is_finite(const arb_mat_t mat);
	132	int arb_mat_is_triu(const arb_mat_t mat);
	133	int arb_mat_is_tril(const arb_mat_t mat);
	134
	135	ARB_MAT_INLINE int
	136	arb_mat_is_diag(const arb_mat_t mat)
	137	{
	138	return arb_mat_is_tril(mat) && arb_mat_is_triu(mat);
	139	}
	140
128	141	/* Radius and interval operations */
129	142
130	143	ARB_MAT_INLINE void

154	167	void arb_mat_one(arb_mat_t mat);
155	168
156	169	void arb_mat_ones(arb_mat_t mat);
	170
	171	void arb_mat_indeterminate(arb_mat_t mat);
157	172
158	173	void arb_mat_hilbert(arb_mat_t mat, slong prec);
159	174

356	371	int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec);
357	372	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);
358	373	int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec);
	374	int arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec);
359	375
360	376	int arb_mat_inv(arb_mat_t X, const arb_mat_t A, slong prec);
361	377

393	409
394	410	void arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec);
395	411
396		void _arb_mat_charpoly(arb_ptr cp, const arb_mat_t mat, slong prec);
397
398		void arb_mat_charpoly(arb_poly_t cp, const arb_mat_t mat, slong prec);
	412	void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec);
	413	void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec);
	414	void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec);
	415	void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec);
399	416
400	417	void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec);
	418
	419	void _arb_mat_diag_prod(arb_t res, const arb_mat_t A, slong a, slong b, slong prec);
	420	void arb_mat_diag_prod(arb_t res, const arb_mat_t A, slong prec);
401	421
402	422	/* Sparsity structure */
403	423

-2

configure less more

9	9	# arb => soname
10	10	# 2.7.0 => 0.0.0
11	11	ARB_MAJOR=2
12		ARB_MINOR=6
13		ARB_PATCH=1
	12	ARB_MINOR=7
	13	ARB_PATCH=0
14	14
15	15	PREFIX="/usr/local"
16	16	GMP_DIR="/usr/local"

-0

doc/source/acb.rst less more

142	142
143	143	Adds err to the error bounds of both the real and imaginary
144	144	parts of x, modifying x in-place.
	145
	146	.. function:: void acb_get_mid(acb_t m, const acb_t x)
	147
	148	Sets m to the midpoint of x.
145	149
146	150	Input and output
147	151	-------------------------------------------------------------------------------

+337

-73

doc/source/acb_mat.rst less more

4	4
5	5	An :type:`acb_mat_t` represents a dense matrix over the complex numbers,
6	6	implemented as an array of entries of type :type:`acb_struct`.
7
8	7	The dimension (number of rows and columns) of a matrix is fixed at
9	8	initialization, and the user must ensure that inputs and outputs to
10	9	an operation have compatible dimensions. The number of rows or columns
11	10	in a matrix can be zero.
12	11
	12	.. note::
	13
	14	Methods prefixed with acb_mat_approx treat all input entries
	15	as floating-point numbers (ignoring the radii of the balls) and
	16	compute floating-point output (balls with zero radius) representing
	17	approximate solutions without error bounds.
	18	All other methods compute rigorous error bounds.
	19	The approx methods are typically useful for computing initial
	20	values or preconditioners for rigorous solvers.
	21	Some users may also find approx methods useful
	22	for doing ordinary numerical linear algebra in applications where
	23	error bounds are not needed.
13	24
14	25	Types, macros and constants
15	26	-------------------------------------------------------------------------------

92	103	Sets mat to a random matrix with up to prec bits of precision
93	104	and with exponents of width up to mag_bits.
94	105
	106	.. function:: void acb_mat_randtest_eig(acb_mat_t mat, flint_rand_t state, acb_srcptr E, slong prec)
	107
	108	Sets mat to a random matrix with the prescribed eigenvalues
	109	supplied as the vector E. The output matrix is required to be
	110	square. We generate a random unitary matrix via a matrix
	111	exponential, and then evaluate an inverse Schur decomposition.
	112
95	113	Input and output
96	114	-------------------------------------------------------------------------------
97	115

107	125	Comparisons
108	126	-------------------------------------------------------------------------------
109	127
	128	Predicate methods return 1 if the property certainly holds and 0 otherwise.
	129
110	130	.. function:: int acb_mat_equal(const acb_mat_t mat1, const acb_mat_t mat2)
111	131
112		Returns nonzero iff the matrices have the same dimensions
113		and identical entries.
	132	Returns whether the matrices have the same dimensions and identical
	133	intervals as entries.
114	134
115	135	.. function:: int acb_mat_overlaps(const acb_mat_t mat1, const acb_mat_t mat2)
116	136
117		Returns nonzero iff the matrices have the same dimensions
	137	Returns whether the matrices have the same dimensions
118	138	and each entry in mat1 overlaps with the corresponding entry in mat2.
119	139
120	140	.. function:: int acb_mat_contains(const acb_mat_t mat1, const acb_mat_t mat2)

123	143
124	144	.. function:: int acb_mat_contains_fmpq_mat(const acb_mat_t mat1, const fmpq_mat_t mat2)
125	145
126		Returns nonzero iff the matrices have the same dimensions and each entry
	146	Returns whether the matrices have the same dimensions and each entry
127	147	in mat2 is contained in the corresponding entry in mat1.
128	148
129	149	.. function:: int acb_mat_eq(const acb_mat_t mat1, const acb_mat_t mat2)
130	150
131		Returns nonzero iff mat1 and mat2 certainly represent the same matrix.
	151	Returns whether mat1 and mat2 certainly represent the same matrix.
132	152
133	153	.. function:: int acb_mat_ne(const acb_mat_t mat1, const acb_mat_t mat2)
134	154
135		Returns nonzero iff mat1 and mat2 certainly do not represent the same matrix.
	155	Returns whether mat1 and mat2 certainly do not represent the same matrix.
136	156
137	157	.. function:: int acb_mat_is_real(const acb_mat_t mat)
138	158
139		Returns nonzero iff all entries in mat have zero imaginary part.
	159	Returns whether all entries in mat have zero imaginary part.
140	160
141	161	.. function:: int acb_mat_is_empty(const acb_mat_t mat)
142	162
143		Returns nonzero iff the number of rows or the number of columns in mat is zero.
	163	Returns whether the number of rows or the number of columns in mat is zero.
144	164
145	165	.. function:: int acb_mat_is_square(const acb_mat_t mat)
146	166
147		Returns nonzero iff the number of rows is equal to the number of columns in mat.
148
	167	Returns whether the number of rows is equal to the number of columns in mat.
	168
	169	.. function:: int acb_mat_is_exact(const acb_mat_t mat)
	170
	171	Returns whether all entries in mat have zero radius.
	172
	173	.. function:: int acb_mat_is_zero(const acb_mat_t mat)
	174
	175	Returns whether all entries in mat are exactly zero.
	176
	177	.. function:: int acb_mat_is_finite(const acb_mat_t mat)
	178
	179	Returns whether all entries in mat are finite.
	180
	181	.. function:: int acb_mat_is_triu(const acb_mat_t mat)
	182
	183	Returns whether mat is upper triangular; that is, all entries
	184	below the main diagonal are exactly zero.
	185
	186	.. function:: int acb_mat_is_tril(const acb_mat_t mat)
	187
	188	Returns whether mat is lower triangular; that is, all entries
	189	above the main diagonal are exactly zero.
	190
	191	.. function:: int acb_mat_is_diag(const acb_mat_t mat)
	192
	193	Returns whether mat is a diagonal matrix; that is, all entries
	194	off the main diagonal are exactly zero.
149	195
150	196	Special matrices
151	197	-------------------------------------------------------------------------------

162	208	.. function:: void acb_mat_ones(acb_mat_t mat)
163	209
164	210	Sets all entries in the matrix to ones.
	211
	212	.. function:: void acb_mat_indeterminate(acb_mat_t mat)
	213
	214	Sets all entries in the matrix to indeterminate (NaN).
165	215
166	216	.. function:: void acb_mat_dft(acb_mat_t mat, int type, slong prec)
167	217

267	317	Sets res to mat raised to the power exp. Requires that mat
268	318	is a square matrix.
269	319
	320	.. function:: void acb_mat_approx_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
	321
	322	Approximate matrix multiplication. The input radii are ignored and
	323	the output matrix is set to an approximate floating-point result.
	324	For performance reasons, the radii in the output matrix will not
	325	necessarily be written (zeroed), but will remain zero if they
	326	are already zeroed in res before calling this function.
270	327
271	328	Scalar arithmetic
272	329	-------------------------------------------------------------------------------

434	491	versions and additionally handles small or triangular matrices
435	492	by direct formulas.
436	493
437		Characteristic polynomial
438		-------------------------------------------------------------------------------
439
440		.. function:: void _acb_mat_charpoly(acb_ptr cp, const acb_mat_t mat, slong prec)
441
442		.. function:: void acb_mat_charpoly(acb_poly_t cp, const acb_mat_t mat, slong prec)
443
444		Sets cp to the characteristic polynomial of mat which must be
445		a square matrix. If the matrix has n rows, the underscore method
446		requires space for `n + 1` output coefficients.
447		Employs a division-free algorithm using `O(n^4)` operations.
448
449		Special functions
450		-------------------------------------------------------------------------------
451
452		.. function:: void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec)
453
454		Sets S to the truncated exponential Taylor series `S = \sum_{k=0}^{N-1} A^k / k!`.
455		See :func:`arb_mat_exp_taylor_sum` for implementation notes.
456
457		.. function:: void acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
458
459		Sets B to the exponential of the matrix A, defined by the Taylor series
460
461		.. math ::
462
463		\exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}.
464
465		The function is evaluated as `\exp(A/2^r)^{2^r}`, where `r` is chosen
466		to give rapid convergence of the Taylor series.
467		Error bounds are computed as for :func:`arb_mat_exp`.
468
469		.. function:: void acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec)
470
471		Sets trace to the trace of the matrix, i.e. the sum of entries on the
472		main diagonal of mat. The matrix is required to be square.
473
474
475		Component and error operations
476		-------------------------------------------------------------------------------
477
478		.. function:: void acb_mat_get_mid(acb_mat_t B, const acb_mat_t A)
479
480		Sets the entries of B to the exact midpoints of the entries of A.
481
482		.. function:: void acb_mat_add_error_mag(acb_mat_t mat, const mag_t err)
483
484		Adds err in-place to the radii of the entries of mat.
485
486		Approximate solving
487		-------------------------------------------------------------------------------
488
489		.. function:: void acb_mat_approx_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
490
491		Approximate matrix multiplication. The input radii are ignored and
492		the output matrix is set to an approximate floating-point result.
493		The radii in the output matrix will not necessarily be zeroed.
494
495	494	.. function:: void acb_mat_approx_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)
496	495
497	496	.. function:: void acb_mat_approx_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)

501	500	.. function:: void acb_mat_approx_solve_lu_precomp(acb_mat_t X, const slong * perm, const acb_mat_t A, const acb_mat_t B, slong prec)
502	501
503	502	.. function:: int acb_mat_approx_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)
	503
	504	.. function:: int acb_mat_approx_inv(acb_mat_t X, const acb_mat_t A, slong prec)
504	505
505	506	These methods perform approximate solving without any error control.
506	507	The radii in the input matrices are ignored, the computations are done

510	511	output matrices are set to the approximate floating-point results with
511	512	zeroed error bounds.
512	513
513		Approximate solutions are useful for computing preconditioning matrices
514		for certified solutions. Some users may also find these methods useful
515		for doing ordinary numerical linear algebra in applications where
516		error bounds are not needed.
	514	Characteristic polynomial and companion matrix
	515	-------------------------------------------------------------------------------
	516
	517	.. function:: void _acb_mat_charpoly(acb_ptr poly, const acb_mat_t mat, slong prec)
	518
	519	.. function:: void acb_mat_charpoly(acb_poly_t poly, const acb_mat_t mat, slong prec)
	520
	521	Sets poly to the characteristic polynomial of mat which must be
	522	a square matrix. If the matrix has n rows, the underscore method
	523	requires space for `n + 1` output coefficients.
	524	Employs a division-free algorithm using `O(n^4)` operations.
	525
	526	.. function:: void _acb_mat_companion(acb_mat_t mat, acb_srcptr poly, slong prec)
	527
	528	.. function:: void acb_mat_companion(acb_mat_t mat, const acb_poly_t poly, slong prec)
	529
	530	Sets the n by n matrix mat to the companion matrix of the polynomial
	531	poly which must have degree n.
	532	The underscore method reads `n + 1` input coefficients.
	533
	534	Special functions
	535	-------------------------------------------------------------------------------
	536
	537	.. function:: void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec)
	538
	539	Sets S to the truncated exponential Taylor series `S = \sum_{k=0}^{N-1} A^k / k!`.
	540	See :func:`arb_mat_exp_taylor_sum` for implementation notes.
	541
	542	.. function:: void acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
	543
	544	Sets B to the exponential of the matrix A, defined by the Taylor series
	545
	546	.. math ::
	547
	548	\exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}.
	549
	550	The function is evaluated as `\exp(A/2^r)^{2^r}`, where `r` is chosen
	551	to give rapid convergence of the Taylor series.
	552	Error bounds are computed as for :func:`arb_mat_exp`.
	553
	554	.. function:: void acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec)
	555
	556	Sets trace to the trace of the matrix, i.e. the sum of entries on the
	557	main diagonal of mat. The matrix is required to be square.
	558
	559	.. function:: void _acb_mat_diag_prod(acb_t res, const acb_mat_t mat, slong a, slong b, slong prec)
	560
	561	.. function:: void acb_mat_diag_prod(acb_t res, const acb_mat_t mat, slong prec)
	562
	563	Sets res to the product of the entries on the main diagonal of mat.
	564	The underscore method computes the product of the entries between
	565	index a inclusive and b exclusive (the indices must be in range).
	566
	567	Component and error operations
	568	-------------------------------------------------------------------------------
	569
	570	.. function:: void acb_mat_get_mid(acb_mat_t B, const acb_mat_t A)
	571
	572	Sets the entries of B to the exact midpoints of the entries of A.
	573
	574	.. function:: void acb_mat_add_error_mag(acb_mat_t mat, const mag_t err)
	575
	576	Adds err in-place to the radii of the entries of mat.
	577
	578	.. _acb-mat-eigenvalues:
	579
	580	Eigenvalues and eigenvectors
	581	-------------------------------------------------------------------------------
	582
	583	The functions in this section are experimental. There are classes
	584	of matrices where the algorithms fail to converge even as
	585	prec is increased, or for which the error bounds are much worse
	586	than necessary. In some cases, it can help to manually precondition
	587	the matrix A by applying a similarity transformation `T^{-1} A T`.
	588
	589	* If A is badly scaled, take `T` to be a matrix such that the entries
	590	of `T^{-1} A T` are more uniform (this is known as balancing).
	591	* Simply taking `T` to be a random invertible matrix can help if an
	592	algorithm fails to converge despite `A` being well-scaled. (This
	593	can be the case when dealing with multiple eigenvalues.)
	594
	595	.. function:: int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)
	596
	597	Computes floating-point approximations of all the n eigenvalues
	598	(and optionally eigenvectors) of the
	599	given n by n matrix A. The approximations of the
	600	eigenvalues are written to the vector E, in no particular order.
	601	If L is not NULL, approximations of the corresponding left
	602	eigenvectors are written to the rows of L. If R is not NULL,
	603	approximations of the corresponding
	604	right eigenvectors are written to the columns of R.
	605
	606	The parameters tol and maxiter can be used to control the
	607	target numerical error and the maximum number of iterations
	608	allowed before giving up. Passing NULL and 0 respectively results
	609	in default values being used.
	610
	611	Uses the implicitly shifted QR algorithm with reduction
	612	to Hessenberg form.
	613	No guarantees are made about the accuracy of the output. A nonzero
	614	return value indicates that the QR iteration converged numerically,
	615	but this is only a heuristic termination test and does not imply
	616	any statement whatsoever about error bounds.
	617	The output may also be accurate even if this function returns zero.
	618
	619	.. function:: void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)
	620
	621	Given an n by n matrix A, a length-n vector E
	622	containing approximations of the eigenvalues of A,
	623	and an n by n matrix R containing approximations of
	624	the corresponding right eigenvectors, computes a rigorous bound
	625	`\varepsilon` such that every eigenvalue `\lambda` of A satisfies
	626	`\|\lambda - \hat \lambda_k\| \le \varepsilon`
	627	for some `\hat \lambda_k` in E.
	628	In other words, the union of the balls
	629	`B_k = \{z : \|z - \hat \lambda_k\| \le \varepsilon\}` is guaranteed to
	630	be an enclosure of all eigenvalues of A.
	631
	632	Note that there is no guarantee that each ball `B_k` can be
	633	identified with a single eigenvalue: it is possible that some
	634	balls contain several eigenvalues while other balls contain
	635	no eigenvalues. In other words, this method is not powerful enough
	636	to compute isolating balls for the individual eigenvalues (or even
	637	for clusters of eigenvalues other than the whole spectrum).
	638	Nevertheless, in practice the balls `B_k` will represent
	639	eigenvalues one-to-one with high probability if the
	640	given approximations are good.
	641
	642	The output can be used to certify
	643	that all eigenvalues of A lie in some region of the complex
	644	plane (such as a specific half-plane, strip, disk, or annulus)
	645	without the need to certify the individual eigenvalues.
	646	The output is easily converted into lower or upper bounds
	647	for the absolute values or real or imaginary parts
	648	of the spectrum, and with high probability these bounds will be tight.
	649	Using :func:`acb_add_error_mag` and :func:`acb_union`, the output
	650	can also be converted to a single :type:`acb_t` enclosing
	651	the whole spectrum of A in a rectangle, but note that to
	652	test whether a condition holds for all eigenvalues of A, it
	653	is typically better to iterate over the individual balls `B_k`.
	654
	655	This function implements the fast algorithm in Theorem 1 in
	656	[Miy2010]_ which extends the Bauer-Fike theorem. Approximations
	657	E and R can, for instance, be computed using
	658	:func:`acb_mat_approx_eig_qr`.
	659	No assumptions are made about the structure of A or the
	660	quality of the given approximations.
	661
	662	.. function:: void acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t R, const acb_mat_t A, const acb_t lambda_approx, const acb_mat_t R_approx, slong prec)
	663
	664	Given an n by n matrix A and an approximate
	665	eigenvalue-eigenvector pair lambda_approx and R_approx (where
	666	R_approx is an n by 1 matrix), computes an enclosure
	667	lambda guaranteed to contain at least one of the eigenvalues of A,
	668	along with an enclosure R for a corresponding right eigenvector.
	669
	670	More generally, this function can handle clustered (or repeated)
	671	eigenvalues. If R_approx is an n by k
	672	matrix containing approximate eigenvectors for a presumed cluster
	673	of k eigenvalues near lambda_approx,
	674	this function computes an enclosure lambda
	675	guaranteed to contain at
	676	least k eigenvalues of A along with a matrix R guaranteed to
	677	contain a basis for the k-dimensional invariant subspace
	678	associated with these eigenvalues.
	679	Note that for multiple eigenvalues, determining the individual eigenvectors is
	680	an ill-posed problem; describing an enclosure of the invariant subspace
	681	is the best we can hope for.
	682
	683	For `k = 1`, it is guaranteed that `AR - R \lambda` contains
	684	the zero matrix. For `k > 2`, this cannot generally be guaranteed
	685	(in particular, A might not diagonalizable).
	686	In this case, we can still compute an approximately diagonal
	687	k by k interval matrix `J \approx \lambda I` such that `AR - RJ`
	688	is guaranteed to contain the zero matrix.
	689	This matrix has the property that the Jordan canonical form of
	690	(any exact matrix contained in) A has a k by k submatrix
	691	equal to the Jordan canonical form of
	692	(some exact matrix contained in) J.
	693	The output J is optional (the user can pass NULL to omit it).
	694
	695	The algorithm follows section 13.4 in [Rum2010]_, corresponding
	696	to the ``verifyeig()`` routine in INTLAB.
	697	The initial approximations can, for instance, be computed using
	698	:func:`acb_mat_approx_eig_qr`.
	699	No assumptions are made about the structure of A or the
	700	quality of the given approximations.
	701
	702	.. function:: int acb_mat_eig_simple_rump(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	703
	704	.. function:: int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	705
	706	.. function:: int acb_mat_eig_simple(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	707
	708	Computes all the eigenvalues (and optionally corresponding
	709	eigenvectors) of the given n by n matrix A.
	710
	711	Attempts to prove that A has n simple (isolated)
	712	eigenvalues, returning 1 if successful
	713	and 0 otherwise. On success, isolating complex intervals for the
	714	eigenvalues are written to the vector E, in no particular order.
	715	If L is not NULL, enclosures of the corresponding left
	716	eigenvectors are written to the rows of L. If R is not NULL,
	717	enclosures of the corresponding
	718	right eigenvectors are written to the columns of R.
	719
	720	The left eigenvectors are normalized so that `L = R^{-1}`.
	721	This produces a diagonalization `LAR = D` where D is the
	722	diagonal matrix with the entries in E on the diagonal.
	723
	724	The user supplies approximations E_approx and R_approx
	725	of the eigenvalues and the right eigenvectors.
	726	The initial approximations can, for instance, be computed using
	727	:func:`acb_mat_approx_eig_qr`.
	728	No assumptions are made about the structure of A or the
	729	quality of the given approximations.
	730
	731	Two algorithms are implemented:
	732
	733	* The rump version calls :func:`acb_mat_eig_enclosure_rump` repeatedly
	734	to certify eigenvalue-eigenvector pairs one by one. The iteration is
	735	stopped to return non-success if a new eigenvalue overlaps with
	736	previously computed one. Finally, L is computed by a matrix inversion.
	737	This has complexity `O(n^4)`.
	738
	739	* The vdhoeven_mourrain version uses the algorithm in [HM2017]_ to
	740	certify all eigenvalues and eigenvectors in one step. This has
	741	complexity `O(n^3)`.
	742
	743	The default version currently uses vdhoeven_mourrain.
	744
	745	By design, these functions terminate instead of attempting to
	746	compute eigenvalue clusters if some eigenvalues cannot be isolated.
	747	To compute all eigenvalues of a matrix allowing for overlap,
	748	:func:`acb_mat_eig_multiple_rump` may be used as a fallback,
	749	or :func:`acb_mat_eig_multiple` may be used in the first place.
	750
	751	.. function:: int acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	752
	753	.. function:: int acb_mat_eig_multiple(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
	754
	755	Computes all the eigenvalues of the given n by n matrix A.
	756	On success, the output vector E contains n complex intervals,
	757	each representing one eigenvalue of A with the correct
	758	multiplicities in case of overlap.
	759	The output intervals are either disjoint or identical, and
	760	identical intervals are guaranteed to be grouped consecutively.
	761	Each complete run of k identical intervals thus represents a cluster of
	762	exactly k eigenvalues which could not be separated from each
	763	other at the current precision, but which could be isolated
	764	from the other `n - k` eigenvalues of the matrix.
	765
	766	The user supplies approximations E_approx and R_approx
	767	of the eigenvalues and the right eigenvectors.
	768	The initial approximations can, for instance, be computed using
	769	:func:`acb_mat_approx_eig_qr`.
	770	No assumptions are made about the structure of A or the
	771	quality of the given approximations.
	772
	773	The rump algorithm groups approximate eigenvalues that are close
	774	and calls :func:`acb_mat_eig_enclosure_rump` repeatedly to validate
	775	each cluster. The complexity is `O(m n^3)` for m clusters.
	776
	777	The default version, as currently implemented, first attempts to
	778	call :func:`acb_mat_eig_simple_vdhoeven_mourrain` hoping that the
	779	eigenvalues are actually simple. It then uses the rump algorithm as
	780	a fallback.

+114

-48

doc/source/arb_mat.rst less more

4	4
5	5	An :type:`arb_mat_t` represents a dense matrix over the real numbers,
6	6	implemented as an array of entries of type :type:`arb_struct`.
7
8	7	The dimension (number of rows and columns) of a matrix is fixed at
9	8	initialization, and the user must ensure that inputs and outputs to
10	9	an operation have compatible dimensions. The number of rows or columns
11	10	in a matrix can be zero.
12	11
	12	.. note::
	13
	14	Methods prefixed with arb_mat_approx treat all input entries
	15	as floating-point numbers (ignoring the radii of the balls) and
	16	compute floating-point output (balls with zero radius) representing
	17	approximate solutions without error bounds.
	18	All other methods compute rigorous error bounds.
	19	The approx methods are typically useful for computing initial
	20	values or preconditioners for rigorous solvers.
	21	Some users may also find approx methods useful
	22	for doing ordinary numerical linear algebra in applications where
	23	error bounds are not needed.
13	24
14	25	Types, macros and constants
15	26	-------------------------------------------------------------------------------

103	114	Comparisons
104	115	-------------------------------------------------------------------------------
105	116
	117	Predicate methods return 1 if the property certainly holds and 0 otherwise.
	118
106	119	.. function:: int arb_mat_equal(const arb_mat_t mat1, const arb_mat_t mat2)
107	120
108		Returns nonzero iff the matrices have the same dimensions
109		and identical entries.
	121	Returns whether the matrices have the same dimensions
	122	and identical intervals as entries.
110	123
111	124	.. function:: int arb_mat_overlaps(const arb_mat_t mat1, const arb_mat_t mat2)
112	125
113		Returns nonzero iff the matrices have the same dimensions
	126	Returns whether the matrices have the same dimensions
114	127	and each entry in mat1 overlaps with the corresponding entry in mat2.
115	128
116	129	.. function:: int arb_mat_contains(const arb_mat_t mat1, const arb_mat_t mat2)

119	132
120	133	.. function:: int arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2)
121	134
122		Returns nonzero iff the matrices have the same dimensions and each entry
	135	Returns whether the matrices have the same dimensions and each entry
123	136	in mat2 is contained in the corresponding entry in mat1.
124	137
125	138	.. function:: int arb_mat_eq(const arb_mat_t mat1, const arb_mat_t mat2)
126	139
127		Returns nonzero iff mat1 and mat2 certainly represent the same matrix.
	140	Returns whether mat1 and mat2 certainly represent the same matrix.
128	141
129	142	.. function:: int arb_mat_ne(const arb_mat_t mat1, const arb_mat_t mat2)
130	143
131		Returns nonzero iff mat1 and mat2 certainly do not represent the same matrix.
	144	Returns whether mat1 and mat2 certainly do not represent the same matrix.
132	145
133	146	.. function:: int arb_mat_is_empty(const arb_mat_t mat)
134	147
135		Returns nonzero iff the number of rows or the number of columns in mat is zero.
	148	Returns whether the number of rows or the number of columns in mat is zero.
136	149
137	150	.. function:: int arb_mat_is_square(const arb_mat_t mat)
138	151
139		Returns nonzero iff the number of rows is equal to the number of columns in mat.
	152	Returns whether the number of rows is equal to the number of columns in mat.
	153
	154	.. function:: int arb_mat_is_exact(const arb_mat_t mat)
	155
	156	Returns whether all entries in mat have zero radius.
	157
	158	.. function:: int arb_mat_is_zero(const arb_mat_t mat)
	159
	160	Returns whether all entries in mat are exactly zero.
	161
	162	.. function:: int arb_mat_is_finite(const arb_mat_t mat)
	163
	164	Returns whether all entries in mat are finite.
	165
	166	.. function:: int arb_mat_is_triu(const arb_mat_t mat)
	167
	168	Returns whether mat is upper triangular; that is, all entries
	169	below the main diagonal are exactly zero.
	170
	171	.. function:: int arb_mat_is_tril(const arb_mat_t mat)
	172
	173	Returns whether mat is lower triangular; that is, all entries
	174	above the main diagonal are exactly zero.
	175
	176	.. function:: int arb_mat_is_diag(const arb_mat_t mat)
	177
	178	Returns whether mat is a diagonal matrix; that is, all entries
	179	off the main diagonal are exactly zero.
140	180
141	181	Special matrices
142	182	-------------------------------------------------------------------------------

153	193	.. function:: void arb_mat_ones(arb_mat_t mat)
154	194
155	195	Sets all entries in the matrix to ones.
	196
	197	.. function:: void arb_mat_indeterminate(arb_mat_t mat)
	198
	199	Sets all entries in the matrix to indeterminate (NaN).
156	200
157	201	.. function:: void arb_mat_hilbert(arb_mat_t mat)
158	202

292	336	order and B is a linear array of coefficients in column-major order.
293	337	This function assumes that all exponents are small and is unsafe
294	338	for general use.
	339
	340	.. function:: void arb_mat_approx_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)
	341
	342	Approximate matrix multiplication. The input radii are ignored and
	343	the output matrix is set to an approximate floating-point result.
	344	The radii in the output matrix will not necessarily be zeroed.
295	345
296	346	Scalar arithmetic
297	347	-------------------------------------------------------------------------------

467	517	versions and additionally handles small or triangular matrices
468	518	by direct formulas.
469	519
	520	.. 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)
	521
	522	.. 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)
	523
	524	.. function:: int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)
	525
	526	.. 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)
	527
	528	.. function:: int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)
	529
	530	.. function:: int arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec)
	531
	532	These methods perform approximate solving without any error control.
	533	The radii in the input matrices are ignored, the computations are done
	534	numerically with floating-point arithmetic (using ordinary
	535	Gaussian elimination and triangular solving, accelerated through
	536	the use of block recursive strategies for large matrices), and the
	537	output matrices are set to the approximate floating-point results with
	538	zeroed error bounds.
	539
	540	Approximate solutions are useful for computing preconditioning matrices
	541	for certified solutions. Some users may also find these methods useful
	542	for doing ordinary numerical linear algebra in applications where
	543	error bounds are not needed.
	544
470	545	Cholesky decomposition and solving
471	546	-------------------------------------------------------------------------------
472	547

567	642	The inverse is calculated using the method of [Kri2013]_ which is more
568	643	efficient than solving `AX = I` with :func:`arb_mat_solve_ldl_precomp`.
569	644
570		Characteristic polynomial
571		-------------------------------------------------------------------------------
572
573		.. function:: void _arb_mat_charpoly(arb_ptr cp, const arb_mat_t mat, slong prec)
574
575		.. function:: void arb_mat_charpoly(arb_poly_t cp, const arb_mat_t mat, slong prec)
576
577		Sets cp to the characteristic polynomial of mat which must be
	645	Characteristic polynomial and companion matrix
	646	-------------------------------------------------------------------------------
	647
	648	.. function:: void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec)
	649
	650	.. function:: void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec)
	651
	652	Sets poly to the characteristic polynomial of mat which must be
578	653	a square matrix. If the matrix has n rows, the underscore method
579	654	requires space for `n + 1` output coefficients.
580	655	Employs a division-free algorithm using `O(n^4)` operations.
	656
	657	.. function:: void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec)
	658
	659	.. function:: void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec)
	660
	661	Sets the n by n matrix mat to the companion matrix of the polynomial
	662	poly which must have degree n.
	663	The underscore method reads `n + 1` input coefficients.
581	664
582	665	Special functions
583	666	-------------------------------------------------------------------------------

622	705	Sets trace to the trace of the matrix, i.e. the sum of entries on the
623	706	main diagonal of mat. The matrix is required to be square.
624	707
	708	.. function:: void _arb_mat_diag_prod(arb_t res, const arb_mat_t mat, slong a, slong b, slong prec)
	709
	710	.. function:: void arb_mat_diag_prod(arb_t res, const arb_mat_t mat, slong prec)
	711
	712	Sets res to the product of the entries on the main diagonal of mat.
	713	The underscore method computes the product of the entries between
	714	index a inclusive and b exclusive (the indices must be in range).
	715
625	716	Sparsity structure
626	717	-------------------------------------------------------------------------------
627	718

659	750
660	751	Adds err in-place to the radii of the entries of mat.
661	752
662		Approximate solving
663		-------------------------------------------------------------------------------
664
665		.. function:: void arb_mat_approx_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)
666
667		Approximate matrix multiplication. The input radii are ignored and
668		the output matrix is set to an approximate floating-point result.
669		The radii in the output matrix will not necessarily be zeroed.
670
671		.. 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)
672
673		.. 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)
674
675		.. function:: int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)
676
677		.. 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)
678
679		.. function:: int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)
680
681		These methods perform approximate solving without any error control.
682		The radii in the input matrices are ignored, the computations are done
683		numerically with floating-point arithmetic (using ordinary
684		Gaussian elimination and triangular solving, accelerated through
685		the use of block recursive strategies for large matrices), and the
686		output matrices are set to the approximate floating-point results with
687		zeroed error bounds.
688
689		Approximate solutions are useful for computing preconditioning matrices
690		for certified solutions. Some users may also find these methods useful
691		for doing ordinary numerical linear algebra in applications where
692		error bounds are not needed.
	753	Eigenvalues and eigenvectors
	754	-------------------------------------------------------------------------------
	755
	756	To compute eigenvalues and eigenvectors, one can convert to an
	757	:type:`acb_mat_t` and use the functions in :ref:`acb_mat.h: Eigenvalues and eigenvectors<acb-mat-eigenvalues>`.
	758	In the future dedicated methods for real matrices will be added here.

-1

doc/source/credits.rst less more

194	194
195	195	.. [Hoe2001] \J. van der Hoeven. "Fast evaluation of holonomic functions near and in regular singularities", Journal of Symbolic Computation, 31(6):717-743 (2001).
196	196
	197	.. [HM2017] \J. van der Hoeven and B. Mourrain. "Efficient certication of numeric solutions to eigenproblems", MACIS 2017, 81-94, (2017), https://hal.archives-ouvertes.fr/hal-01579079
	198
197	199	.. [JB2018] \F. Johansson and I. Blagouchine. "Computing Stieltjes constants using complex integration", preprint (2018), https://arxiv.org/abs/1804.01679
198	200
199	201	.. [Joh2012] \F. Johansson, "Efficient implementation of the Hardy-Ramanujan-Rademacher formula", LMS Journal of Computation and Mathematics, Volume 15 (2012), 341-359, http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=8710297

224	226
225	227	.. [Kri2013] \A. Krishnamoorthy and D. Menon, "Matrix Inversion Using Cholesky Decomposition" Proc. of the International Conference on Signal Processing Algorithms, Architectures, Arrangements, and Applications (SPA-2013), pp. 70-72, 2013.
226	228
	229	.. [Miy2010] \S. Miyajima, "Fast enclosure for all eigenvalues in generalized eigenvalue problems", Journal of Computational and Applied Mathematics, 233 (2010), 2994-3004, https://dx.doi.org/10.1016/j.cam.2009.11.048
	230
227	231	.. [MPFR2012] The MPFR team, "MPFR Algorithms" (2012), http://www.mpfr.org/algo.html
228	232
229	233	.. [NIST2012] National Institute of Standards and Technology, Digital Library of Mathematical Functions (2012), http://dlmf.nist.gov/

245	249	.. [Tak2000] \D. Takahashi, "A fast algorithm for computing large Fibonacci numbers", Information Processing Letters 75 (2000) 243-246, http://www.ii.uni.wroc.pl/~lorys/IPL/article75-6-1.pdf
246	250
247	251	.. [Tre2008] \L. N. Trefethen, "Is Gauss Quadrature Better than Clenshaw-Curtis?", SIAM Review, 50:1 (2008), 67-87, https://doi.org/10.1137/060659831
248

+25

-11

doc/source/examples.rst less more

49	49
50	50	Sample output::
51	51
52		> build/examples/hilbert_matrix 200
53		prec=20: 0 +/- 5.5777e-330
54		prec=40: 0 +/- 2.5785e-542
55		prec=80: 0 +/- 8.1169e-926
56		prec=160: 0 +/- 2.8538e-1924
57		prec=320: 0 +/- 6.3868e-4129
58		prec=640: 0 +/- 1.7529e-8826
59		prec=1280: 0 +/- 1.8545e-17758
60		prec=2560: 2.955454297e-23924 +/- 6.4586e-24044
	52	$ build/examples/hilbert_matrix 200
	53	prec=20: [+/- 1.32e-335]
	54	prec=40: [+/- 1.63e-545]
	55	prec=80: [+/- 1.30e-933]
	56	prec=160: [+/- 3.62e-1926]
	57	prec=320: [+/- 1.81e-4129]
	58	prec=640: [+/- 3.84e-8838]
	59	prec=1280: [2.955454297e-23924 +/- 8.29e-23935]
61	60	success!
62		cpu/wall(s): 9.06 9.095
63		virt/peak/res/peak(MB): 55.52 55.52 35.50 35.50
	61	cpu/wall(s): 8.494 8.513
	62	virt/peak/res/peak(MB): 134.98 134.98 111.57 111.57
	63
	64	Called with ``-eig n``, instead of computing the determinant,
	65	the program computes the smallest eigenvalue of the Hilbert matrix
	66	(in fact, it isolates all eigenvalues and prints the smallest eigenvalue)::
	67
	68	$ build/examples/hilbert_matrix -eig 50
	69	prec=20: nan
	70	prec=40: nan
	71	prec=80: nan
	72	prec=160: nan
	73	prec=320: nan
	74	prec=640: [1.459157797e-74 +/- 2.49e-84]
	75	success!
	76	cpu/wall(s): 1.84 1.841
	77	virt/peak/res/peak(MB): 33.97 33.97 10.51 10.51
64	78
65	79	keiper_li.c
66	80	-------------------------------------------------------------------------------

+28

-0

doc/source/history.rst less more

11	11	Old versions of the documentation
12	12	-------------------------------------------------------------------------------
13	13
	14	* http://arblib.org/arb-2.16.0.pdf
14	15	* http://arblib.org/arb-2.15.0.pdf
15	16	* http://arblib.org/arb-2.14.0.pdf
16	17	* http://arblib.org/arb-2.13.0.pdf

26	27	* http://arblib.org/arb-2.5.0.pdf
27	28	* http://arblib.org/arb-2.4.0.pdf
28	29	* http://arblib.org/arb-2.3.0.pdf
	30
	31	2018-12-07 -- version 2.16.0
	32	-------------------------------------------------------------------------------
	33
	34	* Linear algebra and arithmetic
	35
	36	* Added acb_mat_approx_eig_qr for approximate computation of eigenvalues
	37	and eigenvectors of complex matrices.
	38	* Added acb_mat_eig_enclosure_rump implementing Rump's algorithm for
	39	certification of eigenvalue-eigenvector pairs as well as clusters.
	40	* Added acb_mat_eig_simple_rump for certified diagonalization of matrices
	41	with simple eigenvalues.
	42	* Added acb_mat_eig_simple_vdhoeven_mourrain, acb_mat_eig_simple for fast
	43	certified diagonalization of matrices with simple eigenvalues.
	44	* Added acb_mat_eig_multiple_rump, acb_mat_eig_multiple for certified
	45	computation of eigenvalues with possible overlap.
	46	* Added acb_mat_eig_global_enclosure for fast global inclusion of eigenvalues
	47	without isolation.
	48	* Added arb_mat_companion, acb_mat_companion for constructing companion
	49	matrices.
	50	* Added several arb_mat and acb_mat helper functions: indeterminate,
	51	is_exact, is_zero, is_finite, is_triu, is_tril, is_diag, diag_prod.
	52	* Added arb_mat_approx_inv, acb_mat_approx_inv.
	53	* Optimized arb_mat_mul_block by using arb_dot when the blocks are small.
	54	* Added acb_get_mid.
	55	* Updated hilbert_matrix example program.
	56
29	57
30	58	2018-10-25 -- version 2.15.1
31	59	-------------------------------------------------------------------------------

+16

-18

doc/source/setup.rst less more

167	167	Computer algebra systems and wrappers
168	168	-------------------------------------------------------------------------------
169	169
170		SageMath (http://sagemath.org/) includes Arb as a standard package and
171		contains a high-level Python interface. Refer to the SageMath documentation:
	170	* Python-FLINT (https://github.com/fredrik-johansson/python-flint) is a
	171	convenient Python interface to both FLINT and Arb.
172	172
173		* RealBallField: http://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/real_arb.html
174		* ComplexBallField: http://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/complex_arb.html
	173	* SageMath (http://sagemath.org/) includes Arb as a standard package and
	174	contains a high-level Python interface. Refer to the SageMath documentation:
175	175
176		Nemo (http://nemocas.org/) is a computer algebra package for the
177		Julia programming language which includes a high-level Julia interface to Arb.
178		The Nemo installation script will create a local installation of
179		Arb along with other dependencies.
	176	* RealBallField: http://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/real_arb.html
	177	* ComplexBallField: http://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/complex_arb.html
180	178
181		* Real balls: http://nemocas.github.io/Nemo.jl/latest/arb.html
182		* Complex balls: http://nemocas.github.io/Nemo.jl/latest/acb.html
	179	* Nemo (http://nemocas.org/) is a computer algebra package for the
	180	Julia programming language which includes a high-level Julia interface to Arb.
	181	The Nemo installation script will create a local installation of
	182	Arb along with other dependencies.
183	183
184		Other wrappers are also available:
	184	* Real balls: http://nemocas.github.io/Nemo.jl/latest/arb.html
	185	* Complex balls: http://nemocas.github.io/Nemo.jl/latest/acb.html
185	186
186		* An experimental standalone Python interface to FLINT and Arb (not requiring SageMath): https://github.com/fredrik-johansson/python-flint
187		* A Java wrapper using JNA: https://github.com/crowlogic/arb/
188		* Another Julia interface: https://github.com/JuliaArbTypes/ArbFloats.jl
	187	* Other wrappers include:
189	188
190		Since the various wrappers add some overhead and do not expose all the
191		functionality in Arb, you may consider using them to try out Arb and prototype
192		algorithms while doing a final implementation in C.
193
	189	* ArbNumerics (Julia): https://github.com/JeffreySarnoff/ArbNumerics.jl
	190	* ArbFloats (Julia): https://github.com/JuliaArbTypes/ArbFloats.jl
	191	* A Java wrapper using JNA: https://github.com/crowlogic/arb/

+53

-8

examples/hilbert_matrix.c less more

0	0	/* This file is public domain. Author: Fredrik Johansson. */
1	1
	2	#include <string.h>
2	3	#include "arb_mat.h"
	4	#include "acb_mat.h"
3	5	#include "flint/profiler.h"
4	6
5	7	int main(int argc, char *argv[])
6	8	{
7	9	arb_mat_t A;
8	10	arb_t det;
9		slong prec, n;
	11	slong i, prec, n;
	12	int eig;
10	13
11		if (argc < 2)
	14	if (argc < 2 \|\| (argc == 3 && strcmp(argv[1], "-eig")))
12	15	{
13		flint_printf("usage: build/examples/hilbert_matrix n\n");
	16	flint_printf("usage: build/examples/hilbert_matrix [-eig] n\n");
14	17	return 1;
15	18	}
16	19
17		n = atol(argv[1]);
	20	if (argc == 2)
	21	eig = 0;
	22	else
	23	eig = 1;
	24
	25	n = atol(argv[argc-1]);
	26
	27	if (eig && (n == 0))
	28	{
	29	flint_printf("smallest eigenvalue: undefined for n == 0\n");
	30	return 1;
	31	}
18	32
19	33	arb_mat_init(A, n, n);
20	34	arb_init(det);

24	38	for (prec = 20; ; prec *= 2)
25	39	{
26	40	arb_mat_hilbert(A, prec);
27
28	41	flint_printf("prec=%wd: ", prec);
29	42
30		arb_mat_det(det, A, prec);
	43	if (eig == 0)
	44	{
	45	arb_mat_det(det, A, prec);
	46	}
	47	else
	48	{
	49	acb_mat_t C, R;
	50	acb_ptr E;
31	51
32		arb_printd(det, 10);
	52	acb_mat_init(R, n, n);
	53	acb_mat_init(C, n, n);
	54	E = _acb_vec_init(n);
	55
	56	acb_mat_set_arb_mat(C, A);
	57	acb_mat_approx_eig_qr(E, NULL, R, C, NULL, 0, prec);
	58	if (acb_mat_eig_simple(E, NULL, NULL, C, E, R, prec))
	59	{
	60	/* A is symmetric so the eigenvalues are real */
	61	for (i = 0; i < n; i++)
	62	arb_zero(acb_imagref(E + i));
	63	/* With isolated eigenvalues, sorting midpoints gives the
	64	right result. */
	65	_acb_vec_sort_pretty(E, n);
	66	acb_get_real(det, E + 0);
	67	}
	68	else
	69	{
	70	arb_indeterminate(det);
	71	}
	72
	73	acb_mat_clear(R);
	74	acb_mat_clear(C);
	75	_acb_vec_clear(E, n);
	76	}
	77
	78	arb_printn(det, 10, 0);
33	79	flint_printf("\n");
34	80
35	81	if (!arb_contains_zero(det))

48	94	flint_cleanup();
49	95	return 0;
50	96	}
51