Anson W K Ma

Malcolm R Mackley

Francisco Chinesta

The microstructure and rheology of carbon nanotube suspensions

Keywords: Microstructure, Rheology, Carbon nanotube suspensions

This paper gives a brief overview of microstructure and rheology that have been observed for a range of carbon nanotube (CNT) suspensions. In general, untreated CNT suspensions show a much higher level of observable optical microstructure reflecting their preference to aggregate; they also show higher levels of viscoelasticity over treated CNT suspensions. An unexpected Helical Band texture for untreated CNTs is reported together with a series of parallel plate optical observations showing a broad spectrum of behaviour for different shear conditions. Both steady shear and linear viscoelastic data are presented for treated and untreated systems.

Introduction

The availability of carbon nanotubes (CNTs) in grams and increasingly kilogram quantities now means that optical and rheological characterisation of CNT suspensions can be carried out in a systematic way. There have already been studies on both the optical [START_REF] Hobbie | Optical measurements of structure and orientation in sheared carbon-nanotube suspensions[END_REF][START_REF] Hobbie | Orientation of carbon nanotubes in a sheared polymer melt[END_REF][START_REF] Lin-Gibson | Elastic flow instability in nanotube suspensions[END_REF][START_REF] Rahatekar | Optical microstructure and viscosity enhancement for an epoxy resin matrix containing multiwall carbon nanotubes[END_REF] and rheological response [START_REF] Pötschke | Rheological behavior of multiwalled carbon nanotube/polycarbonate composites[END_REF][START_REF] Song | Influence of dispersion states of carbon nanotubes on physical properties of epoxy nanocomposites[END_REF][START_REF] Saib | Carbon nanotube composites for broadband microwave absorbing materials[END_REF][START_REF] Fan | Rheology of multiwall carbon nanotube suspensions[END_REF][START_REF] Ma | Filament stretching of carbon nanotube suspensions[END_REF] of different systems and a scientific picture of the state of organisation is slowly emerging such that physical modelling of systems is now possible [10,11]. If the CNTs are suspended in a polymer melt, there is the complexity of rheological and optical effects from both the CNT and matrix [START_REF] Hobbie | Orientation of carbon nanotubes in a sheared polymer melt[END_REF][START_REF] Saib | Carbon nanotube composites for broadband microwave absorbing materials[END_REF]12]. In this paper, a series of experiments are reported where CNTs are suspended in a Newtonian matrix thereby simplifying the problem. In addition, in order to assist mixing and making rheological measurements, a relatively high viscosity epoxy matrix (Araldite LY556, Huntsman LLC) was chosen.

Optical microstructure of CNT suspensions

Helical band formation in untreated CNT suspensions A series of optical observations were carried out using the Cambridge Shear System (CSS450, Linkam Scientific Instruments). As shown in Fig. 1, the CSS450 system consists of a flow cell coupled with an optical microscope which allows for the in situ visualisation of CNT response towards shear. By applying appropriate steady shear to untreated CNTs suspended within a low viscosity epoxy resin, an unusual type of highly anisotropic aggregate structures named "Helical Bands (HBs)" was observed. Detailed experimental methods and results were reported in [13]. Figure 2 shows the time evolution of optical microstructures of the untreated CNT suspension with a constant shear rate of 0.5 s -1 but different gap sizes were used. The formation of HBs is best viewed in video mode, but it can be summarised as follows. In the presence of low shear, initially isotropic aggregates of CNT at different depths within the field of view translated at different speeds due to a velocity gradient across the gap. This led to collision of aggregates and subsequent formation of larger aggregates. These aggregates further grew in mass and length as they rotated in the simple shear flow and captured aggregates that came into contact. The resulting HB structures have a relatively uniform diameter and a preferential alignment in the vorticity direction as shown in Fig. 2c,f and i. In terms of the diameter of these structures, HBs obtained using 55 μm and 130 μm gap sizes showed similar filament diameter (∼60 μm) whereas for a gap size of 180 μm, HBs with an average diameter of about 120 μm were obtained. For a gap size larger than 180 μm, no HB structures were formed within the experimental shear time of 600 s [13]. The diameter of HBs showed a certain non-linear dependence on the gap size used and it is clear that wall confinement is essential to the formation of the HB structures. Figure 3 is a schematic diagram showing the mechanism by which HB structures were formed and the mechanism is analogous to the hand-spinning of cotton into fibres. Intuitively, shear flow at low rates induced aggregation-a well-known phenomenon in the context of colloidal dispersions [START_REF] Smoluchowski | Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen[END_REF], but due to wall confinement at small gap sizes, CNT aggregates can only grow as a filament rather a spherical entity. In addition to gap size, formation of HB structures was found to be sensitive to the magnitude of shear rate applied and the HB structures became unstable and broke up into isotropic aggregates at high shear [13]. In essence, the formation of HBs was believed to be a consequence of both mechanical aggregation and wall confinement. Similar banding structures have been reported for other suspensions (see review [START_REF] Vermant | Flow-induced structure in colloidal suspensions[END_REF]) and although at present there is no universal consensus to describe all the results, the physics behind the formation of HB can be of a more general origin.

Post-curing optical microstructures observed across parallel plate shear domain

Recently, Mykhaylyk and co-workers have exploited the effect of variable shear rate with radial distance from the centre within parallel plate rotary shear [START_REF] Mykhaylyk | The specific work of flow as a criterion for orientation in polymer crystallisation[END_REF]. By freezing optical microstructures within a disc of sheared material they were able to establish the effect of shear rate on certain polymer crystallisation. Using a similar idea, a series of experiments were carried out where CNT suspensions were sheared within two parallel plates of a standard straincontrolled rheometer (TA Instruments Inc.). The upper plate was made of quartz, which allows the passage of UV light whereas shear was applied via rotating the bottom plate. Photo-curing within a rheometer allowed for a precise gap size control (+/-10 μm) and a gap size of 300 μm was used in all the experimental results reported in this paper.

