* ======================================================================
* NIST Guide to Available Math Software.
* Fullsource for module 599 from package TOMS.
* Retrieved from NETLIB on Wed Jun  4 20:21:26 1997.
* ======================================================================
C     ALGORITHM 599, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2,
C     JUN., 1983, P. 255-257.
C**********************************************************************
C**********************************************************************
C**********************************************************************
C                                                                      
C                                                                      
C                                                                      
C     F O R T R A N  SOFTWARE PACKAGE FOR RANDOM NUMBER GENERATION     
C                                                                      
C                                                                      
C                                                                      
C**********************************************************************
C**********************************************************************
C**********************************************************************
C                                                                      
C                                                                      
C                                                                      
C     CONTENTS:                                                        
C                                                                      
C     1) SUNIF  -  0,1 -UNIFORM DISTRIBUTION                           
C                                                                      
C     2) SEXPO  - (STANDARD-) EXPONENTIAL DISTRIBUTION                 
C                                                                      
C     3) SNORM  - (STANDARD-) NORMAL DISTRIBUTION                      
C                                                                      
C     4) SGAMMA - (STANDARD-) GAMMA DISTRIBUTION                       
C                                                                      
C     5) KPOISS - POISSON DISTRIBUTION                                 
C                                                                      
C                                                                      
C     THIS PACKAGE CONSTITUTES A FORTRAN-77 DOCUMENTATION OF A SET OF  
C     ASSEMBLER FUNCTIONS FOR SAMPLING FROM THE ABOVE DISTRIBUTIONS.   
C     ALL ROUTINES MAKE AMPLE USE OF BINARY REPRESENTATIONS OF NUMBERS,
C     THEY ARE AMONG THE MOST ACCURATE AND FAST SAMPLING FUNCTIONS     
C     KNOWN. THE FORTRAN PROGRAMS BELOW YIELD THE SAME RANDOM NUMBER   
C     SEQUENCES AS THE ONES FROM OUR ASSEMBLER PACKAGE, BUT THEY ARE   
C     OF COURSE MUCH SLOWER (BY FACTORS 5-8 ON OUR SIEMENS 7760        
C     COMPUTER.)                                                       
C     THE SET OF ROUTINES WILL ALSO BE ACCEPTABLE TO FORTRAN IV        
C     COMPILERS WHICH ALLOW DATA STATEMENTS FOR ARRAYS WITHOUT         
C     IMPLICIT DO-LOOPS.                                               
C                                                                      
C                                                                      
C     REMARKS:                                                         
C                                                                      
C     -  NO CARE IS TAKEN TO ENSURE THAT THE PARAMETER VALUES LIE      
C        IN THE ALLOWED RANGE (E.G. A/MU > 0.0 FOR SGAMMA/KPOISS).     
C                                                                      
C     -  THE PARAMETER 'IR' MUST BE SET TO SOME  4*K+1 > 0  BEFORE     
C        THE FIRST CALL OF ANY OF THE GENERATORS. THEREAFTER IR        
C        MUST NOT BE ALTERED UNTIL A NEW INITIALIZATION IS DESIRED.    
C                                                                      
C     -  THE PACKAGE PROVIDES RANDOM DEVIATES OF 6-7 DIGITS ACCURACY.  
C        ON MORE ACCURATE COMPUTERS THE CONSTANTS IN SEXPO, SNORM,     
C        SGAMMA AND KPOISS OUGHT TO BE ADJUSTED ACCORDING TO LOCAL     
C        COMMENTS OR WITH THE AID OF THE TABLES IN THE LITERATURE      
C        QUOTED AT THE BEGINNING OF EACH FUNCTION.                     
C                                                                      
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                     
C                                                                      
C       0 , 1   - U N I F O R M  DISTRIBUTION                          
C                                                                      
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C     FOR DETAILS SEE:                                                 
C                                                                      
C               AHRENS, J.H., DIETER, U. AND GRUBE, A.                
C               PSEUDO-RANDOM NUMBERS:  A NEW PROPOSAL                
C                     FOR THE CHOICE OF MULTIPLICATORS                 
C               COMPUTING, 6 (1970), 121 - 138                         
C                                                                      
C**********************************************************************
C                                                                      
      REAL FUNCTION SUNIF(IR)                                          
      DOUBLE PRECISION R,FACTOR,TWO28
      INTEGER IR
C
C     FACTOR - INTEGER OF THE FORM 8*K+5 AS CLOSE AS POSSIBLE
C              TO  2**26 * (SQRT(5)-1)/2     (GOLDEN SECTION)
C     TWO28  = 2**28  (I.E. 28 SIGNIFICANT BITS FOR DEVIATES)
C
      DATA FACTOR /41475557.0D0/, TWO28 /268435456.0D0/
C
C     RETURNS SAMPLE U FROM THE  0,1 -UNIFORM DISTRIBUTION
C     BY A MULTIPLICATIVE CONGRUENTIAL GENERATOR OF THE FORM
C        R := R * FACTOR (MOD 1) .
C     IN THE FIRST CALL R IS INITIALIZED TO
C        R := IR / 2**28 ,
C     WHERE IR MUST BE OF THE FORM  IR = 4*K+1.
C     THEN R ASSUMES ALL VALUES  0 < (4*K+1)/2**28 < 1 DURING
C     A FULL PERIOD 2**26 OF SUNIF.
C     THE PARAMETER IR IS USED ONLY IN THE FIRST CALL FOR
C     INITIALIZATION OF SUNIF. THEREAFTER (WHEN NEGATIVE)
C     IR BECOMES A DUMMY VARIABLE.
C
      IF (IR .GE. 0) GO TO 1
C
C     STANDARD CASE:  SAMPLING
C
      R=DMOD(R*FACTOR,1.0D0)
      SUNIF=SNGL(R)
      RETURN
C
C     FIRST CALL: INITIALIZATION
C
1     R=DBLE(FLOAT(IR))/TWO28
      R=DMOD(R*FACTOR,1.0D0)
      SUNIF=SNGL(R)
      IR=-1
      RETURN
      END
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C                                                                     
C     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                
C                                                                      
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C     FOR DETAILS SEE:                                                 
C                                                                      
C               AHRENS, J.H. AND DIETER, U.                            
C               COMPUTER METHODS FOR SAMPLING FROM THE                
C               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  
C               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               
C                                                                     
C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       
C     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)      
C                                                                      
C**********************************************************************
C                                                                      
      REAL FUNCTION SEXPO(IR)                                          
      REAL Q(8)
      EQUIVALENCE (Q(1),Q1)
      REAL Q1,A,U,SUNIF,USTAR,UMIN
      INTEGER IR,I
C
C     Q(N) = SUM(ALOG(2.0)**K/K])    K=1,..,N ,      THE HIGHEST N
C     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
C
      DATA Q/.6931472,.9333737,.9888778,.9984959,
     ,.9998293,.9999833,.9999986,.9999999/
C
   1  A=0.0
      U=SUNIF(IR)
      GO TO 2
   3  A=A+Q1
   2  U=U+U
      IF (U.LE.1.0) GO TO 3
   4  U=U-1.0
      IF (U.GT.Q1) GO TO 6
   5  SEXPO=A+U
      RETURN
   6  I=1
      USTAR=SUNIF(IR)
      UMIN=USTAR
   7  USTAR=SUNIF(IR)
      IF (USTAR.LT.UMIN) UMIN=USTAR
   8  I=I+1
      IF (U.GT.Q(I)) GO TO 7
   9  SEXPO=A+UMIN*Q1
      RETURN
      END
C                                                                       
C**********************************************************************
C**********************************************************************
C                                                                      
C                                                                      
C     (STANDARD-)  N O R M A L  DISTRIBUTION                           
C                                                                     
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C     FOR DETAILS SEE:                                                 
C                                                                     
C               AHRENS, J.H. AND DIETER, U.                           
C               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM            
C               SAMPLING FROM THE NORMAL DISTRIBUTION.                
C               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.         
C                                                                     
C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  
C     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  
C                                                                      
C**********************************************************************
C                                                                      
      REAL FUNCTION SNORM(IR)                                          
      REAL A(32),D(31),T(31),H(31),SUNIF,S,USTAR,AA,W,U,TT,Y
      INTEGER IR,I
C
C     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
C     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
C
      DATA A/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,
     ,.1970991,.2372021,.2776904,.3186394,.3601299,.4022501,
     ,.4450965,.4887764,.5334097,.5791322,.6260990,.6744898,
     ,.7245144,.7764218,.8305109,.8871466,.9467818,1.009990,
     ,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,
     ,1.675940,1.862732,2.153875/
      DATA D/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
     ,.1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
     ,.1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
     ,.1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
     ,.1134023,.1114027,.1095039/
      DATA T/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
     ,.7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
     ,.1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
     ,.2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
     ,.4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
     ,.9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
     ,.5847031/
      DATA H/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
     ,.4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
     ,.4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
     ,.4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
     ,.5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
     ,.8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
     ,.7010474/
C
   1  U=SUNIF(IR)
      S=0.0
      IF (U.GE.0.5) S=1.0
      U=U+U-S
   2  U=32.0*U
      I=INT(U)
      IF (I.EQ.0) GO TO 9
C
C                                START CENTER
C
   3  USTAR=U-FLOAT(I)
      AA=A(I)
   4  IF (USTAR.LE.T(I)) GO TO 5
      W=(USTAR-T(I))*H(I)
C
C                                EXIT   (BOTH CASES)
C
  17  Y=AA+W
      SNORM=Y
      IF (S.EQ.1.0) SNORM=-Y
      RETURN
C
C                                CENTER CONTINUED
C
   5  U=SUNIF(IR)
      W=U*(A(I+1)-AA)
      TT=(0.5*W+AA)*W
      GO TO 6
   8  TT=U
      USTAR=SUNIF(IR)
   6  IF (USTAR.GT.TT) GO TO 17
   7  U=SUNIF(IR)
      IF (USTAR.GE.U) GO TO 8
      USTAR=SUNIF(IR)
      GO TO 4
C
C                                START TAIL
C
   9  I=6
      AA=A(32)
      GO TO 10
  11  AA=AA+D(I)
      I=I+1
  10  U=U+U
      IF (U.LT.1.0) GO TO 11
  12  U=U-1.0
  13  W=U*D(I)
      TT=(0.5*W+AA)*W
      GO TO 14
  16  TT=U
  14  USTAR=SUNIF(IR)
      IF (USTAR.GT.TT) GO TO 17
  15  U=SUNIF(IR)
      IF (USTAR.GE.U) GO TO 16
      U=SUNIF(IR)
      GO TO 13
      END
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C                                                                      
C     (STANDARD-)  G A M M A  DISTRIBUTION                             
C                                                                     
C                                                                     
C**********************************************************************
C**********************************************************************
C                                                                      
C               PARAMETER  A >= 1.0  ]                                
C                                                                      
C**********************************************************************
C                                                                      
C     FOR DETAILS SEE:                                                 
C                                                                      
C               AHRENS, J.H. AND DIETER, U.                           
C               GENERATING GAMMA VARIATES BY A                        
C               MODIFIED REJECTION TECHNIQUE.                         
C               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                 
C                                                                     
C     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     
C                                 (STRAIGHTFORWARD IMPLEMENTATION)     
C                                                                      
C**********************************************************************
C                                                                     
C               PARAMETER  0.0 < A < 1.0  ]                           
C                                                                      
C**********************************************************************
C                                                                      
C     FOR DETAILS SEE:                                                
C                                                                     
C               AHRENS, J.H. AND DIETER, U.                           
C               COMPUTER METHODS FOR SAMPLING FROM GAMMA,             
C               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              
C               COMPUTING, 12 (1974), 223 - 246.                      
C                                                                     
C     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)   
C                                                                      
C**********************************************************************
C                                                                      
      REAL FUNCTION SGAMMA(IR,A)                                      

      REAL A,Q1,Q2,Q3,Q4,Q5,Q6,Q7,A1,A2,A3,A4,A5,A6,A7
      REAL E1,E2,E3,E4,E5,AA,AAA,SQRT32,S2,S,D,T,SNORM,X,U
      REAL SUNIF,R,Q0,B,SI,C,V,Q,E,SEXPO,W,P
      INTEGER IR
C
C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
C             A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
C     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
C
C     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
C     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
C     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
C
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7 /.04166669,.02083148,
     ,.00801191,.00144121,-.00007388,.00024511,.00024240/
      DATA A1,A2,A3,A4,A5,A6,A7 /.3333333,-.2500030,
     ,.2000062,-.1662921,.1423657,-.1367177,.1233795/
      DATA E1,E2,E3,E4,E5 /1.,.4999897,.1668290,.0407753,.0102930/
C
C     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
C     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
C
      DATA AA /0.0/, AAA /0.0/, SQRT32 /5.656854/
C
      IF (A .EQ. AA) GO TO 1
      IF (A .LT. 1.0) GO TO 12
C
C     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
C
      AA=A
      S2=A-0.5
      S=SQRT(S2)
      D=SQRT32-12.0*S
C
C     STEP  2:  T=STANDARD NORMAL DEVIATE,
C               X=(S,1/2)-NORMAL DEVIATE.
C               IMMEDIATE ACCEPTANCE (I)
C
   1  T=SNORM(IR)
      X=S+0.5*T
      SGAMMA=X*X
      IF (T .GE. 0.0) RETURN
C
C     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
C
      U=SUNIF(IR)
      IF (D*U .LE. T*T*T) RETURN
C
C     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
C
      IF (A .EQ. AAA) GO TO 4
      AAA=A
      R=1.0/A
      Q0=((((((Q7*R+Q6)*R+Q5)*R+Q4)*R+Q3)*R+Q2)*R+Q1)*R
C
C               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
C               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
C               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
C
      IF (A .LE. 3.686) GO TO 3
      IF (A .LE. 13.022) GO TO 2
C
C               CASE 3:  A .GT. 13.022
C
      B=1.77
      SI=.75
      C=.1515/S
      GO TO 4
