
      Subroutine Meshdet(X, X1, Nmsh)

      Implicit Double Precision (A-H,O-Z)

C.... This Subroutine Is Used For Saving And Re-Inserting Meshes

      Dimension X(*), X1(*)

      Do 10, J = 1,Nmsh
         X(J) = X1(J)
 10   Continue

      Return
      End

      Subroutine Soldet(U, U1, Ncomp, Nmsh, Idim1, Idim2)

      Implicit Double Precision (A-H,O-Z)

C.... This Subroutine Is Used For Saving And Re-Inserting Solutions.

      Dimension U(Idim1,*), U1(Idim2,*)

      Do 14, J = 1,Ncomp
         Do 12 K = 1,Nmsh
            U(J,K) = U1(J,K)
 12      Continue
 14   Continue

      Return
      End












      Subroutine Dblmsh (Nmsh, Nmax, Xx, Nmold, Xxold, Maxmsh)
      Implicit Double Precision (A-H,O-Z)
      Dimension Xx(*), Xxold(*)
      Logical Maxmsh
      LOGICAL use_c, comp_c
      Common/Algprs/ Nminit, Iprint, Maxcon, Itsaim, Uval0,use_c,comp_c
      Common /Flags/ Ifinal,Iatt,Iback,Iprec,Iucond
      character (len=150) msg
*  Blas: Dcopy

      Parameter (Half = 0.5d+0)

*  This Routine Is Used To Double The Mesh, I.e., Produce A Mesh
*  With Twice As Many Intervals In Which Each New Interval Is
*  Half The Corresponding Old Interval.

*  On Entry To Dblmsh, The Integer Nmsh And The Array Xx
*  Specify A Set Of Mesh Points Xx(1),..., Xx(Nmsh) (Assumed
*  To Be In Ascending Order).

*  If The Number Of Mesh Points In The Doubled Mesh Would
*  Exceed The Maximum Allowed Number Nmax, The Flag Maxmsh Is
*  Set To True, And We Exit Without Changing Any Other Parameters.

*  Otherwise, Nmold Is Set To The Old Number Of Mesh Points,
*  Xxold Is Set To The Old Mesh, Nmsh Is The New Number Of Mesh
*  Points, And Xx Contains The New Mesh Points.

      Nmold = Nmsh
      Iprec = Min(Iprec,1)
      Call Dcopy(Nmold, Xx, 1, Xxold, 1)

      Ninnew = 2*(Nmsh-1)
      Nmnew = Ninnew + 1
      If(Nmnew .ge. Nmax) Then
         If (Iprint .ge. 0)  then
           Write(msg,901) Nmnew
           call rprint(msg)
         end if
         Maxmsh = .true.
         Return
      Endif
      Maxmsh = .false.

*  Loop Backwards Through The Old Mesh Points To Create The New Ones.

      Xx(Nmnew) = Xx(Nmsh)
      Do 100 I = Ninnew, 4, -2
         Id2 = I/2
         Xx(I) = Half*(Xx(I+1) + Xx(Id2))
         Xx(I-1) = Xx(Id2)
         If (Xx(I) .eq. Xx(I+1) .or. Xx(I) .eq. Xx(I-1)) Then
            Iprec = 2
            Return
         Endif
  100 Continue

*  Calculate The New Xx(2). Xx(1) Remains Unchanged.

      Xx(2) = Half*(Xx(3) + Xx(1))
      If (Xx(2) .eq. Xx(3) .or. Xx(2) .eq. Xx(1)) Then
         Iprec = 2
         Return
      Endif
      Nmsh = Nmnew
      Return
  901 Format (1h , ' Dblmsh.  Maximum Mesh Exceeded, Nmnew  = ', I8)
      End


      Subroutine Unimsh(Nmsh, Aleft, Aright, Nfxpnt, Fixpnt, Xx)
      Implicit Double Precision (A-H,O-Z)
      Integer  Nmsh, Nfxpnt
      Dimension Fixpnt(*), Xx(Nmsh)

      Intrinsic Max

*  Given A Left Endpoint Aleft, A Right Endpoint Aright,
*  A Set Of Nfxpnt Fixed Points Fixpnt(I), I = 1,...,Nfxpnt,
*  (Where Fixpnt(I) Is Different From Aleft And Aright For All I),
*  And An Initial Target Number Nmsh Of Mesh Points,
*  The Subroutine Unimsh Generates A Piecewise Uniform Mesh
*  Beginning At Aleft, Ending At Aright, And With Equally
*  Spaced Points Between Aleft And Fixpnt(1), Then Between
*  Fixpnt(1) And Fixpnt(2), ..., And Finally Between
*  Fixpnt(Nfxpnt) And Aright.  The Final Number Of Intervals
*  Is The Maximum Of Nfxpnt+2 And The Initial Value Of Nmsh.

*  In The Simplest Case When Nfxpnt = 0, Unimsh Generates A
*  Uniform Mesh With Nmsh Intervals In The Closed Interval
*  (Aleft, Aright).

*  On Exit, The Integer Nmsh Contains The Number Of Mesh Points
*  (Which Is The Maximum Of The Initial Nmsh And Nfxpnt).
*  The Array Xx (Of Dimension Nmsh) Contains The Mesh Points.

      If (Nfxpnt .eq. 0) Then

*  If There Are No Interior Fixed Points, The Spacing Is Uniform
*  Throughout The Interval.  Calculate The Spacing Dx
*  And Set Up The Xx Array.

        Ninter = Nmsh - 1

         Dx = (Aright - Aleft)/Ninter
         Do 10 I = 1, Ninter
            Xx(I) = Aleft + (I-1)*Dx
   10    Continue
         Xx(Nmsh) = Aright
         Return
      Endif

*  We Know That There Is At Least One Fixed Point Strictly Between
*  The Endpoints.

      If (Nmsh .lt. Nfxpnt+2)  Nmsh = Nfxpnt + 2
      Ninter = Nmsh - 1
      Xx(1) = Aleft
      Ileft = 1
      Xleft = Aleft
      Totint = Aright - Aleft
      Ndif = Ninter - Nfxpnt
      Do 50 J = 1, Nfxpnt + 1

*  Deal In Turn With The Subintervals Defined By The Interval
*  Boundaries And The Fixed  Points.

         If (J .lt. Nfxpnt+1) Then

*  The J-Th Fixed Point Is Xright.  Calculate Where It Should
*  Fall In The Mesh.

            Xright = Fixpnt(J)
            Nmin = Int(Ninter*(Xright-Aleft)/Totint + 1.5d+0)
            If (Nmin .gt. Ndif+J) Nmin = Ndif + J
            Iright = Max(Ileft+1, Nmin)
         Else
            Xright = Aright
            Iright = Nmsh
         Endif

*  Npt Is The Number Of Equally Spaced Points That Should
*  Lie Strictly Between The (J-1)-Th And J-Th Fixed Points.

         Xx(Iright) = Xright
         Npt = Iright - Ileft - 1
         Dx = (Xright - Xleft)/(Npt + 1)
         Do 30 I = 1, Npt
            Xx(Ileft+I) = Xleft + I*Dx
   30    Continue
         Ileft = Iright
         Xleft = Xright
   50 Continue

      Return
C
      End


* File Rest.f

      Subroutine Errest (Ncomp, Nmsh, Ntol, Ltol, Tol,
     *   Nudim, U, Uold, Etest, Errok)

      Implicit Double Precision (A-H,O-Z)
      Dimension Ltol(Ntol), Tol(Ntol), U(Nudim,Nmsh),
     *              Uold(Ncomp,Nmsh), Etest(Ntol)
      Logical Errok

      Intrinsic Abs
      Intrinsic Max

      Parameter ( One = 1.0d+0 )


*  Given Current And Previous Solutions U And Uold On The Same
*  Mesh, Errest Calculates An Error Measure For Each
*  Component For Which A Tolerance Is Specified.
*  The Error Measure Is The Usual Relative Error Normalized
*  By Dividing By The Tolerance.

*  The Array Etest Specifies The Error Test To Be Applied To Each
*  Error Measure.

*  On Exit, The Logical Flag Errok
*   -- Is .false. If Any Of The Error Measures Exceeds The
*      Corresponding Value In The Array Etest
*   -- Is .true. If All Error Measures Are Less Than The
*      Corresponding Values Of Etest.

      Errok = .true.

      Do 10 Im = 1, Nmsh
      Do 10 It = 1, Ntol
         Icmp = Ltol(It)
         Er = U(Icmp,Im) - Uold(Icmp,Im)
         Denom = Max(One, Abs(Uold(Icmp,Im)))
         Errel = Abs(Er/(Tol(It)*Denom))
         If (Errel .gt. Etest(It)) Errok = .false.
   10 Continue

      Return
      End


      Double Precision Function D1mach(I)
