KRONECKER ALGEBRA-BASED DEADLOCK ANALYSIS FOR RAILWAY SYSTEMS

Deadlock analysis for railway systems differs in several aspects from deadlock analysis in computer science. While the problem of deadlock analysis for standard computer systems is well-understood, multi-threaded embedded computer systems pose new challenges. A novel approach in this area can easily be applied to deadlock analysis in the domain of railway systems. The approach is based on Kronecker algebra. A lazy implementation of the matrix operations even allows analysing exponentially sized systems in a very efficient manner. The running time of the algorithm does not depend on the problem size but on the size of the solution. While other approaches suffer from the fact that additional constraints make the problem and its solution harder, our approach delivers its results faster if constraints are added. In addition, our approach is complete and sound for railway systems, i.e., it generates neither false positives nor false negatives.


INTRODUCTION
The deadlock problem and its solutions were studied in the earliest days of computer science.Although computer science borrowed several concepts and namings from the railway domain such as semaphores or tokens, this was not the case for the deadlock problem.Deadlocks must have been impending in railway systems from the very beginning and railwaymen certainly were aware of them.In contrast to these facts, no solutions were commonly known in the middle of the 20 th century.So computer scientists were the first to study the problem intensively and to deliver adequate solutions.
The following definition of a deadlock is given in [1]: Deadlock: An impasse that occurs when multiple processes are waiting for the availability of a resource that will not become available because it is being held by another process that is in a similar wait state.
For a deadlock to occur, four necessary conditions were found [2]: 1.Mutual exclusion: a resource that cannot be used by more than one process at a time.2. Hold and Wait: processes already holding resources may request new resources held by other processes.3.No preemption: No resource can be forcibly removed from a process holding it; resources can be released only by the explicit action of the process.4. Circular wait: two or more processes form a circular chain where each process waits for a resource that the next process in the chain holds.If one of these conditions cannot be met, a deadlock cannot occur.Clearly all four conditions apply to railway systems.
To handle deadlock problems one distinguishes between deadlock prevention, deadlock avoidance, and deadlock detection.
-Deadlock prevention tries to remove one of the four conditions above in order to make it impossible for a deadlock to occur.-Deadlock avoidance utilizes certain information about processes known in advance to decide whether or not to allocate resources to processes.
Examples for deadlock avoidance algorithms include the Banker's algorithm [3], Wait/Die, Wound/ Wait algorithms, and [4] to name a few.
-When a deadlock is detected at runtime, a process is terminated and restarted later on.Deadlocks in railway systems attracted attention when railway simulation systems became available.Since simulations were supposed to run without human interaction, deadlocks became an obstacle.Obviously, autonomous dispatching systems also suffer from deadlock problems if no adequate solutions can be found.
The outline of our paper is as follows.In Section 2 we contrast deadlocks in computer systems with deadlocks in railway systems.Section 3 introduces Kronecker algebra.Section 4 gives our standard system model.In Section 5 we present a simple example.Section 6 shows how our resulting graphs can be employed for finding optimal solutions.In Section 7 we describe our lazy implementation of Kronecker operations.In Section 8 we extend our standard model in order to solve several important practical problems.In Section 9 we present experimental data.In Section 10 we relate our approach to other work, before we conclude our paper in Section 11.

DEADLOCKS IN COMPUTERS VS. DEADLOCKS IN RAILWAY SYSTEMS
In railway systems trains, routes, and track sections correspond to processes, program code, and shared resources in computer systems, respectively.
The table below shows major differences between the deadlock problem in computer and railway systems.
From the table below it becomes obvious that deadlock analysis methods and algorithms well-suited for standard computer systems cannot be simply transferred to railway systems.However, safety-related embedded systems show more resemblance to railway system.This is the reason why we try to apply the approach we have developed for embedded systems [4] to railway systems.We will show in this paper that this approach is in fact well-suited for railway systems.
In the following we will argue that only deadlock avoidance makes sense in the domain of railway systems; deadlock prevention and deadlock detection cannot be applied in this domain.
Deadlock prevention is done by ensuring that one of the four necessary conditions stated above is not fulfilled.It is easy to see that this cannot be done for railway systems: 1.Mutual exclusion: Track sections have to be seized by trains exclusively.2. Hold and Wait: It is typical that a train occupies two or more track sections at the same time.3.No preemption: A track section cannot be removed from a train "holding" it.4. Circular wait: It is impossible to impose a strict order on the track sections without posing the problem that some trains would seize all track sections of their route from start to destination before they start running.Deadlock detection is of no interest for railway systems because "terminating a train" is no option.
The only remaining option is deadlock avoidance which will be discussed in the rest of this paper.
Before we give an introduction to Kronecker algebra, we state a few additional differences between computer and railway systems: -Computer systems usually consist of a number of processes (threads) and a number of shared resources.In most cases the number of resources is much smaller than the number of processes.In railway systems the number of track sections is usually much greater than the number of trains.-Besides the track sections being shared resources in railway systems there are additional constraints that have to be met.For example, connections and overtaking have to be taken care of, tardy trains may cause penalty, ... While previous approaches perform poorly if additional constraints are added to the problem domain,