C
C               CASE 2:  3.686 .LT. A .LE. 13.022
C
   2  B=1.654+.0076*S2
      SI=1.68/S+.275
      C=.062/S+.024
      GO TO 4
C
C               CASE 1:  A .LE. 3.686
C
   3  B=.463+S-.178*S2
      SI=1.235
      C=.195/S-.079+.016*S
C
C     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
C
   4  IF (X .LE. 0.0) GO TO 7
C
C     STEP  6:  CALCULATION OF V AND QUOTIENT Q
C
      V=T/(S+S)
      IF (ABS(V) .LE. 0.25) GO TO 5
      Q=Q0-S*T+0.25*T*T+(S2+S2)*ALOG(1.0+V)
      GO TO 6
   5  Q=Q0+0.5*T*T*((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V
C
C     STEP  7:  QUOTIENT ACCEPTANCE (Q)
C
   6  IF (ALOG(1.0-U) .LE. Q) RETURN
C
C     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
C               U= 0,1 -UNIFORM DEVIATE
C               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
C
   7  E=SEXPO(IR)
      U=SUNIF(IR)
      U=U+U-1.0
      T=B+SIGN(SI*E,U)
C
C     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
C
      IF (T .LT. (-.7187449)) GO TO 7
C
C     STEP 10:  CALCULATION OF V AND QUOTIENT Q
C
      V=T/(S+S)
      IF (ABS(V) .LE. 0.25) GO TO 8
      Q=Q0-S*T+0.25*T*T+(S2+S2)*ALOG(1.0+V)
      GO TO 9
   8  Q=Q0+0.5*T*T*((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V
C
C     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
C
   9  IF (Q .LE. 0.0) GO TO 7
      IF (Q .LE. 0.5) GO TO 10
      W=EXP(Q)-1.0
      GO TO 11
  10  W=((((E5*Q+E4)*Q+E3)*Q+E2)*Q+E1)*Q
C
C               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
C
  11  IF (C*ABS(U) .GT. W*EXP(E-0.5*T*T)) GO TO 7
      X=S+0.5*T
      SGAMMA=X*X
      RETURN
C
C     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
C
  12  AA=0.0
      B=1.0+.3678794*A
  13  P=B*SUNIF(IR)
      IF (P .GE. 1.0) GO TO 14
      SGAMMA=EXP(ALOG(P)/A)
      IF (SEXPO(IR) .LT. SGAMMA) GO TO 13
      RETURN
  14  SGAMMA=-ALOG((B-P)/A)
      IF (SEXPO(IR) .LT. (1.0-A)*ALOG(SGAMMA)) GO TO 13
      RETURN
      END
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C                                                                      
C     P O I S S O N  DISTRIBUTION                                      
C                                                                     
C                                                                      
C**********************************************************************
C**********************************************************************
C                                                                      
C     FOR DETAILS SEE:                                                 
C                                                                     
C               AHRENS, J.H. AND DIETER, U.                            
C               COMPUTER GENERATION OF POISSON DEVIATES                
C               FROM MODIFIED NORMAL DISTRIBUTIONS.                    
C               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. 
C                                                                      
C     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE) 
C                                                                      
C**********************************************************************
C                                                                      
      INTEGER FUNCTION KPOISS(IR,MU)                                   

      INTEGER IR,L,KFLAG,M,J,K
      REAL A0,A1,A2,A3,A4,A5,A6,A7,FACT,S,D,G,SNORM,FK,DIFMUK
      REAL U,SUNIF,OMEGA,B1,B2,C,C0,C1,C2,C3,FY,PY,PX,FX,E
      REAL SEXPO,T,DEL,V,X,XX,P,Q,P0,PP

C
C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
C             MU=MEAN MU OF THE POISSON DISTRIBUTION
C     OUTPUT: KPOISS=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
C
      REAL MU, MUPREV, MUOLD
C
C     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
C     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
C     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
C
      DIMENSION FACT(10), PP(35)
      DATA MUPREV,MUOLD /0.,0./
      DATA A0,A1,A2,A3,A4,A5,A6,A7 /-.5,.3333333,-.2500068,
     ,.2000118,-.1661269,.1421878,-.1384794,.1250060/
      DATA FACT /1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
C
C     SEPARATION OF CASES A AND B
C
      IF (MU .EQ. MUPREV) GO TO 1
      IF (MU .LT. 10.0) GO TO 12
C
C     C A S E  A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
C
      MUPREV=MU
      S=SQRT(MU)
      D=6.0*MU*MU
C
C             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
C             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
C             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
C
      L=IFIX(MU-1.1484)
C
C     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
C
   1  G=MU+S*SNORM(IR)
      IF (G .LT. 0.0) GO TO 2
      KPOISS=IFIX(G)
C
C     STEP I. IMMEDIATE ACCEPTANCE IF KPOISS IS LARGE ENOUGH
C
      IF (KPOISS .GE. L) RETURN
C
C     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
C
      FK=FLOAT(KPOISS)
      DIFMUK=MU-FK
      U=SUNIF(IR)
      IF (D*U .GE. DIFMUK*DIFMUK*DIFMUK) RETURN
C
C     STEP P. PREPARATIONS FOR STEPS Q AND H.
C             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
C             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
C             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
C             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
C             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
C
   2  IF (MU .EQ. MUOLD) GO TO 3
      MUOLD=MU
      OMEGA=.3989423/S
      B1=.4166667E-1/MU
      B2=.3*B1*B1
      C3=.1428571*B1*B2
      C2=B2-15.*C3
      C1=B1-6.*B2+45.*C3
      C0=1.-B1+3.*B2-15.*C3
      C=.1069/MU
   3  IF (G .LT. 0.0) GO TO 5
C
C             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
C
      KFLAG=0
      GO TO 7
C
C     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
C
   4  IF (FY-U*FY .LE. PY*EXP(PX-FX)) RETURN
C
C     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
C             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
C             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
C
   5  E=SEXPO(IR)
      U=SUNIF(IR)
      U=U+U-1.0
      T=1.8+SIGN(E,U)
      IF (T .LE. (-.6744)) GO TO 5
      KPOISS=IFIX(MU+S*T)
      FK=FLOAT(KPOISS)
      DIFMUK=MU-FK
C
C             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
C
      KFLAG=1
      GO TO 7
C
C     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
C
   6  IF (C*ABS(U) .GT. PY*EXP(PX+E)-FY*EXP(FX+E)) GO TO 5
      RETURN
C
C     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
C             CASE KPOISS .LT. 10 USES FACTORIALS FROM TABLE FACT
C
   7  IF (KPOISS .GE. 10) GO TO 8
      PX=-MU
      PY=MU**KPOISS/FACT(KPOISS+1)
      GO TO 11
C
C             CASE KPOISS .GE. 10 USES POLYNOMIAL APPROXIMATION
C             A0-A7 FOR ACCURACY WHEN ADVISABLE
C             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
C
   8  DEL=.8333333E-1/FK
      DEL=DEL-4.8*DEL*DEL*DEL
      V=DIFMUK/FK
      IF (ABS(V) .LE. 0.25) GO TO 9
      PX=FK*ALOG(1.0+V)-DIFMUK-DEL
      GO TO 10
   9  PX=FK*V*V*(((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V+A0)-DEL
  10  PY=.3989423/SQRT(FK)
  11  X=(0.5-DIFMUK)/S
      XX=X*X
      FX=-0.5*XX
      FY=OMEGA*(((C3*XX+C2)*XX+C1)*XX+C0)
      IF (KFLAG) 4,4,6
C
C     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
C
  12  MUPREV=0.0
      IF (MU .EQ. MUOLD) GO TO 13
      MUOLD=MU
      M=MAX0(1,IFIX(MU))
      L=0
      P=EXP(-MU)
      Q=P
      P0=P
C
C     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
C
  13  U=SUNIF(IR)
      KPOISS=0
      IF (U .LE. P0) RETURN
C
C     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
C             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
C             (0.458=PP(9) FOR MU=10)
C
      IF (L .EQ. 0) GO TO 15
      J=1
      IF (U .GT. 0.458) J=MIN0(L,M)
      DO 14 K=J,L
      IF (U .LE. PP(K)) GO TO 18
  14  CONTINUE
      IF (L .EQ. 35) GO TO 13
C
C     STEP C. CREATION OF NEW POISSON PROBABILITIES P
C             AND THEIR CUMULATIVES Q=PP(K)
C
  15  L=L+1
      DO 16 K=L,35
      P=P*MU/FLOAT(K)
      Q=Q+P
      PP(K)=Q
      IF (U .LE. Q) GO TO 17
  16  CONTINUE
      L=35
      GO TO 13
  17  L=K
  18  KPOISS=K
      RETURN
      END

C* ======================================================================
C* NIST Guide to Available Math Software.
C* Fullsource for module RPQR79 from package CMLIB.
C* Retrieved from TIBER on Thu Jun  5 17:40:41 1997.
C* ======================================================================
      SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR)
C***BEGIN PROLOGUE  HQR
C***DATE WRITTEN   760101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***CATEGORY NO.  D4C2B
C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
C***AUTHOR  SMITH, B. T., ET AL.
C***PURPOSE  Computes eigenvalues of a real upper Hessenberg matrix
C            using the QR method.
C***DESCRIPTION
C
C     This subroutine is a translation of the ALGOL procedure HQR,
C     NUM. MATH. 14, 219-231(1970) by Martin, Peters, and Wilkinson.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
C
C     This subroutine finds the eigenvalues of a REAL
C     UPPER Hessenberg matrix by the QR method.
C
C     On INPUT
C
C        NM must be set to the row dimension of two-dimensional
C          array parameters as declared in the calling program
C          dimension statement.
C
C        N is the order of the matrix.
C
C        LOW and IGH are integers determined by the balancing
C          subroutine  BALANC.  If  BALANC  has not been used,
C          set LOW=1, IGH=N.
C
C        H contains the upper Hessenberg matrix.  Information about
C          the transformations used in the reduction to Hessenberg
C          form by  ELMHES  or  ORTHES, if performed, is stored
C          in the remaining triangle under the Hessenberg matrix.
C
C     On OUTPUT
C
C        H has been destroyed.  Therefore, it must be saved
C          before calling  HQR  if subsequent calculation and
C          back transformation of eigenvectors is to be performed.
C
C        WR and WI contain the real and imaginary parts,
C          respectively, of the eigenvalues.  The eigenvalues
C          are unordered except that complex conjugate pairs
C          of values appear consecutively with the eigenvalue
C          having the positive imaginary part first.  If an
C          error exit is made, the eigenvalues should be correct
C          for indices IERR+1,...,N.
C
C        IERR is set to
C          Zero       for normal return,
C          J          if the J-th eigenvalue has not been
C                     determined after a total of 30*N iterations.
C
C     Questions and comments should be directed to B. S. Garbow,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
C                 1976.
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  HQR
C
      INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR
      REAL H(NM,N),WR(N),WI(N)
      REAL P,Q,R,S,T,W,X,Y,ZZ,NORM,S1,S2
      LOGICAL NOTLAS
C
C***FIRST EXECUTABLE STATEMENT  HQR
      IERR = 0
      NORM = 0.0E0
      K = 1
C     .......... STORE ROOTS ISOLATED BY BALANC
C                AND COMPUTE MATRIX NORM ..........
      DO 50 I = 1, N
C
         DO 40 J = K, N
   40    NORM = NORM + ABS(H(I,J))
C
         K = I
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
         WR(I) = H(I,I)
         WI(I) = 0.0E0
   50 CONTINUE
C
      EN = IGH
      T = 0.0E0
      ITN = 30*N
C     .......... SEARCH FOR NEXT EIGENVALUES ..........
   60 IF (EN .LT. LOW) GO TO 1001
      ITS = 0
      NA = EN - 1
      ENM2 = NA - 1
C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
   70 DO 80 LL = LOW, EN
         L = EN + LOW - LL
         IF (L .EQ. LOW) GO TO 100
         S = ABS(H(L-1,L-1)) + ABS(H(L,L))
         IF (S .EQ. 0.0E0) S = NORM
         S2 = S + ABS(H(L,L-1))
         IF (S2 .EQ. S) GO TO 100
   80 CONTINUE
C     .......... FORM SHIFT ..........
  100 X = H(EN,EN)
      IF (L .EQ. EN) GO TO 270
      Y = H(NA,NA)
      W = H(EN,NA) * H(NA,EN)
      IF (L .EQ. NA) GO TO 280
      IF (ITN .EQ. 0) GO TO 1000
      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C     .......... FORM EXCEPTIONAL SHIFT ..........
      T = T + X
C
      DO 120 I = LOW, EN
  120 H(I,I) = H(I,I) - X
C
      S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
      X = 0.75E0 * S
      Y = X
      W = -0.4375E0 * S * S
  130 ITS = ITS + 1
      ITN = ITN - 1
C     .......... LOOK FOR TWO CONSECUTIVE SMALL
C                SUB-DIAGONAL ELEMENTS.
C                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
      DO 140 MM = L, ENM2
         M = ENM2 + L - MM
         ZZ = H(M,M)
         R = X - ZZ
         S = Y - ZZ
         P = (R * S - W) / H(M+1,M) + H(M,M+1)
         Q = H(M+1,M+1) - ZZ - R - S
         R = H(M+2,M+1)
         S = ABS(P) + ABS(Q) + ABS(R)
         P = P / S
         Q = Q / S
         R = R / S
         IF (M .EQ. L) GO TO 150
         S1 = ABS(P) * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))
         S2 = S1 + ABS(H(M,M-1)) * (ABS(Q) + ABS(R))
         IF (S2 .EQ. S1) GO TO 150
  140 CONTINUE
C
  150 MP2 = M + 2
C
      DO 160 I = MP2, EN
         H(I,I-2) = 0.0E0
         IF (I .EQ. MP2) GO TO 160
         H(I,I-3) = 0.0E0
  160 CONTINUE
C     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
C                COLUMNS M TO EN ..........
      DO 260 K = M, NA
         NOTLAS = K .NE. NA
         IF (K .EQ. M) GO TO 170
         P = H(K,K-1)
         Q = H(K+1,K-1)
         R = 0.0E0
         IF (NOTLAS) R = H(K+2,K-1)
         X = ABS(P) + ABS(Q) + ABS(R)
         IF (X .EQ. 0.0E0) GO TO 260
         P = P / X
         Q = Q / X
         R = R / X
  170    S = SIGN(SQRT(P*P+Q*Q+R*R),P)
         IF (K .EQ. M) GO TO 180
         H(K,K-1) = -S * X
         GO TO 190
  180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
  190    P = P + S
         X = P / S
         Y = Q / S
         ZZ = R / S
         Q = Q / P
         R = R / P
C     .......... ROW MODIFICATION ..........
         DO 210 J = K, EN
            P = H(K,J) + Q * H(K+1,J)
            IF (.NOT. NOTLAS) GO TO 200
            P = P + R * H(K+2,J)
            H(K+2,J) = H(K+2,J) - P * ZZ
  200       H(K+1,J) = H(K+1,J) - P * Y
            H(K,J) = H(K,J) - P * X
  210    CONTINUE
C
         J = MIN0(EN,K+3)
C     .......... COLUMN MODIFICATION ..........
         DO 230 I = L, J
            P = X * H(I,K) + Y * H(I,K+1)
            IF (.NOT. NOTLAS) GO TO 220
            P = P + ZZ * H(I,K+2)
            H(I,K+2) = H(I,K+2) - P * R
  220       H(I,K+1) = H(I,K+1) - P * Q
            H(I,K) = H(I,K) - P
  230    CONTINUE
C
  260 CONTINUE
C
      GO TO 70
C     .......... ONE ROOT FOUND ..........
  270 WR(EN) = X + T
      WI(EN) = 0.0E0
      EN = NA
      GO TO 60
C     .......... TWO ROOTS FOUND ..........
  280 P = (Y - X) / 2.0E0
      Q = P * P + W
      ZZ = SQRT(ABS(Q))
      X = X + T
      IF (Q .LT. 0.0E0) GO TO 320
C     .......... REAL PAIR ..........
      ZZ = P + SIGN(ZZ,P)
      WR(NA) = X + ZZ
      WR(EN) = WR(NA)
      IF (ZZ .NE. 0.0E0) WR(EN) = X - W / ZZ
      WI(NA) = 0.0E0
      WI(EN) = 0.0E0
      GO TO 330
C     .......... COMPLEX PAIR ..........
  320 WR(NA) = X + P
      WR(EN) = X + P
      WI(NA) = ZZ
      WI(EN) = -ZZ
  330 EN = ENM2
      GO TO 60
C     .......... SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30*N ITERATIONS ..........
 1000 IERR = EN
 1001 RETURN
      END


      SUBROUTINE RPQR79(NDEG,COEFF,ROOT,IERR,WORK)
C***BEGIN PROLOGUE  RPQR79
C***DATE WRITTEN   800601   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  F1A1A
C***KEYWORDS  POLYNOMIAL ROOTS,REAL,ROOTS,ZEROES,ZEROS
C***AUTHOR  VANDEVENDER, W. H., (SNLA)
C***PURPOSE  To find the zeros of a polynomial with real coefficients.
C***DESCRIPTION
C
C     This routine is an interface to an eigenvalue routine in EISPACK.
C     This interface was written by Walter H. Vandevender.
C
C     Abstract
C
C         This routine computes all roots of a polynomial with real
C         coefficients by computing the eigenvalues of the
C         companion matrix.
C
C     Description of Parameters
C         The user must dimension all arrays appearing in the call list
C              COEFF(NDEG+1),ROOT(NDEG),WORK(NDEG*NDEG+2*NDEG)
C
C      Input -
C       NDEG   degree of polynomial
C
C       COEFF  array of coeffficients in order of descending powers of
C              Z,   i.e. COEFF(1)*(Z**NDEG) + ... + COEFF(NDEG+1)
C
C       WORK   real work vector of length at least NDEG*(NDEG+2)
C
C     Output-
C       ROOT   complex vector of computed roots.
C
C       IERR   output error code
C             - normal code
C              0  means the roots were computed.
C             - abnormal codes
C              1  more than 30 QR iterations on some eigenvalue
C                 of the companion matrix
C              2  COEF(1) = 0.0
C              3  NDEG invalid (less than 0).
C***REFERENCES  (NONE)
C***ROUTINES CALLED  HQR
C***END PROLOGUE  RPQR79
      INTEGER NDEG,IERR,K,KH,KWR,KWI,KCOL,KWEND,KM1
      REAL COEFF(NDEG+1),WORK(NDEG*NDEG+2*NDEG),SCALE
      COMPLEX ROOT(NDEG)
C***FIRST EXECUTABLE STATEMENT  RPQR79
      IERR = 0
      IF(COEFF(1).EQ.0.0)   GOTO 5
      IF(NDEG.LT.0) GOTO 6
      IF(NDEG.EQ.1) GO TO 7
      SCALE=1.E0/COEFF(1)
      KH=1
      KWR=KH+NDEG*NDEG
      KWI=KWR+NDEG
      KWEND=KWI+NDEG-1
      DO 1 K=1,KWEND
1     WORK(K)=0.E0
      DO 2 K=1,NDEG
       KCOL=(K-1)*NDEG+1
       WORK(KCOL)=-COEFF(K+1)*SCALE
       IF(K.EQ.NDEG) GO TO 2
        WORK(KCOL+K)=1.E0
2     CONTINUE
      CALL HQR(NDEG,NDEG,1,NDEG,WORK(KH),WORK(KWR),WORK(KWI),IERR)
      IF(IERR.NE.0) GO TO 4
      DO 3 K=1,NDEG
       KM1=K-1
3     ROOT(K)=CMPLX(WORK(KWR+KM1),WORK(KWI+KM1))
      RETURN
C
C       IERR   output error code
C             - normal code
C              0  means the roots were computed.
C             - abnormal codes
C              1  more than 30 QR iterations on some eigenvalue
C                 of the companion matrix
C              2  COEF(1) = 0.0
C              3  NDEG invalid (less than 0).

4     IERR = 1
      RETURN
5     IERR = 2
      RETURN
6     IERR = 3
      RETURN
7      ROOT(1)=CMPLX(-COEFF(2)/COEFF(1),0.)
      RETURN
      END

* ======================================================================
* NIST GUIDE TO AVAILABLE MATH SOFTWARE.
* FULLSOURCE FOR MODULE SPOSV FROM PACKAGE LAPACK.
* RETRIEVED FROM NETLIB ON TUE JUN 17 13:13:58 1997.
* ======================================================================
      SUBROUTINE SPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
*  -- LAPACK DRIVER ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     MARCH 31, 1993
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, LDB, N, NRHS
*     ..
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * ), B( LDB, * )
*     ..
*
*  PURPOSE
*  =======
*
*  SPOSV COMPUTES THE SOLUTION TO A REAL SYSTEM OF LINEAR EQUATIONS
*     A * X = B,
*  WHERE A IS AN N-BY-N SYMMETRIC POSITIVE DEFINITE MATRIX AND X AND B
*  ARE N-BY-NRHS MATRICES.
*
*  THE CHOLESKY DECOMPOSITION IS USED TO FACTOR A AS
*     A = U**T* U,  IF UPLO = 'U', OR
*     A = L * L**T,  IF UPLO = 'L',
*  WHERE U IS AN UPPER TRIANGULAR MATRIX AND L IS A LOWER TRIANGULAR
*  MATRIX.  THE FACTORED FORM OF A IS THEN USED TO SOLVE THE SYSTEM OF
*  EQUATIONS A * X = B.
*
*  ARGUMENTS
*  =========
*
*  UPLO    (INPUT) CHARACTER*1
*          = 'U':  UPPER TRIANGLE OF A IS STORED;
*          = 'L':  LOWER TRIANGLE OF A IS STORED.
*
*  N       (INPUT) INTEGER
*          THE NUMBER OF LINEAR EQUATIONS, I.E., THE ORDER OF THE
*          MATRIX A.  N >= 0.
*
*  NRHS    (INPUT) INTEGER
*          THE NUMBER OF RIGHT HAND SIDES, I.E., THE NUMBER OF COLUMNS
*          OF THE MATRIX B.  NRHS >= 0.
*
*  A       (INPUT/OUTPUT) REAL ARRAY, DIMENSION (LDA,N)
*          ON ENTRY, THE SYMMETRIC MATRIX A.  IF UPLO = 'U', THE LEADING
*          N-BY-N UPPER TRIANGULAR PART OF A CONTAINS THE UPPER
*          TRIANGULAR PART OF THE MATRIX A, AND THE STRICTLY LOWER
*          TRIANGULAR PART OF A IS NOT REFERENCED.  IF UPLO = 'L', THE
*          LEADING N-BY-N LOWER TRIANGULAR PART OF A CONTAINS THE LOWER
*          TRIANGULAR PART OF THE MATRIX A, AND THE STRICTLY UPPER
*          TRIANGULAR PART OF A IS NOT REFERENCED.
*
*          ON EXIT, IF INFO = 0, THE FACTOR U OR L FROM THE CHOLESKY
*          FACTORIZATION A = U**T*U OR A = L*L**T.
*
*  LDA     (INPUT) INTEGER
*          THE LEADING DIMENSION OF THE ARRAY A.  LDA >= MAX(1,N).
*
*  B       (INPUT/OUTPUT) REAL ARRAY, DIMENSION (LDB,NRHS)
*          ON ENTRY, THE N-BY-NRHS RIGHT HAND SIDE MATRIX B.
*          ON EXIT, IF INFO = 0, THE N-BY-NRHS SOLUTION MATRIX X.
*
*  LDB     (INPUT) INTEGER
*          THE LEADING DIMENSION OF THE ARRAY B.  LDB >= MAX(1,N).
*
*  INFO    (OUTPUT) INTEGER
*          = 0:  SUCCESSFUL EXIT
*          < 0:  IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE
*          > 0:  IF INFO = I, THE LEADING MINOR OF ORDER I OF A IS NOT
*                POSITIVE DEFINITE, SO THE FACTORIZATION COULD NOT BE
*                COMPLETED, AND THE SOLUTION HAS NOT BEEN COMPUTED.
*
*  =====================================================================
*
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           SPOTRF, SPOTRS, XERBLA
*     ..
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SPOSV ', -INFO )
         RETURN
      END IF
*
*     COMPUTE THE CHOLESKY FACTORIZATION A = U'*U OR A = L*L'.
*
      CALL SPOTRF( UPLO, N, A, LDA, INFO )
      IF( INFO.EQ.0 ) THEN
*
*        SOLVE THE SYSTEM A*X = B, OVERWRITING B WITH X.
*
         CALL SPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
      END IF
      RETURN
*
*     END OF SPOSV
*
      END
      SUBROUTINE SPOTF2( UPLO, N, A, LDA, INFO )
*
*  -- LAPACK ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     FEBRUARY 29, 1992
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, N
*     ..
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * )
*     ..
*
*  PURPOSE
*  =======
*
*  SPOTF2 COMPUTES THE CHOLESKY FACTORIZATION OF A REAL SYMMETRIC
*  POSITIVE DEFINITE MATRIX A.
*
*  THE FACTORIZATION HAS THE FORM
*     A = U' * U ,  IF UPLO = 'U', OR
*     A = L  * L',  IF UPLO = 'L',
*  WHERE U IS AN UPPER TRIANGULAR MATRIX AND L IS LOWER TRIANGULAR.
*
*  THIS IS THE UNBLOCKED VERSION OF THE ALGORITHM, CALLING LEVEL 2 BLAS.
*
*  ARGUMENTS
*  =========
*
*  UPLO    (INPUT) CHARACTER*1
*          SPECIFIES WHETHER THE UPPER OR LOWER TRIANGULAR PART OF THE
*          SYMMETRIC MATRIX A IS STORED.
*          = 'U':  UPPER TRIANGULAR
*          = 'L':  LOWER TRIANGULAR
*
*  N       (INPUT) INTEGER
*          THE ORDER OF THE MATRIX A.  N >= 0.
*
*  A       (INPUT/OUTPUT) REAL ARRAY, DIMENSION (LDA,N)
*          ON ENTRY, THE SYMMETRIC MATRIX A.  IF UPLO = 'U', THE LEADING
*          N BY N UPPER TRIANGULAR PART OF A CONTAINS THE UPPER
*          TRIANGULAR PART OF THE MATRIX A, AND THE STRICTLY LOWER
*          TRIANGULAR PART OF A IS NOT REFERENCED.  IF UPLO = 'L', THE
*          LEADING N BY N LOWER TRIANGULAR PART OF A CONTAINS THE LOWER
*          TRIANGULAR PART OF THE MATRIX A, AND THE STRICTLY UPPER
*          TRIANGULAR PART OF A IS NOT REFERENCED.
*
*          ON EXIT, IF INFO = 0, THE FACTOR U OR L FROM THE CHOLESKY
*          FACTORIZATION A = U'*U  OR A = L*L'.
*
*  LDA     (INPUT) INTEGER
*          THE LEADING DIMENSION OF THE ARRAY A.  LDA >= MAX(1,N).
*
*  INFO    (OUTPUT) INTEGER
*          = 0: SUCCESSFUL EXIT
*          < 0: IF INFO = -K, THE K-TH ARGUMENT HAD AN ILLEGAL VALUE
*          > 0: IF INFO = K, THE LEADING MINOR OF ORDER K IS NOT
*               POSITIVE DEFINITE, AND THE FACTORIZATION COULD NOT BE
*               COMPLETED.
*
*  =====================================================================
*
*     .. PARAMETERS ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. LOCAL SCALARS ..
      LOGICAL            UPPER
      INTEGER            J
      REAL               AJJ