C
C  Double-Precision Machine Constants
C
C  D1mach( 1) = B**(Emin-1), The Smallest Positive Magnitude.
C
C  D1mach( 2) = B**Emax*(1 - B**(-T)), The Largest Magnitude.
C
C  D1mach( 3) = B**(-T), The Smallest Relative Spacing.
C
C  D1mach( 4) = B**(1-T), The Largest Relative Spacing.
C
C  D1mach( 5) = Log10(B)
C
C  To Alter This Function For A Particular Environment,
C  The Desired Set Of Data Statements Should Be Activated By
C  Removing The C From Column 1.
C  On Rare Machines A Static Statement May Need To Be Added.
C  (But Probably More Systems Prohibit It Than Require It.)
C
C  For Ieee-Arithmetic Machines (Binary Standard), One Of The First
C  Two Sets Of Constants Below Should Be Appropriate.  If You Do Not
C  Know Which Set To Use, Try Both And See Which Gives Plausible
C  Values.
C
C  Where Possible, Decimal, Octal Or Hexadecimal Constants Are Used
C  To Specify The Constants Exactly.  Sometimes This Requires Using
C  Equivalent Integer Arrays.  If Your Compiler Uses Half-Word
C  Integers By Default (Sometimes Called Integer*2), You May Need To
C  Change Integer To Integer*4 Or Otherwise Instruct Your Compiler
C  To Use Full-Word Integers In The Next 5 Declarations.
C
C  Comments Just Before The End Statement (Lines Starting With *)
C  Give C Source For D1mach.
C
      Integer Small(2)
      Integer Large(2)
      Integer Right(2)
      Integer Diver(2)
      Integer Log10(2)
      Integer Sc
C
      Double Precision Dmach(5)
C
      Equivalence (Dmach(1),Small(1))
      Equivalence (Dmach(2),Large(1))
      Equivalence (Dmach(3),Right(1))
      Equivalence (Dmach(4),Diver(1))
      Equivalence (Dmach(5),Log10(1))
C
C     Machine Constants For Big-Endian Ieee Arithmetic (Binary Format)
C     Machines In Which The Most Significant Byte Is Stored First,
C     Such As The At&T 3b Series, Motorola 68000 Based Machines (E.g.
C     Sun 3), And Machines That Use Sparc, Hp, Or Ibm Risc Chips.
C
C      Data Small(1),Small(2) /    1048576,          0 /
C      Data Large(1),Large(2) / 2146435071,         -1 /
C      Data Right(1),Right(2) / 1017118720,          0 /
C      Data Diver(1),Diver(2) / 1018167296,          0 /
C      Data Log10(1),Log10(2) / 1070810131, 1352628735 /, Sc/987/
C
C     Machine Constants For Little-Endian (Binary) Ieee Arithmetic
C     Machines In Which The Least Significant Byte Is Stored First,
C     E.g. Ibm Pcs And Other Machines That Use Intel 80x87 Or Dec
C     Alpha Chips.
C
      Data Small(1),Small(2) /          0,    1048576 /
      Data Large(1),Large(2) /         -1, 2146435071 /
      Data Right(1),Right(2) /          0, 1017118720 /
      Data Diver(1),Diver(2) /          0, 1018167296 /
      Data Log10(1),Log10(2) / 1352628735, 1070810131 /, Sc/987/
