C***********************************************************************
C    Module:  xsolve.f
C 
C    Copyright (C) 2000 Mark Drela 
C 
C    This program is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    This program is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with this program; if not, write to the Free Software
C    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C***********************************************************************


      SUBROUTINE GAUSS(NSIZ,NN,Z,R,NRHS)
C     *******************************************************
C     *                                                     *
C     *   Solves general NxN system in NN unknowns          *
C     *    with arbitrary number (NRHS) of righthand sides. *
C     *   Assumes system is invertible...                   *
C     *    ...if it isn't, a divide by zero will result.    *
C     *                                                     *
C     *   Z is the coefficient matrix...                    *
C     *     ...destroyed during solution process.           *
C     *   R is the righthand side(s)...                     *
C     *     ...replaced by the solution vector(s).          *
C     *                                                     *
C     *                              Mark Drela  1984       *
C     *******************************************************
C
      DIMENSION Z(NSIZ,NSIZ), R(NSIZ,NRHS)
C
      DO 1 NP=1, NN-1
        NP1 = NP+1
C
C------ find max pivot index NX
        NX = NP
        DO 11 N=NP1, NN
          IF(ABS(Z(N,NP))-ABS(Z(NX,NP))) 11,11,111
  111      NX = N
   11   CONTINUE
C
        PIVOT = 1.0/Z(NX,NP)
C
C------ switch pivots
        Z(NX,NP) = Z(NP,NP)
C
C------ switch rows & normalize pivot row
        DO 12 L=NP1, NN
          TEMP = Z(NX,L)*PIVOT
          Z(NX,L) = Z(NP,L)
          Z(NP,L) = TEMP
   12   CONTINUE
C
        DO 13 L=1, NRHS
          TEMP = R(NX,L)*PIVOT
          R(NX,L) = R(NP,L)
          R(NP,L) = TEMP
   13   CONTINUE
C
C------ forward eliminate everything
        DO 15 K=NP1, NN
          ZTMP = Z(K,NP)
C
C          IF(ZTMP.EQ.0.0) GO TO 15
C
          DO 151 L=NP1, NN
            Z(K,L) = Z(K,L) - ZTMP*Z(NP,L)
  151     CONTINUE
          DO 152 L=1, NRHS
            R(K,L) = R(K,L) - ZTMP*R(NP,L)
  152     CONTINUE
   15   CONTINUE
C
    1 CONTINUE
C
C---- solve for last row
      DO 2 L=1, NRHS
        R(NN,L) = R(NN,L)/Z(NN,NN)
    2 CONTINUE
C
C---- back substitute everything
      DO 3 NP=NN-1, 1, -1
        NP1 = NP+1
        DO 31 L=1, NRHS
          DO 310 K=NP1, NN
            R(NP,L) = R(NP,L) - Z(NP,K)*R(K,L)
  310     CONTINUE
   31   CONTINUE
    3 CONTINUE
C
      RETURN
      END ! GAUSS


      SUBROUTINE CGAUSS(NSIZ,NN,Z,R,NRHS)
C********************************************
C     Solves general complex linear systems.
C********************************************
      COMPLEX Z(NSIZ,NSIZ), R(NSIZ,NRHS)
      COMPLEX PIVOT, TEMP, ZTMP
C
      DO 1 NP=1, NN-1
        NP1 = NP+1
C
C------ find max pivot index NX
        NX = NP
        DO 11 N=NP1, NN
          IF(ABS(Z(N,NP))-ABS(Z(NX,NP))) 11,11,111
  111      NX = N
   11   CONTINUE
C
        PIVOT = (1.0,0.0)/Z(NX,NP)
C
C------ switch pivots
        Z(NX,NP) = Z(NP,NP)
C
C------ switch rows & normalize pivot row
        DO 12 L=NP1, NN
          TEMP = Z(NX,L)*PIVOT
          Z(NX,L) = Z(NP,L)
          Z(NP,L) = TEMP
   12   CONTINUE
C
        DO 13 L=1, NRHS
          TEMP = R(NX,L)*PIVOT
          R(NX,L) = R(NP,L)
          R(NP,L) = TEMP
   13   CONTINUE
C
C------ forward eliminate everything
        DO 15 K=NP1, NN
          ZTMP = Z(K,NP)
C
C          IF(ZTMP.EQ.0.0) GO TO 15
C
          DO 151 L=NP1, NN
            Z(K,L) = Z(K,L) - ZTMP*Z(NP,L)
  151     CONTINUE
          DO 152 L=1, NRHS
            R(K,L) = R(K,L) - ZTMP*R(NP,L)
  152     CONTINUE
   15   CONTINUE
C
    1 CONTINUE
