[comp.archives] [sci.math.num-analysis] Linpack in ANSI C

schell@llandru.ucdavis.edu (Stephan Schell) (09/27/90)

Archive-name: c-linpack/27-Sep-90
Original-posting-by: schell@llandru.ucdavis.edu (Stephan Schell)
Original-subject: Linpack in ANSI C (some of it is available)
Archive-site: schell@iris.ucdavis.edu
Reposted-by: emv@math.lsa.umich.edu (Edward Vielmetti)

About a month ago I queried the net for its interest in ANSI C versions
of some Linpack routines, and I received a few replies.  I have
finished a very small subset of the routines (for inverting complex matrices,
solving complex linear systems, and computing the singular value decomposition
of a complex matrix) and am willing to distribute the source to whomever wants
it.  Since I don't know how much interest there is in this, I'll start off
by just mailing it out on a piecemeal basis, but will migrate to
anonymous ftp or posting to a newsgroup if there is sufficient interest.

I can be reached at
   schell@llandru.ucdvais.edu
   schell@iris.ucdavis.edu
   ...ucbvax!ucdavis!iris!schell

I have attached the README file that will go along with the source code.
------------------------------------------------------------------------
About c_linpack
---------------
c_linpack contains ANSI C versions of a subset of the LINPACK linear
algebra routines.  Currently, only the single-precision complex versions
of the matrix inverse (cgedi), linear equation solver (cgesl),
and singular value decomposition (csvdc) routines are
included (see list below).  The original Fortran routines were obtained
from netlib@ornl.gov (for more info about that service, send an email
message containing "send index from netlib").  The ANSI C code was
hand-converted and tested, although it is not guaranteed that no flaws
or errors exist.  The simple testing code is included with c_linpack.
If you do detect errors, bugs, etc., please report them to 

                  schell@llandru.ucdavis.edu

If I get around to converting other routines, I will post them, too.
Likewise, if you convert other LINPACK routines from Fortran to ANSI C,
please post them to the net so that we all don't have reinvent this
particular wheel.


DISCLAIMER
----------
The contents of c_linpack are distributed free of charge and are in no
way claimed or guaranteed to work properly by Stephan V. Schell nor by
the Regents of the University of California, nor by anyone else.  No
fault is assumed by said people for any consequences of the use of the
routines in c_linpack.


About matrices in C
-------------------
After trying several different conventions for all of this junk, I
settled on the one described here.  Since matrices in a C function can't
have dimensions passed as arguments to that function (unlike Fortran),
complex matrices in c_linpack are declared to be of type COMPLEX **,
where COMPLEX is the usual struct {float r,i;}.  A supplied routine,
COMPLEX **cmatrix(), allocates space for a vector of pointers, one pointer
per column, and then allocates space for the matrix elements themselves.
Note, then, that the elements in matrix a are addressed as
a[column][row].  This "reverse" ordering allows the structure of the
original Fortran code to be retained (in the sense that the structure is
optimized with respect to the storage ordering of the matrix elements),
easing conversion and validation.  Furthermore, the c_linpack routines 
all assume that matrices and vectors are indexed beginning from 1 
instead of 0.  The details of this indexing are handled by cmatrix() 
(and a companion routine cvector()).  As a result of this convention for
matrices, some of the ANSI C functions have fewer arguments than their
Fortran counterparts, because the leading dimension of an array is not
required (e.g., arguments ldx, ldu, and ldv to csvdc are omitted).


Execution Time
--------------
No comprehensive testing was done.  However, on an HP 9000/835 running
HP-UX 7.0, the c_linpack versions ran slightly faster than the original
Fortran counterparts (both were compiled with -O).  This may either 
support my coding conventions, bode well for HP's C compiler, or bode 
ill for HP's Fortran compiler.  Take your pick.  If anyone tests c_linpack 
against linpack on other machines, I would like to know the results.


Documentation
-------------
The comments (and original authorship) of the original Fortran LINPACK 
code are retained in the modules of c_linpack, although the notation
in the comments has been changed to be consistent with C.  These comments
may suffice for documentation, although you may wish to obtain the useful
references for the original Fortran codes:

J. Dongarra, J.R. Bunch, C.B. Moler, and G.W. Stewart (1978)
_LINPACK User's Guide_
Society for Industrial and Applied Mathematics
Philadelphia, PA

T.F. Coleman and C. Van Loan (1988)
_Handbook for Matrix Computations_
SIAM



List of routines
----------------
BLAS Level 1 routines and some extras:
caxpy.c        complex constant a times vector x plus vector y
ccopy.c        copy vector x to vector y
cdotc.c        dot product of vector x and vector y
cdotu.c        dot product of vector x (unconjugated) and vector y
cmach.c        some code for determining numerical constants of a machine
cmatrix.c      allocate a complex matrix having prescribed column and row 
               index ranges
crotg.c        compute Givens rotation
cscal.c        complex constant a times a vector
csrot.c        apply a Givens rotation
csscal.c       float constant a times a vector
cswap.c        swap two complex vectors
cvector.c      allocate a complex vector having prescribed index range
icamax.c       return the index of the complex vector element having largest
               absolute value
ivector.c      allocate an integer vector having prescribed index range
scasum.c       return the sum of the absolute values of the elements of a 
               complex vector
scnrm2.c       return the unitary norm of a complex vector
snrm2.c        return the Euclidean norm of a float vector
srotg.c        compute a Givens rotation

LINPACK routines:
cgeco.c        factor a matrix by Gaussian elimination and compute condition
cgedi.c        compute inverse and/or determinant of a complex matrix
cgefa.c        factor a matrix by Gaussian elimination
cgesl.c        solve a complex linear system of equations
csvdc.c        compute singular value decomposition of a complex matrix

-----------------------------------------------------------------------
--
-------------------------------------------------------------------------------
Stephan Schell                          schell@llandru.ucdavis.edu
Dept. of Electrical Engineering         {ucbvax,lll-crg}!ucdavis!llandru!schell
      &  Compter Science				
University of California, Davis         (916) 752-1326