*     ..
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      REAL               SDOT
      EXTERNAL           LSAME, SDOT
*     ..
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           SGEMV, SSCAL, XERBLA
*     ..
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX, SQRT
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      UPPER = LSAME( UPLO, 'U' )
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SPOTF2', -INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE
*
      IF( N.EQ.0 )
     $   RETURN
*
      IF( UPPER ) THEN
*
*        COMPUTE THE CHOLESKY FACTORIZATION A = U'*U.
*
         DO 10 J = 1, N
*
*           COMPUTE U(J,J) AND TEST FOR NON-POSITIVE-DEFINITENESS.
*
            AJJ = A( J, J ) - SDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 )
            IF( AJJ.LE.ZERO ) THEN
               A( J, J ) = AJJ
               GO TO 30
            END IF
            AJJ = SQRT( AJJ )
            A( J, J ) = AJJ
*
*           COMPUTE ELEMENTS J+1:N OF ROW J.
*
            IF( J.LT.N ) THEN
               CALL SGEMV( 'TRANSPOSE', J-1, N-J, -ONE, A( 1, J+1 ),
     $                     LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
               CALL SSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
            END IF
   10    CONTINUE
      ELSE
*
*        COMPUTE THE CHOLESKY FACTORIZATION A = L*L'.
*
         DO 20 J = 1, N
*
*           COMPUTE L(J,J) AND TEST FOR NON-POSITIVE-DEFINITENESS.
*
            AJJ = A( J, J ) - SDOT( J-1, A( J, 1 ), LDA, A( J, 1 ),
     $            LDA )
            IF( AJJ.LE.ZERO ) THEN
               A( J, J ) = AJJ
               GO TO 30
            END IF
            AJJ = SQRT( AJJ )
            A( J, J ) = AJJ
*
*           COMPUTE ELEMENTS J+1:N OF COLUMN J.
*
            IF( J.LT.N ) THEN
               CALL SGEMV( 'NO TRANSPOSE', N-J, J-1, -ONE, A( J+1, 1 ),
     $                     LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
               CALL SSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
            END IF
   20    CONTINUE
      END IF
      GO TO 40
*
   30 CONTINUE
      INFO = J
*
   40 CONTINUE
      RETURN
*
*     END OF SPOTF2
*
      END
      SUBROUTINE SPOTRF( UPLO, N, A, LDA, INFO )
*
*  -- LAPACK ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     MARCH 31, 1993
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, N
*     ..
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * )
*     ..
*
*  PURPOSE
*  =======
*
*  SPOTRF COMPUTES THE CHOLESKY FACTORIZATION OF A REAL SYMMETRIC
*  POSITIVE DEFINITE MATRIX A.
*
*  THE FACTORIZATION HAS THE FORM
*     A = U**T * U,  IF UPLO = 'U', OR
*     A = L  * L**T,  IF UPLO = 'L',
*  WHERE U IS AN UPPER TRIANGULAR MATRIX AND L IS LOWER TRIANGULAR.
*
*  THIS IS THE BLOCK VERSION OF THE ALGORITHM, CALLING LEVEL 3 BLAS.
*
*  ARGUMENTS
*  =========
*
*  UPLO    (INPUT) CHARACTER*1
*          = 'U':  UPPER TRIANGLE OF A IS STORED;
*          = 'L':  LOWER TRIANGLE OF A IS STORED.
*
*  N       (INPUT) INTEGER
*          THE ORDER OF THE MATRIX A.  N >= 0.
*
*  A       (INPUT/OUTPUT) REAL ARRAY, DIMENSION (LDA,N)
*          ON ENTRY, THE SYMMETRIC MATRIX A.  IF UPLO = 'U', THE LEADING
*          N-BY-N UPPER TRIANGULAR PART OF A CONTAINS THE UPPER
*          TRIANGULAR PART OF THE MATRIX A, AND THE STRICTLY LOWER
*          TRIANGULAR PART OF A IS NOT REFERENCED.  IF UPLO = 'L', THE
*          LEADING N-BY-N LOWER TRIANGULAR PART OF A CONTAINS THE LOWER
*          TRIANGULAR PART OF THE MATRIX A, AND THE STRICTLY UPPER
*          TRIANGULAR PART OF A IS NOT REFERENCED.
*
*          ON EXIT, IF INFO = 0, THE FACTOR U OR L FROM THE CHOLESKY
*          FACTORIZATION A = U**T*U OR A = L*L**T.
*
*  LDA     (INPUT) INTEGER
*          THE LEADING DIMENSION OF THE ARRAY A.  LDA >= MAX(1,N).
*
*  INFO    (OUTPUT) INTEGER
*          = 0:  SUCCESSFUL EXIT
*          < 0:  IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE
*          > 0:  IF INFO = I, THE LEADING MINOR OF ORDER I IS NOT
*                POSITIVE DEFINITE, AND THE FACTORIZATION COULD NOT BE
*                COMPLETED.
*
*  =====================================================================
*
*     .. PARAMETERS ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. LOCAL SCALARS ..
      LOGICAL            UPPER
      INTEGER            J, JB, NB
*     ..
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           LSAME, ILAENV
*     ..
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           SGEMM, SPOTF2, SSYRK, STRSM, XERBLA
*     ..
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX, MIN
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      UPPER = LSAME( UPLO, 'U' )
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SPOTRF', -INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE
*
      IF( N.EQ.0 )
     $   RETURN
*
*     DETERMINE THE BLOCK SIZE FOR THIS ENVIRONMENT.
*
      NB = ILAENV( 1, 'SPOTRF', UPLO, N, -1, -1, -1 )
      IF( NB.LE.1 .OR. NB.GE.N ) THEN
*
*        USE UNBLOCKED CODE.
*
         CALL SPOTF2( UPLO, N, A, LDA, INFO )
      ELSE
*
*        USE BLOCKED CODE.
*
         IF( UPPER ) THEN
*
*           COMPUTE THE CHOLESKY FACTORIZATION A = U'*U.
*
            DO 10 J = 1, N, NB
*
*              UPDATE AND FACTORIZE THE CURRENT DIAGONAL BLOCK AND TEST
*              FOR NON-POSITIVE-DEFINITENESS.
*
               JB = MIN( NB, N-J+1 )
               CALL SSYRK( 'UPPER', 'TRANSPOSE', JB, J-1, -ONE,
     $                     A( 1, J ), LDA, ONE, A( J, J ), LDA )
               CALL SPOTF2( 'UPPER', JB, A( J, J ), LDA, INFO )
               IF( INFO.NE.0 )
     $            GO TO 30
               IF( J+JB.LE.N ) THEN
*
*                 COMPUTE THE CURRENT BLOCK ROW.
*
                  CALL SGEMM( 'TRANSPOSE', 'NO TRANSPOSE', JB, N-J-JB+1,
     $                        J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ),
     $                        LDA, ONE, A( J, J+JB ), LDA )
                  CALL STRSM( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT',
     $                        JB, N-J-JB+1, ONE, A( J, J ), LDA,
     $                        A( J, J+JB ), LDA )
               END IF
   10       CONTINUE