C
C---- solve for last row
      DO 2 L=1, NRHS
        R(NN,L) = R(NN,L)/Z(NN,NN)
    2 CONTINUE
C
C---- back substitute everything
      DO 3 NP=NN-1, 1, -1
        NP1 = NP+1
        DO 31 L=1, NRHS
          DO 310 K=NP1, NN
            R(NP,L) = R(NP,L) - Z(NP,K)*R(K,L)
  310     CONTINUE
   31   CONTINUE
    3 CONTINUE
C
      RETURN
      END ! CGAUSS



      SUBROUTINE LUDCMP(NSIZ,N,A,INDX)
C     *******************************************************
C     *                                                     *
C     *   Factors a full NxN matrix into an LU form.        *
C     *   Subr. BAKSUB can back-substitute it with some RHS.*
C     *   Assumes matrix is non-singular...                 *
C     *    ...if it isn't, a divide by zero will result.    *
C     *                                                     *
C     *   A is the matrix...                                *
C     *     ...replaced with its LU factors.                *
C     *                                                     *
C     *                              Mark Drela  1988       *
C     *******************************************************
C
      DIMENSION A(NSIZ,NSIZ), INDX(NSIZ)
C
      PARAMETER (NVX=300)
      DIMENSION VV(NVX)
C
      IF(N.GT.NVX) STOP 'LUDCMP: Array overflow. Increase NVX.'
C
      DO 12 I=1, N
        AAMAX = 0.
        DO 11 J=1, N
          AAMAX = MAX( ABS(A(I,J)) , AAMAX )
   11   CONTINUE
        VV(I) = 1.0/AAMAX
   12 CONTINUE
C
      DO 19 J=1, N
        DO 14 I=1, J-1
          SUM = A(I,J)
          DO 13 K=1, I-1
            SUM = SUM - A(I,K)*A(K,J)
   13     CONTINUE
          A(I,J) = SUM
   14   CONTINUE
C
        AAMAX = 0.
        DO 16 I=J, N
          SUM = A(I,J)
          DO 15 K=1, J-1
            SUM = SUM - A(I,K)*A(K,J)
   15     CONTINUE
          A(I,J) = SUM
C
          DUM = VV(I)*ABS(SUM)
          IF(DUM.GE.AAMAX) THEN
           IMAX = I
           AAMAX = DUM
          ENDIF
   16   CONTINUE
C
        IF(J.NE.IMAX) THEN
         DO 17 K=1, N
           DUM = A(IMAX,K)
           A(IMAX,K) = A(J,K)
           A(J,K) = DUM
   17    CONTINUE
         VV(IMAX) = VV(J)
        ENDIF
C
        INDX(J) = IMAX
        IF(J.NE.N) THEN
         DUM = 1.0/A(J,J)
         DO 18 I=J+1, N
           A(I,J) = A(I,J)*DUM
   18    CONTINUE
        ENDIF
C
   19 CONTINUE
C
      RETURN
      END ! LUDCMP


      SUBROUTINE BAKSUB(NSIZ,N,A,INDX,B)
      DIMENSION A(NSIZ,NSIZ), B(NSIZ), INDX(NSIZ)
C
      II = 0
      DO 12 I=1, N
        LL = INDX(I)
        SUM = B(LL)
        B(LL) = B(I)
        IF(II.NE.0) THEN
         DO 11 J=II, I-1
           SUM = SUM - A(I,J)*B(J)
   11    CONTINUE
        ELSE IF(SUM.NE.0.0) THEN
         II = I
        ENDIF
        B(I) = SUM
   12 CONTINUE
C
      DO 14 I=N, 1, -1
        SUM = B(I)
        IF(I.LT.N) THEN
         DO 13 J=I+1, N
           SUM = SUM - A(I,J)*B(J)
   13    CONTINUE
        ENDIF
        B(I) = SUM/A(I,I)
   14 CONTINUE
C
      RETURN
      END ! BAKSUB



      SUBROUTINE BLSOLV
C-----------------------------------------------------------------
C      Custom solver for coupled viscous-inviscid Newton system:
C
C        A  |  |  .  |  |  .  |    d       R       S
C        B  A  |  .  |  |  .  |    d       R       S
C        |  B  A  .  |  |  .  |    d       R       S
C        .  .  .  .  |  |  .  |    d   =   R - dRe S
C        |  |  |  B  A  |  .  |    d       R       S
C        |  Z  |  |  B  A  .  |    d       R       S
C        .  .  .  .  .  .  .  |    d       R       S
C        |  |  |  |  |  |  B  A    d       R       S
C
C       A, B, Z  3x3  blocks containing linearized BL equation coefficients
C       |        3x1  vectors containing mass defect influence 
C                     coefficients on Ue
C       d        3x1  unknown vectors (Newton deltas for Ctau, Theta, m)
C       R        3x1  residual vectors
C       S        3x1  Re influence vectors
C-----------------------------------------------------------------
      INCLUDE 'XFOIL.INC'
