Scientific Computing and Associated Fields Resource Guide

This is a summary of Internet-related resources for a handful of fields related to Scientific Computing, primarily:

This FAQ is available at: www.mathcom.com

Table of Contents:

Specialized Subfields Within Numerical Analysis:

Associated Fields:

Copyright 1995-2002 S. J. Sullivan. This document(s) may be copied and/or reproduced providing that the use is for non-commercial purposes only, and all copies contain this paragraph and copyright notice. The information contained in this document(s) is believed to be true but no guarantees of accuracy are made, and there is no liability of any sort for any consequences of its use.

Steve Sullivan      Mathcom, Inc.      info@mathcom.com

Introduction

Where to find this FAQ: On the web: www.mathcom.com

"[]" Reviews are associated with the name of the reviewer in brackets. Those reviews marked [SJS] are by myself. Unmarked text is generally taken from the package documentation.

Instead of the normal question/answer form, this FAQ is organized as an outline ... hopefully, you'll find your questions answered here.

Overview of Recent Additions

There have been many updates since the last version; some are:

Paul Black's Computational Dictionary is an extensive dictionary of algorithms, data structures, and problems.

deal.II is a C++ class library targeted at adaptive finite elements. See Deal.II

Clemens Valens wrote a wavelet tutorial and associated info. See Wavelets

SBmethod is a C++ code for large scale eigenvalue optimization. See Optimization, Linear and Non-Linear Programming: General Resources

The NEOS Optimization Server solvers offer online solution of optimization problems. Optimization problems are solved automatically with minimal input from the user. Users only need a definition of the optimization problem. See Optimization, Linear and Non-Linear Programming: General Resources

The ZunZun server offers online interactive fitting and graphing for 2D and 3D data sets. See ZunZun Interactive 2D and 3D Data Modeling

GiNaC is an open framework for symbolic computation within the C++ programming language. See Ginac

Ch is an interpretor for C++ numerics. See Ch Interpretor for C++ numerics

UMFPACK is a package written in C for sparse LU factorization See Sparse Linear Algebra Systems

Mondriaan is a sparse matrix partitioning package. See Sparse Linear Algebra Systems

EigTool combines the Pseudospectra GUI with MATLAB's eigs command to provide visualisation of the Arnoldi iteration during an eigs computation. See Sparse Linear Algebra Systems

Acknowledgements

Many thanks to all those who've given their time and advice in creating this FAQ. In order to save contributors from unwanted spam, I do not include their email addresses.

Many thanks also to the organizers of the many services listed herein.

What is Numerical Analysis?

NA is the union of theoretical and computational investigation into the computer solution of mathematical problems. NA generally includes those problems involving continuous functions of real or complex variables, as opposed to solely discrete variables and functions.

The mixing of theoretical and computational concerns leads to a strong emphasis on algorithms: what are the time and memory usage properties of a certain algorithm? What errors are introduced by an algorithm?

The computational aspects of NA usually take place within the scope of floating-point arithmetic, and are implemented on machines ranging from super-computers through PCs to hand-calculators. The theoretical aspects extend into fields such as Calculus, Differential Equations, and Analysis. The field of Linear Algebra is so often used to model physical systems that the theoretical study of Linear Algebra is in itself often considered to be NA at work.

Primary areas of theoretical concern in NA are:

Primary areas of computational concern in NA are:

Indices of NA Software on the Net

For indices of packages oriented towards symbolic algebra, see Symbolic Algebra .

One of the largest is the NIST Guide to Available Mathematical Software (Formerly called GAMS). [SJS]: The guide is maintained by the US National Institute of Standards and Technology (NIST) It is an index and server for a wide variety of mathematical software, including most of netlib (see Netlib ). Much of the software is in Fortran. If you prefer to speak C++ or C, see C++ Resources , and Fortran, C, and f2c .

[Ronald Boisvert]: The main focus is on fine-grained software components, e.g. subroutines, although information about some larger packages are included. As of November 1995, nearly 10,000 components from more than 90 packages have been cross-indexed using a detailed tree-structured problem classification system. Both freely available software (from netlib or developed at NIST) and commercial packages (used by NIST) are indexed, although source code is available only for non-commercial software.

Mathtools.net is a large tech computing portal.

SAL (Scientific applications on Linux) is a large index of software available for linux.

Jack Dongarra's survey is a list of freely available linear algebra software. It compares many free linear algebra packages. It does not discuss the many packages derived from these core solvers.

NA Software Libraries on the Net

Libraries are collections of source code, and source code packages. Much of the code is in Fortran. If you prefer to speak C++ or C, see C++ Resources , and Fortran, C, and f2c . The main library by far is Netlib .

For statistical software, the best resource is Statlib .

Other libraries are

See also the sections on:

Netlib, including LAPACK

Netlib is probably the world's largest repository of numerical methods programs. It is located at Oak Ridge National Laboratory, Knoxville, Tennessee, and at AT&T Bell Laboratories, Murray Hill, NJ.

Some gems of netlib: LAPACK provides a wide variety of linear algebra functions:

LAPACK++ is a C++ version of, sadly, only a subset of LAPACK. LAPACK++ is work in progress, and hopefully the full functionality of LAPACK will be supported soon.

ScaLAPACK is for distributed memory machines.

Machine/architecture dependent Basic Linear Algebra Subroutines (BLAS) are the keystone of Netlib.

Fortran, C, and the f2c Translator

For C++ and C resources, see C++ Resources .

f2c: Most of the programs in netlib are in Fortran. However, netlib contains an excellent Fortran-to-C conversion utility, f2c. While f2c produces working C code, it is visually complex and ugly. Using f2c on a large package like LAPACK can require a good deal of time to get all the options correct. Fortunately, LAPACK has already be converted to C: see CLAPACK.

The utility f2c can also be invoked by email. Send email to netlib@research.att.com, with the subject "execute f2c", and body containing the non-confidential Fortran program to be converted. But the email option is of use only for very small, simple programs, since a resulting C program of any size must be linked with the f2c libraries. Usually one will have to download the f2c package anyway to generate the libraries. Generally it's easier to download the f2c package, build the libraries and the f2c conversion program, and do the conversion locally.

CAUTION: Programs created by f2c conversion use parameter passing conventions different from most C or C++ programs. Their callers must create the appropriate parameters before using them. See the file f2c.ps in the f2c distribution. A good description of this issue may also be found in the "readme" file for clapack in netlib.

Statlib

Statlib is a huge repository of statistics related software and info. Probability, statistics, random variables, distribution functions.

UCAR's Mathematical and Statistical Libraries

UCAR math libraries contain some overlap with netlib.

Hensa Unix Parallel Archive

The Hensa archive contains general info, software, articles, etc., on parallel computing.

Modula-3 NA Library