On the cessation of rotary shear, UV radiation (RX Firefly; Phoseon Technology) was applied as shown in Fig. 4 and this enabled cured acrylic disc to be recovered and subsequently examined using an optical microscope. The local shear rate between the parallel plates is a function of the radial position and the photo-curing method reported herein thus offers the advantage of producing a spectrum of microstructures that corresponds to a wide range of local shear rates. In the absence of CNTs, the acrylic mixture of dipropylene glycol diacrylate (DPGDA) and Additol® BCPK (Cytec Industries) photo-cured within 2 s. Addition of CNTs prolonged the time required for curing, but for low concentration samples (< 0.05% CNT), the curing process generally completed within 2 minutes. Figure 5 shows the optical microstructure of sheared and subsequently cured samples for untreated/treated single/ multi-walled CNT in acrylic. As shown in Fig. 5f-j, treated single-walled CNTs were found to be well-dispersed and steady shear did not result in the coarsening of optical texture or formation of larger aggregates. This suggests that electrostatic repulsion introduced between CNTs by specific chemical treatment [START_REF] Dyke | Overcoming the insolubility of carbon nanotubes through high degrees of sidewall functionalization[END_REF][START_REF] Dyke | Unbundled and highly functionalized carbon nanotubes from aqueous reactions[END_REF] was sufficient to prevent the aggregation of CNTs even in the presence of shear. In contrast, untreated single-walled CNT suspensions (Fig. 5a-e) exhibited a richer optical aggregate structure that depended on the magnitude of shear rate applied. Prior to the application of steady shear, untreated single-walled CNT suspensions showed an optical texture of dense and isotropic aggregates, as shown in Fig. 5a, after 5 hours of high shear mixing (Silverson L4R). Application of low shear (1 s -1 ) led to the formation of Helical Band (HB) structures that aligned perpendicular to the flow direction, resulting in a pattern where anisotropic aggregates radiated from the centre of rotation (Fig. 5b). The area within which HB structures were formed contracted towards the centre as the rim shear rate was increased from 1 to 100 s -1 (Fig. 5c-d), consistent with our earlier in situ observations in [13] that high shear led to the breakup of HBs into more isotropic aggregates. At a (rim) shear rate of 1,000 s -1 , only dense isotropic aggregates remained (Fig. 5e), resembling the microstructures obtained from high shear mixing (Fig. 5a). The setup was based on a standard ARES strain-controlled rheometer which allows easy application of different steady shear conditions with a precise gap size control (+/-10 μm). UV light was introduced immediately after cessation of shear to initiate the photo-curing process of the acrylic mixture, which typically takes less than 2 min to complete in the presence of CNT Similar evolution of aggregate structure as a function of the magnitude of shear rate applied was observed for untreated multi-walled CNTs (Fig. 5k-o). Low shear resulted in the formation of HB structures which became unstable at high shear. In the case of treated multi-walled CNT suspensions (Fig. 5p), the degree of aggregation after mixing is less than that of the untreated multi-walled CNT suspensions (Fig. 5k), but the application of steady shear was found to result in the formation of mesostructure patterns (Fig. 5q-t). As shown in Fig. 5q, CNT aggregates with circumferential alignment instead of HB type alignment were obtained at low shear. At a shear rate of 10 s -1 , the aggregates broke up forming isotropic aggregates (Fig. 5r) and clear circumferential alignment of these isotropic aggregates was observed at a shear rate of 100 s -1 (Fig. 5s). At a shear rate of 1,000 s -1 , larger isotropic aggregates were found close to the centre of rotation, together with a concentration gradient of smaller aggregates in the radial direction (Fig. 5t). The local shear rate close to the centre is comparatively lower than the outer region and this could explain the presence of larger aggregates at the centre; however, the exact physical reason for increasing concentration of small CNT aggregates with increasing radial distance is yet to be investigated.

UV Light

From these observations, it is clear that optical micro-and meso-structures formed within CNT suspensions depended on the type of CNT used as well as the magnitude of shear rate applied. Untreated CNTs tended to exhibit a richer optical texture and low shear resulted in mechanical aggregation of CNTs. A number of micro/meso-structures of CNTs were observed and although a uniform dispersion as in the case of treated single-walled CNTs (Fig. 5f-j) is generally preferred, some of the resulting meso-structures and patterns might offer strategic benefits for some CNT applications.

Experimental rheology of CNT suspensions

Steady shear viscosity

It has been reported by a number of authors (see for example [START_REF] Rahatekar | Optical microstructure and viscosity enhancement for an epoxy resin matrix containing multiwall carbon nanotubes[END_REF][START_REF] Ma | The rheology, microstructure and film processing of carbon nanotube suspensions[END_REF]) that addition of CNTs to a Newtonian suspending medium increased the apparent shear viscosity at low shear and as shear rate increased, the viscosity decreased asymptotically to the suspending medium viscosity. This type of behaviour is known as shear-thinning and is commonly observed in many suspensions [START_REF] Cross | Rheology of non-Newtonian fluids: a new flow equation for pseudo-plastic systems[END_REF]. The extent of shear-thinning was found to depend on a number Fig. 5 UV-cured optical microstructures for four different types of CNT suspensions (conc.=0.05%) after application of steady shear. Optical images were presented as a function of increasing shear from the left to the right with a shear time in excess of 100 s and gap size=300 μm of factors and treated CNT suspensions generally show a smaller shear-thinning effect compared with untreated CNT suspensions as shown in Fig. 6. Figure 6 shows apparent shear viscosity as a function of shear rate with matching optical micrographs for epoxy, treated and untreated CNT suspensions. In terms of the physical understanding, the shear-thinning effect observed in treated CNT suspensions can be attributed to the progressive alignment of CNT in flow direction as the magnitude of shear rate increases. Recent Fokker-Planck (FP) based numerical modelling [11] suggested that treated CNTs behaved essentially as short and rigid fibres that rotate and align in the shear flow and that their rheological behaviour can be described using an FP orientation model with an appropriate rotary diffusion coefficient. Untreated CNTs, however, showed a much larger shear-thinning effect coupled with a hierarchy of CNT aggregate structures (Fig. 6) and CNT orientation alone failed to explain the experimentally observed shear-thinning effect [10]. In this regard, a new model named the Aggregation/ Orientation (AO) model has been formulated and was reported in [10]. The model considered both the effects of orientation and aggregation and was capable of describing the experimental shear data of untreated CNT suspensions. The success of the model confirms the belief that the shearthinning effect of untreated CNT suspensions was contributed by both CNT orientation as well as aggregation.