C
C     Machine Constants For Amdahl Machines.
C
C      Data Small(1),Small(2) /    1048576,          0 /
C      Data Large(1),Large(2) / 2147483647,         -1 /
C      Data Right(1),Right(2) /  856686592,          0 /
C      Data Diver(1),Diver(2) /  873463808,          0 /
C      Data Log10(1),Log10(2) / 1091781651, 1352628735 /, Sc/987/
C
C     Machine Constants For The Burroughs 1700 System.
C
C      Data Small(1) / Zc00800000 /
C      Data Small(2) / Z000000000 /
C
C      Data Large(1) / Zdffffffff /
C      Data Large(2) / Zfffffffff /
C
C      Data Right(1) / Zcc5800000 /
C      Data Right(2) / Z000000000 /
C
C      Data Diver(1) / Zcc6800000 /
C      Data Diver(2) / Z000000000 /
C
C      Data Log10(1) / Zd00e730e7 /
C      Data Log10(2) / Zc77800dc0 /, Sc/987/
C
C     Machine Constants For The Burroughs 5700 System.
C
C      Data Small(1) / O1771000000000000 /
C      Data Small(2) / O0000000000000000 /
C
C      Data Large(1) / O0777777777777777 /
C      Data Large(2) / O0007777777777777 /
C
C      Data Right(1) / O1461000000000000 /
C      Data Right(2) / O0000000000000000 /
C
C      Data Diver(1) / O1451000000000000 /
C      Data Diver(2) / O0000000000000000 /
C
C      Data Log10(1) / O1157163034761674 /
C      Data Log10(2) / O0006677466732724 /, Sc/987/
C
C     Machine Constants For The Burroughs 6700/7700 Systems.
C
C      Data Small(1) / O1771000000000000 /
C      Data Small(2) / O7770000000000000 /
C
C      Data Large(1) / O0777777777777777 /
C      Data Large(2) / O7777777777777777 /
C
C      Data Right(1) / O1461000000000000 /
C      Data Right(2) / O0000000000000000 /
C
C      Data Diver(1) / O1451000000000000 /
C      Data Diver(2) / O0000000000000000 /
C
C      Data Log10(1) / O1157163034761674 /
C      Data Log10(2) / O0006677466732724 /, Sc/987/
C
C     Machine Constants For Ftn4 On The Cdc 6000/7000 Series.
C
C      Data Small(1) / 00564000000000000000b /
C      Data Small(2) / 00000000000000000000b /
C
C      Data Large(1) / 37757777777777777777b /
C      Data Large(2) / 37157777777777777774b /
C
C      Data Right(1) / 15624000000000000000b /
C      Data Right(2) / 00000000000000000000b /
C
C      Data Diver(1) / 15634000000000000000b /
C      Data Diver(2) / 00000000000000000000b /
C
C      Data Log10(1) / 17164642023241175717b /
C      Data Log10(2) / 16367571421742254654b /, Sc/987/
C
C     Machine Constants For Ftn5 On The Cdc 6000/7000 Series.
C
C      Data Small(1) / O"00564000000000000000" /
C      Data Small(2) / O"00000000000000000000" /
C
C      Data Large(1) / O"37757777777777777777" /
C      Data Large(2) / O"37157777777777777774" /
C
C      Data Right(1) / O"15624000000000000000" /
C      Data Right(2) / O"00000000000000000000" /
C
C      Data Diver(1) / O"15634000000000000000" /
C      Data Diver(2) / O"00000000000000000000" /
C
C      Data Log10(1) / O"17164642023241175717" /
C      Data Log10(2) / O"16367571421742254654" /, Sc/987/
C
C     Machine Constants For Convex C-1
C
C      Data Small(1),Small(2) / '00100000'X, '00000000'X /
C      Data Large(1),Large(2) / '7fffffff'X, 'Ffffffff'X /
C      Data Right(1),Right(2) / '3cc00000'X, '00000000'X /
C      Data Diver(1),Diver(2) / '3cd00000'X, '00000000'X /
C      Data Log10(1),Log10(2) / '3ff34413'X, '509f79ff'X /, Sc/987/
C
C     Machine Constants For The Cray 1, Xmp, 2, And 3.
C
C      Data Small(1) / 201354000000000000000b /
C      Data Small(2) / 000000000000000000000b /
C
C      Data Large(1) / 577767777777777777777b /
C      Data Large(2) / 000007777777777777776b /
C
C      Data Right(1) / 376434000000000000000b /
C      Data Right(2) / 000000000000000000000b /
C
C      Data Diver(1) / 376444000000000000000b /
C      Data Diver(2) / 000000000000000000000b /
C
C      Data Log10(1) / 377774642023241175717b /
C      Data Log10(2) / 000007571421742254654b /, Sc/987/
C
C     Machine Constants For The Data General Eclipse S/200
C
C     Small, Large, Right, Diver, Log10 Should Be Declared
C     Integer Small(4), Large(4), Right(4), Diver(4), Log10(4)
C
C     Note - It May Be Appropriate To Include The Following Line -
C     Static Dmach(5)
C
C      Data Small/20k,3*0/,Large/77777k,3*177777k/
C      Data Right/31420k,3*0/,Diver/32020k,3*0/
C      Data Log10/40423k,42023k,50237k,74776k/, Sc/987/
C
C     Machine Constants For The Harris Slash 6 And Slash 7
C
C      Data Small(1),Small(2) / '20000000, '00000201 /
C      Data Large(1),Large(2) / '37777777, '37777577 /
C      Data Right(1),Right(2) / '20000000, '00000333 /
C      Data Diver(1),Diver(2) / '20000000, '00000334 /
C      Data Log10(1),Log10(2) / '23210115, '10237777 /, Sc/987/
C
C     Machine Constants For The Honeywell Dps 8/70 Series.
C
C      Data Small(1),Small(2) / O402400000000, O000000000000 /
C      Data Large(1),Large(2) / O376777777777, O777777777777 /
C      Data Right(1),Right(2) / O604400000000, O000000000000 /
C      Data Diver(1),Diver(2) / O606400000000, O000000000000 /
C      Data Log10(1),Log10(2) / O776464202324, O117571775714 /, Sc/987/
C
C     Machine Constants For The Ibm 360/370 Series,
C     The Xerox Sigma 5/7/9 And The Sel Systems 85/86.
C
C      Data Small(1),Small(2) / Z00100000, Z00000000 /
C      Data Large(1),Large(2) / Z7fffffff, Zffffffff /
C      Data Right(1),Right(2) / Z33100000, Z00000000 /
C      Data Diver(1),Diver(2) / Z34100000, Z00000000 /
C      Data Log10(1),Log10(2) / Z41134413, Z509f79ff /, Sc/987/
C
C     Machine Constants For The Interdata 8/32
C     With The Unix System Fortran 77 Compiler.
C
C     For The Interdata Fortran Vii Compiler Replace
C     The Z'S Specifying Hex Constants With Y'S.
C
C      Data Small(1),Small(2) / Z'00100000', Z'00000000' /
C      Data Large(1),Large(2) / Z'7effffff', Z'Ffffffff' /
C      Data Right(1),Right(2) / Z'33100000', Z'00000000' /
C      Data Diver(1),Diver(2) / Z'34100000', Z'00000000' /
C      Data Log10(1),Log10(2) / Z'41134413', Z'509f79ff' /, Sc/987/
C
C     Machine Constants For The Pdp-10 (Ka Processor).
C
C      Data Small(1),Small(2) / "033400000000, "000000000000 /
C      Data Large(1),Large(2) / "377777777777, "344777777777 /
C      Data Right(1),Right(2) / "113400000000, "000000000000 /
C      Data Diver(1),Diver(2) / "114400000000, "000000000000 /
C      Data Log10(1),Log10(2) / "177464202324, "144117571776 /, Sc/987/
C
C     Machine Constants For The Pdp-10 (Ki Processor).
C
C      Data Small(1),Small(2) / "000400000000, "000000000000 /
C      Data Large(1),Large(2) / "377777777777, "377777777777 /
C      Data Right(1),Right(2) / "103400000000, "000000000000 /
C      Data Diver(1),Diver(2) / "104400000000, "000000000000 /
C      Data Log10(1),Log10(2) / "177464202324, "047674776746 /, Sc/987/
C
C     Machine Constants For Pdp-11 Fortrans Supporting
C     32-Bit Integers (Expressed In Integer And Octal).
C
C      Data Small(1),Small(2) /    8388608,           0 /
C      Data Large(1),Large(2) / 2147483647,          -1 /
C      Data Right(1),Right(2) /  612368384,           0 /
C      Data Diver(1),Diver(2) /  620756992,           0 /
C      Data Log10(1),Log10(2) / 1067065498, -2063872008 /, Sc/987/
C
C      Data Small(1),Small(2) / O00040000000, O00000000000 /
C      Data Large(1),Large(2) / O17777777777, O37777777777 /
C      Data Right(1),Right(2) / O04440000000, O00000000000 /
C      Data Diver(1),Diver(2) / O04500000000, O00000000000 /
C      Data Log10(1),Log10(2) / O07746420232, O20476747770 /, Sc/987/
C
C     Machine Constants For Pdp-11 Fortrans Supporting
C     16-Bit Integers (Expressed In Integer And Octal).
C
C     Small, Large, Right, Diver, Log10 Should Be Declared
C     Integer Small(4), Large(4), Right(4), Diver(4), Log10(4)
C
C      Data Small(1),Small(2) /    128,      0 /
C      Data Small(3),Small(4) /      0,      0 /
C
C      Data Large(1),Large(2) /  32767,     -1 /
C      Data Large(3),Large(4) /     -1,     -1 /
C
C      Data Right(1),Right(2) /   9344,      0 /
C      Data Right(3),Right(4) /      0,      0 /
C
C      Data Diver(1),Diver(2) /   9472,      0 /
C      Data Diver(3),Diver(4) /      0,      0 /
C
C      Data Log10(1),Log10(2) /  16282,   8346 /
C      Data Log10(3),Log10(4) / -31493, -12296 /, Sc/987/
C
C      Data Small(1),Small(2) / O000200, O000000 /
C      Data Small(3),Small(4) / O000000, O000000 /
C
C      Data Large(1),Large(2) / O077777, O177777 /
C      Data Large(3),Large(4) / O177777, O177777 /
C
C      Data Right(1),Right(2) / O022200, O000000 /
C      Data Right(3),Right(4) / O000000, O000000 /
C
C      Data Diver(1),Diver(2) / O022400, O000000 /
C      Data Diver(3),Diver(4) / O000000, O000000 /
C
C      Data Log10(1),Log10(2) / O037632, O020232 /
C      Data Log10(3),Log10(4) / O102373, O147770 /, Sc/987/
C
C     Machine Constants For The Prime 50 Series Systems
C     With 32-Bit Integers And 64v Mode Instructions,
C     Supplied By Igor Bray.
C
C      Data Small(1),Small(2) / :10000000000, :00000100001 /
C      Data Large(1),Large(2) / :17777777777, :37777677775 /
C      Data Right(1),Right(2) / :10000000000, :00000000122 /
C      Data Diver(1),Diver(2) / :10000000000, :00000000123 /
C      Data Log10(1),Log10(2) / :11504046501, :07674600177 /, Sc/987/
C
C     Machine Constants For The Sequent Balance 8000
C
C      Data Small(1),Small(2) / $00000000,  $00100000 /
C      Data Large(1),Large(2) / $ffffffff,  $7fefffff /
C      Data Right(1),Right(2) / $00000000,  $3ca00000 /
C      Data Diver(1),Diver(2) / $00000000,  $3cb00000 /
C      Data Log10(1),Log10(2) / $509f79ff,  $3fd34413 /, Sc/987/
C
C     Machine Constants For The Univac 1100 Series.
C
C      Data Small(1),Small(2) / O000040000000, O000000000000 /
C      Data Large(1),Large(2) / O377777777777, O777777777777 /
C      Data Right(1),Right(2) / O170540000000, O000000000000 /
C      Data Diver(1),Diver(2) / O170640000000, O000000000000 /
C      Data Log10(1),Log10(2) / O177746420232, O411757177572 /, Sc/987/
C
C     Machine Constants For The Vax Unix F77 Compiler
C
C      Data Small(1),Small(2) /        128,           0 /
C      Data Large(1),Large(2) /     -32769,          -1 /
C      Data Right(1),Right(2) /       9344,           0 /
C      Data Diver(1),Diver(2) /       9472,           0 /
C      Data Log10(1),Log10(2) /  546979738,  -805796613 /, Sc/987/
C
C     Machine Constants For The Vax-11 With
C     Fortran Iv-Plus Compiler
C
C      Data Small(1),Small(2) / Z00000080, Z00000000 /
C      Data Large(1),Large(2) / Zffff7fff, Zffffffff /
C      Data Right(1),Right(2) / Z00002480, Z00000000 /
C      Data Diver(1),Diver(2) / Z00002500, Z00000000 /
C      Data Log10(1),Log10(2) / Z209a3f9a, Zcff884fb /, Sc/987/
C
C     Machine Constants For Vax/Vms Version 2.2
C
C      Data Small(1),Small(2) /       '80'X,        '0'X /
C      Data Large(1),Large(2) / 'Ffff7fff'X, 'Ffffffff'X /
C      Data Right(1),Right(2) /     '2480'X,        '0'X /
C      Data Diver(1),Diver(2) /     '2500'X,        '0'X /
C      Data Log10(1),Log10(2) / '209a3f9a'X, 'Cff884fb'X /, Sc/987/
C
C  ***  Issue Stop 779 If All Data Statements Are Commented...
      If (Sc .ne. 987) Stop 779
C  ***  Issue Stop 778 If All Data Statements Are Obviously Wrong...
      If (Dmach(4) .ge. 1.0d0) Stop 778
      If (I .lt. 1  .or.  I .gt. 5) Goto 999
      D1mach = Dmach(I)
      Return
  999 Write(*,1999) I
 1999 Format(' D1mach - I Out Of Bounds',I10)
      Stop
      End





C
C
C-----------------------------------------------------------------------
C
C     THE AUGUST 27 1992 VERSION OF COLROW IN WHICH X IS NO LONGER
C     REQUIRED, WITH THE SOLUTION BEING RETURNED IN B, THE RIGHT
C     HAND SIDE.  IN ADDITION, ALL VARIABLES ARE EXPLICITLY DECLARED.
C     A PARAMETER "JOB" IS INCLUDED, TO SPECIFY WHICH OF A.X = B OR
C     TRANSPOSE(A).X = B IS TO BE SOLVED.

      SUBROUTINE COLROW(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,IFLAG,JOB)
