{'fraction_non_alphanumeric': 0.054982817869415807, 'fraction_numerical': 0.020618556701030927, 'mean_word_length': 3.887885786269158, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Use of Solvable Directed Graphs in a Jacobi-like Algorithm

In this paper, we introduce a Jacobi-like algorithm (we call D-NJLA) to reduce a real nonsymmetric n × n matrix to a real upper triangular form by the help of solvable directed graphs. This method uses only real arithmetic and a sequence of orthogonal similarity transformations and achieves ultimate quadratic convergence. A theoretical analysis is constructed and some experimental results are given.


I. INTRODUCTION
The numerical computation of matrices' eigenvalues is a problem of major importance in many scientific and engineering applications. The applications areas, where eigenvalue calculations arise, are structural dynamics, electrical networks, macro-economics, quantum chemistry, chemical reactions, control theory, etc [26]. It is generally difficult to determine precisely eigenvalues of a matrix by solving the characteristic equation. So, the fact that similarity transformations on a matrix leave its eigenvalues invariant, i.e., a n × n matrix can be reduced to the simpler form through unitary similarity transformations where Q is unitary and B is the simpler matrix form (diagonal, upper or lower triangular matrix). So two matrix A and B are called similar. However, as the size of the matrix A increases, finding the unitary similarity matrix Q becomes more difficult. Therefore, the unitary transformation matrix Q is calculated iteratively, i.e., where Q r+1 is the unitary matrix and r is the number of steps [3], [9]. The most well-known of the iterative methods developed for symmetric matrices is the Jacobi method [5], [7], [12], [25]. Because it is successful and easy to implement for the computation of eigenvalues of symmetric matrices, now it has been generalized for nonsymmetric matrices which is known Jacobi-like algorithm [6], [10], [13], [18], [20], [21]. The convergence rate of these algorithm depends on the sorting of the diagonal elements of A and choosing elements to be annihilated. Therefore, in some algorithms, the diagonal elements, are sorted by the choosing of orthogonal transformation matrices and elements to be annihilated [4], [10], [20].
In this paper, we propose a Jacobi-like algorithm, which computes the Schur form of nonsymmetric square matrices A, with the help of solvable directed graphs, over the real field. The proposed algorithm was designed in two-stages: Stage1 and Stage2. To compute eigenvalues of A with solvable directed graphs, the strategies for choosing the element to be eliminated and the proper orthogonal transformation matrices were determined separately for these stages. Stage1 is composed of reducing the matrix to nearly 2 × 2 upper block triangular form by a solvable directed graph. In each step in Stage1, we choose an entry to be eliminated, which has the property that the strictly lower triangular part element has the largest magnitude and the sign of 2 × 2 subdeterminant is negative, i.e., strictly lower triangular elements whose one diagonal element corresponding to a positive eigenvalue, and another diagonal element corresponding to a negative eigenvalue of the original matrix A are annihilated. Every step of the annihilation, we choose the orthogonal transformation matrix as close as possible to identity and which corresponds to the orthogonal transformation matrix whose rotation parameter falls into [−3, −2] as much as possible.
Stage2 consists of reducing the matrix that is obtained by the result of Stage1 to 2 × 2 partitioned form whose diagonal blocks are upper triangular matrices whose diagonal elements correspond to greater than or equal to zero and negative eigenvalues of the original matrix A with the help of a solvable directed graph. In each step in Stage2, we choose an entry to be eliminated, which has the property that the strictly lower triangular part element has the largest magnitude and the sign of 2 × 2 subdeterminant is positive or 2 × 2 subdeterminant equals to zero, i.e., strictly lower triangular elements whose both diagonal elements corresponding same signed eigenvalues of the original matrix A are annihilated. At every step of the annihilation, we choose the transformation matrix as close as possible to identity and which corresponds to the orthogonal transformation matrix whose rotation parameter falls into [0, 1] as much as possible. The matrix obtained by the result of Stage2 is reduced in Stage1 in the next cycle. Stage1 and Stage2 have been alternatively applied until reducing the matrix A to the upper triangular form. So, processing in this manner eigenvalues of the matrix A are computed by the help of solvable directed graphs. The remainder of the paper is organized as follows: the connection between the annihilation process and directed graphs is given in Section 2. In Section 3 the algorithm for calculating eigenvalues is described. We give some numerical experiments in Section 4. The conclusions are given in Section 5.
II. GRAPH THEORETICAL APPROACH TO THE ALGORITHM The algorithm was built based on some definitions in graph theory. So it would be useful to give some necessary basic graph theory definitions. We consider only finite directed graphs.
Definition 2.1: A directed graph (digraph) is a pair D = (V, E), where V is a set of vertices and E ⊆ V × V is a binary relation representing the directed edges of D [19]. For a vertex x ∈ D, we denote E(x) = {y ∈ D|E(x, y)} the set successors of x.
Definition 2.2: Let D = (V, E) be a directed graph. A subset of vertices K is a kernel of D when: • K is stable: it contains no pair of adjacent vertices • K is absorbing: every vertex v / ∈ K has a successor in K, i.e., ∀ / ∈ K, there is a vertex k ∈ K such that (v, k) ∈ D [22]. A digraph possessing a kernel is called solvable [17]. Now, we associate (2) with Definition 2.1. Each Q r+1 in (2) is calculated by annihilating an element of strictly lower triangular part of A r . Let the indices of the elements annihilated after performing r steps be (j 1 , i 1 ), (j 2 , i 2 ), . . . , (j r , i r ).
Here, each step of (2) makes an element of strictly lower triangular part of A r zero, and computation to do this can destroy pairs of zeros already created. So, let consider any distinct points in the plane, which we shall call vertex. For every sequential annihilation, we connect the vertex P (jm,im) to the vertex P (jm+1,im+1) by means of a directed edge P (jm,im) P (jm+1,im+1) direct from P (jm,im) to P (jm+1,im+1) . In this way, the annihilation process can be associated with a finite directed graph G(A r+1 ) = D(E, V ).
To reduce a (n×n) nonsymmetric matrix by means of solvable finite directed graphs the elements to be annihilated and orthogonal transformation matrices are chosen strategically at each cycle of Stage1 and Stage2 of the D-NJLA.
The key problem is how to choose the elements to be eliminated and orthogonal matrices during the computation. These are the main topic that we discuss in the next section.

III. COMPUTATION OF EIGENVALUES OF A SQUARE
NONSYMMETRIC MATRIX BY D-NJLA ALGORITHM In this section, a Jacobi-like algorithm (which we call D-NJLA) is proposed. Before depict D-NJLA, we first give a necessary definition of matrices.
Definition 3.1: Let a non-symmetric n × n matrix A be divided into two blocks as follows: where, the size of block A 1,1 ,A 2,2 ,A 1,2 and A 2,1 is p × p, n − p × n − p, p × n − p and n − p × p, respectively. If A 2,1 2 = O( ), then A will be nearly 2 × 2 block upper triangular. If A 1,1 = U 1,1 and A 2,2 = U 2,2 are upper triangular matrices, then A will be 2 × 2 partitioned form [16]. In D-NJLA, at each cycle consist of two stages: Stage1 and Stage2. Let A = (a i,j ) = A 1 0 be a nonsymmetric real n × n matrix. In Stage1, a sequence of matrices A k r − = (a r − ,k i,j ), is obtained by a convergent sequence of successive transformations of the form : where k ≥ 1 is the cycle number, r − ≥ 0 is the step number in k-th cycle and Q k r − +1 is an orthogonal matrix in k-th cycle such that where c = y + 1 2 and (10) and (11) correspond to [0, π 3 ] and [ 2π 3 , π], respectively.
Case I: for ∆ k r − > 0, where ∆ k r − is the discriminant of (51), yields and Case II: for ∆ k r − = 0, yields , the roots are out of the scope in our modification. It is detailed in Appendix A and Appendix B how the orthogonal transformation matrix Q k r − +1 and y roots are obtained. At the k-th cycle of Stage1, after performing r − steps we which is obtained at the (k − 1)-th cycle after performing r + steps and is in form like (6).
On the other hand, in the Stage2 a sequence of matrices A k r + = (a r + ,k i,j ), is obtained by a convergent sequence of successive transformations of the form : where k ≥ 1 is the cycle number, r + ≥ 0 is the step number in k-th cycle and Q k r + +1 is an orthogonal matrix given by (9) with (10) and (11). y is computed by (12) and (13), where K k r + , L k r + and δ r + are given by , which is in form like (7).
IV. CHOOSING THE ELEMENT TO BE ANNIHILATED AND PROPER Q k r+1 MATRICES DURING THE COMPUTATION The characteristic equation of a n × n matrix A = (a i,j ) is, in general, an equation of the nth degree polynomial as given If µ 1 , µ 2 , . . . , µ n are eigenvalues of A, we have where T r(A) is the trace of A, and where D is the determinant of A.
Then the coefficients c 1 , c 2 , . . . , c n of the characteristic equation are . . .
Determinant of A can be calculated by where c n is given above [15].
From (19), the sign of D, which is the product of eigenvalues of A, is related to the coefficients of the characteristic equation. Therefore, we use the sign of the subdeterminant to choose the elements to be annihilated in the Stage1 and Stage2. So, before explaining the strategy of choosing the element to be annihilated, we give the following definition. (2,4), . . . , (2, n), . . . , (n− 1, n)) be a finite sequence. Then a 2 × 2 submatrix of A is a matrix obtained from A by retaining only the elements belonging to rows with suffixes i, j columns with same suffixes i, j where (i, j) ∈ S. A 2 × 2 submatrix of A will be denoted by Also, determinant of the submatrix A[i, j] is 2 × 2 subdeterminant that will be denoted by where ((i, j) ∈ S) [16].
The annihilation strategy and choosing proper orthogonal transformation matrices of Stage1 and Stage2 is as follows with the help of Definition 4.1.
In the r − -th step of the k-th cycle of Stage1, while annihilating the strict lower triangular elements A k r − , to every r − , k there corresponds a pair (j * , i * ) of two indices i * and j * (satisfying 1 ≤ i * < j * ≤ n) which are indices of strictly lower triangular part element largest in magnitude and satisfy From (12) and (13), there are at least two and at most four Q k r − +1 matrices that will be used to annihilate a r − ,k j * ,i * . Also, for Q k r − +1 matrices there are 0 ≤ y ≤ 1 and −3 ≤ y ≤ −2 intervals (given by (10) and (11)). So, in order to annihilate a r − ,k j * ,i * , one must choose from among the Q k r − +1 matrices. The procedure used to find the proper Q k r − +1 matrices during the Stage1 computations is given as follows: transformation matrices by using (12) or (13) roots. b) Exclude both extraneous transformation matrices and the transformation matrices which correspond to the roots not fall into 0 ≤ y ≤ 1 or −3 ≤ y ≤ −2 interval. c) If the number of remaining roots is one then, the transformation matrix which corresponds to this root is used for annihilation. d) If the number of remaining roots is two and fall into same then, choose the transformation matrix that is the closest identity. e) If the number of remaining roots is two and fall into different interval then, choose the transformation matrix corresponding to the root falling in −3 ≤ y ≤ −2 interval. f) If ∆ k r − < 0 then, skip without annihilation in the next steps because the roots are out of the scope. 2) If δ k r − < 0 then, skip without annihilation in the next steps because the roots are complex numbers. On the other hand, in the r + -th step of the k-th cycle of Stage2, while annihilating the strict lower triangular elements of A k r + , to every r + , k there corresponds a pair (j * * , i * * ) of two indices i * * and j * * (satisfying 1 ≤ i * * < j * * ≤ n) which are indices of strictly lower triangular part element largest in magnitude and satisfy D(i * * , j * * ) ≥ 0 Again we can say here, from (12) and (13), there are at least two and at most four Q k r + +1 matrices that will be used to annihilate a r + ,k j * * ,i * * . Also, for the Q k r + +1 matrices there are 0 ≤ y ≤ 1 and −3 ≤ y ≤ −2 intervals (given by (10) and (11)). So, in order to annihilate a r + ,k j * * ,i * * , one must choose from among the Q k r + +1 matrices. The procedure used to find the proper Q k r + +1 matrices during the Stage2 computations is the same except that • If the number of remaining roots is two and fall into different interval then, choose the transformation matrix corresponding to the root falling in 0 ≤ y ≤ 1. item is used instead of (1e)-th item and Corollary 1: At the end of each of the Stage1 and Stage2 of the D-NJLA, the matrix forms obtained are similar to (5) and (7), respectively, since the eigenvalues are not arranged with the same signs consecutively on the diagonal. Stage1 and Stage2 continue sequentially until the matrix A is reduced to the upper triangular form.   (12) or (13) and find M by extracting both extraneous transformations and those corresponding to y roots not fall into (10) or (11) among them if |M | = 1 then use Q k r − +1 which corresponds to this root if |M | = 2 then if both roots fall into (10) or (11) then use Q k r − +1 which is close to identity matrix if one root falls into (10) and the other falls into (11) then use Q k r − +1 which corresponds y root falls into (11) if (|M | = 0 or δ k r − < 0) or ∆ k r − < 0 then skip without annihilation (12) or (13) and find M by extracting both extraneous transformations and those corresponding to y roots not fall into (10) or (11) among them if |M | = 1 then use Q k r + +1 which corresponds to this root if |M | = 2 then if both roots fall into (10) or (11) then use Q k r + +1 which corresponds to be close to identity if one root falls into (10) and the other falls into (11) then use Q k r + +1 which corresponds y root falls into (10) if (|M | = 0 or δ k r + < 0) or ∆ k r + < 0 then skip without annihilation The proof of asymptotic quadratic convergence of the classic symmetric Jacobi method uses monotonicity of the offnorm, i.e., offnorm converges to a limit [14], [23], [24]. For the nonsymmetric case, the offnorm is no longer monotone. So, the possible changes in the moduli of entries in the strictly lower triangular part of the matrix in a single step of k-th cycle of Stage1 will be investigated in detail. Let A k 0 − be the starting matrix, A k r − = (a r − ,k ij ) be the matrix obtained after performing r − step of k-th cycle of the Stage1. Additionally, we denote such that 2 × 2 subdeterminant is negative. Here, ξ k r − is the modulus of the largest element in modulus of A k r − such that 2 × 2 subdeterminant is negative. The modulus of the largest element in modulus of the strictly lower triangular part of A k r − is ζ k r − such that 2×2 subdeterminant is negative. The smallest distance between to two diagonal elements of the matrix A k r − is k r − such that 2 × 2 subdeterminant is negative. And, the minimum difference between the K k r − and L k r − values of the matrix A k r − is δ k r − such that 2 × 2 subdeterminant is negative. In the following, let r − be fixed and suppose we have Assume that the (j * , i * )-element of A k r − is the pivot element of the current ((r − + 1)st) step of k-th cycle of Stage1. We compute the orthogonal matrix that is as close as possible to identity among all matrices satisfying (27) and which corresponds to the rotation matrix whose rotation angle falls into [−3, −2] as much as possible.
To obtain an upper bound for the modulus of 3 − 2y − y 2 2 , we will use the following lemma, that is similar to Theorem V.2.1 in [11].
where 4 Using this lemma and (25), we find that Using the above lemma and the following theorem, we will show that k-th cycle of the Stage1 is quadratically convergent.
Theorem 1: after perfoming r − Stage1 steps of D-NJLA. Let ζ k 0 − and ζ k µ − be the largest modulus of an element in the strictly lower triangular part of A k 0 − and A k µ − that is obtained from A k 0 − after having completed µ − sweeps such that 2 × 2 subdeterminant is negative, i.e., and let φ k 0 − and φ k µ − be the largest modulus of an element in such that 2 × 2 subdeterminant is negative. Morover, set N = n 1 n 2 , where n 1 and n 2 are the number of negative and positive signed eigenvalues of A, respectively and If δ k − > 0 and if ψ k 0 − is sufficiently small such that then for all µ − ∈ IN ∪ {0}, we have that i.e., k-th cycle of the Stage1 of D-NJLA sweep converges quadratically over the number of sweeps.
Similarly, it can be said k-th cycle of the Stage2 of D-NJLA sweep converges quadratically over the number of sweeps. So, we can say cycles of the Stage1 and Stage2 converge quadratically, i.e., the D-NJLA converges quadratically.
a) Complexity of D-NJLA: D-NJLA calculates the eigenvalues of nonsymmetric matrices with solvable directed graphs. At each step of the Stage 1 of the D-NJLA, an entry, which has the property that the strictly lower triangular part element has the largest magnitude and the sign of 2 × 2 subdeterminant is negative, is chosen. At each step of the Stage2, an entry, which has the property that the strictly lower triangular part element has the largest magnitude and the sign of 2 × 2 subdeterminant is zero or positive, is chosen. Therefore, n(n − 1) 2 elements must be compared in each step of both Stage1 and Stage2. So, D-NJLA's complexity is O(n 2 ).