Computer systems Railway systems
In standard computer applications all resources are held until the process terminates (not true for embedded systems where processes may not terminate at all).
A train seizes a track section shortly before it enters it, and holds it until it leaves it (which may be long before it reaches its destination).
Program code consists of straight line code, if, and loop statements.
Routes compare to straight line code; there may be alternative routes which compare to if statements; loops may be useful for model railroad applications.
If a multithreaded computer program contains a deadlock, the program is erroneous.
If in a railway system it is possible to bypass a deadlock situation, then the system is correct.Deadlocks are always possible, but can usually be bypassed.
A computer program is "correct" if we can prove it deadlock-free.
A set of trains is schedulable if all possible deadlocks can be bypassed.
If a deadlock is detected, it is viable to terminate a program and restart it later on (not true for safety-related embedded systems).
Trains cannot simply be taken out of the overall system.
R. Mittermayr, J. Blieberger, A. Schöbel: Kronecker Algebra-based Deadlock Analysis for Railway Systems our approach delivers its results faster if constraints are added.This makes it a very promising approach.Details will be given in the following sections.

KRONECKER ALGEBRA -AN OVERVIEW
Kronecker product and Kronecker sum form Kronecker algebra.In the following we define both operations, state properties, and give examples.
We define the set of matrices , where L denotes a set of labels such that , ,0 L + is a commutative monoid and , ,1 L $ is a monoid.In the remaining parts of this paper only matrices M M ! will be used, except where stated explicitly.Let o M ^h refer to the order 1 of matrix M M ! .In addition we will use n-by-n zero matrices Z z, n ij = ^h, where , : i j z 0 , i j

6
= .Definition 1 (Kronecker product).Given an m-by-n matrix A and a p-by-q matrix B, their Kronecker product -denoted by A B 7 -is an mp-by-nq block matrix defined by As stated in [5] the Kronecker product is also being referred to as Zehfuss product or direct product of matrices or matrix direct product. 2 Next we list some basic properties of the Kronecker product.Proofs and additional properties can be found in [6], [7], [8], [9] This operation must not be confused with the direct sum of matrices, group direct product or direct product of modules for which the symbol 5 is used, too.
In the following we list basic properties of the Kronecker sum of matrices A, B, and C. Additional properties can be found in [10].The Kronecker sum is non-commutative because for element-wise comparison in general A B B A 5 5 !
. Anyway, it essentially commutes because from a graph point of view, the graphs represented by matrices A B 5 and B A 5 are structurally isomorphic.In addition, we state a property of the Kronecker sum which we call Mixed Sum Rule.
Let the matrices A and C have order m and B and D have order n.Then we have A proof of the Mixed Sum Rule can be found in [4].
In addition, the Kronecker sum is associative (cf.[4]), i.e., A B C A B C 5 5 5 5 = ^ĥ h .The associativity properties of the operations 7 and 5 imply that the k-fold operations The Kronecker sum calculates all possible interleavings (see e.g.[11] for a proof).The following example illustrates interleaving and how the Kronecker sum handles it.
Example 3: Let matrices C and D be defined as follows: As already mentioned above, there is a correspondence between matrices and graphs.In general, a directed labelled graph , , G V E ne consists of a set of labelled nodes V, a set of labelled directed edges , and an entry node ne.In this paper we label the graph nodes simply by positive integers.The correspondence between graphs and matrices -frequently called adjacency matrices -is as follows.If there exists an edge labelled a from node i to node j, then the corresponding adjacency matrix M has m a , i j = .If there is no edge between node i and node j, then m 0 , i j = .The graph corresponding to matrix C from above is depicted in Figure 1, whereas the graph of matrix D is shown in All possible temporal interleavings of these actions are shown in Table 1.In Figure 3 the graph represented by the adjacency matrix C D 5 is depicted.It is easy to see that all possible temporal interleavings are generated correctly.
In general, by calculating the Kronecker sum of the adjacency matrices of two graphs the adjacency matrix of the Cartesian product graph [12] is computed (cf.[13]).Now assume that Ta and Tb denote the same track section.It is clear that in this case the temporal interleavings of Table 1 are no more valid.The trains have to synchronise in order to perform their actions correctly.This can be modelled by Kronecker product and an additional matrix of the form where p denotes the action "Enter the track section" and v means "Train has left the track section".The correct system behaviour can be described by matrix . The corresponding graph is shown in Figure 4.It is interesting that the graph became decomposed into seven parts (sub-graphs).Clearly, only the part reachable from the entry node is responsible for the system behaviour, the other six parts can safely be ignored.Thus we concentrate on the sub-graph with entry node 1. Studying this sub-graph we see that now the trains enter the track section one after the other.Note that the two paths in the subgraph correctly mirror the two cases where C enters the track section before D and vice versa.A proof that Kronecker product models synchronisation correctly can be found in [4].It is also worth noting, that parts unreachable from the entry node always occur if synchronisation via Kronecker product takes place.This observation gives rise to the lazy implementation described in Section 7. Table 1 Interleavings