Small amplitude oscillatory measurement: linear viscoelasticity A number of authors have reported Linear Viscoelasticity (LVE) effects of CNT suspensions [START_REF] Pötschke | Rheological behavior of multiwalled carbon nanotube/polycarbonate composites[END_REF][START_REF] Song | Influence of dispersion states of carbon nanotubes on physical properties of epoxy nanocomposites[END_REF][START_REF] Fan | Rheology of multiwall carbon nanotube suspensions[END_REF][START_REF] Xu | Shear and extensional rheology of carbon nanofibers suspensions[END_REF][START_REF] Huang | Dispersion rheology of carbon nanotubes in a polymer matrix[END_REF][START_REF] Hough | Viscoelasticity of single wall carbon nanotube suspensions[END_REF]. The level of elasticity in samples tested within this paper was very dependent of the type of CNT used. Untreated CNT suspensions tended to show a larger storage modulus (G') compared with treated CNT suspensions at the same concentration level (Fig. 7). In untreated CNT suspensions, G' plateau at low frequencies is commonly observed and normally this is a rheological signature of a physical network [START_REF] Fan | Rheology of multiwall carbon nanotube suspensions[END_REF]12]. The presence of a network of CNT aggregates was confirmed by optical observations and step strain experiments [10]. In a step strain experiment, a finite strain is applied to the suspension and the stress relaxation as a function of time is then followed. Figure 8 shows the stress relaxation of 0.1% concentration with different strain magnitudes. At low strain, the untreated CNT suspension essentially behaved as a viscoelastic solid with a very long relaxation time. At high strain, the stress relaxation time reduced and this strain-softening behaviour is probably an indication of the yielding of the physical network [START_REF] Amari | Stress relaxation of carbon blacklinseed oil suspensions[END_REF][START_REF] Mewis | Yielding in weakly flocculated systems[END_REF][START_REF] Winter | Soft polymeric materials near the transition from liquid to solid state[END_REF]. In the case of treated CNT suspensions, mild elasticity (G') was observed. In the literature, macroscopic fibres suspended in a Newtonian matrix were reported to demonstrate no elasticity (see review paper [START_REF] Petrie | The rheology of fibre suspensions[END_REF]), although elasticity was reported for sub-micron rod-like systems where Brownian motion could be a source of elasticity [START_REF] Larson | The Structure and Rheology of Complex Fluids[END_REF]. Understanding the viscoelasticity of treated CNT suspensions remains an important part of ongoing research and the findings will be reported in a future paper.

Conclusions

In this paper, the microstructure and rheology of both treated and untreated carbon nanotubes (CNTs) suspended in a Newtonian fluid have been summarised. In terms of microstructure, untreated CNT suspensions showed richer optical texture than suspensions with chemically treated CNTs. In situ microstructure of CNT suspensions was followed using the Cambridge Shear System. Helical Band (HB) structures were formed when low shear was applied to untreated CNTs suspended in a low viscosity matrix and within small confinement. HB structures are highly anisotropic mesostructure that align perpendicular to the flow direction and they form in an intriguing way that has similarity to tumbleweed formation in desert. Apart from in situ microstructure studies, a photo-curing protocol was developed to preserve flow-induced structures and the technique has been applied to different types of CNT suspensions. A number of "aboriginal" mesostructure patterns were obtained. Although a uniform dispersion of CNT is desirable for most applications, these flow-induced structures might offer strategic benefits for certain applications and the exact physical origin for each individual pattern is yet to be investigated.

In terms of rheology, both treated and untreated CNT suspensions exhibited steady shear-thinning characteristic. The treated CNT suspension with little optical microstructure, however, showed a less pronounced shear-thinning behaviour and this behaviour can be explained by considering treated CNTs as short rigid fibres that align in shear flows. Mild elasticity was recorded when small-amplitude oscillatory measurements were carried out on the treated CNT suspensions and the exact origin for this elasticity is unclear. In the case of the untreated CNT suspension, step strain experiments confirmed the presence of a physical CNT (aggregate) network and it is believed that this network is responsible for the large elasticity observed in small-amplitude oscillatory measurements as well as the more significant shear-thinning behaviour of untreated CNT suspensions in steady shear.

Fig. 1 Fig. 2

 12 Fig.1The Cambridge Shear System (CSS450, Linkam Scientific Instruments) used to capture in situ optical microstructures of CNT suspensions.[13]

Fig. 4

 4 Fig.4Experimental setup used in UV-cured microstructure studies. The setup was based on a standard ARES strain-controlled rheometer which allows easy application of different steady shear conditions with a precise gap size control (+/-10 μm). UV light was introduced immediately after cessation of shear to initiate the photo-curing process of the acrylic mixture, which typically takes less than 2 min to complete in the presence of CNT

Fig. 3

 3 Fig. 3 Schematic diagram of Helical Band (HB) formation: A HB nucleus rotates within the steady shear and captures initially isotropic aggregates of nanotubes. The nanotubes are then wound helically to form a nanotube HB with an axis perpendicular to the direction of flow. [13]

Fig. 6

 6 Fig. 6 Apparent steady shear viscosity (h a ) of (a) epoxy, (b) 0.3% treated CNT and (c) 0.3% untreated CNT with matching quiescent optical micrographs (optical depth=130 μm)

Fig. 8

 8 Fig.8Stress relaxation of 0.1% untreated CNT-epoxy suspension after step strain with different magnitudes (Gap size=0.5 mm; Temp.= 25°C; 50mm parallel plates; ARES-strain-controlled rheometer).[START_REF] Ma | The rheology, microstructure and film processing of carbon nanotube suspensions[END_REF]

Fig. 7

 7 Fig.7Linear Viscoelasticity (LVE) data (G': storage modulus, G": loss modulus and h à : complex viscosity) of (a) 0.5% untreated CNT in epoxy and (b) 0.5% treated CNT in epoxy. (Strain=1%; Gap size= 0.5 mm; Temp.=25°C; 50mm parallel plates; ARES-strain-controlled rheometer)

Upper Plate Stainless Steel Bottom Plate CNT in photo-curable acrylic UV Light Quartz Upper Plate Stainless Steel Bottom Plate CNT in photo-curable acrylic


ItamarArel Ken Habgood

Itamar Arel

A condensation-based application of Cramer's rule for solving large-scale linear systems

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Systems of linear equations are central to many science and engineering application domains. Fast linear solvers generally use a form of Gaussian elimination [START_REF] Galoppo | LU-GPU: Efficient algorithms for solving dense linear systems on graphics hardware[END_REF], the most common of which is LU-decomposition. The latter involves a computation complexity of W LU ≈ 2 3 N 3 [START_REF] Karniadakis | Parallel Scientific Computing in C ++ and MPI: A Seamless Approach to Parallel Algorithms and Their Implementation[END_REF], where N denotes the number of linearly independent columns in a matrix. The factor 2 accounts for one addition and one multiplication. If only multiplications are considered, then W LU ≈ N 3 3 , which is the operation count sometimes quoted in the literature.

In most implementations of Gaussian elimination, the row with the largest lead value and the first row are interchanged during each reduction. This is referred to as partial pivoting and facilitates the minimization of truncation errors. Historically, Cramer's rule has been considered inaccurate when compared to these methods. As this paper will discuss, the perceived inaccuracy does not originate from Cramer's rule but rather from the method utilized for obtaining determinants.

This paper revisits Cramer's rule [START_REF] Dianant | Advanced Linear Algebra for Engineers with Matlab[END_REF] and introduces an alternative framework to the traditional LU-decomposition methods offering similar computational complexity and storage requirements. To the best of the authors' knowledge, this is the first work to demonstrate such characterization. In utilizing a technique similar to partial pivoting to calculate determinant values, the algorithm's stability properties are derived and shown to be comparable to LU-decomposition for asymmetric systems.

The rest of this paper is structured as follows. Section 2 describes the algorithm and presents its computational complexity. Section 3 discusses the stability of the algorithm including forward and backward error analysis. Section 4 presents empirical results to support the proposed complexity assertions. Finally, a brief discussion and drawn conclusions are provided in Section 5.

1

Revisiting Cramer's rule

The proposed algorithm centers on the mathematically elegant Cramer's rule, which states that the components of the solution to a linear system in the form Ax = b (where A is invertible) are given by x i = det A i (b) / det( A), [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] where A i (b) denotes the matrix A with its ith column replaced by b [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. Unfortunately, when working with large matrices, Cramer's rule is generally considered impractical. This stems from the fact that the determinant values are calculated via minors. As the number of variables increases, the determinant computation becomes unwieldy [START_REF] Chapra | Numerical Methods for Engineers[END_REF]. The time complexity is widely quoted as O(N!), which would make it useless for any practical application when compared to a method like LU-decomposition at O(N 3 ). Fortunately, Cramer's rule can be realized in far lower complexity than O(N!). The complexity of Cramer's rule depends exclusively on the determinant calculations. If the determinants are calculated via minors, that complexity holds. In an effort to overcome this limitation, matrix condensation techniques can reduce the size of the original matrix to one that may be solved efficiently and quickly. As a result, Cramer's rule becomes O(N 3 ) similar to LU-decomposition.

The other concern with Cramer's rule pertains to the numerical instability, which is less studied [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. A simple example put forward in [START_REF] Moler | Cramer's rule on 2-by-2 systems[END_REF] suggests that Cramer's rule is unsatisfactory even for 2-by-2 systems, mainly because of round error difficulties. However, that argument heavily depends on the method for obtaining the determinants. If an accurate method for evaluating determinants is used then Cramer's rule can, in fact, be numerically stable. In fact, a later paper [START_REF] Dunham | Cramer's rule reconsidered or equilibration desirable[END_REF] revisited the cited example and provided an example where Cramer's rule yielded a highly accurate answer while Gaussian elimination with pivoting a poor one. We thus provide a detailed treatment of the stability and accuracy aspects of the proposed determinant calculations employed by the proposed algorithm.

Chio's condensation

Chio's condensation [START_REF] Eves | Elementary Matrix Theory[END_REF] method reduces a matrix of order N to order N -1 when evaluating its determinant. As will be shown, repeating the procedure numerous times can reduce a large matrix to a size convenient for the application of Cramer's rule. Chio's pivotal condensation theorem is described as follows. Let A = [a ij ] be an N × N matrix for which a 11 = 0. Let D denote the matrix obtained by replacing elements a ij not in the lead row or lead column by a 11 a 1 j a i1 a ij , then it can be shown that |A| = |D| a n-2 11 [START_REF] Eves | Elementary Matrix Theory[END_REF].

Note that this process replaces each element in the original matrix with a 2 × 2 determinant consisting of the a 11 element, the top value in the element's column, the first value in the element's row and the element being replaced. The calculated value of this 2 × 2 determinant replaces the initial a i, j with a i, j . The first column and first row are discarded, thereby reducing the original N × N matrix to an (N -1) × (N -1) matrix with an equivalent determinant. As an example, we consider the following 3 × Obtaining each 2 × 2 determinant requires two multiplications and one subtraction. However, if the value of a 1,1 is one, then only a single multiplication is required. In the example above we note that a 1,1 is used in each element as a multiplier to the matrix element, for example, the equation for the matrix element in position (2, 2) is a 11 a 22a 21 a 12 . If in this situation a 11 = 1, then the equation changes to a 22a 21 a 12 . This holds true for every element in the matrix. Therefore for each condensation step k, if a kk = 1 then (Nk) 2 multiplications are removed.

In order to guarantee a kk = 1, an entire row or column must be divided by a kk . This value would need to be stored because the determinant value calculated by Chio's condensation would be reduced by this factor. To find the true value at the end of the condensation, the calculated answer would need to be multiplied by each a kk that was factored out. Multiplying all of these values over numerous condensation steps would result in an extremely large number that would exceed the floating point range of most computers. This is where the elegance of Cramer's rule is exploited. Cramer's rule determines each variable by a ratio of determinants,

x i = det( A i (b))/ det( A).
Given that both determinants are from the same condensation line, they both are reduced by the same a kk values. The a kk values factored out during Chio's condensation cancel during the application of Cramer's rule. This allows the algorithm to simply discard the a kk values in the final computations. The actual determinant values are not correct, however the ratio evaluated at the core of Cramer's rule remains correct.

The cost of using Chio's condensation is equivalent to computing (Nk) 2 2 × 2 determinants and (Nk) divisions to create a kk = 1. Hence, the computational effort required to reduce an N × N matrix to a 1

× 1 matrix is O(N 3 ), since N-1 k=1 2(N -k) 2 + (N -k) = 2N 3 3 - N 2 2 - N 6 -∼ O N 3 . (2)
Combined with Cramer's rule, this process can yield the determinant to find a single variable. In order to find all N variables, this would need to be repeated N times, suggesting that the resulting work would amount to O(N 4 ). We next describe an efficient way of overcoming the need for O(N 4 ) computations by sharing information pertaining to each variable.

Matrix mirroring

The overarching goal of the proposed approach is to obtain an algorithm with O(N 3 ) complexity and low storage requirement overhead. As discussed above, Chio's condensation and Cramer's rule provide an elegant solution with O(N 4 ) complexity. In order to retain O(N 3 ) computational complexity, it is necessary to reuse some of the intermediate calculations performed by prior condensation steps. This is achieved by constructing a binary, tree-based data flow in which the algorithm mirrors the matrix at critical points during the condensation process, as detailed next.

A matrix A and a vector of constants b are passed as arguments to the algorithm. The latter begins by appending b to A creating an augmented matrix. All calculations performed on this matrix are also performed on b. Normal utilization of Cramer's rule would involve substitution of the column corresponding to a variable with the vector b, however the proposed algorithm introduces a delay in such substitution such that multiple variables can be solved utilizing one line of Chio's condensation. In order to delay the column replacement, b must be subject to the same condensation manipulations that would occur had it already been in place. This serves as the motivation for appending b to the matrix during condensation.

The condensation method removes information associated with discarded columns, which suggests that the variables associated with those columns cannot be computed once condensed. For this reason, a mirror of the matrix is created each time the matrix size is halved. The mirrored matrix is identical to the original except the order of its columns is reversed. For example, the first and last column are swapped, the second and second to last column are swapped, and so on. A simple 3 × 3 matrix mirroring operation would be: a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 → mirrored → a 13 a 12 -a 11 a 23 a 22 -a 21 a 33 a 32 -a 31 .

In the example above, the third column of the mirror is negated, which is performed in order to retain the correct value of the determinant. Any exchange of columns or rows requires the negation of one to preserve the correct determinant value. As discussed in more detail below, this negation is not necessary to arrive at the correct answer, but is applied for consistency.

Following the mirroring, each matrix is assigned half of the variables. The original matrix can solve for the latter half while the mirrored matrix solves for the first half of the variables. In the example above, there are three variables: x 1 , x 2 , x 3 . The original matrix could solve for x 2 , x 3 , and the mirrored matrix would provide x 1 . Each matrix uses condensation to yield a reduced matrix with size at least equal to the number of variables it's responsible for. For the case of a 3 × 3 matrix, we have the pair of matrices:

× × × × a 22 a 23 × a 32 a 33 original: x 2 ,x 3 × × × × a 22 a 21 × a 32 a 31 . mirrored: x 1
Once this stage is reduced, the algorithm either solves for the variables using Cramer's rule or mirrors the matrix and continues with further condensation. A process flow depicting the proposed framework is illustrated in Fig. 1.

Extended condensation

Chio's condensation reduces the matrix by one order per repetition. Such an operation is referred to here as a condensation step of size one. It's possible to reduce the matrix by more than one order during each step. Carrying out the condensation one stage further, with leading pivot a 11 a 12 a 21 a 22 (assumed to be non-zero), we have [START_REF] Aitken | Determinants and Matrices[END_REF] Fig. 1. Tree architecture applied to solving a linear system using the proposed algorithm.

× × × × × × × × × × × a 33 a 34 × × a 43 a 44 . ( 4
)
In this case, each of the matrix elements {a 33 , a 34 , a 43 , a 44 } is replaced by a 3 × 3 determinant instead of a 2 × 2 determinant. This delivers a drastic reduction in the number of repetitions needed to condense a matrix. Moreover, a portion of each minor is repeated, namely the a 11 a 12 a 21 a 22 component. Such calculation can be performed once and then reused multiple times during the condensation process. Letting M denote the size of the condensation, and using the example above, M = 2 while for the basic Chio's condensation technique M = 1. As an example, a 6 × 6 matrix could be reduced to a 2 × 2 matrix in only two condensation steps, whereas it would take four traversals of the matrix to arrive at a 2 × 2 matrix if M = 1. The trade-off is a larger determinant to calculate for each condensed matrix element. Instead of having 2 × 2 determinants to calculate, 3 × 3 determinants are needed, suggesting a net gain of zero.

However, the advantage of this formulation is that larger determinants in the same row or column share a larger number of the minors. The determinant minors can be calculated at the beginning of each row and then reused for every element in that row. This reduces the number of operations required for each element with a small penalty at the outset of each row. In practical computer implementation, this also involves fewer memory access operations, thus resulting in higher overall execution speed.

Pivoting in the case of M > 1 requires identifying a lead determinant that is not small. As with pivoting for LUdecomposition, ideally the largest possible lead determinant would be moved into the top left portion of the matrix. Unfortunately, this severely compromises the computational complexity, since an exhaustive search for the largest lead determinant is impractical. Instead, a heuristic method should be employed to select a relatively large lead determinant when compared to the alternatives.

Algorithm 1 Extended condensation.

{A[][] = current matrix, N = size of current matrix} {mirrorsize = unknowns for this matrix to solve for} while

(N -M) > mirrorsize do lead_determinant = CalculateMinor(A[][], M + 1, M + 1) if lead_determinant = 0 then return(error) end if A[1:N][1] = A[1:N][1] / lead_determinant;
[i] = leadminor[i] + (-1) j-1 A[row][j] × reusableminor[i][j]
end for end for {Core Loop; find the M × M determinant for each A[][] item} for col = (M + 1) to (N + 1) step by 1 do for i = 1 to M step by 1 do {calculate MxM determinant} The pseudocode in Algorithm 1 details a basic implementation of the proposed algorithm using extended condensation. The number of rows and columns condensed during each pass of the algorithm is represented by the variable M, referred to earlier as the condensation step size. The original matrix is passed to the algorithm along with the right-hand side vector in A, which is an N × (N + 1) matrix. The mirrorsize variable represents the number of variables a particular matrix solves for. In other words, it reflects the smallest size that a matrix can be condensed to before it must be mirrored. The original matrix, A, has a mirrorsize = N since it solves for all N variables. It should be noted that this will result in bypassing the while loop completely at the initial call of the algorithm. Prior to performing any condensations, the original matrix is mirrored. The first mirror will have mirrorsize = N 2 , since it only has to solve for half of the variables. After this mirror has been created, the algorithm will begin the condensation identified within the while loop. , since an M + 1 value excludes a column beyond the M rows and M columns used to calculate the determinant. This is used in the algorithm to find the lead determinant at each condensation step. When M = 1, the method simply returns the value at A [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF].

A[row][col] = A[row][col] + (-1) j leadminor[i] × A[i][col]
The CramersRule method solves for the variables associated with that particular matrix using Cramer's rule. The method replaces a particular column with the condensed solution vector, b, finds the determinant, and divides it by the determinant of the condensed matrix to find each variable. The method then replaces the other columns until all variables for that particular phase are calculated. The Mirror function creates a mirror of the given matrix as described in Section 2.2.

The arrays labeled reusable_minor and lead_minor maintain the pre-calculated values discussed earlier. The array reusable_minor is populated once per condensation step and holds M 2 minors that will be used at the outset of each row. The latter will then populate the lead_minor array from those reusable minors. The lead_minor array holds the M minors needed for each matrix element in a row. Hence, the lead_minor array will be repopulated N -M times per condensation step.

The algorithm divides the entire first column by the top left M × M determinant value. This is performed for the same reason that the first column of the matrix was divided by the a 11 element in the original Chio's condensation formulation. Dividing the first row by that determinant value causes the lead determinant calculation to retain a value of one during the condensation. This results in a '1' being multiplied by the a i, j element at the beginning of every calculation, thus saving one multiplication operation per each element.

Computation complexity

As illustrated in the pseudocode, the core loop of the algorithm involves the calculation of the M × M determinants for each element of the matrix during condensation. Within the algorithm, each M × M determinant requires M multiplications and M additions/subtractions. Normally, this would necessitate the standard computational workload to calculate a determinant, i.e. 2 3 M 3 , using a method such as Gaussian elimination. However, the reuse of the determinant minors described earlier reduces the effort to 2M operations within the core loop.

An additional workload outside the core loop is required since M 2 minors must be pre-calculated before Chio's condensation commences. Assuming the same 2 3 M 3 workload using Gaussian elimination to find a determinant and repetition of this at each of the N M condensation steps, yields an overhead of

N M × M 2 × 2 3 M 3 = 2M 4 N 3 .
(

) 5
In situations where M N, this effort is insignificant, although with larger values of M relative to N, it becomes nonnegligible.

The optimal size of M occurs when there's a balance between pre-calculation done outside the core loop and the savings during each iteration. In order to reuse intermediate calculation results, a number of determinant minors must be evaluated in advance of each condensation step. These reusable minors save time during the core loops of the algorithm, but do not utilize the most efficient method. If a large number of minors are required prior to each condensation, their additional computation annuls the savings obtained within the core iterations.

The optimal size of M is thus calculated by setting the derivative of the full complexity formula, 2

3 N 3 -M N 2 + N 2 M + 2 3 M 4 N + M 2 N
, with respect to M to zero. This reduces to 8 3 M 5 + 2M 3 = N M 2 + N, suggesting an optimal value of 3 3 8 N

for M. As an example, the optimal point for a 1000 × 1000 matrix is ≈ 7.22. Empirical results indicate that the shortest execution time for a 1000 × 1000 matrix was achieved when M = 8, supporting the theoretical result.

In order to condense a matrix from N × N to (N -M) × (N -M), the core calculation is repeated (N -M) 2 times. The algorithm requires N/M condensation steps to reduce the matrix completely and solve using Cramer's rule. In terms of operations, this equates to

γ = N/M k=1 2M(N -kM) 2 = 2M N/M k=1 N 2 -2N Mk + M 2 k 2 = 2M N M N 2 -2N M N M ( N M + 1) 2 + M 2 N M ( N M + 1)( 2N M + 1) 6 = 2 3 N 3 -M N 2 + M 2 N 3 , ( 6
)
resulting in a computational complexity, γ , of 2 3 N 3 to obtain a single variable solution.

Mirroring occurs with the initial matrix and then each time a matrix is reduced in half. An N × N matrix is mirrored when it reaches the size of

N 2 × N 2 .
Once the matrix is mirrored, there is double the work. In other words, two N 2 × N 2 matrices each require a condensation, where previously there was only one. However, the amount of work for two matrices of half the size is much lower than that of one N × N matrix, which avoids the O (N 4 ) growth pattern in computations. This is due to the O (N 3 ) nature of the condensation process.

Since mirroring occurs each time the matrix is reduced to half, log 2 N matrices remain when the algorithm concludes. The work associated with each of these mirrored matrices needs to be included in the overall computation load estimate. The addition of the mirrors follows a geometric series resulting in roughly 2.5 times the original workload, which leads to a computational complexity of 5 3 N 3 when ignoring the lower order terms.

The full computational complexity is the combination of the work involved in reducing the original matrix, γ , and that of reducing the mirrors generated by the algorithm. Hence, the total complexity can be expressed as follows:

γ + log 2 N k=0 2 k 2 3 N 2 k 3 -M N 2 k 2 + M 2 3 N 2 k . ( 7
)
The latter summation can be simplified using the geometric series equivalence n-1 k=0 ar k = a 1-r n 1-r [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF], which when ignoring the lower order terms reduces to:

γ + 8 9 N 3 = 14N 3 9 ≈ 5 3 N 3 . (8)

Mirroring considerations and related memory requirements

Eq. ( 8) expresses the computational complexity assuming a split into two matrices each time the algorithm performs mirroring. One may consider a scenario in which the algorithm creates more than two matrices during each mirroring step. In the general case, the computational complexity is given by

γ + log S N k=0 (S -1)S k 2 3 N S k 3 -M N S k 2 + M 2 3 N S k ( 9
)
where S denotes the number of splits. When lower order terms are ignored, this yields

γ + 2S 2 3(S + 1) N 3 . ( 10
)
As S increases the complexity clearly grows. The optimal number of splits is thus two, since that represents the smallest value of S that can still solve for all variables. Additional splits could facilitate more work in parallel, however they would generate significantly greater overall workload.

The memory requirement of the algorithm is 2 × (N + 1) 2 , reflecting the need for sufficient space for the original matrix and the first mirror. The rest of the algorithm can reuse that memory space. Since the memory requirement is double the amount required by typical LU-decomposition implementations and similar to LU-decomposition, the original matrix is overwritten during calculations.

Algorithm stability

The stability properties of the proposed algorithm are very similar to those of Gaussian Elimination techniques. Both schemes are mathematically accurate yet subject to truncation and rounding errors. As with LU-decomposition, if these errors are not accounted for, the algorithm returns poor accuracy. LU-decomposition utilizes partial or complete pivoting to minimize truncation errors. As will be shown, the proposed algorithm employs a similar technique.

Each element during a condensation step is affected by the lead determinant and the 'lead minors' discussed in the earlier section. In order to avoid truncation errors, these values should go from largest on the first condensation step to smallest on the last condensation step. This avoids a situation where matrix values are drastically reduced, causing truncation, and then significantly enlarged later, magnifying the truncation error. The easiest method to avoid this problem is by moving the rows that would generate the largest determinant value to the lead rows before each condensation step. This ensures the largest determinant values available are used.

Once the matrix is rearranged with the largest determinant in the lead, normalization occurs. Each lead value is divided by the determinant value, resulting in the lead determinant equaling unity. This not only reduces the number of floating point calculations but serves to normalize the matrix.

Backward error analysis

The backward stability analysis presented yields results similar to LU-decomposition. This coupled with empirical findings provides evidence that the algorithm yields accuracy comparable to that of LU-decomposition. As with any computer calculations, rounding errors affect the accuracy of the algorithm's solution. Backward stability analysis shows that the solution provided is the exact solution to a slightly perturbed problem. The typical notation for this concept is

(A + F )x = b + δb. ( 11
)
Here A denotes the original matrix, b gives the constant values, and x gives the solution calculated using the algorithm. F represents the adjustments to A, and δb the adjustment to b that provides a problem that would result in the calculated solution using exact arithmetic. In this analysis the simplest case is given, namely where the algorithm uses M = 1.

In the first stage of condensation, A (2) is computed from A (1) , which is the original matrix. It should be noted that each condensation step also incurs error on the right-hand side due to the algorithm carrying those values along during reduction. This error must also be accounted for in each step, so in the first stage of condensation, b (2) is computed from b (1) just as A. Before Chio's pivotal condensation occurs, the largest determinant is moved into the lead position. Since M = 1, the determinant is simply the value of the element. This greatly simplifies the conceptual nature for conveying this analysis.

Normally M × M determinants would need to be calculated, and then all the rows comprising the largest determinant moved into the lead. Here, the row with the largest lead value a i1 is moved to row one, followed by each element in column one being divided by a kk . This normalizes the matrix so that the absolute values of all row leads are smaller or equal to one, such that a kj = a kj a kk (1 + η kj ), where η kj β -t+1 . [START_REF] Moler | Cramer's rule on 2-by-2 systems[END_REF] In this case, β -t+1 is the base number system used for calculation with t digits. This is equivalent to machine epsilon, .

The computed elements a [START_REF] Aitken | Determinants and Matrices[END_REF] ij are derived from this basic equation a

(2) ij = a ij -a ik ×a kj .
When the error analysis is then added, we get: (1) ij , γ (1) ij β -t+1 and α (1)

a (2) ij = a (1) ij -a ik × a kj × 1 + γ (1) ij 1 + α
ij β -t+1 , (13) a (2) ij = a (1) ij -a ik × a kj + e (1) ij (14) where e (1) ij = a (2) ij × α (1) ij (1 + α (1) ij )
-a ik × a kj × γ (1) ij , ij = 2, . . . ,n.

(15)

This then provides the elements for E (1) such that A + E (1) provides the matrix that would condense to A (2) with exact arithmetic. The lead column is given by e (1) ij = a ik × η kj . This follows for each condensation step

A (k+1) = A (k) + E (k) and similarly for the right-hand side, b (k+1) = b (k) + E (k) b
, where E includes an additional column to capture the error incurred on b. In this case, the E matrix will capture the variability represented by δb found in Eq. [START_REF] Karniadakis | Parallel Scientific Computing in C ++ and MPI: A Seamless Approach to Parallel Algorithms and Their Implementation[END_REF]. If taken through all steps of condensation, then E = E (1)

+ • • • + E (n-1) , giving (A + E)x = b. ( 16
)
Bounds on E need evaluation, since this controls how different the matrix used for computation is from the original matrix. It's important to note that, due to the use of Cramer's rule, the algorithm can normalize a matrix and simply discard the value used to normalize. Cramer's rule is a ratio of values so as long as both values are divided by the same number the magnitude of that number is unimportant. This is a crucial trait, since needing to retain and later use these values would cause instability.

Consider a = max |a ij |, g = 1 a max |a (k)
ij | and Eq. ( 15). If these are combined along with the knowledge that |a

1 j | 1, the following is obtained e (k) ij β -t+1 1 -β -t+1 a (k+1) ij + β -t+1 × a (k) ij 2 1 -β -t+1 agβ -t+1 . ( 17
)
In essence the E matrix yields the following

|E| ag Υ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 . . . . . . 0 β -t+1 2 . . . . . . 2 : : : : : : β -t+1 2 . . . . . . 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 . . . 0 0 0 0 . . . 0 0 β -t+1 2 . . . 2 : : : : 0 β -t+1 2 . . . 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + • • • ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ , ( 18
)
where Υ = β -t+1

(1-β -t+1 ) , such that

|E| ag Υ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 . . . . . . 0 β -t+1 2 . . . . . . 2 : 2 + β -t+1 4
. . . 4 : : : :

β -t+1 2 + β -t+1 4 + β -t+1 . . . 2n ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 19
)
The bottom row of the matrix clearly provides the largest possible value, whereby the summation is roughly n 2 . When combined with the other factors, it yields the equality E ∞ 2n 2 ag

β -t+1
1-β -t+1 . If it's assumed that 1β -t+1 ≈ 1 and a is simply a scaling factor of the matrix, two values of interest are left: n 2 and the growth factor g. The growth factor is the element that has the greatest impact on the overall value, since it provides a measure of the increase in value

) ij = max a (k) ij -a ik × a kj 2 × max a (k) ij . ( 20
)
The value of a (k) ij is at most the largest value in the matrix. The value of a ik × a kj is also at most the largest value in the matrix. Since a ik is the row lead, it's guaranteed to be one or less. The value of a kj could possibly be the largest value in the matrix. Therefore, the greatest value that could result from the equation a (k) ija ik × a kj is twice the maximum value in the matrix or 2 max |a (k) ij |. This can then repeat at most n times, which results in a growth factor of 2 n . The growth factor given for LU-decomposition with partial pivoting in the literature is g 2 n-1 [START_REF] Ortega | Numerical Analysis: A Second Course[END_REF]. The slight difference being that this algorithm computes a solution directly, whereas LU-decomposition analysis must still employ forward and backward substitution to compute a solution vector. As with LU-decomposition, it can be seen that generally the growth rate will less than double at each step. In fact, the values tend to cancel each other leaving the growth rate around 1 in actual usage.