C
C***************************************************************
C
C  THIS PROGRAM SOLVES ONE OF THE LINEAR SYSTEMS  A*X = B OR
C  TRANSPOSE(A)*X = B, WHERE  A IS AN ALMOST BLOCK DIAGONAL
C  MATRIX OF THE FORM
C
C               TOPBLK
C               ARRAY(1)
C                     ARRAY(2)
C                          .
C                             .
C                                .
C                                   .
C                                    ARRAY(NBLOKS)
C                                           BOTBLK
C
C  WHERE
C           TOPBLK IS  NRWTOP  BY NOVRLP
C           ARRAY(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP
C           BOTBLK IS NRWBOT BY NOVRLP,
C  AND
C           NOVRLP = NRWTOP + NRWBOT
C  WITH
C           NOVRLP.LE.NRWBLK .
C
C  THE LINEAR SYSTEM IS OF ORDER  N = NBLOKS*NRWBLK + NOVRLP.
C
C  THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH
C  ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING,
C  WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX  A
C  WITHOUT INTRODUCING FILL-IN.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE, REMOVE
C  ALL DOUBLE PRECISION STATEMENTS.  THERE IS ONE SUCH STATEMENT
C  IN C O L R O W, THREE IN C R D C M P, AND TWO IN C R S O L V.
C  IN ADDITION, REFERENCES TO BUILT-IN FUNCTIONS DABS AND DMAX1
C  MUST BE REPLACED BY ABS AND AMAX1, RESPECTIVELY.  DABS OCCURS
C  NINE TIMES, IN C R D C M P.  DMAX1 OCCURS FOUR TIMES, IN
C  C R D C M P.  FINALLY, ZERO IS INITIALISED TO 0.D0 IN A
C  DATA STATEMENT IN C R D C M P.  THIS MUST BE REPLACED BY:
C               DATA ZERO/0.0/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C               JOB    - INTEGER, INDICATING:
C                      = 0: SOLVE A*X = B;
C                      NON-ZERO: SOLVE TRANSPOSE(A)*X = B.
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                    B - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR (IF IFLAG = 0)
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  AUXILIARY PROGRAMS  *****
C
C       CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,IFLAG)
C            - DECOMPOSES THE MATRIX  A  USING MODIFIED
C              ALTERNATE ROW AND COLUMN ELIMINATON WITH
C              PARTIAL PIVOTING, AND IS USED FOR THIS
C              PURPOSE IN C O L R O W.
C              THE ARGUMENTS ARE AS IN C O L R O W.
C
C       CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,B,JOB)
C            - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED.
C              THE ARGUMENTS ARE ALL AS IN C O L R O W.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       THE SUBROUTINE  C O L R O W  AUTOMATICALLY SOLVES THE
C  INPUT SYSTEM WHEN IFLAG=0.  C O L R O W  IS CALLED ONLY ONCE
C  FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C  HAND SIDES CAN BE OBTAINED BY ONE CALL TO  C O L R O W  AND
C  P-1 CALLS TO CRSLVE ONLY. SINCE THE ARRAYS TOPBLK,ARRAY,
C  BOTBLK AND PIVOT CONTAIN THE DECOMPOSITION OF THE GIVEN
C  COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C  C O L R O W , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C  CALLS TO CRSLVE WITH THE SAME LEFT HAND SIDES. FOR THE
C  SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C  MATRIX, THE ARRAYS TOPBLK,ARRAY,BOTBLK MUST BE COPIED
C  BEFORE A CALL TO  C O L R O W .
C
C*************************************************************************
        IMPLICIT NONE
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,B
        INTEGER N,NRWTOP,NOVRLP,NRWBLK,NCLBLK,NBLOKS,NRWBOT,PIVOT(*),
     *          IFLAG,JOB, IDAMAX,i,IFAIL,j,jj
        DIMENSION TOPBLK(NRWTOP,*),ARRAY(NRWBLK,NCLBLK,*),
     *          BOTBLK(NRWBOT,*),B(*)

c        NOUT = 6

C       DO THE FACTORIZATION USING CRDCMP:
C
        CALL CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
     *          BOTBLK,NRWBOT,PIVOT,IFLAG)
        IF(IFLAG.NE.0) RETURN
C
c     *****************solving the linear system********************
        job=0

        CALL CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
     *          BOTBLK,NRWBOT,PIVOT,B,JOB)
        RETURN
        END

C******************************************************************



      SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,
     *   NBLOKS,BOTBLK,NRWBOT,PIVOT,IFLAG)
C
C***************************************************************
C
C  C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A
C  USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH
C  PARTIAL PIVOTING.  THE MATRIX  A  IS STORED IN THE ARRAYS
C  TOPBLK, ARRAY, AND BOTBLK.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A TO BE DECOMPOSED
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER PIVOT(*)
      DIMENSION TOPBLK(NRWTOP,*),ARRAY(NRWBLK,NCLBLK,*),BOTBLK(NRWBOT,*)
      DATA ZERO / 0.0D+0 /
      Common/Mchprs/ Flmin, Flmax, Epsmch
C
C***************************************************************
C
C          ****  DEFINE THE CONSDTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
      IFLAG = 0
      PIVMAX = ZERO
      NRWTP1 = NRWTOP+1
      NROWEL = NRWBLK-NRWTOP
      NRWEL1 = NROWEL+1
      NVRLP0 = NOVRLP-1
      Pivtol = 10.0d+0*Epsmch
C
C***************************************************************
C
C          ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C               IF PARAMETERS ARE INVALID THEN TERMINATE AT 10;
C                                         ELSE CONTINUE AT 100.
C
C***************************************************************
C
      IF (N.NE.NBLOKS*NRWBLK+NOVRLP) GO TO 10
      IF (NOVRLP.NE.NRWTOP+NRWBOT) GO TO 10
      IF (NCLBLK.NE.NOVRLP+NRWBLK) GO TO 10
      IF (NOVRLP.GT.NRWBLK) GO TO 10
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      GO TO 20
 10    CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = 1
      RETURN
 20    CONTINUE
C
C***************************************************************
C
C               ****  FIRST, IN TOPBLK....
C
C***************************************************************
C
C          ***  APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN
C                 PIVOTING ....
C
C***************************************************************
C
      DO 110 I = 1, NRWTOP
         IPLUS1 = I+1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         IPVT = I
         COLMAX = DABS(TOPBLK(I,I))
         DO 30 J = IPLUS1, NOVRLP
            TEMPIV = DABS(TOPBLK(I,J))
            IF (TEMPIV.LE.COLMAX) GO TO 30
            IPVT = J
            COLMAX = TEMPIV
 30             CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C

         If (Colmax .le. Pivtol) Then
            Iflag = -1
            Return
         Endif

c         IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410
         PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         PIVOT(I) = IPVT
         IF (IPVT.EQ.I) GO TO 60
         DO 40 L = I, NRWTOP
            SWAP = TOPBLK(L,IPVT)
            TOPBLK(L,IPVT) = TOPBLK(L,I)
            TOPBLK(L,I) = SWAP
 40             CONTINUE
         DO 50 L = 1, NRWBLK
            SWAP = ARRAY(L,IPVT,1)
            ARRAY(L,IPVT,1) = ARRAY(L,I,1)
            ARRAY(L,I,1) = SWAP
 50             CONTINUE
 60                 CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         COLPIV = TOPBLK(I,I)
         DO 100 J = IPLUS1, NOVRLP
            COLMLT = TOPBLK(I,J)/COLPIV
            TOPBLK(I,J) = COLMLT
            IF (IPLUS1.GT.NRWTOP) GO TO 80
            DO 70 L = IPLUS1, NRWTOP
               TOPBLK(L,J) = TOPBLK(L,J)-COLMLT*TOPBLK(L,I)
 70                   CONTINUE
 80                          CONTINUE
            DO 90 L = 1, NRWBLK
               ARRAY(L,J,1) = ARRAY(L,J,1)-COLMLT*ARRAY(L,I,1)
 90                   CONTINUE
 100                      CONTINUE
 110                       CONTINUE
C
C***************************************************************
C
C          ****  IN EACH BLOCK ARRAY(,,K)....
C
C***************************************************************
C
      INCR = 0
      DO 320 K = 1, NBLOKS
         KPLUS1 = K+1
C
C          *****************************************************
C
C          ***  FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH
C                       ROW PIVOTING....
C
C          *****************************************************
C
         DO 180 J = NRWTP1, NRWBLK
            JPLUS1 = J+1
            JMINN = J-NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            IPVT = JMINN
            ROWMAX = DABS(ARRAY(JMINN,J,K))
            LOOP = JMINN+1
            DO 120 I = LOOP, NRWBLK
               TEMPIV = DABS(ARRAY(I,J,K))
               IF (TEMPIV.LE.ROWMAX) GO TO 120
               IPVT = I
               ROWMAX = TEMPIV
 120                  CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         If (Rowmax .le. Pivtol) Then
            Iflag = -1
            Return
         Endif

c            IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410
            PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            INCRJ = INCR+J
            PIVOT(INCRJ) = INCR+IPVT+NRWTOP
            IF (IPVT.EQ.JMINN) GO TO 140
            DO 130 L = J, NCLBLK
               SWAP = ARRAY(IPVT,L,K)
               ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K)
               ARRAY(JMINN,L,K) = SWAP
 130                  CONTINUE
 140                         CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            ROWPIV = ARRAY(JMINN,J,K)
            DO 150 I = LOOP, NRWBLK
               ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV
 150                  CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            DO 170 L = JPLUS1, NCLBLK
               ROWMLT = ARRAY(JMINN,L,K)
               DO 160 I = LOOP, NRWBLK
                  ARRAY(I,L,K) = ARRAY(I,L,K)-ROWMLT*ARRAY(I,J,K)
 160                        CONTINUE
 170                               CONTINUE
 180                                   CONTINUE
C
C          *****************************************************
C
C          ***  NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH
C                      COLUMN PIVOTING....
C
C          *****************************************************
C
         DO 310 I = NRWEL1, NRWBLK
            IPLUSN = I+NRWTOP
            IPLUS1 = I+1
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            IPVT = IPLUSN
            COLMAX = DABS(ARRAY(I,IPVT,K))
            LOOP = IPLUSN+1
            DO 190 J = LOOP, NCLBLK
               TEMPIV = DABS(ARRAY(I,J,K))
               IF (TEMPIV.LE.COLMAX) GO TO 190
               IPVT = J
               COLMAX = TEMPIV
 190                  CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C

         If (Colmax .le. Pivtol) Then
            Iflag = -1
            Return
         Endif