The Modula-3 library (old link: http://m3.polymtl.ca/m3/pkg/contrib/m3na/) is a collection of numerical analysis routines written in Modula-3. Includes linear algebra, roots, ffts, and a bit of statistics.

Forth Numerical/Scientific Library

One library of Forth scientific software is at Skip Carter's site at Taygeta

Eiffel Numerical/Scientific Library

Some commercial Eiffel scientific software is at eiffel.com

Lisp Numerical/Scientific Libraries

Some widely used lisp scientific libraries are:

Lisp-Stat

U. Waterloo Quail

U Mass instrumentation and statistical analysis packages

CMU's library

Java Numerical/Scientific Libraries

The Visual Numerics' Java package, JNL, is a Numerical Library for Java, is a set of classes for the most important numerical functions missing in Java. The library is comprised of one numerical type class, Complex, and three categories of numerical functions classes: the special functions class, the linear algebra classes, and the statistics class. All classes use double precision floating point as the underlying float type.

The JLAPACK project provides the LAPACK and BLAS numerical subroutines translated from their Fortran 77 source into class files, executable by the Java Virtual Machine (JVM) and suitable for use by Java programmers. This makes it possible for a Java application or applet distributed on the web to use established legacy numerical code that was originally written in Fortran. The translation was accomplished using a special purpose Fortran-to-Java (source-to-source) compiler.

Linear Algebra for Statistics Java Package

JAMA is a basic linear algebra package for Java. It provides user-level classes for constructing and manipulating real, dense matrices. Five fundamental matrix decompositions are provided:

JAMA is by no means a complete linear algebra environment. For example, there are no provisions for matrices with particular structure (e.g., banded, sparse) or for more specialized decompositions (e.g. Shur, generalized eigenvalue). Complex matrices are not included. It is not our intention to ignore these important problems. We expect that some of these (e.g. complex) will be addressed in future versions. It is our intent that the design of JAMA not preclude extension to some of these additional areas.

The SciMark Java numerical benchmark consists of various kernels (FFT, Monte Carlo, sparse matrix computation, finite-difference stencils, and LU factorization) and is meant to provide an indication of how well Java environments perform on numeric and scientific applications. SciMark can be run interactively within your browser, or downloaded to run in other Java environments. The web site includes bar-graph comparisons between various computer/Java platforms, as well as an archive of previous results.

Lucent Libraries

The Lucent Software Libraries include the following items:

PLAPACK: Parallel Linear Algebra

PLAPACK is an MPI based Parallel Linear Algebra Package (PLAPACK) designed to provide a user friendly infrastructure for building parallel dense linear algebra libraries. The Users' Guide, "Using PLAPACK: Parallel Linear Algebra Package" is available from The MIT Press. WHAT IS DIFFERENT: PLAPACK provides three features not currently found in other publically available parallel dense linear algebra libraries:

WGS: SLICOT Control Theory Libraries

The SLICOT library objectives are to bring together the existing numerical software for control and systems theory in a widely available library, called SLICOT, and to extend this library to cover as far as possible the area of industrial applications.

GSL: GNU Scientific Library

The GNU Scientific Library (GSL) is a collection of routines for numerical computing in C. The routines have been written from scratch over a five year period by the GSL team using modern coding conventions. The subject areas covered by the library include:

ATLAS: Automatically Tuned Linear Algebra Software

The ATLAS (Automatically Tuned Linear Algebra Software) provides a complete, high performance implementation of the BLAS library, and a small subsection of the LAPACK library. In addition, the associated developer release (ATLAS 3.3) possesses support for Intel's SSE2, allowing for maximal DGEMM performance of around 2 Gflop/s on a 1.5Ghz Pentium 4 (SSE1 provides a roughly 4 Gflop/s peak SGEMM on the same machine). Prebuilt archives are available for many architectures, including well-tested version of the developer release. In particular, SSE2-enabled Pentium 4 libraries are available for both Linux and Windows.

ATLAS is a software package that will automatically generate highly optimized numerical kernels for our commodity processors. As the underlying computing hardware doubles its speed every eighteen months, it often takes more than a year for software to be optimized or "tuned" for performance on a newly released CPU. Users tend to see only a fraction of the power available from any new processor until it is well on the way to obsolescence.

NA Software Packages on the Net

Packages generally include an NA library and an interpretive language for a front end. Also see Symbolic Algebra , for free symbolic algebra packages.

Octave

Octave is considered the closest-to-Matlab of the Matlab clones.

Octave is a high-level language, primarily intended for numerical computations. It provides a convenient command line interface for solving linear and nonlinear problems numerically.

Octave can do arithmetic for real and complex scalars and matrices, solve sets of nonlinear algebraic equations, integrate functions over finite and infinite intervals, and integrate systems of ordinary differential and differential-algebraic equations.

The Octave distribution includes a 200+ page Texinfo manual. Two and three dimensional plotting is fully supported using gnuplot.

The underlying numerical solvers are currently standard Fortran ones like Lapack, Linpack, Odepack, the Blas, etc., packaged in a library of C++ classes.

RLaB

RLab is an interactive, interpreted scientific programming environment. Rlab is a very high level language intended to provide fast prototyping and program development, as well as easy data-visualization, and processing.

Rlab is not a clone of languages such as those used by tools like Matlab or Matrix_X/Xmath. However, as Rlab focuses on creating a good experimental environment (or laboratory) in which to do matrix math, it can be called "MATLAB-like" since the programming language possesses similar operators and concepts.

Scilab

Scilab is another good Matlab clone.

Scilab is a matrix-based scientific software package resembling Matlab-Simulink and Xmath-SystemBuild. Scilab contains hundreds of built-in mathematical functions, a rich set of data structures which includes polynomials, rationals, linear systems, lists, sparse matrices and comes with a number of specific toolboxes for control, signal processing, ...

It features: Elaborate data structures (polynomial, rational and string matrices, lists, multivariable linear systems,...). Sophisticated interpreter and programming language with Matlab-like syntax. Hundreds of built-in math functions (new primitives can easily be added). Stunning graphics (2d, 3d, animation). Open structure (easy interfacing with Fortran and C via online dynamic link).

Many built-in libraries:

Tela

Tela is a general NA package with graphics, linear algebra, FFT, etc. Is this another Matlab clone?

It is mainly targeted for prototyping large-scale numerical simulations and doing pre- and postprocessing for them, and it replaces a compiled language like C++ or Fortran in this respect. The feature set is therefore biased to operations needed in partial differential equation solvers.

Euler

Euler started as a MatLab clone. It is now a program, which can handle real, complex and interval numbers and matrices, has a 2D/3D graphics, a builtin modern programming language (extension of MatLab's), an exact scalar product, and the Windows 95 version can call functions in an external DLL. The OS/2 and Windows versions interact nicely with the GUI, and have a notebook style interface.

The Unix version is free, the OS/2 version free for educational use, and the Windows version cheap shareware.

These features make EULER an ideal tool for the tasks such as

Prophet

Prophet is an NIH-sponsored Unix workstation software package for life science computing. Prophet includes tools for data management, statistical analysis, curve fitting, data graphing, mathematical modeling, and genetic sequence analysis.

One of PROPHET's greatest assets is its new graphical user interface . Employing the latest advances in software technology, PROPHET lets you store, analyze and present Data Tables, Graphs, Statistical Analyses and Mathematical Modeling, and Sequence Analyses with high-resolution graphics and multiple windows. Anyone, from the computer-naive to the computer-sophisticate, can learn to use it quickly and effectively. PROPHET is a National Computing Resource for Life Science Research sponsored by the National Center for Research Resources of the National Institutes of Health.

Unfortunately, prophet is distributed in binary form only.

Yorick

Yorick is an interpreted language. It has:

Because Yorick can read either text or binary files, it can be used "out of the box" as a pre- and post-processor for most existing physics simulation programs.

As a pre-processor, you can write a Yorick program that produces complicated input files for a simulation. These might be based on output from other programs, or might require evaluation of complicated functions or involve a lot of repetition.

As a post-processor, Yorick allows you to compare the results of several simulations or to analyze results of a single simulation in ways you did not foresee when you ran it.

SpYorick

SpYorick is a plug-in for Yorick that adds support for sparse matrices, operations among matrices and vectors and solution of systems of linear equations. See Yorick .

PETSc

The Portable, Extensible Toolkit for Scientific Computation (PETSc) provides many tools for the parallel (and uniprocessor), numerical solution of PDEs that require solving large-scale, sparse nonlinear systems of equations. PETSc includes nonlinear and linear equation solvers that employ a variety of Newton techniques and Krylov subspace methods. In addition, PETSc provides several parallel sparse matrix formats, including compressed row, block compressed row, and

PETSc is a set of parallel software libraries for the implicit solution of PDEs and related problems. New features include:

The Blitz++ Numerical Library Project

Blitz++ is a C++ template class library for scientific computing. It offers a high level of abstraction, but performance which appears to be rivalling that of Fortran. The current alpha version supports arrays and vectors. Matrices are only partially implemented and undocumented; use at your own peril.

Ascend

ASCEND IV is a free, large-scale, equation-based, environment featuring a strongly-typed, object-oriented model-description language. ASCEND is designed to reduce the time needed for creating, debugging, and solving mathematical models by orders of magnitude in comparison with C++-like and FORTRAN-like languages.

ASCEND includes interactive support tools for modeling, debugging, and solving systems with tens of thousands of nonlinear algebraic or differential equations. Including:

Algae, formerly called Alki

Algae is a high-level interpreted language for numerical analysis. Algae borrows from languages like MATLAB and APL, but it was developed because we needed a free, efficient, and versatile language capable of handling large problems.

Algae is fast. It's generally faster than MATLAB, RLaB, and Octave, often by a significant margin. Of course, words like "generally" and "often" (and benchmarks themselves, for that matter) don't mean much when it comes time to run your particular problem, and I'm sure that you can find cases where Algae gets beat. Still, Algae is fast.

Algae's arrays may be stored in sparse form; only the non-zero elements and their locations are stored. This type of storage is required for practical analysis in many fields. In my own field (structural dynamics), a matrix with 20,000 rows and columns is not considered large.

Commercial NA Libraries and Packages

Commercial libraries and packages tend to merge, so I've combined them in one category. Typically a commercial product contains:

Many symbolic algebra packages also contain NA packages. For info on these packages, see Symbolic Algebra .

NAG

NAG provides numerical, symbolic, statistical, and visualization libraries in Fortran 77, Fortran 90, C, Pascal, Ada, and parallel machine versions. High performance Fortran 90 and Fortran 77 compilers.

IMSL and PV-WAVE

IMSL from Visual Numerics, Inc. is a set of routines in C, C++, and Fortran libraries for general NA, statistics and graphics. PV-WAVE is a visual programming environment that includes IMSL as a "plug-in".

PV-WAVE is a software environment for solving problems requiring the application of graphics, mathematics, numerics and statistics to data and equations.

Matlab and Simulink

Matlab from Mathhworks is an interactive general NA package, including graphics. A huge variety of "toolboxes" are available, both from the vendor and on the net, for various specialized NA areas: control systems, neural nets, optimization, symbolic math, and on and on.

Simulink is modeling, simulation, and system analysis tool.

For a comparison of Matlab and IDL, see Comparison of IDL and Matlab .

MATLAB is a technical computing environment for high-performance numeric computation and visualization. MATLAB integrates numerical analysis, matrix computation, signal processing, and graphics in an easy-to-use environment where problems and solutions are expressed just as they are written mathematically - without traditional programming.

MATLAB has evolved over a period of years with input from many users. In university environments, it has become the standard instructional tool for introductory courses in applied linear algebra, as well as advanced courses in other areas. In industrial settings, MATLAB is used for research and to solve practical engineering and mathematical problems. Typical uses include general purpose numeric computation, algorithm prototyping, and special purpose problem solving with matrix formulations that arise in disciplines such as automatic control theory, statistics, and digital signal processing (time-series analysis).

MATLAB also features a family of application-specific solutions that we call toolboxes. Very important to most users of MATLAB, toolboxes are comprehensive collections of MATLAB functions (M-files) that extend the MATLAB environment in order to solve particular classes of problems. Areas in which toolboxes are available include signal processing, control systems design, dynamic systems simulation, systems identification, neural networks, and others.

SIMULINK is a tool for modeling, analyzing, and simulating an extraordinarily wide variety of physical and mathematical systems, including those with nonlinear elements and those which make use of continuous and discrete time.

WavBox

Wavbox is a wavelet Toolbox for Matlab. A software toolbox for wavelet transforms and adaptive wavelet packet decompositions with new search algorithms. Requires Matlab.

IDL

Research Systems Inc. IDL, the Interactive Data Language, is software for data analysis, visualization, and cross-platform application development. IDL combines all of the tools you need for any type of project, from "quick-look," interactive analysis and display to large-scale commercial programming projects. From Research Systems, Inc.

For a comparison of IDL and Matlab, see Comparison of IDL and Matlab .

Following are two sets of comments on IDL:

Comments by Pierre Maxted

I find that IDL is good for "playing" with data. This works well for astronomers who seem to end up always wanting to do something a little different to last time to data that always has slightly different quirks every time. I also find that it is a rather easy language in which to write my own routines. This is probably because I can start with interactive IDL to get the feel for what the data is like and what I want to do with it - this then becomes a simple batch file which can be turned into a routine if the procedure is useful - this seems to be a natural way to develop things. These libraries of routines are what makes IDL really powerful in my opinion. I found that adding the astronomy user's library to IDL was like adding wheels to a car. I would recommend to anyone considering using IDL to find out what libraries are out there (e.g. starting at the IDL WWW home page).

Whatever you add to the FAQ, make one point clear - calling IDL a fancy plotting package is like calling a Formula 1 racing car good for picking up the kids from school - IDL can do plotting, but that is not its strength.

Well, I agree that the hard copy manuals are rather opaque but Version 4 of IDL has online help (Hyperhelp) that is rather good - especially since it had text searching capabilities so that you can go straight to the bit you need (usually).

Comments by Amara Graps

Following is an excerpt of comments by [Amara Graps]: For the full text of her review, please see: Amara Graps' Papers

If you install IDL without a valid license, you will get IDL's 7 minute demo mode. This mode is designed for users who are considering buying the package.

IDL is a vector-based language that makes it easy to manipulate arrays and matrices. I've done testing comparing IDL speed to Fortran in various actions, and IDL was as fast as a Fortran program for the IDL array computations where loops were removed (i.e., when using implicit loops in IDL instead of explicit FOR statements).

The scientific functions and procedures that come with IDL are often all that scientists need. In addition, there are net archives containing contributed routines. The archives at John Hopkins and at Goddard are especially good (see below).

The language, for the most part is "open", i.e. you can see the text of any particular procedure or function, in case you doubt the technique, or want to modify it. Some functions and procedures are black-box, intrinsic functions or procedures, but not nearly as many as Matlab (see below) are.

Most work in IDL is done at the command line level. However, IDL supplies rudimentary "widgets" to wrap a GUI around your procedures and functions. You can create buttons, menus, scrollboxes etc.

Three-d plotting is currently not very well documented, and the way that IDL does it is very convoluted. Other users and I have complained about it, and I think RSI are taking steps to better document how to do it.

Image processing and animation is pretty slick. If you need to do "slicing and dicing" of a volume, in a way like Spyglass Dicer, IDL has a really great widget routine to do it. The IDL plots are high quality enough to use in initial journal submissions.

RSI's support (writing to support@rsinc.com) is pretty good, I usually get responses within 24-34 hours. You have to pay yearly technical support costs, though- about $200 year (don't remember exactly how much). The Usenet group: comp.lang.idl-pvwave has some smart programmers giving answers if you don't want to pay for the IDL technical support. RSI usually doesn't answer questions on that newsgroup (they have a company policy against promoting IDL there because it's shared by two products: IDL and PV-WAVE).

I've never liked the IDL documentation very much. The information that you need probably *is* in the manuals, but it's somewhat hard to find (the manuals are organized in a weird way). [Note, however, the comments by Pierre Maxted above].

Public IDL code

The anonymous ftp sites below contain public domain IDL code. JHU/APL/S1R IDL library NASA IDL Astro Library IUE RDAF library at NASA IUE RDAF library at U. of Colorado ICUR Spectral Analysis Software IDL ROSAT software IDLmeteo library ESRG library Hal Mueller has a Digital U.S. Map browser based on images created by Ray Sterner at Johns Hopkins University using IDL.

Comparison of IDL and Matlab

Following is an excerpt of a paper by [Amara Graps]: For the full text of her review, please see: Amara Graps' paper.

IDL is a package that began life as an image-processing utility that has grown to be a general-purpose numerical analysis tool. Matlab started as a numerical analysis package that now includes [at extra cost] image processing tools. Now the two have a similar scientific data-analysis environment, with capabililties to build GUI programs and do very robust data analysis.

(Note: all prices are approximate October 1995 prices - SJS) They each cost about the same: ~$1500 for Mac and PC versions and more for Unix (~4000 -- single user to $15,000 -- unlimited number of users).

Matlab is popular among education institutions because it has exceptional educational discounts. If you are an academic, Matlab can be had for $495 and each toolbox only $195. My NASA colleagues thought that MathWorks "nickled- and-dimed" them with the costs of the Toolkits (like the signal processing toolkit), but given what you get, it probably isn't that unreasonable.

IDL seems to be more widespread in the NASA communities probably because the original developer used several spacecraft teams (Pioneer Venus and Voyager) as test beds for the IDL software.

IDL is more of a true programming language. Matlab has scripts and functions and no way to explicitly type a variable. IDL has programs, procedures, and functions and a language syntax sort of like a cross between Fortran, Pascal, and APL. If you have programmed in Fortran before, then the syntax will be a snap to learn. Matlab's syntax is much more compact than IDL's. For example: x = transpose(y) in IDL is x=y' in Matlab.

Matlab has many more built-in, intrinsic functions than IDL. MatLab has many optional Toolkits, such as a Signal Processing Toolkit and an Image Processing Toolkit, which are libraries of more intrinsic functions.

Reading and writing files, and handling formats such as GIF, PICT, GDF, and custom formats, seems much easier in IDL than MATLAB. Handling directories is difficult in MATLAB when run on non-unix machines.

Matlab has more types of graph types than IDL, and handling colors is simpler than IDL. However, I found most other Matlab graphical programming non-intuitive. It uses a system where each element in a graph is an "object." These objects can have sub-objects. So to change an element in a graph, say the axis color, you have to first find the object (a "get" function), and then set it to the color you want. IDL has system variables storing all graphics elements which can be easily changed. One can also customize a graph upon making the graph, with a keyword.

IDL's technical support is pretty good, but Matlab's is better. Post a question on comp.soft-sys.matlab and either a developer, the company president, or a tech support person will respond that day. You can call them, too, but it's not a toll-free call.

MLAB

MLAB (for Modeling LABoratory) is a program for interactive mathematical and statistical modeling. MLAB was originally developed at the National Institutes of Health. It includes curve-fitting, differential equations, statistics and graphics as some of its major capabilities.

Gauss

GAUSS, from Econotron Software, is available for IBM PCs and compatibles as well as UNIX workstations

As a complete programming language, the GAUSS system is both flexible and powerful. Immediately available to the GAUSS user is a wide variety of statistical, mathematical and matrix handling routines. Powerful data handling capabilities including a data loop allow transformations in a data set by directly using variable names in expressions. This greatly simplifies data transformations and makes for shorter more readable programs. GAUSS can be used in either command mode(interactively) or in edit mode. In command mode; one-line commands, or small screen-resident programs, can be issued and the results of calculations seen immediately. In edit mode you can write complex programs and store them in files.

GAUSS has over 400 functions built in, including LINPACK and EISPACK routines.

MathViews

MathViews for Windows is matlab look-alike. It has a full set of linear algebra and signal processing functionality. It provides easy access to: matrix and linear algebra, digital signal processing, instrument control, image processing, time series analysis, data visualization and waveform display and editing. MathViews is highly compatible with the matlab syntax and will execute most matlab m-files with no changes. We also have WaveTool. WaveTool is an interactive software tool for creating, editing and analyzing captured waveshapes. Waveforms can be created using any combination of drawing, math expressions (matlab syntax), insertion from a library of waveforms or data values pasted from other applications such as Microsoft Excel.

Matlab to C++ Compiler and C++ Matrix Class Library

MATCOM V2 is a Matlab(R) to C++ compiler. MATCOM creates MEX files and standalone C++ applications, with royalty free distribution. MATCOM translates Matlab code to C++, which is compiled by your optimizing C++ compiler. The resulting code runs significantly faster than the original interpreted source. Prior knowledge of C++ is not necessary to use MATCOM. The compilation is fully automated by a smart project manager. Fully functional, time limited evaluation version of MATCOM V2 can be downloaded freely from the MathTools web site.

MATLIB, a Matlab Compatible C++ Matrix Class Library, is designed for development of advanced scientific high-level C++ code. Evalution version of the MATLIB can be downloaded from the home page noted above. The library includes Complex math, Binary and unary operators, Powerful indexing capabilities, Signal processing, File I/O, Linear algebra, String operations and Graphics. Over 300 mathematical functions are included in MATLIB. MATLIB supports matrices of doubles, floats, ints and chars mixed in the program. Images can be stored in matrices of chars, using 1/8 memory storage. On many applications, where 8 digits of precision are sufficient, float-precision matrices can save half the memory usage. Memory allocation and de-allocation is managed automatically.

O-Matrix

O-Matrix is a complete solution for analysis, visualization, and rapid application development, (RAD). The interactive environment enables both point-and-shoot oriented analysis for rapid solutions and a powerful language and debugger for large demanding applications. The integrated matrix language provides performance that far exceeds typical integrated environments and often rivals compiled code for numerically-intensive applications.

SPT - The Signal Processing Toolbox for O-Matrix: SPT for O-Matrix provides a broad, extensible set of functions and utilities that expands the capabilities of O-Matrix in the area of digital signal processing. Built on the high-performance capabilities and flexibility of O-Matrix, the toolbox contains tools for algorithm development, data analysis, and visual presentation. SPT enables the rapid development of robust, high-performance applications in areas including digital filtering, spectral estimation, digital communication system simulation, time series analysis, real time control, geophysics, and more. The package comes with complete source code, online tutorials and examples, and unlimited free technical support.

KBF - The Kalman Filter Interface Pack for O-Matrix: The KBF interface simplifies the design of a filter or smoother by dividing the design process into simple steps with each step corresponding to an input window or dialog. Extensive plotting of results, including residuals and correlations is also automated. In addition, a simulation feature is included to aid in testing the validity of the design. KBF is used in a wide range of applications such as tracking, weather and earth process modeling, economic forecasting and bioengineering.

UCALC

UCALC for Windows includes many features useful for scientific calculations and graphing:

In addition to the above, UCALC introduces an innovative feature - User Solution Modules (USM). USMs are created by simply providing a math formula, and giving it a name. UCALC parses the formula and creates a template form which allows the user to solve for any item by simply filling in the blanks. Useful USMs are included to get you started.

Laipe

Laipe includes fast parallel solvers for Windows NT running on an SMP computer.

VisualMath

VisualMath is an online symbolic math and computer algebra system. It can perform exact, numeric, symbolic and graphic computation, e.g. arbitrary-precision calculation, solve equation, plot data and user-defined functions, linear regression, symbolic differentation and integration, pattern-match. It is a programming language, in which you can define conditional, case, piecewise, recursive, multi-value functions and procedures, derivatives, integrals and rules. It runs in any computer that supports Java.

Harwell Subroutine Library

HSL 2000 contains codes for Automatic differentiation, Differential Equations, Eigenvalues and Eigenvectors, Mathematical Functions, Sorting, Linear Programming, Linear Algebra, Nonlinear Equations, Polynomials, Optimization and Non-linear data fitting. It contains many new routines, which have been developed as a result of recent research into advanced numerical algorithms.

Ch Interpretor for C++ numerics

Ch is a superset of C interpreter. It is embeddable. C/Ch/C++ allows software developers to use one language, anywhere and everywhere, for any programming and numerical computing tasks.

By incorporating external libraries, it provides access to linear algebra, ODE, integration, graphics, etc.

Newsgroups for NA

Newsgroups related to numerical analysis are:

Professional Societies for NA

The (AMS) American Mathematical Society

The AMS provides:

General organization information, preprint titles, pointers to other preprint servers and net resources. Also includes:

MathDoc, the document delivery service offered by the AMS, provides copies of original journal, collection and conference proceedings articles from publications covered by Mathematical Reviews, Current Mathematical Publications, and the MathSci database. This costs roughly US$14. per ten pages, as of October 1995.

MathSciNet is a searchable database available on the World Wide Web. It is based on the data in Mathematical Reviews and Current Mathematical Publications, leading publications that catalog and review research literature in mathematics. This costs roughly US$5500. per year, as of October 1995.

(SIAM) The Society for Industrial and Applied Mathematics

SIAM provides: General organizational information, tables of contents of SIAM journals, and recently accepted articles. Society for Industrial and Applied Mathematics 3600 University City Science Center Philadelphia, PA 19104-2688 (215) 382-9800

Journals include:

ACM, Inc. (Association for Computing Machinery)

The ACM provides: General organizational information, info on journals and conferences. Of particular interest are: ACM Transactions on Mathematical Software (TOMS)

The collected Algorithms of the ACM (CALGO) software can be found in netlib (toms directory); see Netlib .

IEEE The Institute of Electrical and Electronic Engineers

The IEEE provides general organizational information, journals, books, conferences.

Electronic Newsletters for NA

NA-Net is a system developed to serve the community of numerical analysts and other researchers. The Na-Net provides two independent databases and a weekly digest to its members. The Email Database is the electronic mail address of each of its members, and is capable of forwarding mail to them. In addition, this database serves as the distribution list for the NA Digest (see below). The White Pages Database is basically a directory service. It provides a way to exchange personal information among its members. Contained in the database are phone numbers, postal mailing addresses, research interests, affiliations, etc. The NA Digest is a way to provide its members with a weekly collection of articles on topics related to numerical analysis and those who practice it. To get on the mailing list for the digest and enter yourself in their database, you can use email or the World Wide Web. Take advantage of this very useful service!

Electronic Journals for NA

Indices of Journals See the large list of online journals, and other resources, at Penn State Mathematics. See also:

The Southwest Journal of Pure and Applied Mathematics

BIT emphasizes numerical methods in approximation, linear algebra, and ordinary and partial differential equations, but also publishes papers in areas such as numerical functional analysis and numerical optimization.

ETNA (Electronic Transactions in Numerical Analysis). is an online peer reviewed journal. Keyword searching. Most documents are postscript and may be downloaded. Started in 1993. Indexes only itself.

EJDE, The Electronic Journal of Differential Equations, is a purely on-line peer-reviewed journal.

Electronic J. of Linear Algebra

NY J. of Mathematics Appears more abstract than applied. Started in 1994. Indexes only itself.

Online Papers for NA

General See the listings for professional societies, especially the AMS, in Professional Societies for NA .

Springer Verlag preprints. Not free.

Virtual Lib Math Jnls Preprints links to numerous organizations offering preprints.

The NCSTRL CoRR (Networked Computer Science Technical Reference Library) is a partnership of ACM, the Los Alamos e-Print archive, and NCSTRL (Networked Computer Science Technical Reference Library), an online Computing Research Repository (CoRR) is being established. The Repository has been integrated into the collection of over 20,000 computer science research reports and other material available.

Miscellaneous Web Sites for NA

See also:

Caltech Guide to Math Resources has a good collection of pointers to quality sites.

The huge Penn State Mathematics site has information on many fields within mathematics.

The U.C. Berkeley site includes courseware, pointers to inet libraries, lecture notes, seminars, and software.

The World-Wide Web Virtual Library of Mathematics has info on specialized fields (topology, cryptography, optimization, etc), academic departments, miscellaneous math societies and institutes, pointers to commercial software, newsgroups, nice collection of preprints pointers, electronic journals

U. Tennessee Knoxville Mathematics Archives WWW Server

Mathematical Constants, by Steven Finch, provides constants and algorithms for their generation.

The reclusive Dr. TP's list of www sites for numerical methods:

Amara Graps' list of science links

Catalog of Math Resources, By M. Maheswaran, University of Wisconsin Marathon Center, is another large site, with a sizeable section on Applied Mathematics.

Alan Miller's page provides a collection of code for least squares, random number generation, quadruple precision, and for various statistical purposes. All of this code is written in Fortran 90 and most of it is in the subset of Fortran 90 chosen for the free ELF90 compiler.

Altug Koker's list of simulation software (old link: http://piranha.eng.buffalo.edu/simulation/comp.simulation/FAQ.html) contains pointers to software for: aerospace, automotive, chemistry / biotechnology, graphics and imaging, electronics / electrical engineering, petroleum

Google's numerical analysis index

Scientific Applications on Linux

S. Baum's site has a wide variety of links to mathematical and scientific software.

C++ Resources for NA

Indices:

The Object-Oriented Numerics page has pointers to C++ libraries and classes.

The The Blitz++ Numerical Library Project is a C++ template class library for scientific computing.

Trumphurst FAQ on C++ Libraries

See also Skip Carter's page.

See also Roldan Pozo's page.

Individual C libraries in Netlib :

Also see Joerg Arndt's web page on FFT code, at Transforms (FFT, etc) and digital signal processing (DSP)

Fortran Resources for NA

The Fortran Library is a comprehensive guide to online Fortran resources, including compiler vendors and resellers, benchmarks, programming tools, books and articles on Fortran and numerical methods, and commerical and public domain Fortran software. The Fortran Market provides many links to Fortran resources and is a reseller for several Fortran 90 and 95 compilers, benchmarking and test suites, and books on Fortran.

Math FAQ

Alex Lopez-Ortiz's sci.math FAQ is a good reference for many mathematical questions. It is more oriented towards pure mathematics than NA.

Math Dictionaries and Glossaries

Harvey Greenberg's Mathematical Programming Glossary

Eric Weisstein's World of Mathematics

Weisstein's World of Mathematics developed with the support of Wolfram Research, is one of the web's most extensive mathematics resources.

Paul Black's Computational Dictionary

Paul Black's site contains an extensive dictionary of algorithms, data structures, and problems.

Books, With and Without Software, for NA

For an on-line mathematics dictionary, see: Math Dictionary .

Acton, Forman S. 1990, Numerical methods that [usually] work, Harper & Row, Publishers, ISBN 0883854503. [Daniel Pick] This book is almost worth its price just for the cathartic interlude in the middle of the book on what not to compute. You should require your students to read it, learn it, live it. You may find just giving them the railroad problem found at the beginning of the book a worthwhile exercise. [Bill Frensley] Amen, brother! The only complaint that I have about Acton's interlude is that after demolishing the notion of "fitting exponential data," he fails to point out that this is the inverse Laplace transform problem. Perhaps if everyone read this and made the connection, we would be spared the monthly "is there any good algorithm for the inverse Laplace transform?"

Arge, Bruaset and Langtangen, editors, Modern Software Tools for Scientific Computing. Contains selected contributions from the SciTools'96 conference. Ssee SciTools96 This book consists of papers discussing topics related to computational differential equations, computational geometry and software development. For further information, see birkhauser

Baker, Louis. 1991, C Mathematical Function Handbook

Daehlen and Tveito, editors, Numerical Methods and Software Tools in Industrial Mathematics, Birkhauser It includes an introduction to object-oriented numerics and a presentation of the software tools Diffpack and Siscat. See software tools Also a series of challenging applications are discussed. For more information, see birkhauser

Dahlquist, Germund; Bjorck, Ake, 1974, Numerical Methods translated by Ned Anderson, Prentice-Hall, 1974. A nice mix of theory and practice. Used as a text at Stanford, among other places.[John Chandler]

Dubois, Paul. 1996, Object Technology for Scientific Computing, Object-Oriented Numerical Software in Eiffel and C, Prentice Hall, 1996, paper, 350 pages ISBN 0-13-267808-X

Duff, I.S., and G.A. Watson, editors, 1996, The State of the Art in Numerical Analysis. The 1996 conference on the State of the Art in Numerical Analysis was organized to provide the numerical analysis community, and users of numerical methods, with a forum where an account of the important recent developments in the subject could be presented in a coherent and concentrated way in a manner accessible to the non-specialist in the sub-area.

Forsythe, George; Moler, Cleve B. 1967, "Computer Solution of Linear Algebraic Systems", Prentice-Hall I consider this possibly the best textbook I have ever seen in any field. Covers only linear systems, of course.[John Chandler]

Gautschi, Walter, 1997, Numerical analysis : an introduction, Boston, Mass. : Birkhauser, 1997 A very readable, carefully written book with emphasis on nonlinear aspects of numerical analysis. The book is suitable for a first year graduate course and contains a rich set of examples and exercises. [Jorg Peters]

Golub, Gene H.; Van Loan, Charles F., 1996, Matrix Computations, Third edition, Johns Hopkins, Baltimore, ISBN 0-8018-3739-1, [SJS] A classic for handling matrices. Many current programs are based on this text. Good mix of theory and implementation.

Golub, Gene H., 1984, Studies in Numerical Analysis, Mathematical Association of America, ISBN 0883851261. [Daniel Pick] This contains several outstanding essays from several numerical analysts, including Wilkinson's The Perfidious Polynomial, which explains why rootfinding of polynomials numerically is such a tricky problem. It gives an great introduction to the thinking of recent numerical analysts. [Amara Graps] All of the chapters are really good- my favorites are: "Fast Poisson Solvers" and "Multigrid Methods for Partial Differential Equations."

Kahaner, David; Moler, Cleve; Nash, Stephen, 1989, Numerical Methods and Software, Prentice Hall, ISBN 0-13-627258-4. An excellent book which touches on a variety of topics and makes use of the publicly available software like the QUADPACK and SLATEC libraries to illustrate the points. [Vijay Gupta]

Knuth, Donald E. 1981 The Art of Computer Programming: Seminumerical algorithms, 3rd edition, Addison-Wesley. Once was the reference; now a bit dated.

Lau, H. T., 1995, A Numerical Library in C for Scientists and Engineers CRC Press, ISBN 0-8493-7376-X. This book is basically a compilation of program listings, with a diskette containing source code. The listings are accompanied by brief overviews of the algorithms involved, and generally include references. There is no discussion of theory. While the text by Stoer & Bulirsch is at the theoretical end of the NA spectrum, this text is at the application end. Although the program calling parameters are well described, as far as I could see the programs contain no internal documentation whatsoever. Although this book is copyright 1995, the references contain one source dated 1992 (Press et al's volume), one source dated 1981 (NUMAL in Fortran), and one source dated 1980 (NUMAL in Algol). The remainder of the references are dated 1976 and earlier. It's not clear to me that this book offers anything over Press et al's text. Lau has far less discussion of theory and methodology, and while Press's internal documentation of programs may leave something to be desired, Lau's book has none whatsoever. [SJS]

Mathews, John H. 1992, NUMERICAL METHODS: for Mathematics, Science & Engineering, Prentice Hall Source for the programs is available in several languages: Matlab, C, Fortran, Pascal, Mathematica.

Muller, Jean-Michel, 1997, Elementary functions, algorithms and Implementation, Birkhauser Boston. See: Birkhauser and Muller.

Nash, John C., 1987, 1996, Nonlinear Parameter Estimation: an integrated system in BASIC.

Peters, Martin, 1997, Solving Problems in Scientific Computing Using Maple and MATLAB. All Maple and MATLAB programs can be obtained from: Peters text.

Petkovsek, Marko; Wilf, Herbert; Zeilberger, Doron, 1995, "A=B", AK PETERS, Ltd., ISBN 1-56881-063-6. [Donald Knuth] I'm especially pleased to see the appearance of this book, because its authors have not only played key roles in the new developments, they are also master expositors of mathematics. It is always a treat to read their publications, especially when they are discussing really important stuff. Science advances whenever an Art becomes a Science. And the state of the Art advances too, because people always leap into new territory once they have understood more about the old. This book will help you reach new frontiers.

Press, William H.; Teukolsky, Saul A.; Vetterling, William A.; Flannery, Brian P. 1992 Numerical Recipes in C, Second edition A long time standard practical text. A compendium of a wide variety of NA areas. Contains some good introductions to theory and overviews of algorithms. The bridge from algorithm overview to implementation is often missing. Some of the programs should be viewed with skepticism. Some are poorly documented, and some users have reported numerical problems with the various programs. [SJS] Their more recent C++ book appears to be simply a copy of the C book, with minimal C++ classes wrapped around the code.

The Numerical Recipes in C version may be downloaded for free.

Sedgewick, Robert. 1988, Algorithms, Second edition, Addison-Wesley.

Stewart, Ian. 1999, Notes on Matrix Algorithms I: Basic Decomposition, Notes on Matrix Algorithms II: Eigensystems. Siam Press. Parts may be available online at Ian Stewart's page.

Stoer, Josef; Bulirsch, Roland, Introduction to Numerical Analysis Springer-Verlag, New York, 1980, ISBN 0-387-90420-4. The classic NA text. A standard for graduate and upper undergraduate courses in NA. Covers a wide set of fields in depth. Strong theoretical orientation. [SJS]

Strang, Gilbert, 1988, Linear Algebra and It's Applications, Third Edition, Harcourt Brace, ISBN 0-15-551005-3. A well-written introduction to theory.

Strang, Gilbert and Borre, Kai, 1997, Linear Algebra, Geodesy, and GPS, Wellesley-Cambridge Press. The book has three closely connected parts. The first is Linear Algebra, which is the fundamental tool in positioning calculations. The second is Geodesy, which includes the correct weightings for observation errors and the development of least squares. The third is GPS, which we describe in detail at several levels : first the basic structure of the GPS system, then the algorithms that yield accurate positions from inaccurate pseudoranges, and finally the Kalman filtering (and "Bayes" filtering) that give superior accuracy in postprocessing a long series of observations.

Trefethen, Nick, and Bau, David, 1997, Numerical Linear Algebra, Siam.

Tyrtyshnikov, Eugene E, 1997, A Brief Introduction to Numerical Analysis. Birkhauser, Boston. Tyrtyshnikov writes: Without any doubt, there are many quite good and excellent books on the subject. But I know definitely that I did not realize this when I was a student. In this book, my first desire was to present those lectures that I wished I would have heard when I was a student.

Watkins, David S., 1991, Fundamentals of Matrix Computations John Wiley, ISBN 0-471-61414-9 More readable than Stoer & Bulirsch or Golub & Van Loan, and contains some implementation techniques not present in Golub & Van Loan. Has good descriptions of theory and implementations, and many implementations are covered as straightforward exercises. Not as wide a variety of fields as either Golub & Van Loan or Stoer & Bulirsch. [SJS]

Dense Linear Algebra Systems

The best resources for linear algebra are:

The texts by Golub & Van Loan and by Watkins.

On the net, library package LAPACK, with BLAS, from netlib (see Netlib ) is a widely recommended replacement for EISPACK and LINPACK. It handles dense matrices only.

As far as complete packages, nearly all NA packages, both on the net and commercial, include linear algebra as part of their core. See the listings above.

See also: Jack Dongarra's list of Linear Algebra Software

National HPCC Software Exchange Catalog

Simpler or smaller systems for linear algebra include:

C.R. Birchenhall's MatClass is a C++ class for numerical computation Offers a general purpose dense, real matrix class. Has a family of decomposition classes based on LU, Cholesky, Householder QR and SVD. Has a family of OLS regression classes based on above decompositons. Has a family of special function classes. Random number class. Has a simplified I/O structure. Documents : Very thorough tex manual, with discussion of design philosophy. Currently the manual does not cover all the features of the I/O.

One user comments that: From my experience I can say that is very well done. It covers all the possible statistical distributions, as well as, random numbers generation from these distributions.

newmat0.9 is a matrix library. Simple, well documented, not templated but easy to use. TNT (The Template Numerical Toolkit) is a collection of mathematical libraries for numeric computation in C++. Its fundamental classes include vectors, matrices, and multidimensional arrays. The basic goal is to allow one to express mathematical computation at a higher level of abstraction, while still retaining some control over performance and optimization issues. Doing so requires a careful analysis to balance these tradeoffs.

The toolkit provides an integrated collection of generic matrix/vector classes based on components of the Standard Template Library (STL) and ANSI C++, together with specialization of generic algorithms for maximal efficiency. The TNT component for linear algebra is a successor to the Lapack++, Sparselib++, IML++, and MV++ packages. Its goal is to formally integrate these ideas into a generic algorithmic library, supporting user-defined data types, and increasing its functionality. It also takes advantage of the latest features in the ANSI C++ specification which were not available when we designed our earlier packages.

[SJS] TNT appears to have a good base of basic linear algebra, but is nowhere near as complete as LAPACK. My cursory check found no eigenvalue facilities. I suspect the project was abandoned half done. Correct me if I'm wrong!

Sparse Linear Algebra Systems

Netlib (see Netlib ) has a number of libraries for handling sparse systems. See the netlib directories: sparse, sparse-blas, sparsepak.

See also: Jack Dongarra's list of Linear Algebra Software

National HPCC Software Exchange Catalog

Victor Eijkhout's Overview of Iterative Linear System Solver Packages. [N. Sukumar]: Recommends the TEMPLATES book ["TEMPLATES for the soln of linear systems: building blocks for iterative methods"] that is availabe in postscript format at netlib (see Netlib ). Also in the same dir is code (in C and Fortran) for algorithms discussed in the book. "The book is awesome (easy to read, good pointers on selection of solvers etc.), especially for one who is not aware of how and why sparse solvers work."

SPARSKIT is a linear algebra system for sparse matrices.

PSPASES PSPASES: A Scalable Parallel Direct Solver for Sparse Symmetric Positive Definite Systems is a stand-alone MPI-based parallel library for solving linear systems of equations involving sparse symmetric positive definite matrices. The library efficiently implements the scalable parallel algorithms developed by authors, for each of the four phases of direct solution method; viz. ordering, symbolic factorization, numerical Cholesky factorization, and solution of triangular systems.

AZTEC: A Parallel Iterative Package for Linear Systems, Aztec is an iterative library that greatly simplifies the parallelization process when solving the linear systems of equations Ax = b where A is a user supplied nxn sparse matrix, b is a user supplied vector of length n and x is a vector of length n to be computed. Aztec is intended as a software tool for users who want to avoid cumbersome parallel programming details but who have large sparse linear systems which require an efficiently utilized parallel processing system. A collection of data transformation tools are provided that allow for easy creation of distributed sparse unstructured matrices for parallel solution. Once the distributed matrix is created, computation can be performed on any of the parallel machines running Aztec: workstation clusters (DEC, SGI, SUN, LINUX, etc.), Cray T3E, Intel TeraFlop, Intel Paragon, IBM SP2, nCUBE 2 as well as other MPI platforms, vector machines or serial machines.

Aztec includes a number of Krylov iterative methods such as conjugate gradient (CG), generalized minimum residual (GMRES) and stabilized biconjugate gradient (BiCGSTAB) to solve systems of equations. These Krylov methods are used in conjunction with various preconditioners such as polynomial or domain decomposition methods using LU or incomplete LU factorizations within subdomains. Although the matrix A can be general, the package has been designed for matrices arising from the approximation of partial differential equations (PDEs).

S+ provides sparse LU factorization with partial pivoting is important for many scientific applications and delivering high performance for this problem is difficult on distributed memory machines. This project studies the properties of elimination forests and uses them to guide supernode partitioning/amalgamation and execution scheduling. This design with 2D mapping effectively identifies dense structures without introducing too many zeros in the BLAS computation and exploits asynchronous parallelism with low buffer space cost. The implementation of this code, called S+, uses supernodal matrix multiplication which retains the BLAS-3 level efficiency and avoids unnecessary arithmetic operations. This project also studies two space optimization techniques which can greatly improve the worst-case performance of static symbolic factorization.The experiments show that S+ can achieve up to 10.85GFLOPS on 128 Cray T3E 450MHz nodes, which is the highest performance reported in the literature.

WSMP, the Watson Sparse Matrix Package, (WSMP Part 1, WSMP Part 2) is a robust, high-performance, easy- to-use serial and parallel software for the direct solution of general and symmetric sparse linear systems on IBM RS6000 and SP platforms.

UMFPACK is a package written in C for sparse LU factorization. The current version (V4.0beta) can now handle real or complex matrices, square or rectangular matrices, and non-singular or singular matrices. A simple MATLAB interface is provided. For example, x = umfpack (A, '\', b) is like the MATLAB statement x = A\b.

Mondriaan is a sparse matrix partitioning package. It is a two-dimensional, multilevel, hypergraph-based package for the partitioning of a rectangular sparse matrix, as a sequential preprocessing step for parallel sparse matrix-vector multiplication. It is written in C. It is distributed under the GNU license.

EigTool combines the Pseudospectra GUI with MATLAB's eigs command to provide visualisation of the Arnoldi iteration during an eigs computation. It includes all of the existing features of the Pseudospectra GUI along with several new ones, notably norms of matrix powers and exponentials and a collection of over 30 demos.

Random Number Generators (RNGs)

Web Sites for Random Number Generators

RANDLIB (C and Fortran versions) provides 32 virtual random number generators. Each generator can provide 1,048,576 blocks of numbers, and each block is of length 1,073,741,824. Any generator can be set to the beginning or end of the current block or to its starting value. Packaging is provided so that if these capabilities are not needed, there is a single generator with period 2.3 X 10^18. Using this base, routines are provided that return:

pLab is an excellent web site created by the Department of Mathematics at the Salzburg University. It covers RNGs and tests for randomness, and contains publications, software, and pointers to other resources.

Additional services may be available from pLab -- see: pLab Overview

See also the Diehard page, described at: Tests for Randomness

Another web page with info on RNGs is Skip Carter's page at Taygeta

Also see Pierre L'Ecuyer's papers:

Info on 1/f noise, called "flicker noise" or "pink noise", is at Wentian Li's site.

Some code in Fortran 90 for generating random numbers from a range of discrete and continuous distributions is at Alan Miller's page.

SPRNG is a scalable parallel random number generator library. (It can be used for the usual serial random number generation too.)

Types of Random Number Generators

Random number generators typically consist of three components:

There are four common types of basic generators.

Linear Congruential Generators (LCG)

Linear congruential generators are of the form

x(n) = a * x(n-1) + b mod M

A good discussion of this type is in the text by Press et al. These are by far the most common form in use, and have periods in the range 10^6 to 10^9 when using 32 bit words. Unfortunately, they have some unpleasant properties. When taken in pairs, triplets, or n-tuples, the random values often fall on only a few planes in n-space. Using a shuffle box can break up this uniformity. For more information on the problems associated with linear congruential generators, see Marsaglia 1968, Marsaglia 1993 and Sullivan 1993. See also Marsaglia's web site, described under Tests for Randomness .

Add-with-carry and Subtract-with-borrow Generators (AWCG, SWBG)

Add-with-carry generators have the form:

x(n) = x(n-s) + x(n-r) + carry mod M

Subtract-with-borrow generators have the form:

x(n) = x(n-s) - x(n-r) - carry mod M

See Marsaglia 1991. These generators have much longer periods, from 10^200 to 10^500, than LCGs, and they are almost as fast as LCGs. Unfortunately they too have the problem of falling mainly on the planes. See Tezuka 1993. Again, a shuffle box, as described in Press et al, can eliminate this problem. The implementation by Marsaglia and Zaman may be found at: Marsaglia's generators.

Multiply-with-carry Generators (MWCG)

Multiply-with-carry generators have the form:

x(n) = a*x(n-1) + carry mod M

[George Marsaglia]: Here is a simple version, one so simple and good I predict it will be the system generator for many future computers:

x(n)=a*x(n-1)+carry mod 2^32

With multiplier 'a' chosen from a large set of easily found integers, the period is a*2^31-1, around 2^60, and I have yet to find a test it will not pass!
[Charles Knechtel]: George Marsaglia has some interesting new ideas for pseudorandom number generation using MWCG and MWCG combined with other pseudorandom number generators (RNGs); these new ideas are implemented as part of the Diehard suite of RNG tests. The best of these seem to pass all the Diehard tests, even tests that many common pseudorandom number generators fail. Source is at the Diehard page, described at: Tests for Randomness .

Inversive Conguential Generators (ICG)

Inversive congruential generators are of the form

x(n) = a * ~(x(n-1)) + b mod M

where:
M is prime
"~y" denotes the multiplicative inverse of y in the field over {0, 1, ..., M-1} using arithmetic modulo M, and ~0 is 0.

ICGs do not have the "mainly in the planes" problem. See the pLab, described above at Web Sites for Random Number Generators , for more info on ICGs.

[Charles Knechtel]: Marsaglia's Diehard suite of RNG tests, at Tests for Randomness also contains a subroutine named makeinvc which implements an ICG and generates a file of pseudorandom numbers named INVCONG.32; however, examination of the FORTRAN source code for makeinvc shows that a modulus M equal to 2^32 was used (see file makewhat.f. Eichenauer-Herrmann (1992, p. 173, 175) notes that certain ICGs with power of two modulus have distributional disadvantages compared to the "excellent quality" of ICGs with prime modulus. Indeed, the documentation for INVCONG.32 remarks that makeinvc fails to pass many of the Diehard tests (see the last section of file makef.txt from the Diehard archive fortran.tar.gz. An ICG with prime modulus should have better uniform properties than makeinvc, and Otmar Lendl's ICG implementation might be speedier than makeinvc.

Tests for Randomness

The Diehard suite of RNG generators and tests, by George Marsaglia, is available from: Diehard overview, source, & binaries For snapshots of some of the graphs, see: Diehard snapshots

A GUI version of Diehard is in preparation by Dr. Balasubramanian Narasimhan at: Diehard GUI

Also see pLab, described above at Web Sites for Random Number Generators , for descriptions of some tests.

Commercial RNGs and Hardware RNGs

Robert Davies' page has a comparison of hardware and software RNGs.

Protego Hardware Generator is a hardware RNG.

Function Evaluation

If you'd like to write a paragraph for the FAQ on function evaluation, please contact me ...

Finding Roots

This is a standard topic in nearly all general texts on NA, and many routines exist in netlib (see Netlib ) and other libraries.

See also Alan Miller's page on an algorithm for finding the roots of polynomial equations. It is a translation into Fortran 90 of the version of Aberth's method by D. Bini in Numerical Algorithms, v.13 (1996), 179-200. The Fortran 77 code for this algorithm can be obtained from the numeralgo directory of netlib (see Netlib ).

Curve Fitting, Data Modelling, Interpolation, Extrapolation

This vague heading covers a multitude of closely related endeavors. There are many approaches, described in the following sections.

Commercial packages:

Other approaches to data modelling are:

Interpolation

For one dimensional work, see Stoer & Bulirsch 1980 or Press et al 1992. No doubt some of the techniques described there can be found in netlib (see Netlib ) by using the NIST guide (see Indices of NA Software on the Net ).

For greater than or equal one dimensions, some references are:

Dave Watson's home page: Dave Watson's Software

Watson, Dave, 1992, Contouring: A guide to the analysis and display of spatial data, Pergamon Press, ISBN 0 08 040286 0. Text includes a floppy with source code in Basic (sigh). Subject areas: isoline maps, isochor maps, methodology, data sorting, subset selection, local coordinates, binary local coordinates, barycentric coordinates, rectangular local coordinates, natural neighbour coordinates, gradient estimation, least squares gradients, least squares quadratics, spline gradients, cross product gradients, neighbourhood-based gradients, interpolation, inverse distance weighting, inverse distance weighted gradients, minimum curvature splines, neighbourhood-based interpolation, selected references to contouring literature

[Dave Watson] writes: Linear triangle-based interpolation is quick and dirty - probably the best compromise available if you are in a hurry. Natural neighbor interpolation is much slower but does provide a result that most people consider to be closest to that which an experienced manual draftsman would generate. The 'kriging' method is a variation on the splines technique - the main difference being that a subjectively selected probabilistic function is used rather than one of the analytically derived radial basis functions. This class of method generates discontinuities at points where the interpolation subset changes, and kriging attempts to ameliorate this effect be recomputing the coefficients for each interpolation point - this makes it very slow.

Interpolation on spheres: Dave Rusin's FAQ See also the comp.graphics.algorithms FAQ and the graphics FAQ.

Software:

[Albrecht Preusser] wrote: You could take routines from ACM Alg. 624 "Triangulation and Interpolation of Arbitrarily Distributed Points in the Plane" by Robert J. Renka. This Algorithm is available from netlib (see Netlib ) in the toms directory. Another algorithm in FORTRAN is contained in Akima's ACM Alg. 526, it is faster, but takes much more memory. By the way, I have developed methods based on the Delaunay Triangulation from both algorithms for computing nonlinear contours and nonlinear interpolation, which are published as ACM Alg. 626 and 684.

Spatial and Geometric Analysis toolbox (SaGA) - a MATLAB package for various aspects of spatial data interpolation and analysis and geometric modeling. It is available at: Glenn Kirill OAX is an spacial objective analysis program. OAX applies the method of optimal interpolation (objective analysis) to estimate the values of variables at specified points in a multidimensional space. For each point, the computation is based on a weighted average of the values of a specified number of data points (i.e. "nearest neighbours") which are closest to that point. For the sake of efficiency, several dependent variables can be interpolated simultaneously using the same weights if the underlying statistical model and relative noise level in the input data are assumed to be same for each dependent variable.

In this version, OAX extends the capability of earlier implementations (OPTIMAL_ESTIMATOR directive, oax3.0, oax4.0) by allowing the parameters controlling the calculation of the weights to vary as a function of the independent variables, and by allowing the specification of a relative uncertainty (NOISE) for each input data record separately.

Statistical Curve Fitting

See also: Probability and Statistics .
See also: Levenberg-Marquardt algorithm

Books

Jobson, J.D., Applied Multivariate Data Analysis, volumes 1 and 2, Springer-Verlag, New York. 1991.

Johnson, Richard A.; Wichern, Dean W, Applied Multivariate Statistical Analysis, Third Edition. [Ramin Samadani] writes: I like the book by Lawson and Hanson, solving least squares problems. I think its supposed to be back in print now or soon. Lawson, Charles L.; Hanson, Richard J. Prentice-Hall, Englewood Cliffs, N.J. 1974

Press et al 1991 also has a description of the field.

Spaeth, Helmuth: "Mathematical Algorithms for Linear Regression" 1987, Academic Press.

Software

ODRPACK, from Netlib , is a portable collection of Fortran subprograms for fitting a model to data. It is designed primarily for instances when the explanatory as well as the response variables have significant errors, implementing a highly efficient algorithm for solving the weighted orthogonal distance regression problem, i.e., for minimizing the sum of the squares of the weighted orthogonal distances between each data point and the curve described by the model equation.

Levenberg-Marquardt algorithm

[Alan Miller]: The non-linear least-squares problem is to find a vector of parameter values, a, to minimize the sum of squares of differences between a given vector, y, of observed values and a vector of fitted values f(a,x) where the values of x are known. e.g. a particularly notorious problem is to fit:

y = a(0) + a(1).exp(-a(2).x) + a(3).exp(-a(4).x)

sometimes with more exponential terms, to a set of (x,y)-pairs.

Writing S = sum [y(i) - f(a,x(i))]^2, where i denotes the i-th observation pair, we have that the first and second derivatives of S with respect to the parameters, a(j), are:

S' = -2.sum [r(i).(df/da(j))]
S'' = 2.sum [(df/da(j)).(df/da(k))] - 2.sum [r(i).d2f/(da(j)da(k))]

where r(i) is the i-th residual = y(i) - f(a,x(i)).

The Gauss-Newton (G-N) method uses the linearization of f about the latest estimates of a, and solves a linear least sqaures problem to obtain a new vector a.

If J denotes the Jacobian = the n x k matrix of values of (df/da(j)) at the current estimates of the parameter values, where n is the number of observations and k is the number of parameters, then the G-N method ignores the second derivatives of f in S'' and uses

a(new) = a(old) + inv(J'J).J'r

This guarantees a downhill search direction but often overshoots unless good starting values are available.

The Levenberg-Marquardt (L-M) method uses a trust-region approach, also known as ridge regression or regularization, to shrink the step sizes to ensure a reduction in S at each iteration. In the L-M method

a(new) = a(old) + inv(J'J + D).J'r

where D is a diagonal matrix. Levenberg suggested that D = d.I, that is that the same scalar value be added to each diagonal element. In practice, the diagonal elements of J'J are often of vastly different magnitudes, and Marquardt suggested

D = lambda.diag(J'J).

For the first iteration, lambda is given a fairly large value (say 1.0). If the first attempt at an iteration reduces S then lambda is reduced for the next iteration by a factor of say 5-10. If the first step increases S then lambda is increased by a factor of perhaps 3-5 until S is reduced. The final value of lambda for the iteration is then used for the next iteration.

The L-M method converges upon a stationary point (or sometimes a sub-space) which may be the minimum which is being sought, but it could be a saddle-point, a local minimum, or one or more of the parameter values may tend to infinity.

The Jacobian of first derivatives is often ill-conditioned so that QR factorizations are usually employed to solve the LS equations.

The L-M algorithm is often very efficient. Methods which use second derivatives require substantially more computation and then it is common to find that the Hessian matrix is not positive definite very close to the minimum.

Software:

MINPACK includes L-M algorithms using either analytic derivatives or derivatives estimated from differences. Obtain from the minpack directory of netlib (see Netlib ).

For many problems, the function f is linear in some of the parameters. In the sum of exponentials example mentioned, f is linear in a(0), a(1) and a(3). The VARPRO code uses an L-M type of algorithm but separates the linear from the non-linear parameters. This code is available from the opt directory of netlib (see Netlib ).

The package NL2SOL (TOMS algorithm 573) uses a hybrid combination of G-N and estimation of a positive-definite approximation to the Hessian from differences of first derivatives. It is available from the toms directory of netlib (see Netlib ). (Fortran 90 version from Alan Miller ).

Updates which allow for bounds on parameter values can be found in netlib directory port with file names starting n2f, n2g, rn2f or rn2g. (See Netlib ).

There is a very easy-to-use shareware package for PCs including graphics (NLREG) which is based upon the NL2SOL package, and includes graphics. It is more successful than many of the non-linear least-squares solvers found in commercial statistical packages. It can be downloaded from the math directory of simtelnet.

Time series analysis

Time series analysis may be used for some types of data. See:

Brockwell, Peter J.; Davis, Richard A., 1991, Time Series: Theory and Methods, Second Edition, Springer-Verlag, New York. ISBN 0-387-97429-6. [SJS]: Beautiful approach with emphasis on theory. Emphasis on the linear model.

Hamilton, James D., 1994, Time Series Analysis, Princeton U. Press, Princeton, NJ. ISBN 0-691-04289-6. [SJS]: Applied approach, with heavy emphasis on applications from economics. Not having an economics background is a detriment in reading this text. Broad mix of models.

Montgomery, Douglas C. / Forecasting and time series analysis (1990)

Pankratz, Alan / Forecasting with Univariate Time Series Models: Concepts and Cases / (1983) / good, clear

Vandaele, Walter / Applied Time Series and Box-Jenkins Models (1983). Good; pedagogical and practical

ARfit is a Matlab package for the estimation and spectral decomposition of multivariate autoregressive models.

ARfit is a collection of Matlab routines for

ARfit may be freely distributed. Suggestions for improvement are greatly appreciated. M-files, documentation, and a paper describing the algorithms used are available.

Also see Statlib .

Also see Probability and Statistics .

Splines

Useful for modelling single and multiple dimensional data.

Books:

de Boor, Carl, 1978, A Practical Guide to Splines, Springer-Verlag, New York, ISBN 0-540-90356-9. [Tim Strotman]: Provides an excellent blending of theory and practical implementation.

Dierckx, Paul, 1993, Curve and surface fitting with splines, Clarendon Press, ISBN 0198534418

Farin, Gerald E., 1995, NURB curves and surfaces : from projective geometry to practical use, A.K. Peters, ISBN 1568810385

Farin, Gerald E., 1993, Curves and surfaces for computer aided geometric design : a practical guide, third edition, Academic Press

Nurnberger, Gunther, 1989, Approximation by Spline Functions, Springer-Verlag.

Piegl, Les; Tiller, Wayne, 1995, The NURBS Book, Springer-Verlag, New York, ISBN 0-387-55069-0. [Tim Strotman]: Provides an excellent blending of theory and practical implementation.

Spaeth, Helmuth, 1995, One Dimensional Spline Interpolation Algorithms, A.K. Peters, ISBN 1-56881-016-4

Spaeth, Helmuth, and A.K. Peters, 1995, Two Dimensional Spline Interpolation Algorithms, A.K. Peters, ISBN 1-56881-017-2

Wahba, Grace, 1990, Spline Models for Observational Data

Software:

See Indices of NA Software on the Net .

Nearest Neighbor algorithms

Books:

Dasarathy, Belur V., 1991, Nearest Neighbor (NN) Norms: NN Pattern Classification Techniques, IEEE Computer Society Press, Los Alamitos, CA. [SJS]: Excellent.

Fitting a circle to a set of points

Problem: Given N Points:

(x1,y1); (x2,y2).. (xN,yN),

we want to find for best-fit circle:

(X0,Y0), R.

(Note: for fitting an *ellipse*, substitute the equation for an ellipse for the equation for a circle in the "Brute Force Approach").

Brute Force Approach (leads to a non-linear system): [Amara Graps]

Idea: Minimize by least squares the root-mean-squared error of the radius in the equation for a circle. In this method, one minimizes the sum of the squares of the *length* differences, which gives you a non-linear problem with all the associated pitfalls.

(x-X0)^2 + (y-Y0)^2 = R^2 Equation for a circle
R = SQRT[ (x-X0)^2 + (y-Y0)^2) ] Radius of the circle

where:

(X0,Y0) = center of circle
(x,y) = point coordinates
R = radius

1) Get first estimate for (X0,Y0,R).

(Details: Find some points to make first estimates- either solve the circle exactly (3 equations, 3 unknowns) to get a first estimate of the center and radius, or just do a kind of centroid calculation on all points- to get a rough estimate of the center and radius.)

2) Calculate r (r1, r2,.. rN) for each of your N points from the equation for a radius of a circle.

3) Calculate the root-mean-squared error For example, for 5 given points on the circle:

RMS error = SQRT[ [ (r1-R)^2 + (r2-R)^2 + (r3-R)^2 + (r4-R)^2
+ (r5-R)^2] / 3 ]

(dividing by "3" because we have 3 free parameters.)

4) Change (X0,Y0,R) slightly in your minimization algorithm to try a new (X0,Y0,R).

Minimization details:

Because minimization algorithms can get very computationally intensive, if one's circle-fitting problem is simple, I would look for a "canned" minimization routine. Some commercial computer programs for plotting and spreadsheets do this sort of thing. For example, the Excel spreadsheet has the built-in "Solver" that will perform minimzation.

Other possibilties for minimization software:

Matlab with the optimization toolbox.

IDL with Markwardt's MPFIT or MACSYMA.

ODRPACK from Netlib (see Netlib ) is an excellent library for a production mode approach to fitting these data, with a graphical user interface at Plato Guide At the very beginning of the User's Guide there is a description of how to fit data to an ellipse or circle. [Don Wells, Hans D. Mittelmann]. See Plato: Least Squares

Gaussfit is a superb response to the whole family of problems of adjusting parameters to fit functions to data; it is especially effective and easy to use for one-time problems and for exploratory work [Don Wells].

The L-BFGS-B Nonlinear Optimization package allows one to specify bounds on the variables, from: L-BFGS-B [Helmut Jarausch].

The Fortran programs fcircle.f, fellipse.f and the Lawson & Hanson least-squares routines ls.f shows how to implement these least-squares problems at: fcircle.f, fellipse.f, etc.

Alan Miller has updated the above for Fortran 90. The least-squares software is at: Alan Miller's Sofware with a working example at: Alan Miller's Example

Some Web Guides [Amara Graps]

GraphPad Guide to Nonlinear Regression

User's Manual -- GaussFit: {A} system for least squares and robust estimation", William H. Jefferys and Michael J. Fitzpatrick and Barbara E. McArthur and James E. McCartney", [Don Wells]

Other Papers

I. Kasa, "A Curve Fitting Procedure and Its Error Analysis", IEEE Trans. on Inst. Meas., 3/96, describes an algorithm that is a modified least squares analysis that uses the difference of the square of the radii rather than the absolute difference. [Rod Loewen]

M. Pilu, A. Fitzgibbon, R.Fisher, "Ellipse-specific Direct least-square Fitting", IEEE International Conference on Image Processing, Lausanne, September 1996. Pile et al A Java demo is at: Java Demo [Amara Graps]

Moura, L. and Kitney, R. I. (1991). A direct method for least -squares circle fitting. Computer Physics Communications 64, 57-63. [Colin B. Desmond]

"A Buyers Guide to Conic Fitting", A.W. Fitzgibbon, British Machine Vision Conference, 1995.

Cunningharn, Robert W.: Comparison of three methods for determining fit parameter uncertainties for the Marquardt Compromise, Computers in Physics 7(5), 570, (1993) [Amara Graps]

W. Gander, G. H. Golub and R. Strebel Least Squares Fitting of Circles and Ellipses BIT vol.34, 1994, pgs. 558-578 [Suhail A Islam].

Berman & Somlo: "Efficient procedures for fitting circles..." IEEE Instrum & Meast., IM-35, no.1, pp. 31-35, 1986. [Peter Somlo]

Paul T. Boggs and Richard H. Byrd and Robert B. Schnabel, "A Stable and efficient algorithm for nonlinear orthogonal distance regression", SIAM Journal on Scientific and Statistical Computing, 1987, 8, 6,1052-1078 [Don Wells].

P. T. Boggs and J. R. Donaldson and R. H. Byrd and R. B. Snabel, "Algorithm 676, {ODRPACK} -- Software for Weighted orthogonal distance regression", ACM Transactions On Mathematical Software, 1989, 15, 348-364 [Don Wells].

5) Calculate r (r1, r2 etc.) again from new (X0,Y0,R) above.

6) Calculate RMS error again.