Mirroring does not affect the stability analysis of the algorithm. The matrices that are used to calculate the answers may have been mirrored numerous times. Since no calculations take place during the mirroring, and it does not introduce an additional condensation step, the mirroring has no bearing on the accuracy.

Backward error measurements

One of the most beneficial traits of LU-decomposition is that, although it has a growth rate of 2 n-1 , in practice it generally remains stable. This can be demonstrated by relating the relative residual to the relative change in the matrix, giving the following inequality:

b -A x A • x E A . ( 21
)
The symbol x represents the norm of the calculated solution vector. When the residual found from this solution set is divided by the norm of the original matrix multiplied by the norm of the solution set, an estimate is produced of how close the solved problem is to the original problem. If an algorithm produces a solution to a problem that is very close to the original problem then the algorithm is considered stable. A reasonable expectation for how close the solved and given problems should be is expressed as [START_REF] Heath | Scientific Computing: An Introductory Survey[END_REF].

E A ≈ n machine . ( 22
)
A pragmatic value of machine ≈ 2.2 × 10 -16 reflects the smallest value the hardware can accurately support, and n represents the size of the linear system. Table 1 shows this relative residual calculations when using Cramer's rule in comparison to those obtained with Matlab and for the target values given by Eq. ( 22). The infinite norm is used for all norm calculations and Cramer's rule used a condensation step size (M) of 8. As shown in Table 1, both Matlab's implementation of LU-decomposition and Cramer's rule deliver results below the target level to suggest a stable algorithm for the test matrices considered. The latter were created by populating the matrices with random values between -5 and 5 using the standard C language random number generator. Results produced by the proposed algorithm for these matrices were measured over numerous trials.