SYSTEM MODEL
We model a general railway system S by a set of track sections T T i r , .Each track section Ti is modelled by matrix In addition, a railway system consists of a set of trains L L j t , .The route Rj of train Lj is a sequence of track sections , , T T l l s , .The behaviour of railway system , , , where during the evaluation of the Kronecker product we let p p p . The different paths in the graph corresponding to matrix S mirror all possible behaviours of the railway system in terms of temporal interleavings of the actions of trains, namely entering and leaving track sections.
From the discussion above and from [4] it is clear that deadlocks appear as purely structural properties of the underlying graph.Deadlocks manifest themselves as non-final nodes 4 with no successors.

A SIMPLE EXAMPLE
In this section we give a small example on how deadlocks can be avoided by our Kronecker algebrabased approach.
In Figure 5 a typical scenario is shown that may lead to a deadlock.The routes of the three involved trains are given in Table 2.If train L3 enters track section T3 before the other trains move, a deadlock is unavoidable.
The matrices for the five track sections have the form The order of the resulting matrix can be computed by multiplying the order of the involved matrices, i.e., the order of matrices Rj and the order of matrices Ti .In our example the resulting matrix will have order 6912 6 6 6 2 5  $ $ $ = .Due to synchronisation only a small part of the resulting graph is reachable from the entry node.Only this part -consisting of 42 nodes -is displayed in Figure 6.To increase readability our implementation distinguishes between the following node types: -Diamond nodes denote deadlocks or nodes from which only deadlocks can be reached.-Solid nodes denote safe states.A state is safe if all trains can perform their actions without having to take into account the moves of other trains in the system, provided that the track section which they are to enter is not occupied by another train. 5-From dotted nodes both diamond and solid nodes can be reached.
If only deadlock avoidance is of interest, graphs like that in Figure 6 can be scaled down.In fact, diamond nodes are of no interest because we want to avoid deadlocks.From the set of solid nodes only those are of interest, which have dotted predecessors only.In addition, we do not have to show the node labels, i.e., the matrix line numbers.If we eliminate the non-interesting nodes from Figure 6, we obtain the graph in Figure 7.We call such graphs NDLS-graphs because they contain no deadlock nodes and only a limited number of safe nodes.
In order to avoid deadlocks the minimal number of actions can be obtained by selecting a path from the entry node to a solid node which is situated nearest to the entry node.This, however, may not be the best alternative.In the next section we elaborate on the problem of finding the "cheapest" path.
In order to evaluate the quality of algorithms or tools it is convenient to introduce false positives and false negatives.Let us assume we have an algorithm at hand that examines a system to find out whether the system deadlocks or not.Then there are four possible outcomes.complete and sound.There is only one exception: In case of alternative routes (cf.Section 8) deadlocks in all possible paths will be reported although only one path can be chosen by the train.This situation may somehow be classified as reporting false positives.
For computer systems our approach can be modified such that it stops as soon as the first deadlock is found 7 .For railway systems the algorithm may be modified to stop as soon as a safe state is found.