7) Compare current RMS error with previous RMS error. If it doesn't vary by more some small amount, say 10^{-3} then you're done, otherwise continue steps 4 -- 7.

Other (more elegant) approaches that reduce to a linear system.

Linear solution

If you choose to minimize the squares of the *area* differences, you get a linear problem, which is a much safer way. [Pawel Gora, Zdislav V. Kovarik, Daniel Pfenniger, Condensed by Amara Graps]

1. Function to minimize is (sum of the area differences):

Q = Sum { [ (xi - X0)^2 + (yi -Y0)^2 - R^2 ]^2 , i=1..N }

2. A necessary condition is the system of 3 equations with 3 unknowns X0, Y0, R. Calculate the partial derivatives of Q, with respect to X0, Y0, R. (all d's are partial)

dQ / dX0 = 0
dQ / dY0 = 0
dQ / dR  = 0

3. Developing we get the linear least-squares problem:

| x1 y1 1 | | a |    | -x1^2-y1^2 |
| x2 y2 1 | | b | =~ | -x2^2-y2^2 |
| x3 y3 1 | | c |    | -x3^2-y3^2 |
| x4 y4 1 |          | -x4^2-y4^2 |
| x5 y5 1 |          | -x5^2-y5^2 |
 . . . . .   . .     . . . . . . .

(for example, for 5 points)

where

a = -2 X0; b = -2 Y0 ; c = X0^2 + Y0^2 - R^2.

Take any good least-squares algorithm to solve it, yielding a,b,c. So the final circle solution will be given with

X0 = -a/2; Y0 = -b/2; R^2 = X0^2+Y0^2 - c.

By the way, with 5 points you can also find the unique quadratic form (ellipse, hyperbola) which passes exactly through 5 points. With more than 5 points one can do a linear least-squares as well. The problem is then to minimize:

| x1^2-y1^2 x1y1 x1 y1 1 | | a |    | -x1^2-y1^2 |
| x2^2-y2^2 x2y2 x2 y2 1 | | b | =~ | -x2^2-y2^2 |
| x3^2-y3^2 x3y3 x3 y3 1 | | c |    | -x3^2-y3^2 |
| x4^2-y4^2 x4y4 x4 y4 1 | | e |    | -x4^2-y4^2 |
| x5^2-y5^2 x5y5 x5 y5 1 | | f |    | -x5^2-y5^2 |
..... .........

The robust or L1 or least-first-power approximation [Zdislav V. Kovarik].

If you try to minimize

W