Matlab includes a number of test matrices in a matrix gallery that were used for further comparison. In particular, a set of dense matrices from this gallery were selected. Each type had four samples of 1000 × 1000 matrices. In many cases, Cramer's rule resorted to a condensation size of one (M = 1) for improved accuracy. The relative residuals were then calculated in the manner discussed above. Table 2 provides a summary of those findings. The results show relative residuals near Matlab for a majority of the specialized matrices.

Forward measurements

The forward error is typically defined as the relative difference between the true values and the calculated ones. Here, the actual answers are generated by Matlab and GSL (GNU scientific library). The solution vector provided by the algorithm was compared to those given by both software packages. Table 3 details the observed relative difference between the software packages and the solutions provided by the proposed algorithm.

Implementation results

The algorithm has been implemented on a single processor platform. Compiled in the C programming environment, it has been compared to LAPACK (Linear Algebra PACKage) via Matlab on a single core Intel Pentium machine, to provide a baseline for comparison. The processor was an Intel Pentium M Banias with a frequency of 1.5 GHz using a 32 KB L1 data cache and 1 MB L2 cache. The manufacturer quotes a maximum GFLOPS rate of 2.25 for the specific processor [START_REF]Intel_Corporation, Intel microprocessor export compliance metrics[END_REF]. The Linpack benchmark MFLOPS for this processor is given as 755.35.