c            IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410
            PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            INCRN = INCR+IPLUSN
            PIVOT(INCRN) = INCR+IPVT
            IRWBLK = IPLUSN-NRWBLK
            IF (IPVT.EQ.IPLUSN) GO TO 240
            DO 200 L = I, NRWBLK
               SWAP = ARRAY(L,IPVT,K)
               ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K)
               ARRAY(L,IPLUSN,K) = SWAP
 200                  CONTINUE
            IPVBLK = IPVT-NRWBLK
            IF (K.EQ.NBLOKS) GO TO 220
            DO 210 L = 1, NRWBLK
               SWAP = ARRAY(L,IPVBLK,KPLUS1)
               ARRAY(L,IPVBLK,KPLUS1) = ARRAY(L,IRWBLK,KPLUS1)
               ARRAY(L,IRWBLK,KPLUS1) = SWAP
 210                  CONTINUE
            GO TO 240
 220               CONTINUE
            DO 230 L = 1, NRWBOT
               SWAP = BOTBLK(L,IPVBLK)
               BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK)
               BOTBLK(L,IRWBLK) = SWAP
 230                  CONTINUE
 240                         CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
            COLPIV = ARRAY(I,IPLUSN,K)
            DO 300 J = LOOP, NCLBLK
               COLMLT = ARRAY(I,J,K)/COLPIV
               ARRAY(I,J,K) = COLMLT
               IF (I.EQ.NRWBLK) GO TO 260
               DO 250 L = IPLUS1, NRWBLK
                  ARRAY(L,J,K) = ARRAY(L,J,K)-COLMLT*ARRAY(L,IPLUSN,K)
 250                        CONTINUE
 260                                  CONTINUE
               JRWBLK = J-NRWBLK
               IF (K.EQ.NBLOKS) GO TO 280
               DO 270 L = 1, NRWBLK
                  ARRAY(L,JRWBLK,KPLUS1) = ARRAY(L,JRWBLK,KPLUS1)-COLMLT
     *               *ARRAY(L,IRWBLK,KPLUS1)
 270                        CONTINUE
               GO TO 300
 280                     CONTINUE
               DO 290 L = 1, NRWBOT
                  BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK)-COLMLT*BOTBLK(L,
     *               IRWBLK)
 290                        CONTINUE
 300                               CONTINUE
 310                                   CONTINUE
         INCR = INCR+NRWBLK
 320      CONTINUE
C
C***************************************************************
C
C          ****  FINALLY, IN BOTBLK....
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          ***  APPLY NRWBOT ROW ELIMINATIONS WITH ROW
C                  PIVOTING....
C
C               IF BOT HAS JUST ONE ROW GO TO 500
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C

      IF (NRWBOT.EQ.1) GO TO 400
      DO 390 J = NRWTP1, NVRLP0
         JPLUS1 = J+1
         JMINN = J-NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         IPVT = JMINN
         ROWMAX = DABS(BOTBLK(JMINN,J))
         LOOP = JMINN+1
         DO 330 I = LOOP, NRWBOT
            TEMPIV = DABS(BOTBLK(I,J))
            IF (TEMPIV.LE.ROWMAX) GO TO 330
            IPVT = I
            ROWMAX = TEMPIV
 330            CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         If (Rowmax .le. Pivtol) Then
            Iflag = -1
            Return
         Endif

c         IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410
         PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         INCRJ = INCR+J
         PIVOT(INCRJ) = INCR+IPVT+NRWTOP
         IF (IPVT.EQ.JMINN) GO TO 350
         DO 340 L = J, NOVRLP
            SWAP = BOTBLK(IPVT,L)
            BOTBLK(IPVT,L) = BOTBLK(JMINN,L)
            BOTBLK(JMINN,L) = SWAP
 340            CONTINUE
 350                CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         ROWPIV = BOTBLK(JMINN,J)
         DO 360 I = LOOP, NRWBOT
            BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV
 360            CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         DO 380 L = JPLUS1, NOVRLP
            ROWMLT = BOTBLK(JMINN,L)
            DO 370 I = LOOP, NRWBOT
               BOTBLK(I,L) = BOTBLK(I,L)-ROWMLT*BOTBLK(I,J)
 370                  CONTINUE
 380                      CONTINUE
 390                       CONTINUE
 400                        CONTINUE
C
C***************************************************************
C
C          DONE PROVIDED THE LAST ELEMENT IS NOT ZERO
C
C***************************************************************
C
       if (DABS(BOTBLK(NRWBOT,NOVRLP)).gt.pivtol) return
c      IF (PIVMAX+DABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN

C
C***************************************************************
C
C       ****  MATRIX IS SINGULAR - SET IFLAG = - 1.
C                                  TERMINATE AT 1000.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
 410   CONTINUE
      IFLAG = -1
      RETURN
      END

        SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,JOB)
C
C***************************************************************
C
C  C R S L V E  SOLVES THE LINEAR SYSTEM
C                       A*X = B
C  USING THE DECOMPOSITION ALREADY GENERATED IN  C R D C M P.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY  ...
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         THE PIVOT VECTOR FROM  C R D C M P
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C               JOB    - INTEGER, INDICATING:
C                      = 0: SOLVE A*X = B;
C                      NON-ZERO: SOLVE TRANSPOSE(A)*X = B.
C
C       *** ON RETURN  ...
C
C                    B - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
        IMPLICIT NONE
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,B
        DOUBLE PRECISION DOTPRD,BJ,XINCRJ,BINCRJ,SWAP,BI
        INTEGER NRWTOP,NOVRLP,NRWBLK,NCLBLK,NBLOKS,NRWBOT,PIVOT(*),
     *          JOB
        DIMENSION TOPBLK(NRWTOP,*),ARRAY(NRWBLK,NCLBLK,*),
     *          BOTBLK(NRWBOT,*),B(*)
        INTEGER NRWTP1,NRWBK1,NVRLP1,NRWBT1,NROWEL,NVRLP0,NBLKS1,
     *          NBKTOP,J,I,LOOP,INCR,INCRJ,INCRI,JPIVOT,JRWTOP,
     *          LL,L1,IPLUSN,INCRN,NRWTP0,NRWEL1,K,INCRTP,NRWBTL,
     *          IPVTN,NRWELL,IPVTI,L
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
        NRWTP1 = NRWTOP+1
        NRWBK1 = NRWBLK+1
        NVRLP1 = NOVRLP+1
        NRWTP0 = NRWTOP-1
        NRWBT1 = NRWBOT+1
        NROWEL = NRWBLK-NRWTOP
        NRWEL1 = NROWEL+1
        NVRLP0 = NOVRLP-1
        NBLKS1 = NBLOKS+1
        NBKTOP = NRWBLK+NRWTOP
C
C       IF JOB IS NON-ZERO, TRANSFER TO THE SECTION DEALING WITH
C       TRANSPOSE(A)*X = B.
C
        IF ( JOB .NE. 0 ) GO TO 530
C
C***************************************************************
C
C               ****  FORWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST, IN TOPBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 130 J = 1,NRWTOP
           B(J) = B(J)/TOPBLK(J,J)
           IF(J.EQ.NRWTOP)GO TO 120
              BJ = -B(J)
              LOOP = J+1
              DO 110 I = LOOP,NRWTOP
                 B(I) = B(I)+TOPBLK(I,J)*BJ
 110          CONTINUE
 120       CONTINUE
 130    CONTINUE
C
C       ********************************************************
C
C          ***  IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCR = 0
        DO 280 K = 1,NBLOKS
           INCRTP = INCR+NRWTOP
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 220 J = 1,NRWTOP
              INCRJ = INCR+J
              XINCRJ = -B(INCRJ)
              DO 210 I = 1,NRWBLK
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
 210          CONTINUE
 220       CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 240 J = NRWTP1,NRWBLK
              INCRJ = INCR+J
              JPIVOT = PIVOT(INCRJ)
              IF(JPIVOT.EQ.INCRJ)GO TO 225
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(JPIVOT)
                 B(JPIVOT) = SWAP
 225          CONTINUE
              BINCRJ = -B(INCRJ)
              LOOP = J-NRWTP0
              DO 230 I = LOOP,NRWBLK
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+ARRAY(I,J,K)*BINCRJ
 230          CONTINUE
 240       CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 270 J = NRWBK1,NBKTOP
              INCRJ = INCR+J
              JRWTOP = J -NRWTOP
              B(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K)
              IF(J.EQ.NBKTOP)GO TO 260
                 XINCRJ = -B(INCRJ)
                 LOOP = J-NRWTP0
                 DO 250 I = LOOP,NRWBLK
                    INCRI = INCRTP+I
                    B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
 250             CONTINUE
 260          CONTINUE
 270       CONTINUE
           INCR = INCR+NRWBLK
 280       CONTINUE
C
C       ********************************************************
C
C          ***  FINALLY, IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRTP = INCR+NRWTOP
        DO 320 J = 1,NRWTOP
           INCRJ = INCR+J
           XINCRJ = -B(INCRJ)
           DO 310 I = 1,NRWBOT
              INCRI = INCRTP+I
              B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ
 310       CONTINUE
 320     CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF(NRWBOT.EQ.1)GO TO 350
           DO 340 J = NRWTP1,NVRLP0
              INCRJ = INCR+J
              JPIVOT = PIVOT(INCRJ)
              IF(JPIVOT.EQ.INCRJ)GO TO 325
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(JPIVOT)
                 B(JPIVOT) = SWAP
 325          CONTINUE
              BINCRJ = -B(INCRJ)
              LOOP = J-NRWTP0
              DO 330 I = LOOP,NRWBOT
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+BOTBLK(I,J)*BINCRJ
 330          CONTINUE
 340          CONTINUE
 350       CONTINUE
