Implementing the Efficient Matrix Solution for Special Classes of Engineering Problems

Yesterday I described how to set up a matrix for efficient solution, today I’ll describe how to actually execute the solution. I plan to go back and revisit the intermediate steps but I wanted to show how this series of techniques pays off. I also wanted to leverage this work to describe my feel for the amount of computing power that might be available for JavaScript execution on various computing platforms, which discussion started on Monday.

Here is an implementation of the symmetric banded matrix solution written in C, which is very easily converted into similar languages (Java, JavaScript, etc.). I chose this example because it was working code from a project I did back in the 90s, and I know it’s correct. A few things can be done to make it more compact (using the += formation, for instance) and there may be other clever rearrangements that take advantage of unique aspects of various languages, compilers, implementations, and so on, but mostly it’s ok to let the compilers do their magic on the code as it is. I also chose C because I can write text out to a file easily. The same would be true if I wrote the code in Java or any other language and ran it directly so it could write to disk. This could also be done in JavaScript but you’d have to figure out what to do with the output (here’s a discussion of one possibility).

This code is also repurposed from the FORTRAN I originally used at Westinghouse while building simulator models. That’s why you see all the capital letters and really short control variable names. The names of the constants were my own twist on Hungarian notation, which has been depopularized in favor of CamelCase names that don’t try to say anything about data types, since refactoring and other operations may override whatever type you were using to begin with. There is also the potential for the default size of standard variables types to vary on different platforms or implementations. The policy on the simulator was that the first letter of a variable represented a specific property (e.g., P=pressure, T=temperature, Y=admittance, X=concentration, V=volume). The next three letters represented the various plant systems (e.g., RWU = Reactor Water Cleanup, OGS = Offgas). That left us four whole letters to differentiate the various examples of each type (say, temperatures) in a system. It was common for systems of the era to limit variable names to eight characters, not to mention 8.3 file names or something similar. I adopted that basic system for my own thermodynamic models and control systems while I worked at Bricmont. The “itc” prefix meant short integer, thermal model, constant. To that I might have added an “a” for array. “u” as the first letter would indicate a structure (or record, Java and JavaScript just use objects), “r” was for floating point numbers, “f” for double-precision floats, “s” for character strings, “l” for logical/boolean, and so on.

There’s one more bizarre thing I did here. C, Java, and the like use zero-based arrays while Pascal and FORTRAN use one based arrays (Pascal’s situation is even more complex but…). When I translated this to C I preserved its “one-basedness.” It shouldn’t be too much trouble to update.

Here’s the solution code in its most basic form. The ZMATX array is the 2-dimensional coefficient matrix and ZSRCE is the 1-dimensional source or driver matrix. This matrix will hold the results when the process completes. Both matrices must be loaded as described yesterday before the function is run. The ZAUX matrix is used on a temporary basis during execution.

Now comes the trick. In the next block of code I’ve embedded statements that write the previous line of code out to a file, with array indices replaced with the constant values for that iteration. I don’t write out any of the loop or control statements or any of the ancillary calculations. I only write out the statements that do actual matrix calculations.

Running this function generates the code you can download via this link. It’s over 11,000 lines and over 500k in storage, so it clearly can’t be displayed directly on a web page, either as a listing or as client-side code (JavaScript). This particular method is appropriate for a dedicated computer, or perhaps a server, if the number of concurrent sessions is limited. Note that I’ve also added a bit of wrapper information, but the beginning of the generated code is indicated with a comment. The automated code generator can easily be expanded to include any desired wrapper code and commentary. Also note that this version of extreme or total loop unrolling can be used to realize maximum efficiency, but it is also possible to use these techniques on a limited basis. You’ll still save a respectable amount of time, but at the cost of far less output code. The balance is up to you.

Here are the first hundred-plus lines of the output as an illustration:

Writing the code out longhand this way saves all of the extra calculations associated with incrementing loop variables, jumping to the beginning of loops, and calculating the indirect indices of the various array variables. This cuts execution time for array solutions by thirty percent. Since matrix calculations make up the bulk of the model system’s running time, these modifications represent the maximum performance improvement possible. That means there was no point in trying to identify major performance improvements elsewhere in the system. It was more important, then, to make sure the rest of the code was clear and maintainable.