Meschach Frequently Asked Questions Version 0.1 2nd March 1994 This is a list of frequently asked questions regarding Meschach. It is maintained by David Stewart and Zbigniew Leyk, the developers and maintainers of Meschach. A copy of this is kept in the anonymous ftp area thrain.anu.edu.au : /pub/meschach/FAQ.meschach ---------------------------------------------------------------------- Q1. What is Meschach? Q2. Why not use C versions of LINPACK, EISPACK or LAPACK? Q2a. What are the differences of Meschach from LAPACK or EISPACK in C? Q2b. How reliable are the routines compared to the above packages? Q3. How do I get documentation? Q4. I've found a bug in Meschach! What do I do? Q5. It won't compile on my BrandX machine running X-OS. What do I do? Q6. Is there a C++ version of Meschach? Q7. Meschach doesn't have an xxxx routine. Why not? What do I do? Q8. Can I use Meschach with MATLAB? Q9. I don't want ALL of Meschach. How do I get just a subset? Q10. What does Meschach mean? Q11. My code that uses Meschach seems to run very slowly/uses incredible amounts of memory. Why? Q12. I have a huge problem and I keep running out of memory. What do I do? Q13. How did you come to write Meschach? Q14. Complex part of Meschach Q15. Can I do an index search in order to know what routines exist? ---------------------------------------------------------------------- A1. What is Meschach? Meschach is a library of functions written in C for performing matrix computations. The operations can be (briefly) summarised as * linear combinations of vectors and matrices; inner products; matrix-vector and matrix-matrix products * solving linear systems of equations * solving least-squares problems * computing eigenvalues (or singular values) * finding "nice" bases for subspaces like ranges of matrices and null spaces. There are, however, a lot of different ways of, for example, solving a system of linear equations depending on whether the matrix is real or complex, symmetric, positive definite or indefinite, sparse or has some other structure. The rule in computational mathematics is: Exploit Structure. So there is the standard Gaussian Elimination method for general real (and general complex) matrices, as well as Cholesky factorisation (for positive definite symmetric matrices), which is about twice as fast, LDL factorisation which just needs symmetric matrices, but is not always well-behaved, BKP (Bunch-Kaufman-Parlett) factorisation for indefinite symmetric matrices which is well-behaved but more complex, QR factorisation which is best used for least squares problems. And then there are sparse matrices, and a variety of different techniques for different types of sparse matrices. A list of the features of Meschach follows. o dynamic allocation, de-allocation, re-sizing and copying of objects o dense complex matrices and vectors as well as real matrices and vectors o input and output routines for these objects, and MATLAB save/load format o error/exception handling o basic numerical linear algebra -- linear combinations, inner products, matrix-vector and matrix-matrix products including transposed and adjoint forms o vector min, max, sorting, componentwise products, quotients o dense matrix factorise and solve -- LU, Cholesky, LDL^T, QR, QR with column pivoting, symmetric indefinite (BKP) o dense matrix factorisation update routines -- LDL^T, QR (real matrix updates only) o eigenvector/eigenvalue routines -- symmetric, real Schur decomposition, SVD, extract eigenvector o sparse matrix "utility" routines o sparse matrix factorise and solve -- Cholesky, LU and BKP (Bunch-Kaufman-Parlett symmetric indefinite factorisation) o sparse incomplete factorisations -- Cholesky and LU o iterative techniques -- pre-conditioned conjugate gradients, CGNE, LSQR, CGS, GMRES, MGCR, Lanczos, Arnoldi o allowance for "procedurally defined" matrices in the iterative techniques o various "torture" routines for checking aspects of Meschach o memory tracking for locating memory leaks Here is a sample piece of code using Meschach to read in a matrix A and a vector b and to solve A.x = b for x and print the result. MAT *A, *Acopy; VEC *b, *x; PERM *pivot; A = m_input(MNULL); /* read in matrix from stdin */ Acopy = m_copy(A,MNULL); /* create & copy matrix */ b = v_input(VNULL); /* read in vector from stdin */ pivot = px_get(A->m); /* get pivot permutation -- size = # rows of A */ LUfactor(A,pivot); /* LU factorisation */ x = LUsolve(A,pivot,b,VNULL); /* ... and solve A.x = b; result is in x */ /* .... and print out solution with a comment */ printf("# x =\n"); v_output(x); Meschach aims to be an easy-to-use, easy-to-debug-your-code library. The basic data objects (vectors, matrices etc.) are represented by self-contained data structures which can be created, destroyed and resized at will. Workspace arrays do not need to be passed to routines (except for some internal routines). Operations are continually checked to see if they are valid so that errors are quickly located. Unlike matrix interpreter/calculators such as MATLAB, the user is given control over memory usage and can efficiently re-use memory which requires garbage collection in a more interpretive environment. ---------------------------------------------------------------------- A2. Why not use C versions of LINPACK, EISPACK or LAPACK? LINPACK, EISPACK and LAPACK are well-known Fortran-language packages of routines for a variety of problems in matrix computations. With the magic of translators like f2c, why not use the translated versions of these libraries? Well, a number of people are already doing this; after all, LINPACK etc. are public domain and they have been incorporated in public domain software such as Octave and RlaB, to name but two. But, as with all translations from Fortran, these do not use any of the facilities that Fortran lacks. In particular, there is no attempt to structure the data, nor is there any use of dynamic memory allocation, C's error handling (using setjmp() and longjmp()). Finally there is no input/output. Meschach uses not only the data structuring facilities in C to provide self-contained data structures, but this is used with the error handling to ensure, for example, that the matrices passed to a matrix multiplication routine have compatible sizes. The data structuring is also used with the dynamic memory allocation to allow efficient resizing -- vectors are only physically reallocated if the requested size is above the current "high water mark". Input and output are both handled using C's I/O facilities to provide output that is both machine and human-readable. Binary MATLAB-compatible output (for dense matrices) is also done. And another useful aspect of Meschach's I/O is that it allows for comments in the data. As well, none of the above packages has anything to do with sparse matrices. There are, of course, sparse matrix packages that are freely available, but even there, they do not seem to handle, in an integrated way, both sparse and dense matrices. Meschach has not only sparse matrix data structures and the routines to make use of them, but as the sparse matrices are self-contained data structures, use of sparse and dense matrices can be intermixed. You may also like to look at answer A2a. ---------------------------------------------------------------------- A2a. What is the difference of meschach from LAPACK or EISPACK in C? First: The biggest difference is that Meschach is written from scratch in C, and works in terms of data structures that are self-contained mathematical objects (matrices, vectors, permutations, sparse matrices etc). These can be created, destroyed and re-sized at will. Second: LAPACK and EISPACK have no sparse matrix routines, and no iterative routines; Meschach has both and includes sparse factorisations, both with fill-in and "incomplete" factorisations. Third: The algorithms are similar; we have **NOT** taken LINPACK/EISPACK/LAPACK and re-implemented this in C. Fourth: Meschach makes comprehensive use of C's facilities (as well as structures a.k.a. records) such as error handling, and dynamic memory allocation. E.g. workspace arrays do not need to be passed. They are created and re-sized as necessary where they are used. With Meschach 1.2 you can destroy any of these local workspace vectors etc if memory is short. Fifth: Meschach is smaller in the amount of memory needed, both for source and executable than e.g. LAPACK. Full source code in shar files is under 1Mbyte for Meschach, more like 4-5Mbyte for LAPACK. There **are** a number of things in LAPACK that are not in Meschach (e.g. banded complex factorisations), but then again, there is a package of sparse matrix operations in Meschach that is not in LAPACK. We suspect that it is partly due to the way the software has been constructed. Meschach has been a 1 or 2 person effort; LAPACK has had something like 10 full-time programmers. We've had to be economical in terms of how the software has been developed. --------------------------------------------------------------------------- A2b. How reliable are the routines compared to the above packages? Well, LINPACK and EISPACK have been around for nearly 20 years. The bugs should be found by now! So we can't honestly claim that we have yet achieved **THAT** level of reliability. But we believe that it probably approaches that level of reliability for the more commonly used routines. Because there is a considerable amount of error/exception checking in Meschach, bugs have been relatively easy to track down. A number of bugs have been found in Meschach 1.2a (and there will be more). But these have been mostly minor and easily fixed. ------------------------------------------------------------------------ A3. How do I get documentation? There is some basic documentation that comes with the source code of Meschach and there will be more on that later. But the comprehensive documentation has been published as a book (both a users' and a reference manual) by Centre for Mathematics and its Applications (CMA) School of Mathematical Sciences Australian National University Canberra, ACT 0200 Australia as "Meschach: Matrix Computations in C", volume 32 in the "Proceedings of the CMA". The cost is about A$30 (approx. US$22) + postage/handling. Orders can be passed to the CMA via David Stewart (david.stewart@anu.edu.au). The manual is nearly 250pp long. If you just want to check out the sort of things that Meschach has before committing yourself to buying the manual, there is the manual for the **OLD** version of Meschach which is available over ftp at ftpmaths.anu.edu.au : /pub/meschach/version1.1b/bookdvi.tar [.Z or .gz] You can also have a look at some of the documentation that comes with the source code. In the DOC directory there is tutorial.txt which is an ASCII file containing a tutorial introduction to Meschach. There is also fnindex.txt which contains an index (in alphabetical order) of the functions in Meschach and a brief description of what they do. The calling sequences can be found from the header files and the actual source code. Of course, this is no substitute for a comprehensive manual, but it does give you a start with using Meschach. The torture.c and *tort.c files can also be useful in showing how Meschach can be used. ---------------------------------------------------------------------- A4. I've found a bug in Meschach! What do I do? Congratulations! Let us know! Not all bugs are serious --- most of the bugs found by users are installation problems. Most bugs in Meschach are very easily fixed. Please try to determine where the bug is occuring (which routine, preferably) and what data results in a problem. Send that information to us and we will aim to fix the problem as quickly as possible. Send the information to either of us. We will pass it on to the other. Sometimes the torture test(s) give a message "Error: ..." indicating that there may be something wrong. Please have a look at the size of the error that is printed out with the "Error: ..." message. If this is of the order of 1e-15 or 1e-14 it is unlikely that there is anything seriously wrong. The reason is that torture and the other testing programs generally use randomly generated problems, which may on occasion be more ill-conditioned than expected and give slightly larger errors than usual. If on the other hand, errors of size 1e-6 (say) or more for the direct methods **DOES** indicate a serious bug. ---------------------------------------------------------------------- A5. It won't compile on my BrandX machine running X-OS. What do I do? Firstly make sure that you run "configure" or "configure --with-..." to get the options you want. Then you use "make". The simplest installation procedure for Unix systems to create the entire library is configure --with-all make all If you want to use gcc (or if that is not the default compiler) then use configgnu --with-all make all The "configure" script is not perfect and sometimes misses some things. For example, when we used configure on Linux, it missed the fact that the type u_int was already declared. This meant that some editing of the generated machine.h file was needed. Users of MS-DOS machines and C compilers, and users of VMS and other non-Unix operating systems will need to manually edit your machine.h and/or makefile for your system. Some guidance in this can come from the examples in the MACHINES/xxxx directories which contain machine.h, makefile (and occasionally machine.c) files for some machines. If configure fails, then the files to modify for your particular machine are "machine.h" and "makefile". (The file "machine.c" can also be modified for speed.) There are a number of macros in machine.h which can be set to control the way the library is compiled. You should consult a local expert (and/or your local manuals!) to determine how things should be set. If this is not successful, send us a message. We are unlikely to know the magic words for your particular machine/compiler/operating system, but we may know someone who does. ---------------------------------------------------------------------- A6. Is there a C++ version of Meschach? There is a C++ version (called Meschach++) which is under development. A beta version can be found in the directory MESCHPP. There is a manual for this version. Note that this manual will not be a substitute for the main manual described above -- the Meschach++ manual will refer to "Meschach: Matrix Computations in C" for many details on the behaviour of the various routines. The C++ version is, basically, a "wrapper" for Meschach. Meschach++ is not being developed by David Stewart and Zbigniew Leyk, but by Stephen Roberts, Brian Austin and Alex Austin also while at the Australian National University, School of Mathematical Sciences. You can contact Stephen Roberts at steve@laplace.anu.edu.au. Some points should be made about the use of C++. The use of operator overloading is one of the most attractive features of C++, especially for tasks such as matrix computations. However, it also has some drawbacks. One drawback is the problem of temporary variables. In a matrix expression of the form: S = A + B + C + D + E; where all the objects are matrices (of the same size) if operator overloading is implemented in the naive way, then this will result in temporary matrices for (A+B), ((A+B)+C), (((A+B)+C)+D) and ((((A+B)+C)+D)+E). In fact, the sum could be accumulated directly into S (provided there is no aliasing) and there do not need to be any temporaries. Instead, a naive implementation would spend most of its time allocating and deallocating temporary variables. This is just one aspect of a problem with optimising evaluation of matrix expressions. There are ways of fixing this problem using delayed evaluation, which can be implemented in C++. . Another is to use 3-term functions, similar in spirit to Meschach's use of a final output parameter. This loses the elegance of operator overloading, but gains greatly in performance. ---------------------------------------------------------------------- A7. Meschach doesn't have an xxxx routine. Why not? What do I do? There is usually one answer to "Why not?" We haven't had time. But some operations should be avoided if possible. It was a long time before we implemented matrix inversion as a single function. Most of the time it is not needed: often only the solution of systems of equations is needed. There is no single "eigenvalue/vector" routine. This is because this is done by computing the Schur decomposition of a matrix, and then extracting information about the eigenvalues and eigenvectors. The Schur decomposition not only provides more information, but in the case of eigenvectors, is much more numerically stable. After all, not all matrices can be transformed into diagonal form -- they have fewer eigenvectors than columns. The other point that should be made is that Meschach is a toolkit. Often it is up to you to make the part you need. You may need to re-consider what you are doing -- is it numerically stable? Is it an efficient way of doing the task? It may be better to build on what Meschach provides, rather than hoping that it will all be there. ---------------------------------------------------------------------- A8. Can I use Meschach with MATLAB? The answer is both "Yes" and "No". The "Yes" part is that Meschach provides input and output of dense matrices and vectors that is compatible with MATLAB. This means that matrices can be transferred between MATLAB and Meschach programs. However, Meschach has not been made into a MEX-compatible library. So Meschach routines are not **yet** ready to be used in MEX routines. ---------------------------------------------------------------------- A9. I don't want ALL of Meschach. How do I get just a subset? Yes you can. If you have the shar files from netlib then the most basic library can be built out of shar files meschach0.shar and meschach1.shar. This contains the core routines and the basic matrix, vector computations. (It is not even enough to get torture to work.) The next shar file needed is meschach2.shar which contains the matrix factorisation routines. With these three shar files the basic library can be built. The shar file meschach3.shar contains the sparse and iterative methods, while meschach4.shar contains the complex vector and matrix operations. If you want the basic library (with the matrix factorisations), get the shar files meschach[012].shar, expand them, and then build the library using configure make basic If you want to add just the sparse routines use configure --with-sparse make sparse If you want to add just the complex routines use configure --with-complex make complex ---------------------------------------------------------------------- A10. What does Meschach mean? Meschach doesn't stand for anything in particular. And it is NOT an acronym. It is, however, someone's name. If you want to find out, look two thirds down page ii of the manual :-) ---------------------------------------------------------------------- A11. My code that uses Meschach seems to run very slowly/uses incredible amounts of memory. Why? It is quite possible to make Meschach run very inefficiently. What it is most likely doing is spending all of its time allocating matrices. (If it is not also deallocating them **explicitly**, then you will use an amazing amount of memory.) Look at where you call Meschach. Is the output parameter NULL? If so, then a new data structure is created to hold the output whenever the routine is called. For example, MAT *A, *B, *C; for ( i = 0; i < 100; i++ ) { .... C = m_add(A,B,MNULL); /* C = A+B */ .... } will create a new C every time through the loop. To avoid this use C as the output parameter: MAT *A, *B, *C; C = m_get(A->m,A->n); for ( i = 0; i < 100; i++ ) { .... C = m_add(A,B,C); .... } See the pages under the index entry "memory management" for a comprehensive discussion of ways of dealing with memory allocation and deallocation. ---------------------------------------------------------------------- A12. I have a huge problem and I keep running out of memory. What do I do? Large problems often require specialised techniques. A large problem is one where just using dense matrices is intolerably slow. For most workstations and most problems this is somewhere around 1000 x 1000, possibly a bit less, possibly a bit more. The first strategy to use is to consider using sparse matrices. Direct methods can be used with sparse matrices, at least for solving systems of equations. The danger here is that fill-in will result in the sparse matrix becoming a bit too dense for the amount of memory available. Re-ordering techniques should be used (unfortunately Meschach does not currently implement these) to reduce the amount of fill-in created by the factorisation process. If this is not enough, then iterative methods should be used. The ones implemented in Meschach are based on "Krylov subspaces". These just require a function to compute the product A.x for any given vector x (or sometimes the transposed-matrix-vector product A^T.x). These do not require the matrix involved to be sparse; it can be dense, but as long as it has the right structure so that A.x can be computed quickly, these methods are appropriate. Preconditioning is desirable for these iterative methods, and can be achieved in a variety of ways. One is to use incomplete factorisations (which simply ignore fill-in). Other methods are problem dependent: for partial differential equations, multigrid preconditioners can work very well. Classical methods (Gauss-Seidel, SOR, SSOR, Chebyshev acceleration) can be used to create pre-conditioners which when combined with a Krylov subspace method can be far better than either alone. ---------------------------------------------------------------------- A13. How did you come to write Meschach? The first of us to start working on Meschach was David Stewart. The first ideas of how to use pointers and dynamic memory allocation for storing matrices was developing while completing an Electrical Engineering degree. The next year was spent as a systems programmer for 5 months and then starting an Honours year in Mathematics (in USA read "Masters", in Germany read "Diplom"). By 1986, when he started his PhD he was ready to start writing code for it. The library was never an end in itself, but a tool for writing the real programs for the thesis (in numerical analysis and differential equations). In 1989 the library was the basis for programs using linear programming and ODE solvers; however, it was still very limited and had no eigenvalue/vector routines, no complex routines, and was not particularly portable. For a while he was waiting for a job on a project at an operations research company on large scale optimisation -- the project eventually fell through, but in the mean time some basic sparse routines were written in anticipation. Later that year David Stewart was hired as a lecturer in Computational Mathematics at the University of Queensland, Australia, and taught Numerical Linear Algebra. As a result, the development of a number of routines took place in parallel with them being taught in the classroom. Meanwhile, Zbigniew Leyk had also completed his PhD and was working in a Post-doc at Cornell. While there he developed his own ideas about a numerical library in C as part of his work on multigrid methods. This had a different flavour to Meschach in some respects, but both approaches emphasised the use of data structures. By 1991 when both David Stewart and Zbigniew Leyk arrived at the Australian National University, the ideas had matured enough; Meschach had been developed in terms of scope and portability to become a really general tool. Meschach was placed in netlib in 1991. A manual had been written over 1990/1991. There was mutual interest in working with Meschach from both of us, and during 1993 we worked intensively on a major new revision, which has become version 1.2. This had complex data structures and operations, a revamped way of dealing with iterative methods, better methods for handling workspace vectors and matrices and a new installation procedure based on GNU's autoconf. ---------------------------------------------------------------------- A14. Complex part of Meschach >> I wrote for myself a file called complex.h >> and this is the one which is included >> automatically in your files. How can this work? >> What type of complex.h is required? complex.h is not for a C++ header file. Rather, the idea is that if a user have a better local complex.h (C) header file then that should be used. However, there may be problems with names being different (e.g. cadd() vs. zadd() -- Meschach uses the latter). It might be necessary to edit machine.h to make it work better for a particular setup. -------------------------------------------------------------------------- A15. Can I do an index search in order to know what routines exist? In meschach0.shar (and in the tar and zip distributions) is a file called fnindex.txt which contains a list of routines and briefly what they do. Use your favourite string search method to find the right routine(s). ------------------------------------------------------------------------ End of Meschach FAQ