```      SUBROUTINE XXSIM(A,IA,M,B,N,X,WKS,ERR,ISTOP)
C-----------------------------------------------------------------------
C
C ROUTINE: XXSIM
C
C PURPOSE: SOLVES THE SYSTEM OF SIMULTANEOUS EQUATIONS AX=B USING THE
C	   NETLIB LINALG ROUTINE LSQR WHICH FOLLOWS THIS ROUTINE ALONG
C          WITH ITS DEPENDENCIES. THIS ROUTINE REPLACES NAG LIBRARY
C          ROUTINE F04ATF. HOWEVER, THE LU DECOMPOSITION IS NOT OUTPUT.
C
C CALLING PROGRAM: GENERAL USE
C
C INPUT:
C
C     (R*8)   A(,)    THE MATRIX A
C     (I*4)   IA      THE FIRST DIMENSION OF THE TWO-DIMENSIONAL
C                     ARRAY A, IA<=M
C     (I*4)   M       NUMBER OF ROWS OF A
C     (R*8)   B()     RIGHT-HAND-SIDE VECTOR B, DIMENSION = M
C     (I*4)   N       NUMBER OF COLUMNS OF A
C     (R*8)   WKS()   WORKSPACE VECTOR, DIMENSION >= M*N+2*N
C
C OUTPUT:
C
C     (R*8)   X()     SOLUTION VECTOR X, DIMENSION = N
C     (R*8)   ERR()   VECTOR OF THE ERROR ESTIMATES OF THE COMPONENTS OF
C                     X. SEE VARIABLE 'SE' IN LSQR.
C     (I*4)   ISTOP   AN ERROR CODE SET TO 0 OR 4 IF NO ERROR. SEE LSQR
C                     FOR A DESCRIPTION OF THE ERROR CODES.
C
C CALLS:
C     (EXT.)  MA      SUBROUTINE TO CALCULATE PRODUCTS OF A WITH GIVEN
C                     VECTORS. GIVEN NEXT IN THE FILE.
C ROUTINES:
C-----------------------------------------------------------------------
C     NAME     SOURCE   PURPOSE
C-----------------------------------------------------------------------
C     LSQR     NETLIB   CALCULATES THE SOLUTION - SEE BELOW
C-----------------------------------------------------------------------
C
C AUTHOR: WILLIAM OSBORN (TESSELLA SUPPORT SERVICES PLC.)
C
C DATE:   10-06-96
C
C VERSION 1.1                             DATE: 10-06-96
C MODIFIED: WILLIAM OSBORN
C           - FIRST VERSION
C VERSION 1.2                             DATE: 13-09-96
C MODIFIED: WILLIAM OSBORN
C           - ADDED IA PARAMETER SO THAT NON-SQUARE MATRICES CAN BE USED
C
C VERSION 1.3                             DATE: 25-09-96
C MODIFIED: WILLIAM OSBORN
C           - CHANGED DIMENSION OF B TO M
C
C VERSION 1.4                             DATE: 14-06-2000
C MODIFIED: Martin O'Mullane
C           - Changed iwk from real*8 to integer.
C
C-----------------------------------------------------------------------

INTEGER             IA,          ISTOP,       M,           N
REAL*8              A(IA,N),     B(M),        ERR(N)
REAL*8              WKS(M*N+2*N),             X(N)
INTEGER             IW(LENIW),   LENIW,       LENRW,       M
INTEGER             MODE,        N
REAL*8              RW(LENRW),   X(N),        Y(M)
DOUBLE PRECISION    ACOND,       ANORM,       ARNORM,      ATOL
DOUBLE PRECISION    BTOL,        CONLIM,      DAMP,        RNORM
DOUBLE PRECISION    RW(LENRW),   SE(N),       U(M),        V(N)
DOUBLE PRECISION    W(N),        X(N),        XNORM
INTEGER             ISTOP,       ITN,         ITNLIM
INTEGER             IW(LENIW),   LENIW,       LENRW,       M
INTEGER             N,           NOUT
DOUBLE PRECISION    X(N),        Y(N)
INTEGER             INCX,        INCY,        N
DOUBLE PRECISION    X(N)
INTEGER             INCX,        N
DOUBLE PRECISION    A,           X(N)
INTEGER             INCX,        N
```