AVOIDING DEADLOCKS UNDER ADDITIONAL CONSTRAINTS
Our approach produces a graph which gives all temporal interleavings of train moves.A short reflection shows that not all of these interleavings may result in an "optimal" behaviour.However, a more sophisticated analysis may choose to assign weights to the edges of the graph.These weights can be chosen such that the term "optimal" from above makes sense.Weights may include temporal constraints originating from physical properties of the trains such as accelerating and braking characteristics depending e.g. on the weight of the train, or penalties due to tardiness.
In any case, well-known algorithms such as those of [14] or [15] can be applied to graphs with nonnegative weights to find the optimal (shortest) path.Dijkstra's algorithm performs in time quadratic in the number of nodes of the graph.The Fredman & Tarjan version uses time proportional to the sum of the number of edges and the number of nodes times the logarithm of the number of nodes.
If only deadlock avoidance is of concern, the optimisations given in the previous section apply.

LAZY IMPLEMENTATION OF KRONECKER ALGEBRA
Until now we primarily focused on a pure mathematical model for deadlock analysis.An alert reader will have noticed that the order of the matrices increases exponentially in the number of trains.On the other hand, we have seen that the Kronecker product results in parts of the resulting matrix that cannot be reached from the entry node of the corresponding graph.This comes solely from the fact that synchronisation excludes some interleavings.In addition, all involved matrices are sparse (cf.[4]), such that space saving data structures suggest them selves.
Choosing a lazy implementation [16] for the matrix operations, however, ensures that, when extracting the reachable parts of the underlying graph, the overall effort is reduced to exactly these parts.By starting from the entry node and calculating all reachable successor nodes, our lazy implementation does exactly this.Thus, for example, if the resulting graph's size is linear in terms of the involved trains, only linear effort will be necessary to generate the resulting graph.
Our implementation distinguishes between two kinds of matrices: Sparse matrices are used for representing trains and track sections.Lazy matrices are employed for representing all the other matrices, i.e. those resulting from the operations of the Kronecker algebra.Besides the employed operation, a lazy matrix simply keeps track of its operands.Whenever an entry of a lazy matrix is retrieved, depending on the operation recorded in the lazy matrix, entries of the operands are retrieved and the recorded operation is performed on these entries to calculate the result.In the course of this computation, even the successors of nodes are evaluated lazily.Retrieving entries of operands is done recursively if the operands are again lazy matrices, or it is done by retrieving the entries from the sparse matrices, where the actual data resides.
In addition, our lazy implementation allows for simple parallelising.For example, retrieving the entries of left and right operands can be done concurrently.Exploiting this, we expect further performance improvements for our implementation if run on multicore architectures.All examples in this paper have  been computed by our lazy implementation.The performance of our implementation is discussed in Section 9.

EXTENSIONS OF THE STANDARD MODEL
In this section we discuss how our standard model defined in Section 4 can be extended in order to solve certain important problems.
The first problem we consider occurs if the length of a train is greater than the length of a track section.Assume a route of a train containing three consecutive track sections t1, t2, and t3, where t2 is shorter than the length of the train.In addition, assume that the train occupies track section t1.If t2 were long enough, the train's moves would be , , , , p v p v 2 1 3 2 f In case of the short track section the train's moves are , , , , p p v v 2 3 1 2 f It is clear that this approach can easily be extended if more than three track sections are used by a train simultaneously.
Furthermore, note that double slip switches can be seen as two normal switches with a zero length track section in between.Thus, routes containing double slip switches represent typical examples of the method above.
Sometimes it is useful to model a railway system from a more abstract view.For example one abstracts away from the track details of a station and instead would like to use only the fact that the station has a capacity of holding c trains at the same time, where each train is able to enter and leave the station independently from the other trains in the station.Such a station can be modelled by a c c Thus, our approach is able to handle abstract views as well.One can envisage even more abstract levels such as abstracting whole railway lines to their capacity.This allows to model nets of railway lines.
Sometimes it is useful when a train is allowed to take alternative routes.This gives more freedom in bypassing a potential deadlock situation.Since our approach was developed for handling embedded computer systems, it is able to handle alternatives which are present as if-statements in computer programs.Generating matrices that model alternative routes is very easy.A line in the matrix representing a fork state of a route has two successors, i.e., two non-zero entries.Later on the different alternatives can be joined and a single (straightforward) route continues.Clearly, this works even for more than two alternative routes.
As already noted, our approach was developed for embedded computer systems.Thus, it also supports loop statements.Matrices modelling loops can be employed to model trains running in loops.This might be of interest for model railroaders.
By adding further matrices we can model synchronisation of trains which can be used to ensure connections and overtaking of trains.For example, assume that we want to ensure that train A shall not leave track section tA before train B leaves track section tB.We introduce two "artificial track sections" d and e for synchronisation with the usual operations , , and ve.In addition, we insert immediately in front of the moves vA and vB the actions , v p d e and , p v d e, respectively.Thus, the two trains are synchronised and we are able to model connections correctly.
We conclude this section with a more elaborate example containing some of the extensions described above.In particular, it contains a double slip switch (track section 7) and it shows how overtaking can be modelled.The system is depicted in Figure 8.The routes are defined as follows: , , , , , , , , , , R p v p p v v p v p v v  , , , , , , , , , , , ,  R , , , , , , , , , , , , , ,  R  = .Note that track section 7 has zero length; it originates from the double-slip switch located between track sections 3, 9, 6, and 8.