C
C***************************************************************
C
C               ****  BACKWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 430 LL = 1,NRWBOT
           J = NVRLP1-LL
           INCRJ = INCR+J
           NRWBTL = NRWBT1-LL
           B(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J)
           IF(LL.EQ.NRWBOT)GO TO 420
              XINCRJ = -B(INCRJ)
              LOOP = NRWBOT-LL
              DO 410 I = 1,LOOP
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ
 410          CONTINUE
 420       CONTINUE
 430    CONTINUE
C
C       ********************************************************
C
C          ***  THEN IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 490 L = 1,NBLOKS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           K = NBLKS1-L
           INCR = INCR-NRWBLK
           DO 450 L1 = NRWEL1,NRWBLK
              I = NRWBLK+NRWEL1-L1
              IPLUSN = I+NRWTOP
              LOOP = IPLUSN+1
              INCRN = INCR+IPLUSN
              DOTPRD = B(INCRN)
              DO 440 J = LOOP,NCLBLK
                 INCRJ = INCR+J
                 DOTPRD = DOTPRD-ARRAY(I,J,K)*B(INCRJ)
 440          CONTINUE
              B(INCRN) = DOTPRD
              IPVTN = PIVOT(INCRN)
              IF(INCRN.EQ.IPVTN)GO TO 445
                 SWAP = B(INCRN)
                 B(INCRN) = B(IPVTN)
                 B(IPVTN) = SWAP
 445          CONTINUE
 450      CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           INCRTP = INCR+NRWTOP
           DO 460 J = NRWBK1,NCLBLK
              INCRJ = INCR+J
              XINCRJ = -B(INCRJ)
              DO 455 I = 1,NROWEL
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
 455          CONTINUE
 460       CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 480 LL = 1,NROWEL
              J = NRWBK1-LL
              INCRJ = INCR+J
              NRWELL = NRWEL1-LL
              B(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K)
              IF(LL.EQ.NROWEL)GO TO 470
                 XINCRJ = -B(INCRJ)
                 LOOP = NROWEL-LL
                 DO 465 I = 1,LOOP
                    INCRI = INCRTP+I
                    B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
 465             CONTINUE
 470          CONTINUE
 480       CONTINUE
 490      CONTINUE
C
C       ********************************************************
C
C          ***  IN TOPBLK FINISH WITH....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 520 L = 1,NRWTOP
           I = NRWTP1-L
           LOOP = I+1
           DOTPRD = B(I)
           DO 510 J = LOOP,NOVRLP
              DOTPRD = DOTPRD-TOPBLK(I,J)*B(J)
 510       CONTINUE
           B(I) = DOTPRD
           IPVTI = PIVOT(I)
           IF(I.EQ.IPVTI)GO TO 515
                 SWAP = B(I)
                 B(I) = B(IPVTI)
                 B(IPVTI) = SWAP
 515       CONTINUE
 520    CONTINUE
C
C       RETURN FROM THE SOLUTION OF A.X = B.
        RETURN
C
C       IF JOB IS NON-ZERO, SOLVE TRANSPOSE(A)*X = B:
C
 530           CONTINUE

C       FIRST, FORWARD ELIMINATION OF RHS USING TRANSPOSE(U).

        DO 540 I = 1,NRWTOP
           IPVTI = PIVOT(I)
           IF ( I .NE. IPVTI ) THEN
              SWAP = B(I)
              B(I) = B(IPVTI)
              B(IPVTI) = SWAP
           ENDIF
           BI = -B(I)
           LOOP = I+1
           DO 535 J = LOOP,NOVRLP
              B(J) = B(J) + BI*TOPBLK(I,J)
 535       CONTINUE
 540    CONTINUE

C       IN EACH BLOCK, K = 1,..,NBLOKS:

        INCR = NRWTOP
        DO 590 K = 1,NBLOKS

C          FIRST, THE FORWARD SOLUTION.

           DO 550 J = 1,NROWEL
              INCRJ = INCR + J
              DO 545 I = 1,J-1
                 B(INCRJ) = B(INCRJ) - ARRAY(I,NRWTOP+J,K)*B(INCR+I)
 545          CONTINUE
              B(INCRJ) = B(INCRJ)/ARRAY(J,NRWTOP+J,K)
 550       CONTINUE

C           FORWARD MODIFICATION.

            DO 570 I = 1,NOVRLP
               INCRI = INCR + NROWEL + I
               LOOP = NRWBLK + I
               DO 560 J = 1,NROWEL
                  INCRJ = INCR + J
                  B(INCRI) = B(INCRI) - ARRAY(J,LOOP,K)*B(INCRJ)
 560           CONTINUE
 570        CONTINUE

C           NOW, FORWARD ELIMINATION OF RHS USING TRANSPOSE(U). THIS
C           CORRESPONDS TO THE LOOP 540 ABOVE.

            INCR = INCR + NROWEL
            DO 580 I = 1,NRWTOP
               INCRI = INCR + I
               IPVTI = PIVOT(INCRI)
               IF ( INCRI .NE. IPVTI ) THEN
                  SWAP = B(INCRI)
                  B(INCRI) = B(IPVTI)
                  B(IPVTI) = SWAP
               ENDIF
               LOOP = NROWEL + I
               BI = -B(INCRI)
               DO 575 J = I+1,NOVRLP
                  INCRJ = INCR+J
                  L = NRWBLK + J


                  B(INCRJ) = B(INCRJ) + BI*ARRAY(LOOP,L,K)
 575           CONTINUE
 580        CONTINUE
            INCR = INCR + NRWTOP
 590      CONTINUE

C       FINALLY, FINISH WITH NRWBOT SOLUTIONS:

        DO 600 J = 1,NRWBOT
           INCRJ = INCR + J
           DO 595 I = 1,J-1
              B(INCRJ) = B(INCRJ) - BOTBLK(I,J+NRWTOP)*B(INCR+I)
 595       CONTINUE
           B(INCRJ) = B(INCRJ)/BOTBLK(J,J+NRWTOP)
 600    CONTINUE


C       NOW, THE BACKWARD PASS:


C       FIRST, BACKWARD SOLUTION IN BOTBLK:

        INCRJ = INCR + NRWBOT
        DO 610 J = 1,NRWBOT-1
           INCRJ = INCRJ - 1
           DO 605 I = NRWBOT-J+1,NRWBOT
              INCRI = INCR + I
              B(INCRJ) = B(INCRJ) - BOTBLK(I,NOVRLP-J)*B(INCRI)
 605       CONTINUE

           IF ( INCRJ .NE. PIVOT(INCRJ) ) THEN
              SWAP = B(INCRJ)
              B(INCRJ) = B(PIVOT(INCRJ))
              B(PIVOT(INCRJ)) = SWAP
           ENDIF
 610    CONTINUE

C       NOW DO THE DEFERRED OPERATIONS IN BOTBLOK:

        DO 620 J = 1,NRWTOP
           INCRJ = INCR - J + 1
           DO 615 I = 1,NRWBOT
              INCRI = INCR + I
              B(INCRJ) = B(INCRJ) - BOTBLK(I,NRWTP1-J)*B(INCRI)
 615       CONTINUE
 620    CONTINUE


C       NOW, IN EACH BLOCK, K = NBLOKS,..,1:
        DO 800 K = NBLOKS,1,-1

C          FIRST, THE BACKSUBSTITUIONS:

           DO 630 J = 1,NRWTOP
              INCRJ = INCR - J + 1
              LOOP = NBKTOP - J + 1
              DO 625 I = 1,J-1
                 INCRI = INCR - I + 1
                 B(INCRJ) = B(INCRJ) - ARRAY(NRWBLK-I+1,LOOP,K)*B(INCRI)
 625          CONTINUE
              B(INCRJ) = B(INCRJ)/ARRAY(NRWBLK-J+1,LOOP,K)
 630       CONTINUE

C          THEN THE BACKWARD SOLUTION IN THE KTH BLOCK:

           DO 650 J = 1,NROWEL
              INCRJ = INCR - NRWTOP -J + 1
              DO 645 I = 1,J+NRWTOP-1
                 INCRI = INCRJ + I
                 B(INCRJ) = B(INCRJ) -
     *           ARRAY(NRWBLK-NRWTOP-J+1+I,NRWBLK-J+1,K)*B(INCRI)
 645          CONTINUE
              IF ( INCRJ .NE. PIVOT(INCRJ) ) THEN
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(PIVOT(INCRJ))
                 B(PIVOT(INCRJ)) = SWAP
              ENDIF
 650       CONTINUE

C          NOW, THE DEFERRED OPERATIONS ON B:

           INCR = INCR - NRWBLK
           DO 660 J = 1,NRWTOP
              INCRJ = INCR + J - NRWTOP
              DO 655 I = 1,NRWBLK
                 INCRI = INCR + I
                 B(INCRJ) = B(INCRJ) - ARRAY(I,J,K)*B(INCRI)
 655          CONTINUE
 660        CONTINUE
 800    CONTINUE