Several optimization techniques have been applied to the implementation, including the SIMD (single instruction multiple data) parallelism that is standard on most modern processors. The program code also employs memory optimizations such as cache blocking to reduce misses. No multi-processor parallelization has been programmed into the implementation such that the algorithm itself could be evaluated against the industry standard prior to any parallelization efforts. The implementation has checks to ensure a certain levels of accuracy. In certain matrices, as mentioned in Section 3.2, the algorithm may drop to a lower execution speed to ensure accuracy of the results.

As can be seen in Table 4, the algorithm runs approximately 2.5 times slower than the execution time of Matlab, independent of matrix size, which closely corresponds to the theoretical complexity analysis presented above. Both pieces of software processed numerous trials on a 1.5 GHz single core processor. The results further show that while the algorithm is slower than the LU-based technique, it is consistent. Even as the matrix sizes grow, the algorithm remains roughly 2.5 times slower than state of the art methodologies. Fig. 2 depicts a comparison between the proposed algorithm and Matlab.

Conclusions

To the best of our knowledge, this is the first paper outlining how Cramer's rule can be applied in a scalable manner. It introduces a novel methodology for solving large-scale linear systems. Unique utilization of Cramer's rule and matrix condensation techniques yield an elegant process that has promise for parallel computing architectures. Implementation results support the theoretical claims that the accuracy and computational complexity of the proposed algorithm are similar to LU-decomposition.

  {calculate the minors that are common} for i = 1 to M step by 1 do for j = 1 to M step by 1 do {each minor will exclude one row & col} reusableminor[i][j] = CalculateMinor(A[][], i, j); end for end for for row = (M + 1) to (N + 1) step by 1 do {find the lead minors for this row} for i = 1 to M step by 1 do Set leadminor[i] = 0; for j = 1 to M step by 1 do leadminor

  {Reduce matrix size by condensation step size} N = N -M; end while if N has reached Cramer's rule size (typically 4) then {solve for the subset of unknowns assigned} x[] = CramersRule(A[][]); else {recursive call to continue condensation} A_mirror[][] = Mirror(A[][]) Recursively call Algorithm (A_mirror[][], N, mirrorsize/2) Recursively call Algorithm (A[][], N, mirrorsize/2) end if

  Three external functions assist the pseudocode: CalculateMinor, CramersRule and Mirror. CalculateMinor finds an M × M determinant from the matrix passed as an argument. The two additional arguments passed identify the row and column that should be excluded from the determinant calculation. For example, CalculateMinor( A[][], 2, 3) would find the M × M determinant from the top left portion of matrix A[][], whereby row 2 and column 3 are excluded. The method would return the top left M × M determinant of A[][] without excluding any rows or columns by calling CalculateMinor( A[][], M + 1, M + 1)