*
         ELSE
*
*           COMPUTE THE CHOLESKY FACTORIZATION A = L*L'.
*
            DO 20 J = 1, N, NB
*
*              UPDATE AND FACTORIZE THE CURRENT DIAGONAL BLOCK AND TEST
*              FOR NON-POSITIVE-DEFINITENESS.
*
               JB = MIN( NB, N-J+1 )
               CALL SSYRK( 'LOWER', 'NO TRANSPOSE', JB, J-1, -ONE,
     $                     A( J, 1 ), LDA, ONE, A( J, J ), LDA )
               CALL SPOTF2( 'LOWER', JB, A( J, J ), LDA, INFO )
               IF( INFO.NE.0 )
     $            GO TO 30
               IF( J+JB.LE.N ) THEN
*
*                 COMPUTE THE CURRENT BLOCK COLUMN.
*
                  CALL SGEMM( 'NO TRANSPOSE', 'TRANSPOSE', N-J-JB+1, JB,
     $                        J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ),
     $                        LDA, ONE, A( J+JB, J ), LDA )
                  CALL STRSM( 'RIGHT', 'LOWER', 'TRANSPOSE', 'NON-UNIT',
     $                        N-J-JB+1, JB, ONE, A( J, J ), LDA,
     $                        A( J+JB, J ), LDA )
               END IF
   20       CONTINUE
         END IF
      END IF
      GO TO 40
*
   30 CONTINUE
      INFO = INFO + J - 1
*
   40 CONTINUE
      RETURN
*
*     END OF SPOTRF
*
      END
      SUBROUTINE SPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
*  -- LAPACK ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     MARCH 31, 1993
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, LDB, N, NRHS
*     ..
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * ), B( LDB, * )
*     ..
*
*  PURPOSE
*  =======
*
*  SPOTRS SOLVES A SYSTEM OF LINEAR EQUATIONS A*X = B WITH A SYMMETRIC
*  POSITIVE DEFINITE MATRIX A USING THE CHOLESKY FACTORIZATION
*  A = U**T*U OR A = L*L**T COMPUTED BY SPOTRF.
*
*  ARGUMENTS
*  =========
*
*  UPLO    (INPUT) CHARACTER*1
*          = 'U':  UPPER TRIANGLE OF A IS STORED;
*          = 'L':  LOWER TRIANGLE OF A IS STORED.
*
*  N       (INPUT) INTEGER
*          THE ORDER OF THE MATRIX A.  N >= 0.
*
*  NRHS    (INPUT) INTEGER
*          THE NUMBER OF RIGHT HAND SIDES, I.E., THE NUMBER OF COLUMNS
*          OF THE MATRIX B.  NRHS >= 0.
*
*  A       (INPUT) REAL ARRAY, DIMENSION (LDA,N)
*          THE TRIANGULAR FACTOR U OR L FROM THE CHOLESKY FACTORIZATION
*          A = U**T*U OR A = L*L**T, AS COMPUTED BY SPOTRF.
*
*  LDA     (INPUT) INTEGER
*          THE LEADING DIMENSION OF THE ARRAY A.  LDA >= MAX(1,N).
*
*  B       (INPUT/OUTPUT) REAL ARRAY, DIMENSION (LDB,NRHS)
*          ON ENTRY, THE RIGHT HAND SIDE MATRIX B.
*          ON EXIT, THE SOLUTION MATRIX X.
*
*  LDB     (INPUT) INTEGER
*          THE LEADING DIMENSION OF THE ARRAY B.  LDB >= MAX(1,N).
*
*  INFO    (OUTPUT) INTEGER
*          = 0:  SUCCESSFUL EXIT
*          < 0:  IF INFO = -I, THE I-TH ARGUMENT HAD AN ILLEGAL VALUE
*
*  =====================================================================
*
*     .. PARAMETERS ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. LOCAL SCALARS ..
      LOGICAL            UPPER
*     ..
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           STRSM, XERBLA
*     ..
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      UPPER = LSAME( UPLO, 'U' )
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SPOTRS', -INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE
*
      IF( N.EQ.0 .OR. NRHS.EQ.0 )
     $   RETURN
*
      IF( UPPER ) THEN
*
*        SOLVE A*X = B WHERE A = U'*U.
*
*        SOLVE U'*X = B, OVERWRITING B WITH X.
*
         CALL STRSM( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT', N, NRHS,
     $               ONE, A, LDA, B, LDB )
*
*        SOLVE U*X = B, OVERWRITING B WITH X.
*
         CALL STRSM( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', N,
     $               NRHS, ONE, A, LDA, B, LDB )
      ELSE
*
*        SOLVE A*X = B WHERE A = L*L'.
*
*        SOLVE L*X = B, OVERWRITING B WITH X.
*
         CALL STRSM( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', N,
     $               NRHS, ONE, A, LDA, B, LDB )
*
*        SOLVE L'*X = B, OVERWRITING B WITH X.
*
         CALL STRSM( 'LEFT', 'LOWER', 'TRANSPOSE', 'NON-UNIT', N, NRHS,
     $               ONE, A, LDA, B, LDB )
      END IF
*
      RETURN
*
*     END OF SPOTRS
*
      END
      INTEGER          FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
     $                 N4 )
*
*  -- LAPACK AUXILIARY ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     SEPTEMBER 30, 1994
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER*( * )    NAME, OPTS
      INTEGER            ISPEC, N1, N2, N3, N4