C       FINALLY, THE LAST SET OF BACK-SUBSTITUTIONS IN TOPBLK:

        DO 900 J = 1,NRWTOP
           INCRJ = NRWTOP -J + 1
           DO 850 I = INCRJ+1,NRWTOP
              B(INCRJ) = B(INCRJ) - TOPBLK(I,INCRJ)*B(I)
 850       CONTINUE
           B(INCRJ) = B(INCRJ)/TOPBLK(INCRJ,INCRJ)
 900    CONTINUE
C
C       RETURN FROM THE SOLUTION OF A-TRANSPOSE.X = B

        RETURN
        END


      double precision function dasum(n,dx,incx)
c#
c#     takes the sum of the absolute values.
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#
      double precision dx(*),dtemp
      integer i,incx,m,mp1,n,nincx
c#
      dasum = 0.0d0
      dtemp = 0.0d0
      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
        dtemp = dtemp + dabs(dx(i))
 10      continue
      dasum = dtemp
      return
c#
c#        code for increment equal to 1
c#
c#
c#        clean-up loop
c#
 20    m = mod(n,6)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dabs(dx(i))
 30      continue
      if( n .lt. 6 ) go to 60
 40    mp1 = m + 1
      do 50 i = mp1,n,6
        dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
     *  + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
 50      continue
 60       dasum = dtemp
      return
      end




      Subroutine Lufac(N, Ndim, A, Ip, Ier)
      Implicit Double Precision (A-H,O-Z)
      Dimension A(Ndim,N), Ip(N)
      Intrinsic Abs
*  Blas: Daxpy, Dscal, Dswap. Idamax

      Parameter ( Zero = 0.0d+0, One = 1.0d+0 )

*  The Subroutine Lufac Is A Very Simple Code To Compute The
*  Lu Decomposition (With Partial Pivoting) Of The N By N Matrix A.

*  The Lu Factors Are Overwritten On A.  The Integer Array Ip
*  Reflects The Pairwise Interchanges Performed.  Note That Ip(K)
*  Therefore Does Not Give The Index In The Original Array Of
*  The K-Th Pivot.

*  On Exit, The Error Flag Ier Is Zero When No Zero Pivots Are
*  Encountered.  Otherwise, Ier Is Equal To The Index Of The
*  Step At Which A Zero Pivot Occurred.

      Ier = 0
      Ip(N) = 0

*  Begin Loop Over Columns 1 Through N-1.  K Is The Current
*  Column Index.

      Do 100 K = 1, N-1

*  Find The Row Index Ipiv Of The Element Of Largest Magnitude In
*  Column K.

         Ipiv = K-1 + Idamax(N-K+1, A(K,K), 1)
         Piv = A(Ipiv,K)
         If (Piv .eq. Zero) Then
            Ier = K
            Return
         Endif
         Ip(K) = Ipiv

*  Perform Interchanges If Necessary.

         If (Ipiv .ne. K) Then
            Call Dswap(N-K+1, A(Ipiv,K), Ndim, A(K,K), Ndim)
         Endif

*  Save The (Negative) Multipliers In The Subdiagonal Elements Of
*  Column K.

         Call Dscal(N-K, (-One/Piv), A(K+1,K), 1)

*  Update The Remaining Matrix.  Note That A(I,K) Now Contains
*  The Negative Multipliers.

         Do 50 J = K+1, N
            Call Daxpy(N-K, A(K,J), A(K+1,K), 1, A(K+1,J), 1)
   50    Continue

*  End Of Loop Over Columns.

  100 Continue
      If (A(N,N) .eq. Zero) Ier = N
      Return
      End



      Subroutine Lusol (N, Ndim, A, Ip, B, X)
      Implicit Double Precision (A-H, O-Z)
      Dimension A(Ndim,N), Ip(N), B(N), X(N)

*  Blas:  Daxpy, Dcopy

*  The Subroutine Lusol Is A Simple-Minded Routine To Solve A
*  Linear System Whose Lu Factors Have Been Computed By Lufac.
*  On Entry, The Matrix A Should Contain The Lu Factors, And
*  Ip Should Contain The Interchange Array Constructed By Lufac.


*  Copy The Right-Hand Side B Into X.

      Call Dcopy (N, B, 1, X, 1)

*  Forward Solution With L (Unit Lower-Triangular Factor), Which
*  Is Stored In The Strict Lower Triangle Of A.

      Do 20 K = 1, N-1
         Ipiv = Ip(K)
         If (Ipiv .ne. K) Then
            Tem = X(Ipiv)
            X(Ipiv) = X(K)
            X(K) = Tem
         Endif
         Call Daxpy ( N-K, X(K), A(K+1,K), 1, X(K+1), 1 )
   20 Continue

*  Backward Solution With U (Upper-Triangular Factor), Which Is Stored
*  In The Upper Triangle Of A.

      Do 40 Kb = N, 1, -1
         X(Kb) = X(Kb)/A(Kb,Kb)
         Call Daxpy(Kb-1, (-X(Kb)), A(1,Kb), 1, X(1), 1)
   40 Continue

      Return
      End

      Subroutine Dcopy ( N, X, Incx, Y, Incy )
      Implicit Double Precision (A-H,O-Z)
      Integer            N, Incx, Incy
      Dimension          X( * ), Y( * )

C  Dcopy  Performs The Operation
C
C     Y : =  X
C
C  Nag Fortran 77 Version Of The Blas Routine Dcopy .
C  Nag Fortran 77 O( N ) Basic Linear Algebra Routine.
C
C  -- Written On 26-November-1982.
C     Sven Hammarling, Nag Central Office.

      Integer            I     , Ix    , Iy

      If( N .lt. 1 )Return
      If( ( Incx .eq. Incy ) .and. ( Incy .gt. 0 ) )Then
         Do 10, Iy = 1, 1 + ( N - 1 )*Incy, Incy
            Y( Iy ) = X( Iy )
   10    Continue
      Else
         If( Incx .ge. 0 )Then
            Ix = 1
         Else
            Ix = 1 - ( N - 1 )*Incx
         End If
         If( Incy .gt. 0 )Then
            Do 20, Iy = 1, 1 + ( N - 1 )*Incy, Incy
               Y( Iy ) = X( Ix )
               Ix      = Ix + Incx
   20       Continue
         Else
            Iy = 1 - ( N - 1 )*Incy
            Do 30, I = 1, N
               Y( Iy ) = X( Ix )
               Iy      = Iy + Incy
               Ix      = Ix + Incx
   30       Continue
         End If
      End If

      Return

*     End Of Dcopy .

      End

      Subroutine Daxpy ( N, Alpha, X, Incx, Y, Incy )
      Implicit Double Precision (A-H,O-Z)
      Integer            N, Incx, Incy
      Dimension   X( * ), Y( * )

C  Daxpy  Performs The Operation
C
C     Y : =  Alpha*X + Y
C
C
C  Modified Nag Fortran 77 Version Of The Blas Routine Daxpy .
C
C  -- Written On 3-September-1982.
C     Sven Hammarling, Nag Central Office.

      Integer            I     , Ix    , Iy
      Parameter        ( Zero  = 0.0d+0 )

      If( N     .lt. 1    )Return
      If( Alpha .eq. Zero )Return

      If( ( Incx .eq. Incy ) .and. ( Incx .gt. 0 ) )Then
         If ( Alpha .eq. 1.0d+0 ) Then
            Do 5, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               Y( Ix ) = X( Ix ) + Y( Ix )
 5          Continue
         Else
            Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               Y( Ix ) = Alpha*X( Ix ) + Y( Ix )
 10         Continue
         Endif
      Else
         If( Incy .ge. 0 )Then
            Iy = 1
         Else
            Iy = 1 - ( N - 1 )*Incy
         End If
         If( Incx .gt. 0 )Then
            Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               Y( Iy ) = Alpha*X( Ix ) + Y( Iy )
               Iy      = Iy + Incy
   20       Continue
         Else
            Ix = 1 - ( N - 1 )*Incx
            Do 30, I = 1, N
               Y( Iy ) = Alpha*X( Ix ) + Y( Iy )
               Ix      = Ix + Incx
               Iy      = Iy + Incy
   30       Continue
         End If
      End If
      Return

*     End Of Daxpy .

      End

      Double Precision Function Ddot  ( N, X, Incx, Y, Incy )
      Implicit Double Precision (A-H,O-Z)
      Integer           N, Incx, Incy
      Dimension         X( * ), Y( * )

C  Ddot   Returns The Value
C
C     Ddot   = X'Y
C
C
C  Modified Nag Fortran 77 Version Of The Blas Routine Ddot  .
C
C  -- Written On 21-September-1982.
C     Sven Hammarling, Nag Central Office.

      Integer             I     , Ix    , Iy
      Parameter         ( Zero  = 0.0d+0 )

      Sum = Zero
      If( N .ge. 1 )Then
         If( ( Incx .eq. Incy ) .and. ( Incx .gt. 0 ) )Then
            Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               Sum = Sum + X( Ix )*Y( Ix )
   10       Continue
         Else
            If( Incy .ge. 0 )Then
               Iy = 1
            Else
               Iy = 1 - ( N - 1 )*Incy
            End If
            If( Incx .gt. 0 )Then
               Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx
                  Sum = Sum + X( Ix )*Y( Iy )
                  Iy  = Iy  + Incy
   20          Continue
            Else
               Ix = 1 - ( N - 1 )*Incx
               Do 30, I = 1, N
                  Sum = Sum + X( Ix )*Y( Iy )
                  Ix  = Ix  + Incx
                  Iy  = Iy  + Incy
   30          Continue
            End If
         End If
      End If

      Ddot   = Sum
      Return

