In-Place Matrix Inversion by Modified Gauss-Jordan Algorithm

The classical Gauss-Jordan method for matrix inversion involves augmenting the matrix with a unit matrix and requires a workspace twice as large as the original matrix as well as computational operations to be performed on both the original and the unit matrix. A modified version of the method for performing the inversion without explicitly generat-ing the unit matrix by replicating its functionality within the original matrix space for more efficient utilization of computational resources is presented in this article. Although the algorithm described here picks the pivots solely from the diagonal which, therefore, may not contain a zero, it did not pose any problem for the author because he used it to invert structural stiffness matrices which met this requirement. Techniques such as row/column swapping to handle off-di-agonal pivots are also applicable to this method but are beyond the scope of this article.


Introduction
Gauss-Jordan [1] is a standard matrix inversion procedure developed in 1887. It requires the original matrix to be appended by a unit (identity) matrix and after the inversion operation is completed the original matrix is transformed into a unit matrix while the appended unit matrix becomes the inverse.
During the early days of his career as a professional engineer and software developer [2], the author created an engineering software library for the Government of India in the mid-1960s. Some of the softwares were in the field of structural analysis which involved computing displacement vectors for multiple load vectors applied to the same stiffness matrix created using the modeling technique suggested by Tezcan [3].
Since his work involved handling large numbers of loading cases applied to the same structure, he opted for matrix inversion to generate flexibility matrices instead of solution of equations to process the load vectors because the latter method would require either regeneration or retrieval of the stiffness matrix to process every load vector, either of which would substantially increase the computational time and effort.

Development
The stiffness matrices encountered in structural analysis have three important characteristics. Their largest elements always fall on the diagonal which can never have a zero element and they are always symmetrical. The author started with the Gauss-Jordan algorithm however, since both computing power and hardware speed were severely limited in those days, he attempted to carry out as much of his work within the computer memory (incore) as possible. Since Gauss-Jordan required an augmenting unit matrix and therefore both memory space and amount of computation for a matrix twice as large as the original matrix, he sought a way around this problem.
Since after a pivotal column of the original matrix is processed with the Gauss-Jordan method, it always contains a unity (1) on the diagonal and zeroes in all other rows and remains unchanged during all subsequent operations without any mathematical utility, the author felt that this memory space could be utilized to replicate the corresponding column of a virtual unit matrix without creating a real unit matrix. He used this idea to develop a modified procedure to simulate a non-existent augmenting unit matrix and called it virtualization of the aug-D. DASGUPTA 1393 menting matrix. It allowed for the inversion procedure to be carried out within the memory space of the original matrix which he called the In-Place Inversion method.
Since it made the entire memory space available to the original matrix, it enabled him to invert matrices twice the size possible with Gauss-Jordan with the available resources and also reduced the computational time and effort because its operations were confined to the original matrix space.
Although the author first used the algorithm to invert relatively small matrices within the computer memory (in-core), he later also used it to invert large matrices such as encountered in Finite-Element analysis which involved disk operation. It increased computational efficiency because not only the computation was confined to a smaller number of elements but the length of the read/ written vectors was the same as the row width of the original matrix and not twice as long as would have been with Gauss-Jordan. It was also particularly useful for PCbased applications.
Since this method uses the same underlying mathematics as Gauss-Jordan and can be enhanced with the same techniques applicable to it, it can be used wherever Gauss-Jordan is used. It produced identical results as Gauss-Jordan as shown in the examples cited in this article and a 5 × 5 Hilbert matrix [4] was also successfully inverted.
The author has used certain terms in the following discussion which are defined as follows: Normalization is dividing an entire row by its pivotal element to transform the pivotal element to unity; Virtualization is replicating an element or a vector of the current augmenting matrix within the original matrix space without creating the real unit matrix; the Complementary of a component of the original matrix is its corresponding component in the virtual augmenting matrix and the Reduction of a row is the modification of the pivotal row by the ratio of the row element on the pivotal column and the pivotal element and then subtracting it from the row, thereby reducing the row element on the pivotal column to zero.