p v p v p v p p v v p v p v v
Kronecker algebra calculations produce a graph with 55,050,240 potential nodes.After the reductions presented in Section 5 only 992 nodes are left.Now we introduce additional constraints.We assume that train L2 overtakes L1 and L4 overtakes L3 within the station.We need one additional "artificial track section" for each train pair.The corresponding operations are denoted by , , , p v p 12 12 13 and v13.Now the routes read as follows: , , , , , , , , , , , R p v p p v v p v p v v v

=
. Kronecker algebra produces 298,721,280 potential nodes.After the usual reductions, however, only 144 nodes are left (cf. Figure 9).This means that the number of potential nodes has grown by a factor of five, but the size of the solution has dropped by factor of approximately 6.8.The resulting NDLS-graph is depicted in Figure 9.
This example shows that if constraints are added to a railway system, our approach delivers its solution faster.This is rather atypical for recent approaches.

EXPERIMENTAL DATA
In this section we give experimental data for ten examples.Let o S ^h and o C ^h refer to the order of the adjacency matrix S, which is not computed by our lazy implementation, and the order of the adjacency matrix C of the resulting graph, respectively.o NDLS ^h denotes the order of the NDLS graph.In addition, k and r refer to the number of trains and the number of track sections, respectively.
Next we refer to the data depicted in Table 3.The runtime of our implementation shows a strong correlation to order o C ^h of the adjacency matrix C of the resulting graph with a Pearson product-moment correlation coefficient of 0.999027713.The correlation coefficient between o NDLS ^h and the runtime is 0.9813607056.In contrast, the values of the purely theoretical order o S ^h of the resulting adjacency matrix S correlate to the runtime only with a correlation coefficient of 0.2370050995.These observations show that the runtime complexity does not depend on the order o S ^h which grows exponentially in the number of trains.We conclude this section by stating that the collected data give strong indication that the runtime complexity of our approach is linear in the number of nodes present in the resulting graph.worked on generating reachability sets in composed automata and for describing networks of synchronized automata.Both approaches employ Kronecker algebra.
In [18] we used similar definitions of Kronecker algebra.The approach is used for analysing the timing behaviour of concurrent programs.
Although not closely related we recognize the work done in the field of stochastic automata networks (SAN) which is based on the work of [19] and in the field of generalized stochastic Petri nets (GSPN) (e.g.[20]) as related work.Compared to ours these fields are completely different.Nevertheless, basic operators are shared and some properties influenced this paper.
In [21] Banker's algorithm [3] is modified such that it can be employed for deadlock analysis in railway systems.Since Banker's algorithm has been designed for standard computer systems it is not well-suited for railway systems.For example, it may prohibit allocating a resource (track section) although a potential deadlock can be bypassed.In contrast to our approach both track sections and switches have to be modelled as resources in [21].
An operations research approach is proposed in [22] to do deadlock analysis in railroad systems.
In [23] Movement Consequence Analysis (MCA) and Dynamic Route Reservation (DRR) are introduced for deadlock analysis.Both are rule-based methods for which correctness cannot be proved.It delivers false positives.
An algorithm which can handle only simple railway networks is given in [24].In addition, it produces false positives.The same restrictions apply to the approach presented in [25].
Deadlock-free algorithms are presented in [26].All of them may generate false positives.
A coloured Petri net model is used in [27] to describe a railway network system and to derive the traffic controller.Safeness and deadlock freedom are guaranteed.A deadlock prevention strategy is defined and expressed by a set of linear inequality constraints.It is shown how collision and deadlock prevention constraints can be expressed as coloured generalised mutual exclusion constraints and that the controller can be implemented by a set of monitor places.It is however not clear how this approach can be extended to find optimal solutions such as our approach suggests (cf.Section 6).
In [28] a coloured Petri net model of a simple railway station operation is constructed and Banker's algorithm is employed for deadlock avoidance.It turned out that Banker's algorithm is not able to handle a large number of processes (trains) and resources (track sections).