Fig. 2 .

 2 Fig.2. Algorithm execution times compared to those obtained using Matlab(TM) .

  11 a 22a 21 a 12 ) (a 11 a 23a 21 a 13 ) 0 (a 11 a 32a 31 a 12 ) (a 11 a 33a 31 a 13 )

			3 matrix:
		a 11 a 12 a 13
	A =	a 21 a 22 a 23	and its condensed form:
		a 31 a 32 a 33
	0	0	0	× × ×
	0 (a =	× a 22 a 23
				× a 32 a 33

a

  a 11 a 12 a 13 a 14 a 21 a 22 a 23 a 24 a 31 a 32 a 33 a 34 a 41 a 42 a 43 a 44 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 a 11 a 12 a 14 a 21 a 22 a 24 a 31 a 32 a 34 a 11 a 12 a 13 a 21 a 22 a 23 a 41 a 42 a 43 a 11 a 12 a 14 a 21 a 22 a 24 a 41 a 42 a 44 a 11 a 12 a 13 a 14 a 21 a 22 a 23 a 24 a 31 a 32 a 33 a 34 a 41 a 42 a 43 a 44

	=	a 11 a 12 a 21 a 22	-1	×	,	(3)
	=	a 21 a 22 a 11 a 12	-1

Table 1

 1 Relative residual measurements for Cramer's algorithm.

	Matrix size	n machine	Cramer's	b-A x A • x	Matlab	b-A x A • x
	1000 × 1000	2.22E-13	5.93E-14		8.05E-16
	2000 × 2000	4.44E-13	5.42E-14		1.04E-15
	3000 × 3000	6.66E-13	9.62E-14		1.95E-15
	4000 × 4000	8.88E-13	3.32E-13		2.49E-15
	5000 × 5000	1.11E-12	8.12E-14		3.05E-15
	6000 × 6000	1.33E-12	5.52E-14		3.35E-15
	7000 × 7000	1.55E-12	7.46E-14		3.55E-15
	8000 × 8000	1.78E-12	8.12E-14		4.28E-15
	over numerous condensation steps. Fortunately, this value is bound because all multipliers are due to the pivoting and the
	division performed before each condensation, such that
	max a	(k+1

Table 2

 2 Comparisons using Matlab matrix gallery.

	Matrix type	Cramer's	b-A x A • x	Matlab	b-A x A • x
	chebspec -Chebyshev spectral differentiation matrix	1.95E-07		1.13E-16
	clement -Tridiagonal matrix with zero diagonal entries	3.66E-16		2.80E-17
	lehmer -Symmetric positive definite matrix	2.67E-15		6.46E-18
	circul -Circulant matrix	6.39E-14		8.35E-16
	lesp -Tridiagonal matrix with real, sensitive eigenvalues	1.22E-16		1.43E-18
	minij -Symmetric positive definite matrix	2.44E-15		2.99E-18
	orthog -Orthogonal and nearly orthogonal matrices	1.41E-17		6.56E-17
	randjorth -Random J-orthogonal matrix	1.40E-09		7.24E-16
	frank -Matrix with ill-conditioned eigenvalues	5.59E-03		1.52E-21

Table 3

 3 Algorithm relative error when compared to Matlab and GSL solution sets.

	Matrix size	κ(A)	Matlab x -x ∞	GSL x -x ∞	Avg Matlab	Avg GSL
	1000 × 1000	506 930	2.39E-9	1.93E-10	1.03E-10	5.38E-12
	2000 × 2000	790 345	4.52E-9	5.36E-9	1 . 0 1 E -10	7.27E-12
	3000 × 3000	1 540 152	1.95E-8	1.84E-8	1.12E-10	2.09E-11
	4000 × 4000	12 760 599	4.81E-8	5.62E-8	1.43E-10	7.91E-11
	5000 × 5000	765 786	2.92E-8	4.39E-8	1.18E-10	3.46E-11
	6000 × 6000	1 499 430	8.67E-8	8 . 7 0 E -8	1.37E-10	6.04E-11
	7000 × 7000	3 488 010	9.92E-8	8.95E-8	1.27E-10	5.15E-11
	8000 × 8000	8 154 020	9.09E-8	9.43E-8	1.86E-10	7.85E-11

Table 4

 4 Execution time comparison to LAPACK.

	Matrix size	Algorithm (s)	Matlab (s)	Ratio
	1000 × 1000	2.06	.91	2.26
	2000 × 2000	16.44	6.32	2.60
	3000 × 3000	52.33	19.92	2.63
	4000 × 4000	115.44	45.10	2.56
	5000 × 5000	220.32	86.90	2.54
	6000 × 6000	380.92	142.05	2.68
	7000 × 7000	583.02	242.61	2.40
	8000 × 8000	872.26	334.68	2.61

Acknowledgements

The authors would like to thank Prof. Jack Dongarra and the Innovative Computing Laboratory (ICL) at the University of Tennessee for their insightful suggestions and feedback.