*     End Of Ddot  .

      End

      Subroutine Dscal ( N, Alpha, X, Incx )
      Implicit Double Precision (A-H,O-Z)
      Integer          N, Incx
      Dimension        X( * )

C  Dscal  Performs The Operation
C
C     X : =  Alpha*X
C
C
C  Modified Nag Fortran 77 Version Of The Blas Routine Dscal .
C
C  -- Written On 26-November-1982.
C     Sven Hammarling, Nag Central Office.

      Integer            Ix
      Parameter        ( One   = 1.0d+0, Zero  = 0.0d+0 )

      If( N .ge. 1 )Then
         If( Alpha .eq. Zero )Then
            Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               X( Ix ) = Zero
   10       Continue
         Else If( Alpha .eq. ( -One ) )Then
            Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               X( Ix ) = -X( Ix )
   20       Continue
         Else If( Alpha .ne. One )Then
            Do 30, Ix = 1, 1 + ( N - 1 )*Incx, Incx
               X( Ix ) = Alpha*X( Ix )
   30       Continue
         End If
      End If

      Return

*     End Of Dscal .

      End

      Subroutine Dswap ( N, X, Incx, Y, Incy )
      Implicit Double Precision (A-H,O-Z)
      Integer      N, Incx, Incy
      Dimension    X( * ), Y( * )

C  Dswap  Performs The Operations
C
C     Temp : =  X,   X : =  Y,   Y : =  Temp.
C
C
C  Modified Nag Fortran 77 Version Of The Blas Routine Dswap .
C
C  -- Written On 26-November-1982.
C     Sven Hammarling, Nag Central Office.

      Integer            I     , Ix    , Iy

      If( N .lt. 1 )Return

      If( ( Incx .eq. Incy ) .and. ( Incy .gt. 0 ) )Then
         Do 10, Iy = 1, 1 + ( N - 1 )*Incy, Incy
            Temp    = X( Iy )
            X( Iy ) = Y( Iy )
            Y( Iy ) = Temp
   10    Continue
      Else
         If( Incx .ge. 0 )Then
            Ix = 1
         Else
            Ix = 1 - ( N - 1 )*Incx
         End If
         If( Incy .gt. 0 )Then
            Do 20, Iy = 1, 1 + ( N - 1 )*Incy, Incy
               Temp    = X( Ix )
               X( Ix ) = Y( Iy )
               Y( Iy ) = Temp
               Ix      = Ix + Incx
   20       Continue
         Else
            Iy = 1 - ( N - 1 )*Incy
            Do 30, I = 1, N
               Temp    = X( Ix )
               X( Ix ) = Y( Iy )
               Y( Iy ) = Temp
               Iy      = Iy + Incy
               Ix      = Ix + Incx
   30       Continue
         End If
      End If

      Return

*     End Of Dswap .

      End

      Integer Function Idamax( N, X, Incx )
      Implicit Double Precision (A-H,O-Z)
      Integer         N, Incx
      Dimension       X( * )

C  Idamax Returns The Smallest Value Of I Such That
C
C     Abs( X( I ) ) = Max( Abs( X( J ) ) )
C                      J
C
C  Nag Fortran 77 Version Of The Blas Routine Idamax.
C  Nag Fortran 77 O( N ) Basic Linear Algebra Routine.
C
C  -- Written On 31-May-1983.
C     Sven Hammarling, Nag Central Office.

      Intrinsic           Abs
      Integer             I     , Imax  , Ix

      If( N .lt. 1 )Then
         Idamax = 0
         Return
      End If

      Imax = 1
      If( N .gt. 1 )Then
         Xmax = Abs( X( 1 ) )
         Ix   = 1
         Do 10, I = 2, N
            Ix = Ix + Incx
            If( Xmax .lt. Abs( X( Ix ) ) )Then
               Xmax = Abs( X( Ix ) )
               Imax = I
            End If
   10    Continue
      End If

      Idamax = Imax
      Return

*     End Of Idamax.

      End
      Subroutine Dload ( N, Const, X, Incx )
      Implicit Double Precision (A-H,O-Z)
      Dimension  X(*)
C
C  Dload  Performs The Operation
C
C     X = Const*E,   E' = ( 1  1 ... 1 ).
C
C
C  Nag Fortran 77 O( N ) Basic Linear Algebra Routine.
C
C  -- Written On 22-September-1983.
C     Sven Hammarling, Nag Central Office.
C
      Parameter        ( Zero = 0.0d+0 )

      If( N .lt. 1 )Return

      If( Const .ne. Zero )Then
         Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx
            X( Ix ) = Const
   10    Continue
      Else
         Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx
            X( Ix ) = Zero
   20    Continue
      End If

      Return

*     End Of Dload .
      End



      Subroutine Maxpy ( Nrow, Ncol, Alpha, Xmat, Nrowy, Ymat )
      Implicit Double Precision (A-H,O-Z)
      Dimension Xmat(Nrow, Ncol), Ymat(Nrowy, Ncol)

*  Subroutine Maxpy Takes As Input The Scalar Alpha And Two Matrices,
*  Xmat And Ymat.  Xmat Has Declared Row Dimension Nrow, And
*  Ymat Has Declared Row Dimension Nrowy, But Both Are
*  Conceptually Nrow By Ncol.
*  On Output, (New Ymat) Is Alpha*Xmat+ (Old Ymat), By Analogy
*  With The Vector Blas Routine Saxpy.

      Do 100 J = 1, Ncol
      Do 100 I = 1, Nrow
         Ymat(I,J) = Ymat(I,J) + Alpha*Xmat(I,J)
  100 Continue
      Return
      End


      Subroutine Matcop( Nrow1, Nrow2, Nrow, Ncol, Xmat1, Xmat2 )
      Implicit Double Precision (A-H,O-Z)
      Dimension Xmat1(Nrow1, Ncol), Xmat2(Nrow2, Ncol)

*  Given 2 Matrices Xmat1 And Xmat2, Where Xmat1 Has Declared
*  Row Dimension Nrow1, Xmat2 Has Declared Row Dimension Nrow2,
*  And Both Have Column Dimension Ncol, The Routine Matcop Copies
*  Rows 1 Through Nrow, And Columns 1 Through Ncol From Xmat1 Into
*  Xmat2.


      If (Nrow .le. 0 .or. Ncol .le. 0) Return
      Do 100 J = 1, Ncol
      Do 100 I = 1, Nrow
           Xmat2(I,J) = Xmat1(I,J)
  100 Continue
      Return
      End

      Subroutine Mtload( Nrow, Ncol, Const, Nrowx, Xmat )
      Implicit Double Precision (A-H,O-Z)
      Dimension Xmat(Nrowx, Ncol)

*  Mtload Sets Elements 1 Through Nrow, 1 Through Ncol, Of The
*  Matrix Xmat (Whose Declared Row Dimension Is Nrowx) To The
*  Scalar Value Const.

      If (Nrow .le. 0 .or. Ncol .le. 0)  Return
      Do 100 J = 1, Ncol
      Do 100 I = 1, Nrow
         Xmat(I,J) = Const
  100 Continue
      Return
      End

      Subroutine Mssq  ( Nrow, Ncol, Xmat, Scale, Sumsq )
      Implicit Double Precision (A-H,O-Z)
      Dimension   Xmat(Nrow, *)

*  Given The Nrow By Ncol Matrix Xmat, Mssq Returns Values
*  Scale And Sumsq Such That
*     (Scale**2) * Sumsq = Sum Of Squares Of Xmat(I,J),
*  Where  Scale = Max  Abs(Xmat(I,J)).

*  Mssq Is A Stripped-Down Matrix Version Of The Blas Routine Sssq.

      Intrinsic    Abs
      Parameter   ( One   = 1.0d+0, Zero  = 0.0d+0 )

      Scale = Zero
      Sumsq = One
      If( Nrow .ge. 1 .and. Ncol .ge. 1) Then
         Do 10 I = 1, Nrow
         Do 10 J = 1, Ncol
            If( Xmat(I,J) .ne. Zero )Then
               Absxij = Abs(Xmat(I,J))
               If( Scale .lt. Absxij ) Then
                  Sumsq = One + Sumsq* (Scale/Absxij)**2
                  Scale = Absxij
               Else
                  Sumsq = Sumsq + (Absxij/Scale)**2
               End If
            End If
   10    Continue
      End If
      Return
      End

      Subroutine Dssq  ( N, X, Incx, Scale, Sumsq )
      Implicit Double Precision (A-H,O-Z)
      Integer            N, Incx
      Dimension   X( * )

*  Given The N-Vector X, Dssq Returns Values Scale And Sumsq Such That
*     (Scale**2) * Sumsq = Sum Of Squares Of X(I),
*  Where  Scale = Max  Abs(X(I)).

      Intrinsic          Abs
      Parameter        ( One   = 1.0d+0, Zero  = 0.0d+0 )

      Scale = Zero
      Sumsq = One
      If( N .ge. 1 )Then
         Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx
            If( X( Ix ) .ne. Zero )Then
               Absxi = Abs( X( Ix ) )
               If( Scale .lt. Absxi )Then
                  Sumsq = One   + Sumsq*( Scale/Absxi )**2
                  Scale = Absxi
               Else
                  Sumsq = Sumsq + ( Absxi/Scale )**2
               End If
            End If
   10    Continue
      End If

      Return

*     End Of Dssq.

      End






