SUBROUTINE GELB

PURPOSE
   TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH A
   COEFFICIENT MATRIX OF BAND STRUCTURE.

USAGE
   CALL GELB(R,A,M,N,MUD,MLD,EPS,IER)

DESCRIPTION OF PARAMETERS
   R	  - M BY N RIGHT HAND SIDE MATRIX (DESTROYED).
	    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
   A	  - M BY M COEFFICIENT MATRIX WITH BAND STRUCTURE
	    (DESTROYED).
   M	  - THE NUMBER OF EQUATIONS IN THE SYSTEM.
   N	  - THE NUMBER OF RIGHT HAND SIDE VECTORS.
   MUD	  - THE NUMBER OF UPPER CODIAGONALS (THAT MEANS
	    CODIAGONALS ABOVE MAIN DIAGONAL).
   MLD	  - THE NUMBER OF LOWER CODIAGONALS (THAT MEANS
	    CODIAGONALS BELOW MAIN DIAGONAL).
   EPS	  - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
	    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
   IER	  - RESULTING ERROR PARAMETER CODED AS FOLLOWS
	    IER=0  - NO ERROR,
	    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
		     TERS M,MUD,MLD OR BECAUSE OF PIVOT ELEMENT
		     AT ANY ELIMINATION STEP EQUAL TO 0,
	    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
		     CANCE INDICATED AT ELIMINATION STEP K+1,
		     WHERE PIVOT ELEMENT WAS LESS THAN OR
		     EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
		     ABSOLUTELY GREATEST ELEMENT OF MATRIX A.

REMARKS
   BAND MATRIX A IS ASSUMED TO BE STORED ROWWISE IN THE FIRST
   ME SUCCESSIVE STORAGE LOCATIONS OF TOTALLY NEEDED MA
   STORAGE LOCATIONS, WHERE
     MA=M*MC-ML*(ML+1)/2    AND    ME=MA-MU*(MU+1)/2	WITH
     MC=MIN(M,1+MUD+MLD),  ML=MC-1-MLD,  MU=MC-1-MUD.
   RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
   IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN SOLUTION
   MATRIX R IS STORED COLUMNWISE TOO.
   INPUT PARAMETERS M, MUD, MLD SHOULD SATISFY THE FOLLOWING
   RESTRICTIONS     MUD NOT LESS THAN ZERO
		    MLD NOT LESS THAN ZERO
		    MUD+MLD NOT GREATER THAN 2*M-2.
   NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
   RESTRICTIONS ARE NOT SATISFIED.
   THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
   PARAMETERS ARE SATISFIED AND IF PIVOT ELEMENTS AT ALL
   ELIMINATION STEPS ARE DIFFERENT FROM 0. HOWEVER WARNING
   IER=K - IF GIVEN - INDICATES POSSIBLE LOSS OF SIGNIFICANCE.
   IN CASE OF A WELL SCALED MATRIX A AND APPROPRIATE TOLERANCE
   EPS, IER=K MAY BE INTERPRETED THAT MATRIX A HAS THE RANK K.
   NO WARNING IS GIVEN IF MATRIX A HAS NO LOWER CODIAGONAL.

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
   NONE

METHOD
   SOLUTION IS DONE BY MEANS OF GAUSS ELIMINATION WITH
   COLUMN PIVOTING ONLY, IN ORDER TO PRESERVE BAND STRUCTURE
   IN REMAINING COEFFICIENT MATRICES.