VI. EXPERIMENTAL RESULTS
We implemented the D-NJLA in MATLAB and run it on an Intel Core i7 processor (2.7 GHz) computer. An off-diagonal element a i,j is considered zero if | a i,j |< eps, where eps is a machine dependent constant. As a stopping criteria we used of f (A) = i<j≤n a 2 i,j = constant.
We illustrate the D-NJLA with the following numerical examples. We give the "effective" number of sweeps of the Stage 1 in k-th cycle, that is the number of rotations divided by n 1 n 2 . "Effective" number of sweeps of the Stage 2 in k-th cycle is the number of rotations divided by n 1 (n 1 − 1) 2 + n 2 (n 2 − 1) 2 . Here, n 1 and n 2 are the number of negative and positive signed eigenvalues of the given matrix, respectively.
We take matrices from [2], which are 7 × 7 'Grund/b1 ss', 9 × 9 'vanHeukelum/cage4', 102 × 102 'MathWorks/pivtol', 135×135 'Rajat/rajat11' and 216×216 'Fidap/ex1' and reduce them to upper triangular form.   The matrices given as an example above have a large number of very close eigenvalues. Therefore, as the number of steps increases in Stage2, r + increases (i.e, δ r + decreases) considerably, so N √ 2 increases considerably. Thus, the rate of convergence is expected to slow down significantly. But, as can be seen from Fig. 2 and Fig. 3, no significant slowdown in the convergence rate in the Stage2 was observed.
In all figures, no stationary phase was observed in both the Stage1 and Stage2. As expected, it can be observed maxof f (A) decreases over the number of sweeps of the Stage1 and Stage2. As predicted by the theory, the convergence rate of the Stage1 and Stage2 becomes asymptotically quadratic.
Also, we take matrices, some of whose eigenvalues are complex and whose the number of eigenvectors is less than the order of the taken matrix. These matrices are reduced nearly upper triangular form by D-NJLA and the algorithm converges linearly for this type of matrices because the roots have not been taken into consideration for ∆ k r − < 0 and ∆ k r + < 0. Additionally, the matrices, which have a zero eigenvalue are reduced upper triangular form by the D-NJLA and the convergence behaviour of the D-NJLA is asymptotic quadratic and no stationary phase is observed.