Mathematics
Gauss-Jordan is a standard matrix inversion procedure as outlined below: For a matrix   A of size n × n, an identity (unit) matrix of size n x n is appended to the matrix. After that the following two operations are iterated on all rows to obtain the inverse.
Operation 1: A pivotal row is selected, usually the one with the largest absolute diagonal element, from the rows whose diagonal elements where j = 1 → 2n. This is a shorter equivalent of the more explicit operation , A  after normalization of the pivotal row. This transforms its pivotal column element , i p A to zero and after it is performed on all non-pivotal rows, their pivotal column elements become 0 (zero), i.e., , After performing these two operations on every row, treating each row once as a pivotal row, the original matrix becomes a unit matrix while the unit matrix becomes the inverse.
Computational details of the inversion of a 3 × 3 matrix [5] are given below in which the original matrix has been multiplied with its inverse in the end to produce a unit matrix as verification of the result. The author's In-Place Inversion algorithm, on the other hand, does not require augmenting with and performing operations on an identity matrix and the procedure is described below:

Gauss-Jordan Inversion-Example
Just as with Gauss-Jordan, the following two operations are iterated on all rows to obtain the inverse.
to virtualize the entire complementary column in the virtual unit matrix for explanatory purposes. In reality, however, this part is performed individually for each non-pivotal row as part of the next operation to preclude the need for saving them all in a vector.
This procedure implicitly duplicates the functionality of the unit matrix of the Gauss-Jordan method within the original matrix. After performing these two operations on every row, treating each row once as a pivotal row, the original matrix is replaced by its inverse.
Notes: In the example below the items titled Saved A p,p and Saved A i,p show the saved pivotal values p P and all i individually although they need only a single memory location because, in reality, they are generated sequentially and not all at the same time as a complete vector.

P
The sequence of operations 1 and 2 can be reversed but in that case the pivotal element , p p A will not be unity during operation 2 and the explicit formula will have to be used, thereby substantially increasing the amount of computation.
Computational details of the same 3 × 3 matrix are given below for comparison. The verification step has been omitted since the result is identical to that obtained by Gauss-Jordan which has already been verified.

Practical Applications
The author's initial application for the method was as the analysis engine for static and dynamic analysis of 2D and 3D frame structures as well as finite-element analysis. He used it for both small and large-scale tasks including plane strain finite-element analysis of a gravity dam under El Centro seismic loading (some disk operation was involved) which produced results closely matching those published by Clough and Chopra [6]. He later incorporated it in a number of IBM-PC/MS-DOS based structural engineering softwares which have been used by both Federal and state governmental agencies as well as many engineering design and production firms in the U.S. [7]. Although this algorithm was used first in India and then in the U.S. for almost five decades, its logic has remained undocumented so far and the author's objective is to present it to other software developers who may find it useful.
The author extensively searched all sources available to him [8][9][10][11] and also corresponded with mathematicians both in the U.S. and abroad for essentially similar inversion methods developed by others but was unable to find anything resembling this particular technique. He, there-fore, concluded that the engineering profession is currently unfamiliar with this procedure and could benefit from its application wherever Gauss-Jordan is used. Although much of the hardware restrictions of the past which prompted him to develop the procedure no longer exist, still there would be no logical reason to use the lengthier method when a shorter but mathematically identical alternative was available.

Conclusions
The In-Place method is a shorter equivalent of the Gauss-Jordan matrix inversion algorithm for making more efficient use of the computational resources by virtualizing the augmenting unit matrix. While the underlying mathematics is the same for both, the In-Place method requires only one-half of the computer resources needed by Gauss-Jordan and less computational time and effort because while its number of steps is the same, its operation is confined to the same number of elements as in the original matrix.
It can be used wherever Gauss-Jordan is used and it has been demonstrated to be equally useful for inverting both small and large matrices. Although the author himself used it exclusively for structural analysis, as a mathematical algorithm it is equally applicable to other engineering, scientific and mathematical tasks which involve matrix inversion since problems encountered across diverse scientific and technological disciplines often exhibit similar characteristics.
The time and resources available to the author did not permit him to further refine the method to utilize the symmetry property of structural stiffness matrices and he hopes that others may pick up where he had left off and make it even more efficient for special matrices with similar properties.
The author wishes to mention that this procedure is meant for those for whom speed and economy of computing power as well as minimization of truncation error are important.
Two areas where he believes his algorithm could also be useful are (a) large-scale real-time applications requiring high-speed inversion and/or high-precision arithmetic and (b) small-scale applications for portable devices with limited memory space.

Acknowledgements
The author gratefully acknowledges the guidance and encouragement he received from his supervisor and mentor late K. Madhavan, former Deputy Director, Central Water and Power Commission, New Delhi, India and the valuable contributions and suggestions of his ex-colleague late M. R. Rao, former Mathematician/Programmer, Computer Center, Planning Commission, New Delhi, India.