Embedded Systems Design - August 2008 - (Page 14)

programmer’s toolbox Trust me, you Mathematically, matrices I have seen homebrewed matrix inversion do not want to see are not commutative. the formulas for You can’t just switch the software that did no pivoting at all, both a order of multiplication, terrible idea and a perfect example of what can higher-order determinants. But as as you can with scalars. happen when embedded programmers go bad. you might surmise, But wait, there’s there exist recurmore. Eagle eyes will have sive algorithms to seen a gotcha in Equation compute them from lower-order determinants. Writing the 3. If the value of x is zero, its inverse doesn’t exist—or, rather, computer code is easier than writing the math equation. is a mathematical inﬁnity. The same is true of the matrix inThe important thing about determinants and singular verse, only more so. Surely the inverse of a matrix like: matrices is that, just as the CPU catches an attempt to divide by zero, we must catch an attempt to invert a singular ⎡0 0 0 ⎤ ⎢ ⎥ matrix. In practice, this means we must calculate the deter0 = ⎢0 0 0 ⎥ (9) minant as we calculate the inverse. If it ever becomes zero, ⎢0 0 0 ⎥ ⎣ ⎦ we must call a halt and report failure. When a CPU detects an attempt to divide by zero, it issues a hardware exception. In a perfect world, we should be cannot exist; there is no matrix we can multiply 0 by to get able to throw a software exception for the singular matrix, as the unit matrix I. But for matrices, 0 is not by any means the well. However, the C++ exception mechanism is complicatonly matrix whose inverse doesn’t exist. ed and usually not allowed in embedded software. We also But if a normal-looking matrix can fail to have an incan’t just send an error message to the console. In general, verse, how can we tell? What can be the rule by which we identify it? To see the answer, you need only think back to the embedded systems don’t have consoles. So instead, I’m rekinds of linear algebra problems that led us to matrices in the turning a ﬂag of sorts: the determinant value. If it’s zero, the inversion failed. ﬁrst place. If I have some set of linear equations like: mInv, the matrix inversion algorithm I’m providing in ax + by = p SimpleMat.cpp, is an oldie but goodie. I cribbed it from (10) the (uncopyrighted) IBM Scientific Subroutine Package, ca. cx + dy = q 1960. I’ve discussed this algorithm at length in earlier issues I can only solve for x and y if the two equations are linof Embedded Systems Design, so I won’t go into it in detail— early independent. That is, they can’t just be scaled versions of perhaps later. The algorithm uses Gauss-Jordan elimination each other. When, in high school, we used to say “n equations with full pivoting. This makes it something of an oddity in in n unknowns,” what we really meant was “n independent the world of numerical analysis. equations . . . .” About pivoting: as the Gauss-Jordan process unfolds, So how can we be sure that the equations are really indethe key elements turn out to be those down the diagonal of pendent? The answer turns out to be: the determinant of the the matrix. We end up dividing by those elements, so needcoefﬁcient matrix is non-zero. A matrix whose determinant less to say, zeros on the diagonal are bad. In fact, near-zero is zero is called singular, and a singular matrix has no inverse. values are also bad. To get the best numerical accuracy, we’d And how, you ask, can I calculate the determinant? Ah, prefer to have large (positive or negative) elements there. that’s the question, and it’s one I won’t answer directly, since The term “pivoting” refers to a process of seeking out the the deﬁnition is horrible. The determinants of small matrices larger elements and moving them to the diagonal, by row, are, however, easy to write down: and perhaps column, transformations. I have seen, in embedded systems, homebrewed matrix inversion software that did no pivoting at all. This is both a ⎡a ⎤ = a ⎣⎦ terrible idea and a perfect example of what can happen ⎡ a11 a12 ⎤ when embedded programmers go bad. It’s perfectly possible ⎢ ⎥ = a11a22 − a21a12 for a matrix to have a zero on the diagonal but still be non⎢a21 a22 ⎥ ⎣ ⎦ (11) singular. By performing no pivoting, the programmers risk ⎡ a11 a12 a13 ⎤ 1 having an algorithm fail when it doesn’t need to. ⎢ ⎥ a21 a22 a23 ⎥ = a11a22a33 + a12a23a31 + a13a21a32 Most implementations of the Gauss-Jordan algorithm use ⎢ ⎢a ⎥ partial pivoting: they seek out the largest element but only in a ⎣ 31 a32 a33 ⎦ given column. Full pivoting looks in the remaining elements of − a11a23a32 − a12a21a33 − a13a22a31 both rows and columns, thereby assuring the largest possible element at each step. Full pivoting takes a little longer. When 14 AUGUST 2008 | embedded systems design | www.embedded.com http://www.embedded.com

Table of Contents Feed for the Digital Edition of Embedded Systems Design - August 2008

Embedded Systems Design - August 2008 Contents Number Include Parity Bit Programmer's Toolbox Cover Feature: Virtualization for Embedded X86 Multiprocessor Applications Wanted: Industry Standards for Benchmarking Embedded VMM Hypervisors Achieving Cache Coherence in a MIPS32 Multicore Design Memory Allocation in C Advertising Index Break Points Marketplace

Embedded Systems Design - August 2008

For optimal viewing of this digital publication, please enable JavaScript and then refresh the page. If you would like to try to load the digital publication without using Flash Player detection, please click here.