CONCLUSION
We have presented a Kronecker algebra-based approach for deadlock analysis in railway systems.It can solve all important practical problems.By employing a lazy implementation for the Kronecker matrix operations the solution is computed very efficiently.
By assigning weights to the edges of the resulting graph various optimisations can be performed.This includes optimisations due to physical properties of the trains such as accelerating and braking characteristics depending e.g. on the weight of the train, or penalties due to tardiness.
The running time of the algorithm does not depend on the problem size but on the size of the solution.While other approaches suffer from the fact that additional constraints make the problem and its solution harder, our approach delivers its results faster if constraints are added.The identity matrix In is an n-by-n matrix with ones on the main diagonal and zeros elsewhere.

4.
A final node corresponds to the destination of a route.
A final node of a railway system corresponds to the state, where all trains have reached their destinations.5.
If a track section is occupied by another train, the movement of the train wanting to enter may be delayed, but no deadlock can occur.6.
This is true for railway systems.It is not true for computer systems because computer programs may e.g.contain dead code, i.e., code that is never executed.Our approach will report a deadlock even if it is contained in dead code.Thus it will deliver false positives for computer systems but no false negatives.7.
This makes sense because a computer system is considered erroneous if at least one deadlock is present.

Example 2 :
We use matrices A and B from Example 1.The Kron ecker sum A B 5 is given by A

Figure 2 .
Now interpret C and D as being trains and , , a b c and d as being actions of the trains with the following meaning: a denotes train C enters track section Ta, b means train C has left track section Ta, c denotes train D enters track section Tb, and d means train D has left track section Tb.

Figure 7
Figure 7 Figure 8 , , , , , , , , , , Rp v p v p p v v p v p v v p

Figure 9
Figure 9 Im and In denote identity matrices3of order m and n, respectively.

Table 3
[17]erning Kronecker algebra probably the closest work to ours was done by Buchholz and Kemper[17].It differs from our work in that we use our approach to avoid deadlocks.Buchholz and Kemper [29]lock-Analyse für Eisenbahn-Dispositionssysteme unterscheidet sich von Deadlock-Analyse in der Informatik in mehreren Punkten.Während Deadlock-Analyse für Standard-Computersysteme im Allgemeinen gut verstanden ist, stellen Eingebettete Systeme eine neue Herausforderung dar.Ein neuer Ansatz aus diesem Bereich kann auf Eisenbahn-Dispositionssysteme angewandt werden.Der Ansatz basiert auf Kronecker Algebra.Eine "Lazy Implementierung" der Matrixoperationen erlaubt sogar die Analyse von exponentiell wachsenden Systemen in sehr effizienter Art und Weise.Die Laufzeit des Algorithmus hängt nicht von der Problemgröße sondern von der Größe des Ergebnisses ab.Während andere Ansätze durch zusätzliche Einschränkungen langsamer werden, liefert unsere Methode das Ergebnis schneller, sobald weitere Einschränkungen hinzukommen.Weiters ist unser Ansatz vollständig (complete) und fehlerfrei (sound), d.h.es treten weder "false positives" noch "false negatives" auf.by-k matrix is known as square matrix of order k.2.Knuth notes in[13]that Kronecker never published anything about it.Zehfuss was actually the first publishing it in the 19 th century[29].He proved that if A and B are matrices of order m and n, respectively, with entries from the domain of real numbers.3.