VII. CONCLUSIONS
In this paper, we have presented the details of D-NJLA, which is a Jacobi-like algorithm. Moreover, we also present numerical experiments with the D-NJLA, that is applied to some examples of real nonsymmetric matrices to reduce them to upper triangular matrices. As seen the given examples, the convergence behavior of the D-NJLA is asymptotic quadratic and no stationary phase is observed, although the matrices given as an example have a large number of very close eigenvalues.
Still, while we encouraged that the D-NJLA works as it does, there are many questions unanswered. Firstly, a proof of global convergence is still not provided. The D-NJLA is not adapted for matrices with complex eigenvalues. Therefore, the cases of ∆ k r − < 0 and ∆ k r + < 0 or the equivalences of (37) and (38) should be take into consideration. So, the algorithm can be adapted for this type matrices. Also, considering the matching between the independence set of graphs and groups [1], a group theory approach can be made to the algorithm.

APPENDIX A ORTHOGONAL TRANSFORMATION MATRIX
The corresponding (n×n) orthogonal transformation matrix Q k r , analogous to the rotation of an axis in the (j, i) plane may be given by: which corresponds to the transformation matrix of the Jacobi method [3]. For more details on the variance of the transformation matrix see [21]. When A is a nonsymmetric matrix, the angle between some eigenvectors can be quite small. In this case, as stated above, the numerical computation of the angle θ, which depends on angles between eigenvectors, can be more delicate [8].
To overcome this problem, we expressed cos θ and sin θ with quadratic polynomials as explained below and the rotation parameter λ 2 corresponding to the rotation angle θ was taken in the following intervals. The characteristic equation of (35) is given by: det(R r − λI) = (1 − λ) n−2 (λ 2 − 2λ cos θ + 1) = 0. (36) To compute the eigenvalues of (36), we set either (1−λ) n−2 or (λ 2 − 2λ cos θ + 1). To avoid complex operations, we consider (λ 2 − 2λ cos θ + 1), hence the rotation angle may be written in terms of λ parameter as: In this case, substituting (37) to sin θ = √ 1 − cos 2 θ, we obtain which is a complex number. To determine the real interval for θ, without taking λ in the denominator of (37) we use which determines rotation angle θ of transformation matrix in terms of the λ 2 parameter is obtained. Hence, for new transformation matrix, we have: where y = λ 2 , c = y + 1 2 and s = 3 − 2y − y 2 2 . So, the angle θ is chosen for appropriate the y parameter, thus, we have 0 ≤ y ≤ 1 and −3 ≤ y ≤ −2.

APPENDIX B CALCULATION OF THE ROTATION PARAMETER y
Let (a r,k j,i ) be any element to be annihilated in Stage1 or Stage2 such that r is any step number in k-th cycle. The element (a r,k j,i ) is annihilated by a new similarity transformation matrix Q k r+1 . This annihilation affects the entries in rows and columns j and i only. The modified a r+1,k j,i element is given by: The y parameter is chosen to annihilate a r+1,k j,i element. So, we obtain a corresponding equation for the y parameter (a r,k j,j − a r,k i,i ) (y + 1) 3 − 2y − y 2