*     ..
*
*  PURPOSE
*  =======
*
*  ILAENV IS CALLED FROM THE LAPACK ROUTINES TO CHOOSE PROBLEM-DEPENDENT
*  PARAMETERS FOR THE LOCAL ENVIRONMENT.  SEE ISPEC FOR A DESCRIPTION OF
*  THE PARAMETERS.
*
*  THIS VERSION PROVIDES A SET OF PARAMETERS WHICH SHOULD GIVE GOOD,
*  BUT NOT OPTIMAL, PERFORMANCE ON MANY OF THE CURRENTLY AVAILABLE
*  COMPUTERS.  USERS ARE ENCOURAGED TO MODIFY THIS SUBROUTINE TO SET
*  THE TUNING PARAMETERS FOR THEIR PARTICULAR MACHINE USING THE OPTION
*  AND PROBLEM SIZE INFORMATION IN THE ARGUMENTS.
*
*  THIS ROUTINE WILL NOT FUNCTION CORRECTLY IF IT IS CONVERTED TO ALL
*  LOWER CASE.  CONVERTING IT TO ALL UPPER CASE IS ALLOWED.
*
*  ARGUMENTS
*  =========
*
*  ISPEC   (INPUT) INTEGER
*          SPECIFIES THE PARAMETER TO BE RETURNED AS THE VALUE OF
*          ILAENV.
*          = 1: THE OPTIMAL BLOCKSIZE; IF THIS VALUE IS 1, AN UNBLOCKED
*               ALGORITHM WILL GIVE THE BEST PERFORMANCE.
*          = 2: THE MINIMUM BLOCK SIZE FOR WHICH THE BLOCK ROUTINE
*               SHOULD BE USED; IF THE USABLE BLOCK SIZE IS LESS THAN
*               THIS VALUE, AN UNBLOCKED ROUTINE SHOULD BE USED.
*          = 3: THE CROSSOVER POINT (IN A BLOCK ROUTINE, FOR N LESS
*               THAN THIS VALUE, AN UNBLOCKED ROUTINE SHOULD BE USED)
*          = 4: THE NUMBER OF SHIFTS, USED IN THE NONSYMMETRIC
*               EIGENVALUE ROUTINES
*          = 5: THE MINIMUM COLUMN DIMENSION FOR BLOCKING TO BE USED;
*               RECTANGULAR BLOCKS MUST HAVE DIMENSION AT LEAST K BY M,
*               WHERE K IS GIVEN BY ILAENV(2,...) AND M BY ILAENV(5,...)
*          = 6: THE CROSSOVER POINT FOR THE SVD (WHEN REDUCING AN M BY N
*               MATRIX TO BIDIAGONAL FORM, IF MAX(M,N)/MIN(M,N) EXCEEDS
*               THIS VALUE, A QR FACTORIZATION IS USED FIRST TO REDUCE
*               THE MATRIX TO A TRIANGULAR FORM.)
*          = 7: THE NUMBER OF PROCESSORS
*          = 8: THE CROSSOVER POINT FOR THE MULTISHIFT QR AND QZ METHODS
*               FOR NONSYMMETRIC EIGENVALUE PROBLEMS.
*
*  NAME    (INPUT) CHARACTER*(*)
*          THE NAME OF THE CALLING SUBROUTINE, IN EITHER UPPER CASE OR
*          LOWER CASE.
*
*  OPTS    (INPUT) CHARACTER*(*)
*          THE CHARACTER OPTIONS TO THE SUBROUTINE NAME, CONCATENATED
*          INTO A SINGLE CHARACTER STRING.  FOR EXAMPLE, UPLO = 'U',
*          TRANS = 'T', AND DIAG = 'N' FOR A TRIANGULAR ROUTINE WOULD
*          BE SPECIFIED AS OPTS = 'UTN'.
*
*  N1      (INPUT) INTEGER
*  N2      (INPUT) INTEGER
*  N3      (INPUT) INTEGER
*  N4      (INPUT) INTEGER
*          PROBLEM DIMENSIONS FOR THE SUBROUTINE NAME; THESE MAY NOT ALL
*          BE REQUIRED.
*
* (ILAENV) (OUTPUT) INTEGER
*          >= 0: THE VALUE OF THE PARAMETER SPECIFIED BY ISPEC
*          < 0:  IF ILAENV = -K, THE K-TH ARGUMENT HAD AN ILLEGAL VALUE.
*
*  FURTHER DETAILS
*  ===============
*
*  THE FOLLOWING CONVENTIONS HAVE BEEN USED WHEN CALLING ILAENV FROM THE
*  LAPACK ROUTINES:
*  1)  OPTS IS A CONCATENATION OF ALL OF THE CHARACTER OPTIONS TO
*      SUBROUTINE NAME, IN THE SAME ORDER THAT THEY APPEAR IN THE
*      ARGUMENT LIST FOR NAME, EVEN IF THEY ARE NOT USED IN DETERMINING
*      THE VALUE OF THE PARAMETER SPECIFIED BY ISPEC.
*  2)  THE PROBLEM DIMENSIONS N1, N2, N3, N4 ARE SPECIFIED IN THE ORDER
*      THAT THEY APPEAR IN THE ARGUMENT LIST FOR NAME.  N1 IS USED
*      FIRST, N2 SECOND, AND SO ON, AND UNUSED PROBLEM DIMENSIONS ARE
*      PASSED A VALUE OF -1.
*  3)  THE PARAMETER VALUE RETURNED BY ILAENV IS CHECKED FOR VALIDITY IN
*      THE CALLING SUBROUTINE.  FOR EXAMPLE, ILAENV IS USED TO RETRIEVE
*      THE OPTIMAL BLOCKSIZE FOR STRTRI AS FOLLOWS:
*
*      NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
*      IF( NB.LE.1 ) NB = MAX( 1, N )
*
*  =====================================================================
*
*     .. LOCAL SCALARS ..
      LOGICAL            CNAME, SNAME
      CHARACTER*1        C1
      CHARACTER*2        C2, C4
      CHARACTER*3        C3
      CHARACTER*6        SUBNAM
      INTEGER            I, IC, IZ, NB, NBMIN, NX
*     ..
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          CHAR, ICHAR, INT, MIN, REAL
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
      GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC
*
*     INVALID VALUE FOR ISPEC
*
      ILAENV = -1
      RETURN
*
  100 CONTINUE
*
*     CONVERT NAME TO UPPER CASE IF THE FIRST CHARACTER IS LOWER CASE.
*
      ILAENV = 1
      SUBNAM = NAME
      IC = ICHAR( SUBNAM( 1:1 ) )
      IZ = ICHAR( 'Z' )
      IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
*
*        ASCII CHARACTER SET
*
         IF( IC.GE.97 .AND. IC.LE.122 ) THEN
            SUBNAM( 1:1 ) = CHAR( IC-32 )
            DO 10 I = 2, 6
               IC = ICHAR( SUBNAM( I:I ) )
               IF( IC.GE.97 .AND. IC.LE.122 )
     $            SUBNAM( I:I ) = CHAR( IC-32 )
   10       CONTINUE
         END IF
*
      ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
*
*        EBCDIC CHARACTER SET
*
         IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
     $       ( IC.GE.145 .AND. IC.LE.153 ) .OR.
     $       ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
            SUBNAM( 1:1 ) = CHAR( IC+64 )
            DO 20 I = 2, 6
               IC = ICHAR( SUBNAM( I:I ) )
               IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
     $             ( IC.GE.145 .AND. IC.LE.153 ) .OR.
     $             ( IC.GE.162 .AND. IC.LE.169 ) )
     $            SUBNAM( I:I ) = CHAR( IC+64 )
   20       CONTINUE
         END IF
*
      ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
*
*        PRIME MACHINES:  ASCII+128
*
         IF( IC.GE.225 .AND. IC.LE.250 ) THEN
            SUBNAM( 1:1 ) = CHAR( IC-32 )
            DO 30 I = 2, 6
               IC = ICHAR( SUBNAM( I:I ) )
               IF( IC.GE.225 .AND. IC.LE.250 )
     $            SUBNAM( I:I ) = CHAR( IC-32 )
   30       CONTINUE
         END IF
      END IF
*
      C1 = SUBNAM( 1:1 )
      SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
      CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
      IF( .NOT.( CNAME .OR. SNAME ) )
     $   RETURN
      C2 = SUBNAM( 2:3 )
      C3 = SUBNAM( 4:6 )
      C4 = C3( 2:3 )
*
      GO TO ( 110, 200, 300 ) ISPEC
*
  110 CONTINUE
*
*     ISPEC = 1:  BLOCK SIZE
*
*     IN THESE EXAMPLES, SEPARATE CODE IS PROVIDED FOR SETTING NB FOR
*     REAL AND COMPLEX.  WE ASSUME THAT NB WILL TAKE THE SAME VALUE IN
*     SINGLE OR DOUBLE PRECISION.
*
      NB = 1