C
      IVTE1 = ISYS(IBLTE(1),1)
C
      VACC1 = VACCEL
      VACC2 = VACCEL * 2.0 / (S(N) - S(1))
      VACC3 = VACCEL * 2.0 / (S(N) - S(1))
C
      DO 1000 IV=1, NSYS
C
        IVP = IV + 1
C
C====== Invert VA(IV) block
C
C------ normalize first row
        PIVOT = 1.0 / VA(1,1,IV)
        VA(1,2,IV) = VA(1,2,IV) * PIVOT
        DO 10 L=IV, NSYS
          VM(1,L,IV) = VM(1,L,IV)*PIVOT
   10   CONTINUE
        VDEL(1,1,IV) = VDEL(1,1,IV)*PIVOT
        VDEL(1,2,IV) = VDEL(1,2,IV)*PIVOT
C
C------ eliminate lower first column in VA block
        DO 15 K=2, 3
          VTMP = VA(K,1,IV)
          VA(K,2,IV) = VA(K,2,IV) - VTMP*VA(1,2,IV)
          DO 150 L=IV, NSYS
            VM(K,L,IV) = VM(K,L,IV) - VTMP*VM(1,L,IV)
  150     CONTINUE
          VDEL(K,1,IV) = VDEL(K,1,IV) - VTMP*VDEL(1,1,IV)
          VDEL(K,2,IV) = VDEL(K,2,IV) - VTMP*VDEL(1,2,IV)
   15   CONTINUE
C
C
C------ normalize second row
        PIVOT = 1.0 / VA(2,2,IV)
        DO 20 L=IV, NSYS
          VM(2,L,IV) = VM(2,L,IV)*PIVOT
   20   CONTINUE
        VDEL(2,1,IV) = VDEL(2,1,IV)*PIVOT
        VDEL(2,2,IV) = VDEL(2,2,IV)*PIVOT
C
C------ eliminate lower second column in VA block
        K = 3
        VTMP = VA(K,2,IV)
        DO 250 L=IV, NSYS
          VM(K,L,IV) = VM(K,L,IV) - VTMP*VM(2,L,IV)
  250   CONTINUE
        VDEL(K,1,IV) = VDEL(K,1,IV) - VTMP*VDEL(2,1,IV)
        VDEL(K,2,IV) = VDEL(K,2,IV) - VTMP*VDEL(2,2,IV)
C
C
C------ normalize third row
        PIVOT = 1.0/VM(3,IV,IV)
        DO 350 L=IVP, NSYS
          VM(3,L,IV) = VM(3,L,IV)*PIVOT
  350   CONTINUE
        VDEL(3,1,IV) = VDEL(3,1,IV)*PIVOT
        VDEL(3,2,IV) = VDEL(3,2,IV)*PIVOT
C
C
C------ eliminate upper third column in VA block
        VTMP1 = VM(1,IV,IV)
        VTMP2 = VM(2,IV,IV)
        DO 450 L=IVP, NSYS
          VM(1,L,IV) = VM(1,L,IV) - VTMP1*VM(3,L,IV)
          VM(2,L,IV) = VM(2,L,IV) - VTMP2*VM(3,L,IV)
  450   CONTINUE
        VDEL(1,1,IV) = VDEL(1,1,IV) - VTMP1*VDEL(3,1,IV)
        VDEL(2,1,IV) = VDEL(2,1,IV) - VTMP2*VDEL(3,1,IV)
        VDEL(1,2,IV) = VDEL(1,2,IV) - VTMP1*VDEL(3,2,IV)
        VDEL(2,2,IV) = VDEL(2,2,IV) - VTMP2*VDEL(3,2,IV)
C
C------ eliminate upper second column in VA block
        VTMP = VA(1,2,IV)
        DO 460 L=IVP, NSYS
          VM(1,L,IV) = VM(1,L,IV) - VTMP*VM(2,L,IV)
  460   CONTINUE
        VDEL(1,1,IV) = VDEL(1,1,IV) - VTMP*VDEL(2,1,IV)
        VDEL(1,2,IV) = VDEL(1,2,IV) - VTMP*VDEL(2,2,IV)
C
C
        IF(IV.EQ.NSYS) GO TO 1000
C
C====== Eliminate VB(IV+1) block, rows  1 -> 3
        DO 50 K=1, 3
          VTMP1 = VB(K, 1,IVP)
          VTMP2 = VB(K, 2,IVP)
          VTMP3 = VM(K,IV,IVP)
          DO 510 L=IVP, NSYS
            VM(K,L,IVP) = VM(K,L,IVP)
     &        - (  VTMP1*VM(1,L,IV)
     &           + VTMP2*VM(2,L,IV)
     &           + VTMP3*VM(3,L,IV) )
  510     CONTINUE
          VDEL(K,1,IVP) = VDEL(K,1,IVP)
     &        - (  VTMP1*VDEL(1,1,IV)
     &           + VTMP2*VDEL(2,1,IV)
     &           + VTMP3*VDEL(3,1,IV) )
          VDEL(K,2,IVP) = VDEL(K,2,IVP)
     &        - (  VTMP1*VDEL(1,2,IV)
     &           + VTMP2*VDEL(2,2,IV)
     &           + VTMP3*VDEL(3,2,IV) )
   50   CONTINUE
C
        IF(IV.EQ.IVTE1) THEN
C------- eliminate VZ block
         IVZ = ISYS(IBLTE(2)+1,2)
C
         DO 55 K=1, 3
           VTMP1 = VZ(K,1)
           VTMP2 = VZ(K,2)
           DO 515 L=IVP, NSYS
             VM(K,L,IVZ) = VM(K,L,IVZ)
     &         - (  VTMP1*VM(1,L,IV)
     &            + VTMP2*VM(2,L,IV) )
  515      CONTINUE
           VDEL(K,1,IVZ) = VDEL(K,1,IVZ)
     &         - (  VTMP1*VDEL(1,1,IV)
     &            + VTMP2*VDEL(2,1,IV) )
           VDEL(K,2,IVZ) = VDEL(K,2,IVZ)
     &         - (  VTMP1*VDEL(1,2,IV)
     &            + VTMP2*VDEL(2,2,IV) )
   55    CONTINUE
        ENDIF
C
        IF(IVP.EQ.NSYS) GO TO 1000
C
C====== Eliminate lower VM column
        DO 60 KV=IV+2, NSYS
          VTMP1 = VM(1,IV,KV)
          VTMP2 = VM(2,IV,KV)
          VTMP3 = VM(3,IV,KV)
C
          IF(ABS(VTMP1).GT.VACC1) THEN
          DO 610 L=IVP, NSYS
            VM(1,L,KV) = VM(1,L,KV) - VTMP1*VM(3,L,IV)
  610     CONTINUE
          VDEL(1,1,KV) = VDEL(1,1,KV) - VTMP1*VDEL(3,1,IV)
          VDEL(1,2,KV) = VDEL(1,2,KV) - VTMP1*VDEL(3,2,IV)
          ENDIF
C
          IF(ABS(VTMP2).GT.VACC2) THEN
          DO 620 L=IVP, NSYS
            VM(2,L,KV) = VM(2,L,KV) - VTMP2*VM(3,L,IV)
  620     CONTINUE
          VDEL(2,1,KV) = VDEL(2,1,KV) - VTMP2*VDEL(3,1,IV)
          VDEL(2,2,KV) = VDEL(2,2,KV) - VTMP2*VDEL(3,2,IV)
          ENDIF
C
          IF(ABS(VTMP3).GT.VACC3) THEN
          DO 630 L=IVP, NSYS
            VM(3,L,KV) = VM(3,L,KV) - VTMP3*VM(3,L,IV)
  630     CONTINUE
          VDEL(3,1,KV) = VDEL(3,1,KV) - VTMP3*VDEL(3,1,IV)
          VDEL(3,2,KV) = VDEL(3,2,KV) - VTMP3*VDEL(3,2,IV)
          ENDIF
C
   60   CONTINUE
C
 1000 CONTINUE
C
C
C
      DO 2000 IV=NSYS, 2, -1
C
C------ eliminate upper VM columns
        VTMP = VDEL(3,1,IV)
        DO 81 KV=IV-1, 1, -1
          VDEL(1,1,KV) = VDEL(1,1,KV) - VM(1,IV,KV)*VTMP
          VDEL(2,1,KV) = VDEL(2,1,KV) - VM(2,IV,KV)*VTMP
          VDEL(3,1,KV) = VDEL(3,1,KV) - VM(3,IV,KV)*VTMP
   81   CONTINUE
C
        VTMP = VDEL(3,2,IV)
        DO 82 KV=IV-1, 1, -1
          VDEL(1,2,KV) = VDEL(1,2,KV) - VM(1,IV,KV)*VTMP
          VDEL(2,2,KV) = VDEL(2,2,KV) - VM(2,IV,KV)*VTMP
          VDEL(3,2,KV) = VDEL(3,2,KV) - VM(3,IV,KV)*VTMP
   82   CONTINUE
C
 2000 CONTINUE
C
      RETURN
      END