*
      IF( C2.EQ.'GE' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
     $            C3.EQ.'QLF' ) THEN
            IF( SNAME ) THEN
               NB = 32
            ELSE
               NB = 32
            END IF
         ELSE IF( C3.EQ.'HRD' ) THEN
            IF( SNAME ) THEN
               NB = 32
            ELSE
               NB = 32
            END IF
         ELSE IF( C3.EQ.'BRD' ) THEN
            IF( SNAME ) THEN
               NB = 32
            ELSE
               NB = 32
            END IF
         ELSE IF( C3.EQ.'TRI' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( C2.EQ.'PO' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( C2.EQ.'SY' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
            NB = 1
         ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
            NB = 64
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            NB = 64
         ELSE IF( C3.EQ.'TRD' ) THEN
            NB = 1
         ELSE IF( C3.EQ.'GST' ) THEN
            NB = 64
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NB = 32
            END IF
         END IF
      ELSE IF( C2.EQ.'GB' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               IF( N4.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            ELSE
               IF( N4.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            END IF
         END IF
      ELSE IF( C2.EQ.'PB' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               IF( N2.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            ELSE
               IF( N2.LE.64 ) THEN
                  NB = 1
               ELSE
                  NB = 32
               END IF
            END IF
         END IF
      ELSE IF( C2.EQ.'TR' ) THEN
         IF( C3.EQ.'TRI' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( C2.EQ.'LA' ) THEN
         IF( C3.EQ.'UUM' ) THEN
            IF( SNAME ) THEN
               NB = 64
            ELSE
               NB = 64
            END IF
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
         IF( C3.EQ.'EBZ' ) THEN
            NB = 1
         END IF
      END IF
      ILAENV = NB
      RETURN
*
  200 CONTINUE
*
*     ISPEC = 2:  MINIMUM BLOCK SIZE
*
      NBMIN = 2
      IF( C2.EQ.'GE' ) THEN
         IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
     $       C3.EQ.'QLF' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         ELSE IF( C3.EQ.'HRD' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         ELSE IF( C3.EQ.'BRD' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         ELSE IF( C3.EQ.'TRI' ) THEN
            IF( SNAME ) THEN
               NBMIN = 2
            ELSE
               NBMIN = 2
            END IF
         END IF
      ELSE IF( C2.EQ.'SY' ) THEN
         IF( C3.EQ.'TRF' ) THEN
            IF( SNAME ) THEN
               NBMIN = 8
            ELSE
               NBMIN = 8
            END IF
         ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
            NBMIN = 2
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
         IF( C3.EQ.'TRD' ) THEN
            NBMIN = 2
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NBMIN = 2
            END IF
         END IF
      END IF
      ILAENV = NBMIN
      RETURN
*
  300 CONTINUE
*
*     ISPEC = 3:  CROSSOVER POINT
*
      NX = 0
      IF( C2.EQ.'GE' ) THEN
         IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
     $       C3.EQ.'QLF' ) THEN
            IF( SNAME ) THEN
               NX = 128
            ELSE
               NX = 128
            END IF
         ELSE IF( C3.EQ.'HRD' ) THEN
            IF( SNAME ) THEN
               NX = 128
            ELSE
               NX = 128
            END IF
         ELSE IF( C3.EQ.'BRD' ) THEN
            IF( SNAME ) THEN
               NX = 128
            ELSE
               NX = 128
            END IF
         END IF
      ELSE IF( C2.EQ.'SY' ) THEN
         IF( SNAME .AND. C3.EQ.'TRD' ) THEN
            NX = 1
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
         IF( C3.EQ.'TRD' ) THEN
            NX = 1
         END IF
      ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NX = 128
            END IF
         END IF
      ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
         IF( C3( 1:1 ).EQ.'G' ) THEN
            IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
     $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
     $          C4.EQ.'BR' ) THEN
               NX = 128
            END IF
         END IF
      END IF
      ILAENV = NX
      RETURN
*
  400 CONTINUE
*
*     ISPEC = 4:  NUMBER OF SHIFTS (USED BY XHSEQR)
*
      ILAENV = 6
      RETURN
*
  500 CONTINUE
*
*     ISPEC = 5:  MINIMUM COLUMN DIMENSION (NOT USED)
*
      ILAENV = 2
      RETURN
*
  600 CONTINUE
*
*     ISPEC = 6:  CROSSOVER POINT FOR SVD (USED BY XGELSS AND XGESVD)
*
      ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
      RETURN
*
  700 CONTINUE
*
*     ISPEC = 7:  NUMBER OF PROCESSORS (NOT USED)
*
      ILAENV = 1
      RETURN
*
  800 CONTINUE
*
*     ISPEC = 8:  CROSSOVER POINT FOR MULTISHIFT (USED BY XHSEQR)
*
      ILAENV = 50
      RETURN
*
*     END OF ILAENV
*
      END
      LOGICAL          FUNCTION LSAME( CA, CB )
*
*  -- LAPACK AUXILIARY ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     SEPTEMBER 30, 1994
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER          CA, CB
*     ..
*
*  PURPOSE
*  =======
*
*  LSAME RETURNS .TRUE. IF CA IS THE SAME LETTER AS CB REGARDLESS OF
*  CASE.
*
*  ARGUMENTS
*  =========
*
*  CA      (INPUT) CHARACTER*1
*  CB      (INPUT) CHARACTER*1
*          CA AND CB SPECIFY THE SINGLE CHARACTERS TO BE COMPARED.
*
* =====================================================================
*
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          ICHAR
*     ..
*     .. LOCAL SCALARS ..
      INTEGER            INTA, INTB, ZCODE
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST IF THE CHARACTERS ARE EQUAL
*
      LSAME = CA.EQ.CB
      IF( LSAME )
     $   RETURN
*
*     NOW TEST FOR EQUIVALENCE IF BOTH CHARACTERS ARE ALPHABETIC.
*
      ZCODE = ICHAR( 'Z' )
*
*     USE 'Z' RATHER THAN 'A' SO THAT ASCII CAN BE DETECTED ON PRIME
*     MACHINES, ON WHICH ICHAR RETURNS A VALUE WITH BIT 8 SET.
*     ICHAR('A') ON PRIME MACHINES RETURNS 193 WHICH IS THE SAME AS
*     ICHAR('A') ON AN EBCDIC MACHINE.
*
      INTA = ICHAR( CA )
      INTB = ICHAR( CB )
*
      IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
*
*        ASCII IS ASSUMED - ZCODE IS THE ASCII CODE OF EITHER LOWER OR
*        UPPER CASE 'Z'.
*
         IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
         IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
*
      ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
*
*        EBCDIC IS ASSUMED - ZCODE IS THE EBCDIC CODE OF EITHER LOWER OR
*        UPPER CASE 'Z'.
*
         IF( INTA.GE.129 .AND. INTA.LE.137 .OR.
     $       INTA.GE.145 .AND. INTA.LE.153 .OR.
     $       INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
         IF( INTB.GE.129 .AND. INTB.LE.137 .OR.
     $       INTB.GE.145 .AND. INTB.LE.153 .OR.
     $       INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64
*
      ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
*
*        ASCII IS ASSUMED, ON PRIME MACHINES - ZCODE IS THE ASCII CODE
*        PLUS 128 OF EITHER LOWER OR UPPER CASE 'Z'.
*
         IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
         IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
      END IF
      LSAME = INTA.EQ.INTB
*
*     RETURN
*
*     END OF LSAME
*
      END
      SUBROUTINE XERBLA( SRNAME, INFO )
*
*  -- LAPACK AUXILIARY ROUTINE (VERSION 2.0) --
*     UNIV. OF TENNESSEE, UNIV. OF CALIFORNIA BERKELEY, NAG LTD.,
*     COURANT INSTITUTE, ARGONNE NATIONAL LAB, AND RICE UNIVERSITY
*     SEPTEMBER 30, 1994
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER*6        SRNAME
      INTEGER            INFO
*     ..
*
*  PURPOSE
*  =======
*
*  XERBLA  IS AN ERROR HANDLER FOR THE LAPACK ROUTINES.
*  IT IS CALLED BY AN LAPACK ROUTINE IF AN INPUT PARAMETER HAS AN
*  INVALID VALUE.  A MESSAGE IS PRINTED AND EXECUTION STOPS.
*
*  INSTALLERS MAY CONSIDER MODIFYING THE STOP STATEMENT IN ORDER TO
*  CALL SYSTEM-SPECIFIC EXCEPTION-HANDLING FACILITIES.
*
*  ARGUMENTS
*  =========
*
*  SRNAME  (INPUT) CHARACTER*6
*          THE NAME OF THE ROUTINE WHICH CALLED XERBLA.
*
*  INFO    (INPUT) INTEGER
*          THE POSITION OF THE INVALID PARAMETER IN THE PARAMETER LIST
*          OF THE CALLING ROUTINE.
*
* =====================================================================
*
*     .. EXECUTABLE STATEMENTS ..
*
      WRITE( *, FMT = 9999 )SRNAME, INFO
*
      STOP
*
 9999 FORMAT( ' ** ON ENTRY TO ', A6, ' PARAMETER NUMBER ', I2, ' HAD ',
     $      'AN ILLEGAL VALUE' )
*
*     END OF XERBLA
*
      END

* ======================================================================
* NIST GUIDE TO AVAILABLE MATH SOFTWARE.
* NEEDEDBLAS FOR MODULE SPOSV FROM PACKAGE LAPACK.
* RETRIEVED FROM CAMSUN ON TUE JUN 17 13:37:44 1997.
* ======================================================================

* BLAS REQUIRED BY LAPACK ROUTINE: SPOSV
* -----------------------------------------------------------
* NOTE: LINK TO BLAS OPTIMIZED FOR YOUR SYSTEM, IF AVAILABLE.
* -----------------------------------------------------------

      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C     MODIFIED 12/3/93, ARRAY(1) DECLARATIONS CHANGED TO ARRAY(*)
C
      REAL SX(*),SY(*),STEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      STEMP = 0.0E0
      SDOT = 0.0E0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        STEMP = STEMP + SX(IX)*SY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      SDOT = STEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        STEMP = STEMP + SX(I)*SY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        STEMP = STEMP + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
     *   SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4)
   50 CONTINUE
   60 SDOT = STEMP
      RETURN
      END
      SUBROUTINE SSCAL(N,SA,SX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1.
C     JACK DONGARRA, LINPACK, 3/11/78.
C     MODIFIED 3/93 TO RETURN IF INCX .LE. 0.
C     MODIFIED 12/3/93, ARRAY(1) DECLARATIONS CHANGED TO ARRAY(*)
C
      REAL SA,SX(*)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF( N.LE.0 .OR. INCX.LE.0 )RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        SX(I) = SA*SX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SX(I) = SA*SX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        SX(I) = SA*SX(I)
        SX(I + 1) = SA*SX(I + 1)
        SX(I + 2) = SA*SX(I + 2)
        SX(I + 3) = SA*SX(I + 3)
        SX(I + 4) = SA*SX(I + 4)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE SGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
     $                   BETA, Y, INCY )
*     .. SCALAR ARGUMENTS ..
      REAL               ALPHA, BETA
      INTEGER            INCX, INCY, LDA, M, N
      CHARACTER*1        TRANS
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * ), X( * ), Y( * )
*     ..
*
*  PURPOSE
*  =======
*
*  SGEMV  PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS
*
*     Y := ALPHA*A*X + BETA*Y,   OR   Y := ALPHA*A'*X + BETA*Y,
*
*  WHERE ALPHA AND BETA ARE SCALARS, X AND Y ARE VECTORS AND A IS AN
*  M BY N MATRIX.
*
*  PARAMETERS
*  ==========
*
*  TRANS  - CHARACTER*1.
*           ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS
*           FOLLOWS:
*
*              TRANS = 'N' OR 'N'   Y := ALPHA*A*X + BETA*Y.
*
*              TRANS = 'T' OR 'T'   Y := ALPHA*A'*X + BETA*Y.
*
*              TRANS = 'C' OR 'C'   Y := ALPHA*A'*X + BETA*Y.
*
*           UNCHANGED ON EXIT.
*
*  M      - INTEGER.
*           ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX A.
*           M MUST BE AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX A.
*           N MUST BE AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - REAL            .
*           ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA.
*           UNCHANGED ON EXIT.
*
*  A      - REAL             ARRAY OF DIMENSION ( LDA, N ).
*           BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY A MUST
*           CONTAIN THE MATRIX OF COEFFICIENTS.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. LDA MUST BE AT LEAST
*           MAX( 1, M ).
*           UNCHANGED ON EXIT.
*
*  X      - REAL             ARRAY OF DIMENSION AT LEAST
*           ( 1 + ( N - 1 )*ABS( INCX ) ) WHEN TRANS = 'N' OR 'N'
*           AND AT LEAST
*           ( 1 + ( M - 1 )*ABS( INCX ) ) OTHERWISE.
*           BEFORE ENTRY, THE INCREMENTED ARRAY X MUST CONTAIN THE
*           VECTOR X.
*           UNCHANGED ON EXIT.
*
*  INCX   - INTEGER.
*           ON ENTRY, INCX SPECIFIES THE INCREMENT FOR THE ELEMENTS OF
*           X. INCX MUST NOT BE ZERO.
*           UNCHANGED ON EXIT.
*
*  BETA   - REAL            .
*           ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS
*           SUPPLIED AS ZERO THEN Y NEED NOT BE SET ON INPUT.
*           UNCHANGED ON EXIT.
*
*  Y      - REAL             ARRAY OF DIMENSION AT LEAST
*           ( 1 + ( M - 1 )*ABS( INCY ) ) WHEN TRANS = 'N' OR 'N'
*           AND AT LEAST
*           ( 1 + ( N - 1 )*ABS( INCY ) ) OTHERWISE.
*           BEFORE ENTRY WITH BETA NON-ZERO, THE INCREMENTED ARRAY Y
*           MUST CONTAIN THE VECTOR Y. ON EXIT, Y IS OVERWRITTEN BY THE
*           UPDATED VECTOR Y.
*
*  INCY   - INTEGER.
*           ON ENTRY, INCY SPECIFIES THE INCREMENT FOR THE ELEMENTS OF
*           Y. INCY MUST NOT BE ZERO.
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 2 BLAS ROUTINE.
*
*  -- WRITTEN ON 22-OCTOBER-1986.
*     JACK DONGARRA, ARGONNE NATIONAL LAB.
*     JEREMY DU CROZ, NAG CENTRAL OFFICE.
*     SVEN HAMMARLING, NAG CENTRAL OFFICE.
*     RICHARD HANSON, SANDIA NATIONAL LABS.
*
*
*     .. PARAMETERS ..
      REAL               ONE         , ZERO
      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     .. LOCAL SCALARS ..
      REAL               TEMP
      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
     $         .NOT.LSAME( TRANS, 'T' ).AND.
     $         .NOT.LSAME( TRANS, 'C' )      )THEN
         INFO = 1
      ELSE IF( M.LT.0 )THEN
         INFO = 2
      ELSE IF( N.LT.0 )THEN
         INFO = 3
      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
         INFO = 6
      ELSE IF( INCX.EQ.0 )THEN
         INFO = 8
      ELSE IF( INCY.EQ.0 )THEN
         INFO = 11
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'SGEMV ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     SET  LENX  AND  LENY, THE LENGTHS OF THE VECTORS X AND Y, AND SET
*     UP THE START POINTS IN  X  AND  Y.
*
      IF( LSAME( TRANS, 'N' ) )THEN
         LENX = N
         LENY = M
      ELSE
         LENX = M
         LENY = N
      END IF
      IF( INCX.GT.0 )THEN
         KX = 1
      ELSE
         KX = 1 - ( LENX - 1 )*INCX
      END IF
      IF( INCY.GT.0 )THEN
         KY = 1
      ELSE
         KY = 1 - ( LENY - 1 )*INCY
      END IF
*
*     START THE OPERATIONS. IN THIS VERSION THE ELEMENTS OF A ARE
*     ACCESSED SEQUENTIALLY WITH ONE PASS THROUGH A.
*
*     FIRST FORM  Y := BETA*Y.
*
      IF( BETA.NE.ONE )THEN
         IF( INCY.EQ.1 )THEN
            IF( BETA.EQ.ZERO )THEN
               DO 10, I = 1, LENY
                  Y( I ) = ZERO
   10          CONTINUE
            ELSE
               DO 20, I = 1, LENY
                  Y( I ) = BETA*Y( I )
   20          CONTINUE
            END IF
         ELSE
            IY = KY
            IF( BETA.EQ.ZERO )THEN
               DO 30, I = 1, LENY
                  Y( IY ) = ZERO
                  IY      = IY   + INCY
   30          CONTINUE
            ELSE
               DO 40, I = 1, LENY
                  Y( IY ) = BETA*Y( IY )
                  IY      = IY           + INCY
   40          CONTINUE
            END IF
         END IF
      END IF
      IF( ALPHA.EQ.ZERO )
     $   RETURN
      IF( LSAME( TRANS, 'N' ) )THEN
*
*        FORM  Y := ALPHA*A*X + Y.
*
         JX = KX
         IF( INCY.EQ.1 )THEN
            DO 60, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  DO 50, I = 1, M
                     Y( I ) = Y( I ) + TEMP*A( I, J )
   50             CONTINUE
               END IF
               JX = JX + INCX
   60       CONTINUE
         ELSE
            DO 80, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  IY   = KY
                  DO 70, I = 1, M
                     Y( IY ) = Y( IY ) + TEMP*A( I, J )
                     IY      = IY      + INCY
   70             CONTINUE
               END IF
               JX = JX + INCX
   80       CONTINUE
         END IF
      ELSE
*
*        FORM  Y := ALPHA*A'*X + Y.
*
         JY = KY
         IF( INCX.EQ.1 )THEN
            DO 100, J = 1, N
               TEMP = ZERO
               DO 90, I = 1, M
                  TEMP = TEMP + A( I, J )*X( I )
   90          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  100       CONTINUE
         ELSE
            DO 120, J = 1, N
               TEMP = ZERO
               IX   = KX
               DO 110, I = 1, M
                  TEMP = TEMP + A( I, J )*X( IX )
                  IX   = IX   + INCX
  110          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  120       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     END OF SGEMV .
*
      END
      SUBROUTINE SGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
     $                   BETA, C, LDC )
*     .. SCALAR ARGUMENTS ..
      CHARACTER*1        TRANSA, TRANSB
      INTEGER            M, N, K, LDA, LDB, LDC
      REAL               ALPHA, BETA
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * ), B( LDB, * ), C( LDC, * )
*     ..
*
*  PURPOSE
*  =======
*
*  SGEMM  PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS
*
*     C := ALPHA*OP( A )*OP( B ) + BETA*C,
*
*  WHERE  OP( X ) IS ONE OF
*
*     OP( X ) = X   OR   OP( X ) = X',
*
*  ALPHA AND BETA ARE SCALARS, AND A, B AND C ARE MATRICES, WITH OP( A )
*  AN M BY K MATRIX,  OP( B )  A  K BY N MATRIX AND  C AN M BY N MATRIX.
*
*  PARAMETERS
*  ==========
*
*  TRANSA - CHARACTER*1.
*           ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSA = 'N' OR 'N',  OP( A ) = A.
*
*              TRANSA = 'T' OR 'T',  OP( A ) = A'.
*
*              TRANSA = 'C' OR 'C',  OP( A ) = A'.
*
*           UNCHANGED ON EXIT.
*
*  TRANSB - CHARACTER*1.
*           ON ENTRY, TRANSB SPECIFIES THE FORM OF OP( B ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSB = 'N' OR 'N',  OP( B ) = B.
*
*              TRANSB = 'T' OR 'T',  OP( B ) = B'.
*
*              TRANSB = 'C' OR 'C',  OP( B ) = B'.
*
*           UNCHANGED ON EXIT.
*
*  M      - INTEGER.
*           ON ENTRY,  M  SPECIFIES  THE NUMBER  OF ROWS  OF THE  MATRIX
*           OP( A )  AND OF THE  MATRIX  C.  M  MUST  BE AT LEAST  ZERO.
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY,  N  SPECIFIES THE NUMBER  OF COLUMNS OF THE MATRIX
*           OP( B ) AND THE NUMBER OF COLUMNS OF THE MATRIX C. N MUST BE
*           AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  K      - INTEGER.
*           ON ENTRY,  K  SPECIFIES  THE NUMBER OF COLUMNS OF THE MATRIX
*           OP( A ) AND THE NUMBER OF ROWS OF THE MATRIX OP( B ). K MUST
*           BE AT LEAST  ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - REAL            .
*           ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA.
*           UNCHANGED ON EXIT.
*
*  A      - REAL             ARRAY OF DIMENSION ( LDA, KA ), WHERE KA IS
*           K  WHEN  TRANSA = 'N' OR 'N',  AND IS  M  OTHERWISE.
*           BEFORE ENTRY WITH  TRANSA = 'N' OR 'N',  THE LEADING  M BY K
*           PART OF THE ARRAY  A  MUST CONTAIN THE MATRIX  A,  OTHERWISE
*           THE LEADING  K BY M  PART OF THE ARRAY  A  MUST CONTAIN  THE
*           MATRIX A.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. WHEN  TRANSA = 'N' OR 'N' THEN
*           LDA MUST BE AT LEAST  MAX( 1, M ), OTHERWISE  LDA MUST BE AT
*           LEAST  MAX( 1, K ).
*           UNCHANGED ON EXIT.
*
*  B      - REAL             ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS
*           N  WHEN  TRANSB = 'N' OR 'N',  AND IS  K  OTHERWISE.
*           BEFORE ENTRY WITH  TRANSB = 'N' OR 'N',  THE LEADING  K BY N
*           PART OF THE ARRAY  B  MUST CONTAIN THE MATRIX  B,  OTHERWISE
*           THE LEADING  N BY K  PART OF THE ARRAY  B  MUST CONTAIN  THE
*           MATRIX B.
*           UNCHANGED ON EXIT.
*
*  LDB    - INTEGER.
*           ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. WHEN  TRANSB = 'N' OR 'N' THEN
*           LDB MUST BE AT LEAST  MAX( 1, K ), OTHERWISE  LDB MUST BE AT
*           LEAST  MAX( 1, N ).
*           UNCHANGED ON EXIT.
*
*  BETA   - REAL            .
*           ON ENTRY,  BETA  SPECIFIES THE SCALAR  BETA.  WHEN  BETA  IS
*           SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT.
*           UNCHANGED ON EXIT.
*
*  C      - REAL             ARRAY OF DIMENSION ( LDC, N ).
*           BEFORE ENTRY, THE LEADING  M BY N  PART OF THE ARRAY  C MUST
*           CONTAIN THE MATRIX  C,  EXCEPT WHEN  BETA  IS ZERO, IN WHICH
*           CASE C NEED NOT BE SET ON ENTRY.
*           ON EXIT, THE ARRAY  C  IS OVERWRITTEN BY THE  M BY N  MATRIX
*           ( ALPHA*OP( A )*OP( B ) + BETA*C ).
*
*  LDC    - INTEGER.
*           ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED
*           IN  THE  CALLING  (SUB)  PROGRAM.   LDC  MUST  BE  AT  LEAST
*           MAX( 1, M ).
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 3 BLAS ROUTINE.
*
*  -- WRITTEN ON 8-FEBRUARY-1989.
*     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
*     IAIN DUFF, AERE HARWELL.
*     JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD.
*     SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD.
*
*
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     .. LOCAL SCALARS ..
      LOGICAL            NOTA, NOTB
      INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
      REAL               TEMP
*     .. PARAMETERS ..
      REAL               ONE         , ZERO
      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     SET  NOTA  AND  NOTB  AS  TRUE IF  A  AND  B  RESPECTIVELY ARE NOT
*     TRANSPOSED AND SET  NROWA, NCOLA AND  NROWB  AS THE NUMBER OF ROWS
*     AND  COLUMNS OF  A  AND THE  NUMBER OF  ROWS  OF  B  RESPECTIVELY.
*
      NOTA  = LSAME( TRANSA, 'N' )
      NOTB  = LSAME( TRANSB, 'N' )
      IF( NOTA )THEN
         NROWA = M
         NCOLA = K
      ELSE
         NROWA = K
         NCOLA = M
      END IF
      IF( NOTB )THEN
         NROWB = K
      ELSE
         NROWB = N
      END IF
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      IF(      ( .NOT.NOTA                 ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
         INFO = 1
      ELSE IF( ( .NOT.NOTB                 ).AND.
     $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
         INFO = 2
      ELSE IF( M  .LT.0               )THEN
         INFO = 3
      ELSE IF( N  .LT.0               )THEN
         INFO = 4
      ELSE IF( K  .LT.0               )THEN
         INFO = 5
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
         INFO = 8
      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
         INFO = 10
      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
         INFO = 13
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'SGEMM ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     AND IF  ALPHA.EQ.ZERO.
*
      IF( ALPHA.EQ.ZERO )THEN
         IF( BETA.EQ.ZERO )THEN
            DO 20, J = 1, N
               DO 10, I = 1, M
                  C( I, J ) = ZERO
   10          CONTINUE
   20       CONTINUE
         ELSE
            DO 40, J = 1, N
               DO 30, I = 1, M
                  C( I, J ) = BETA*C( I, J )
   30          CONTINUE
   40       CONTINUE
         END IF
         RETURN
      END IF
*
*     START THE OPERATIONS.
*
      IF( NOTB )THEN
         IF( NOTA )THEN
*
*           FORM  C := ALPHA*A*B + BETA*C.
*
            DO 90, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 50, I = 1, M
                     C( I, J ) = ZERO
   50             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 60, I = 1, M
                     C( I, J ) = BETA*C( I, J )
   60             CONTINUE
               END IF
               DO 80, L = 1, K
                  IF( B( L, J ).NE.ZERO )THEN
                     TEMP = ALPHA*B( L, J )
                     DO 70, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
   70                CONTINUE
                  END IF
   80          CONTINUE
   90       CONTINUE
         ELSE
*
*           FORM  C := ALPHA*A'*B + BETA*C
*
            DO 120, J = 1, N
               DO 110, I = 1, M
                  TEMP = ZERO
                  DO 100, L = 1, K
                     TEMP = TEMP + A( L, I )*B( L, J )
  100             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  110          CONTINUE
  120       CONTINUE
         END IF
      ELSE
         IF( NOTA )THEN
*
*           FORM  C := ALPHA*A*B' + BETA*C
*
            DO 170, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 130, I = 1, M
                     C( I, J ) = ZERO
  130             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 140, I = 1, M
                     C( I, J ) = BETA*C( I, J )
  140             CONTINUE
               END IF
               DO 160, L = 1, K
                  IF( B( J, L ).NE.ZERO )THEN
                     TEMP = ALPHA*B( J, L )
                     DO 150, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
  150                CONTINUE
                  END IF
  160          CONTINUE
  170       CONTINUE
         ELSE
*
*           FORM  C := ALPHA*A'*B' + BETA*C
*
            DO 200, J = 1, N
               DO 190, I = 1, M
                  TEMP = ZERO
                  DO 180, L = 1, K
                     TEMP = TEMP + A( L, I )*B( J, L )
  180             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  190          CONTINUE
  200       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     END OF SGEMM .
*
      END
      SUBROUTINE SSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
     $                   BETA, C, LDC )
*     .. SCALAR ARGUMENTS ..
      CHARACTER*1        UPLO, TRANS
      INTEGER            N, K, LDA, LDC
      REAL               ALPHA, BETA
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * ), C( LDC, * )
*     ..
*
*  PURPOSE
*  =======
*
*  SSYRK  PERFORMS ONE OF THE SYMMETRIC RANK K OPERATIONS
*
*     C := ALPHA*A*A' + BETA*C,
*
*  OR
*
*     C := ALPHA*A'*A + BETA*C,
*
*  WHERE  ALPHA AND BETA  ARE SCALARS, C IS AN  N BY N  SYMMETRIC MATRIX
*  AND  A  IS AN  N BY K  MATRIX IN THE FIRST CASE AND A  K BY N  MATRIX
*  IN THE SECOND CASE.
*
*  PARAMETERS
*  ==========
*
*  UPLO   - CHARACTER*1.
*           ON  ENTRY,   UPLO  SPECIFIES  WHETHER  THE  UPPER  OR  LOWER
*           TRIANGULAR  PART  OF THE  ARRAY  C  IS TO BE  REFERENCED  AS
*           FOLLOWS:
*
*              UPLO = 'U' OR 'U'   ONLY THE  UPPER TRIANGULAR PART OF  C
*                                  IS TO BE REFERENCED.
*
*              UPLO = 'L' OR 'L'   ONLY THE  LOWER TRIANGULAR PART OF  C
*                                  IS TO BE REFERENCED.
*
*           UNCHANGED ON EXIT.
*
*  TRANS  - CHARACTER*1.
*           ON ENTRY,  TRANS  SPECIFIES THE OPERATION TO BE PERFORMED AS
*           FOLLOWS:
*
*              TRANS = 'N' OR 'N'   C := ALPHA*A*A' + BETA*C.
*
*              TRANS = 'T' OR 'T'   C := ALPHA*A'*A + BETA*C.
*
*              TRANS = 'C' OR 'C'   C := ALPHA*A'*A + BETA*C.
*
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY,  N SPECIFIES THE ORDER OF THE MATRIX C.  N MUST BE
*           AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  K      - INTEGER.
*           ON ENTRY WITH  TRANS = 'N' OR 'N',  K  SPECIFIES  THE NUMBER
*           OF  COLUMNS   OF  THE   MATRIX   A,   AND  ON   ENTRY   WITH
*           TRANS = 'T' OR 'T' OR 'C' OR 'C',  K  SPECIFIES  THE  NUMBER
*           OF ROWS OF THE MATRIX  A.  K MUST BE AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - REAL            .
*           ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA.
*           UNCHANGED ON EXIT.
*
*  A      - REAL             ARRAY OF DIMENSION ( LDA, KA ), WHERE KA IS
*           K  WHEN  TRANS = 'N' OR 'N',  AND IS  N  OTHERWISE.
*           BEFORE ENTRY WITH  TRANS = 'N' OR 'N',  THE  LEADING  N BY K
*           PART OF THE ARRAY  A  MUST CONTAIN THE MATRIX  A,  OTHERWISE
*           THE LEADING  K BY N  PART OF THE ARRAY  A  MUST CONTAIN  THE
*           MATRIX A.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN  THE  CALLING  (SUB)  PROGRAM.   WHEN  TRANS = 'N' OR 'N'
*           THEN  LDA MUST BE AT LEAST  MAX( 1, N ), OTHERWISE  LDA MUST
*           BE AT LEAST  MAX( 1, K ).
*           UNCHANGED ON EXIT.
*
*  BETA   - REAL            .
*           ON ENTRY, BETA SPECIFIES THE SCALAR BETA.
*           UNCHANGED ON EXIT.
*
*  C      - REAL             ARRAY OF DIMENSION ( LDC, N ).
*           BEFORE ENTRY  WITH  UPLO = 'U' OR 'U',  THE LEADING  N BY N
*           UPPER TRIANGULAR PART OF THE ARRAY C MUST CONTAIN THE UPPER
*           TRIANGULAR PART  OF THE  SYMMETRIC MATRIX  AND THE STRICTLY
*           LOWER TRIANGULAR PART OF C IS NOT REFERENCED.  ON EXIT, THE
*           UPPER TRIANGULAR PART OF THE ARRAY  C IS OVERWRITTEN BY THE
*           UPPER TRIANGULAR PART OF THE UPDATED MATRIX.
*           BEFORE ENTRY  WITH  UPLO = 'L' OR 'L',  THE LEADING  N BY N
*           LOWER TRIANGULAR PART OF THE ARRAY C MUST CONTAIN THE LOWER
*           TRIANGULAR PART  OF THE  SYMMETRIC MATRIX  AND THE STRICTLY
*           UPPER TRIANGULAR PART OF C IS NOT REFERENCED.  ON EXIT, THE
*           LOWER TRIANGULAR PART OF THE ARRAY  C IS OVERWRITTEN BY THE
*           LOWER TRIANGULAR PART OF THE UPDATED MATRIX.
*
*  LDC    - INTEGER.
*           ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED
*           IN  THE  CALLING  (SUB)  PROGRAM.   LDC  MUST  BE  AT  LEAST
*           MAX( 1, N ).
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 3 BLAS ROUTINE.
*
*  -- WRITTEN ON 8-FEBRUARY-1989.
*     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
*     IAIN DUFF, AERE HARWELL.
*     JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD.
*     SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD.
*
*
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     .. LOCAL SCALARS ..
      LOGICAL            UPPER
      INTEGER            I, INFO, J, L, NROWA
      REAL               TEMP
*     .. PARAMETERS ..
      REAL               ONE ,         ZERO
      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      IF( LSAME( TRANS, 'N' ) )THEN
         NROWA = N
      ELSE
         NROWA = K
      END IF
      UPPER = LSAME( UPLO, 'U' )
*
      INFO = 0
      IF(      ( .NOT.UPPER               ).AND.
     $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
         INFO = 1
      ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
     $         ( .NOT.LSAME( TRANS, 'T' ) ).AND.
     $         ( .NOT.LSAME( TRANS, 'C' ) )      )THEN
         INFO = 2
      ELSE IF( N  .LT.0               )THEN
         INFO = 3
      ELSE IF( K  .LT.0               )THEN
         INFO = 4
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
         INFO = 7
      ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
         INFO = 10
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'SSYRK ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( ( N.EQ.0 ).OR.
     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     AND WHEN  ALPHA.EQ.ZERO.
*
      IF( ALPHA.EQ.ZERO )THEN
         IF( UPPER )THEN
            IF( BETA.EQ.ZERO )THEN
               DO 20, J = 1, N
                  DO 10, I = 1, J
                     C( I, J ) = ZERO
   10             CONTINUE
   20          CONTINUE
            ELSE
               DO 40, J = 1, N
                  DO 30, I = 1, J
                     C( I, J ) = BETA*C( I, J )
   30             CONTINUE
   40          CONTINUE
            END IF
         ELSE
            IF( BETA.EQ.ZERO )THEN
               DO 60, J = 1, N
                  DO 50, I = J, N
                     C( I, J ) = ZERO
   50             CONTINUE
   60          CONTINUE
            ELSE
               DO 80, J = 1, N
                  DO 70, I = J, N
                     C( I, J ) = BETA*C( I, J )
   70             CONTINUE
   80          CONTINUE
            END IF
         END IF
         RETURN
      END IF
*
*     START THE OPERATIONS.
*
      IF( LSAME( TRANS, 'N' ) )THEN
*
*        FORM  C := ALPHA*A*A' + BETA*C.
*
         IF( UPPER )THEN
            DO 130, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 90, I = 1, J
                     C( I, J ) = ZERO
   90             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 100, I = 1, J
                     C( I, J ) = BETA*C( I, J )
  100             CONTINUE
               END IF
               DO 120, L = 1, K
                  IF( A( J, L ).NE.ZERO )THEN
                     TEMP = ALPHA*A( J, L )
                     DO 110, I = 1, J
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
  110                CONTINUE
                  END IF
  120          CONTINUE
  130       CONTINUE
         ELSE
            DO 180, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 140, I = J, N
                     C( I, J ) = ZERO
  140             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 150, I = J, N
                     C( I, J ) = BETA*C( I, J )
  150             CONTINUE
               END IF
               DO 170, L = 1, K
                  IF( A( J, L ).NE.ZERO )THEN
                     TEMP      = ALPHA*A( J, L )
                     DO 160, I = J, N
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
  160                CONTINUE
                  END IF
  170          CONTINUE
  180       CONTINUE
         END IF
      ELSE
*
*        FORM  C := ALPHA*A'*A + BETA*C.
*
         IF( UPPER )THEN
            DO 210, J = 1, N
               DO 200, I = 1, J
                  TEMP = ZERO
                  DO 190, L = 1, K
                     TEMP = TEMP + A( L, I )*A( L, J )
  190             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  200          CONTINUE
  210       CONTINUE
         ELSE
            DO 240, J = 1, N
               DO 230, I = J, N
                  TEMP = ZERO
                  DO 220, L = 1, K
                     TEMP = TEMP + A( L, I )*A( L, J )
  220             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  230          CONTINUE
  240       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     END OF SSYRK .
*
      END
      SUBROUTINE STRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
     $                   B, LDB )
*     .. SCALAR ARGUMENTS ..
      CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
      INTEGER            M, N, LDA, LDB
      REAL               ALPHA
*     .. ARRAY ARGUMENTS ..
      REAL               A( LDA, * ), B( LDB, * )
*     ..
*
*  PURPOSE
*  =======
*
*  STRSM  SOLVES ONE OF THE MATRIX EQUATIONS
*
*     OP( A )*X = ALPHA*B,   OR   X*OP( A ) = ALPHA*B,
*
*  WHERE ALPHA IS A SCALAR, X AND B ARE M BY N MATRICES, A IS A UNIT, OR
*  NON-UNIT,  UPPER OR LOWER TRIANGULAR MATRIX  AND  OP( A )  IS ONE  OF
*
*     OP( A ) = A   OR   OP( A ) = A'.
*
*  THE MATRIX X IS OVERWRITTEN ON B.
*
*  PARAMETERS
*  ==========
*
*  SIDE   - CHARACTER*1.
*           ON ENTRY, SIDE SPECIFIES WHETHER OP( A ) APPEARS ON THE LEFT
*           OR RIGHT OF X AS FOLLOWS:
*
*              SIDE = 'L' OR 'L'   OP( A )*X = ALPHA*B.
*
*              SIDE = 'R' OR 'R'   X*OP( A ) = ALPHA*B.
*
*           UNCHANGED ON EXIT.
*
*  UPLO   - CHARACTER*1.
*           ON ENTRY, UPLO SPECIFIES WHETHER THE MATRIX A IS AN UPPER OR
*           LOWER TRIANGULAR MATRIX AS FOLLOWS:
*
*              UPLO = 'U' OR 'U'   A IS AN UPPER TRIANGULAR MATRIX.
*
*              UPLO = 'L' OR 'L'   A IS A LOWER TRIANGULAR MATRIX.
*
*           UNCHANGED ON EXIT.
*
*  TRANSA - CHARACTER*1.
*           ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSA = 'N' OR 'N'   OP( A ) = A.
*
*              TRANSA = 'T' OR 'T'   OP( A ) = A'.
*
*              TRANSA = 'C' OR 'C'   OP( A ) = A'.
*
*           UNCHANGED ON EXIT.
*
*  DIAG   - CHARACTER*1.
*           ON ENTRY, DIAG SPECIFIES WHETHER OR NOT A IS UNIT TRIANGULAR
*           AS FOLLOWS:
*
*              DIAG = 'U' OR 'U'   A IS ASSUMED TO BE UNIT TRIANGULAR.
*
*              DIAG = 'N' OR 'N'   A IS NOT ASSUMED TO BE UNIT
*                                  TRIANGULAR.
*
*           UNCHANGED ON EXIT.
*
*  M      - INTEGER.
*           ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF B. M MUST BE AT
*           LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF B.  N MUST BE
*           AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - REAL            .
*           ON ENTRY,  ALPHA SPECIFIES THE SCALAR  ALPHA. WHEN  ALPHA IS
*           ZERO THEN  A IS NOT REFERENCED AND  B NEED NOT BE SET BEFORE
*           ENTRY.
*           UNCHANGED ON EXIT.
*
*  A      - REAL             ARRAY OF DIMENSION ( LDA, K ), WHERE K IS M
*           WHEN  SIDE = 'L' OR 'L'  AND IS  N  WHEN  SIDE = 'R' OR 'R'.
*           BEFORE ENTRY  WITH  UPLO = 'U' OR 'U',  THE  LEADING  K BY K
*           UPPER TRIANGULAR PART OF THE ARRAY  A MUST CONTAIN THE UPPER
*           TRIANGULAR MATRIX  AND THE STRICTLY LOWER TRIANGULAR PART OF
*           A IS NOT REFERENCED.
*           BEFORE ENTRY  WITH  UPLO = 'L' OR 'L',  THE  LEADING  K BY K
*           LOWER TRIANGULAR PART OF THE ARRAY  A MUST CONTAIN THE LOWER
*           TRIANGULAR MATRIX  AND THE STRICTLY UPPER TRIANGULAR PART OF
*           A IS NOT REFERENCED.
*           NOTE THAT WHEN  DIAG = 'U' OR 'U',  THE DIAGONAL ELEMENTS OF
*           A  ARE NOT REFERENCED EITHER,  BUT ARE ASSUMED TO BE  UNITY.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN THE CALLING (SUB) PROGRAM.  WHEN  SIDE = 'L' OR 'L'  THEN
*           LDA  MUST BE AT LEAST  MAX( 1, M ),  WHEN  SIDE = 'R' OR 'R'
*           THEN LDA MUST BE AT LEAST MAX( 1, N ).
*           UNCHANGED ON EXIT.
*
*  B      - REAL             ARRAY OF DIMENSION ( LDB, N ).
*           BEFORE ENTRY,  THE LEADING  M BY N PART OF THE ARRAY  B MUST
*           CONTAIN  THE  RIGHT-HAND  SIDE  MATRIX  B,  AND  ON EXIT  IS
*           OVERWRITTEN BY THE SOLUTION MATRIX  X.
*
*  LDB    - INTEGER.
*           ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED
*           IN  THE  CALLING  (SUB)  PROGRAM.   LDB  MUST  BE  AT  LEAST
*           MAX( 1, M ).
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 3 BLAS ROUTINE.
*
*
*  -- WRITTEN ON 8-FEBRUARY-1989.
*     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
*     IAIN DUFF, AERE HARWELL.
*     JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD.
*     SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD.
*
*
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     .. LOCAL SCALARS ..
      LOGICAL            LSIDE, NOUNIT, UPPER
      INTEGER            I, INFO, J, K, NROWA
      REAL               TEMP
*     .. PARAMETERS ..
      REAL               ONE         , ZERO
      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      LSIDE  = LSAME( SIDE  , 'L' )
      IF( LSIDE )THEN
         NROWA = M
      ELSE
         NROWA = N
      END IF
      NOUNIT = LSAME( DIAG  , 'N' )
      UPPER  = LSAME( UPLO  , 'U' )
*
      INFO   = 0
      IF(      ( .NOT.LSIDE                ).AND.
     $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
         INFO = 1
      ELSE IF( ( .NOT.UPPER                ).AND.
     $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
         INFO = 2
      ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
         INFO = 3
      ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
     $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
         INFO = 4
      ELSE IF( M  .LT.0               )THEN
         INFO = 5
      ELSE IF( N  .LT.0               )THEN
         INFO = 6
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
         INFO = 9
      ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
         INFO = 11
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'STRSM ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( N.EQ.0 )
     $   RETURN
*
*     AND WHEN  ALPHA.EQ.ZERO.
*
      IF( ALPHA.EQ.ZERO )THEN
         DO 20, J = 1, N
            DO 10, I = 1, M
               B( I, J ) = ZERO
   10       CONTINUE
   20    CONTINUE
         RETURN
      END IF
*
*     START THE OPERATIONS.
*
      IF( LSIDE )THEN
         IF( LSAME( TRANSA, 'N' ) )THEN
*
*           FORM  B := ALPHA*INV( A )*B.
*
            IF( UPPER )THEN
               DO 60, J = 1, N
                  IF( ALPHA.NE.ONE )THEN
                     DO 30, I = 1, M
                        B( I, J ) = ALPHA*B( I, J )
   30                CONTINUE
                  END IF
                  DO 50, K = M, 1, -1
                     IF( B( K, J ).NE.ZERO )THEN
                        IF( NOUNIT )
     $                     B( K, J ) = B( K, J )/A( K, K )
                        DO 40, I = 1, K - 1
                           B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
   40                   CONTINUE
                     END IF
   50             CONTINUE
   60          CONTINUE
            ELSE
               DO 100, J = 1, N
                  IF( ALPHA.NE.ONE )THEN
                     DO 70, I = 1, M
                        B( I, J ) = ALPHA*B( I, J )
   70                CONTINUE
                  END IF
                  DO 90 K = 1, M
                     IF( B( K, J ).NE.ZERO )THEN
                        IF( NOUNIT )
     $                     B( K, J ) = B( K, J )/A( K, K )
                        DO 80, I = K + 1, M
                           B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
   80                   CONTINUE
                     END IF
   90             CONTINUE
  100          CONTINUE
            END IF
         ELSE
*
*           FORM  B := ALPHA*INV( A' )*B.
*
            IF( UPPER )THEN
               DO 130, J = 1, N
                  DO 120, I = 1, M
                     TEMP = ALPHA*B( I, J )
                     DO 110, K = 1, I - 1
                        TEMP = TEMP - A( K, I )*B( K, J )
  110                CONTINUE
                     IF( NOUNIT )
     $                  TEMP = TEMP/A( I, I )
                     B( I, J ) = TEMP
  120             CONTINUE
  130          CONTINUE
            ELSE
               DO 160, J = 1, N
                  DO 150, I = M, 1, -1
                     TEMP = ALPHA*B( I, J )
                     DO 140, K = I + 1, M
                        TEMP = TEMP - A( K, I )*B( K, J )
  140                CONTINUE
                     IF( NOUNIT )
     $                  TEMP = TEMP/A( I, I )
                     B( I, J ) = TEMP
  150             CONTINUE
  160          CONTINUE
            END IF
         END IF
      ELSE
         IF( LSAME( TRANSA, 'N' ) )THEN
*
*           FORM  B := ALPHA*B*INV( A ).
*
            IF( UPPER )THEN
               DO 210, J = 1, N
                  IF( ALPHA.NE.ONE )THEN
                     DO 170, I = 1, M
                        B( I, J ) = ALPHA*B( I, J )
  170                CONTINUE
                  END IF
                  DO 190, K = 1, J - 1
                     IF( A( K, J ).NE.ZERO )THEN
                        DO 180, I = 1, M
                           B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  180                   CONTINUE
                     END IF
  190             CONTINUE
                  IF( NOUNIT )THEN
                     TEMP = ONE/A( J, J )
                     DO 200, I = 1, M
                        B( I, J ) = TEMP*B( I, J )
  200                CONTINUE
                  END IF
  210          CONTINUE
            ELSE
               DO 260, J = N, 1, -1
                  IF( ALPHA.NE.ONE )THEN
                     DO 220, I = 1, M
                        B( I, J ) = ALPHA*B( I, J )
  220                CONTINUE
                  END IF
                  DO 240, K = J + 1, N
                     IF( A( K, J ).NE.ZERO )THEN
                        DO 230, I = 1, M
                           B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  230                   CONTINUE
                     END IF
  240             CONTINUE
                  IF( NOUNIT )THEN
                     TEMP = ONE/A( J, J )
                     DO 250, I = 1, M
                       B( I, J ) = TEMP*B( I, J )
  250                CONTINUE
                  END IF
  260          CONTINUE
            END IF
         ELSE
*
*           FORM  B := ALPHA*B*INV( A' ).
*
            IF( UPPER )THEN
               DO 310, K = N, 1, -1
                  IF( NOUNIT )THEN
                     TEMP = ONE/A( K, K )
                     DO 270, I = 1, M
                        B( I, K ) = TEMP*B( I, K )
  270                CONTINUE
                  END IF
                  DO 290, J = 1, K - 1
                     IF( A( J, K ).NE.ZERO )THEN
                        TEMP = A( J, K )
                        DO 280, I = 1, M
                           B( I, J ) = B( I, J ) - TEMP*B( I, K )
  280                   CONTINUE
                     END IF
  290             CONTINUE
                  IF( ALPHA.NE.ONE )THEN
                     DO 300, I = 1, M
                        B( I, K ) = ALPHA*B( I, K )
  300                CONTINUE
                  END IF
  310          CONTINUE
            ELSE
               DO 360, K = 1, N
                  IF( NOUNIT )THEN
                     TEMP = ONE/A( K, K )
                     DO 320, I = 1, M
                        B( I, K ) = TEMP*B( I, K )
  320                CONTINUE
                  END IF
                  DO 340, J = K + 1, N
                     IF( A( J, K ).NE.ZERO )THEN
                        TEMP = A( J, K )
                        DO 330, I = 1, M
                           B( I, J ) = B( I, J ) - TEMP*B( I, K )
  330                   CONTINUE
                     END IF
  340             CONTINUE
                  IF( ALPHA.NE.ONE )THEN
                     DO 350, I = 1, M
                        B( I, K ) = ALPHA*B( I, K )
  350                CONTINUE
                  END IF
  360          CONTINUE
            END IF
         END IF
      END IF
*
      RETURN
*
*     END OF STRSM .
*
      END
