C                 **  VARIABLE DEFINITION LIST  **
C
C     A(3,3)  - STORES THE A-MATRIX FOR A LAMINATE (CLT ANALYSIS).          
C     AI(3,3) - STORES THE INVERSE OF THE A-MATRIX FOR A LAMINATE.
C     CRSH(I) - STORES THE CRITICAL SHEAR LOAD FOR LAMINATE (I).
C     CRSHMS(I) - STORES THE M. OF S. FOR THE CRITICAL SHEAR LOAD
C               - FOR LAMINATE (I).
C     D(3,3)  - STORES THE D-MATRIX FOR A LAMINATE (CLT ANALYSIS).
C     DI(3,3) - STORES THE INVERSE OF THE D-MATRIX FOR A LAMINATE.
C     DIRCR(I) - STORES WHICH STRAIN DIRECTION IS CRITICAL.
C              - (IE: 1-DIR,2-DIR,1/2-DIR) FOR PLIES IN LAMINATE (I).
C     EP1(I,J) - STORES EPSILON-1 FOR PLY (J) IN LAMINATE (I).
C     EP2(I,J) - STORES EPSILON-2 FOR PLY (J) IN LAMINATE (I).
C     EP1MAX(J) - MAX ALLOWABLE STRAIN IN 1-DIRECTION FOR MATERIAL (J).
C     EP2MAX(J) - MAX ALLOWABLE STRAIN IN 2-DIRECTION FOR MATERIAL (J).
C     E1(J),E2(J) - YOUNG'S MODULI FOR EACH MATERIAL (J).
C     FIT(I) - STORES FITNESS FOR LAMINATE (I) IN PARENT POPULATION.
C     GAM12(I,J) - STORES GAMMA-12 FOR PLY (J) IN LAMINATE (I).
C     G12MAX(J) - MAX ALOWABLE SHEAR STRAIN IN 1-2 DIRECTION.
C               - FOR MATERIAL (J)
C     G12(J) - SHEAR MODULUS FOR MATERIAL (J).
C     IS(I,J) - AN ARRAY WHICH STORES THE STACKING SEQUENCE FOR ALL
C       - LAMINATES IN THE POPULATION. 
C       - J = 1 => NPLY;     I = 1 => NPOP.  
C       - PLIES ARE STORED IN SINGLE PLY STACKS AND HALF THE
C       - LAMINATE IS STORED BECAUSE ALL LAMINATES ARE ASSUMED 
C       - SYMMETRIC ABOUT THE MID-PLANE.
C       - 0 = EMPTY PLY; 1 = 0 DEG. PLY; 2 = +/-45 DEG. PLY;
C       - 3 = 90 DEG. PLY.
C     ISD(J) - STORES ENTIRE DECODED STACKING SEQ. FOR A LAMINATE
C        - PLIES ARE STORED IN STACKS OF 1 AND ALL ORIENTATIONS
C        - EXCLUDING EMPTY STACKS ARE STORED.
C        - Ex: IS(*,*) = (0,1,2,3) =>
C        -     ISD(*) = (0,45,90,90,45,0).
C     LAMBDA - A PARMATER USED TO MAKE SURE THE GA DOES NOT CONVERGE
C           - TO VERY THIN LAMINATES.
C     LAMTH(I) - TOTAL LAMINATE THICKNESS FOR EACH LAMINATE (I).
C     M, N -  NUMBER OF TERMS IN SINE SERIES FOR BUCKLING ANALYSIS
C     MATL(I,J) - STORES MATERIAL TYPE FOR EACH PLY (J) IN 
C               - LAMINATE (I).  
C     MATLD(J) - STORES MATERIAL TYPE FOR EACH PLY (J) IN A LAMINATE.
C     NMAT - THE NUMBER OF DIFFERENT MATERIAL TYPES.
C     NPLY - THE NO. OF SINGLE PLIES IN HALF THE LAMINATE 
C           - (LAMINATE IS ASSUMED TO BE SYMMETRIC).
C     NPOP - THE NO. OF LAMINATES IN THE POPULATION.
C     NUM(I) - STORES THE NUMBER OF NON-EMPTY SINGLE PLIES
C             - CONTAINED IN LAMINATE (I). 
C     NU12(J),NU21 - POISSON'S RATIOS FOR EACH MATERIAL (J).
C     NX,NY,NXY - IN-PLANE LOADS APLIED TO THE LAMINATE.
C     N45S(J) - STORES THE NUMBER OF +/- 45 DEGREE PLIES IN A
C             - LAMINATE FOR EACH MATERIAL (J).
C     PLYCR(I) - STORES PLY IN LAMINATE (I) 
C              - THAT HAS LARGEST RATIO: (ALL. STRAIN/PLY STRAIN).
C     PLYMAX - A LARGE INTEGER NUMBER USED FOR DIMENSIONING ARRAYS.
C     POPMAX - THE MAXIMUM NUMBER OF LAMINATES THAT CAN EXIST IN A
C               POPULATION (USED TO DIMENSION ARRAYS).
C     QB(3,3) - STORES THE QBAR MATRIX FOR A LAMINATE (CLT ANALYSIS).
C     STRMS(I) - MIN AMONG PLIES: (ALLOW. PLY STRAIN/PLY STRAIN).
C     THICK(J) - THICKNESS OF EACH INDIVIDUAL PLY MADE OF MATERIAL (J).
C     TOTMAX -  A LARGE INTEGER VALUE USED FOR DIMENSIONING ARRAYS.
C     U(5,J) - STORES THE STIFFNESS INVARIANTS FOR EACH MATERIAL (J).
C     UNBAL - PENALTY MULITPLIER FOR LAMINATES THAT ARE UNBALANCED.
C     VIOL - STORES DESIGN FEASIBILITY.
C     X - LENGTH OF PLATE.
C     Y - WIDTH OF PLATE.
***********************************************************************
C
      SUBROUTINE ANAL(SUBPOP,I,IS,MATL,NPLY,
     $                FIT,X,Y,WRITE)
      IMPLICIT NONE
      INTEGER POPMAX,PLYMAX,TOTMAX
      PARAMETER (POPMAX = 500,PLYMAX = 200,TOTMAX = 400)

      INTEGER SUBPOP
      LOGICAL WRITE

      INTEGER IS(POPMAX,PLYMAX),ISD(TOTMAX),
     $ NMAT,NPLY,VIOL(POPMAX),N15S(10),N30S(10),
     $ N45S(10),N60S(10),N75S(10),MATL(POPMAX,PLYMAX),
     $ MATLD(TOTMAX),M,N,MN,NUM(TOTMAX),DIRCR(POPMAX),
     $ BLCHECK(POPMAX),I,PLYCR(POPMAX),J

      INTEGER COUNT

      DOUBLE PRECISION CRSH(POPMAX),CRSHMS(POPMAX),FIT(POPMAX),
     $ EP1(POPMAX,TOTMAX),EP2(POPMAX,TOTMAX),GAM12(POPMAX,TOTMAX),
     $ U(5,10),STRMS(POPMAX),LAMBDA,UNBAL,E1(10),E2(10),G12(10),
     $ NU12(10),THICK(10),EP1MAX(10),EP2MAX(10),G12MAX(10),NX,NY,
     $ NXY,ALPHA,RHO(10),COST(10),MATCOST(POPMAX),LAMTH(POPMAX),
     $ M1COST(POPMAX),M2COST(POPMAX),LAYCOST(POPMAX),WTMAT1(POPMAX),
     $ WTMAT2(POPMAX),WTTOT(POPMAX),A(3,3),D(3,3),AI(3,3),DI(3,3),
     $ X(POPMAX),Y(POPMAX)


C set maximum number of single plies that can be found in 
C the coded stacking sequence of a laminate
C     this only needs to be called once ever
      COUNT = COUNT + 1
      IF (COUNT .EQ. 1) THEN

      OPEN(UNIT = 1, STATUS='OLD', FILE = 'two_mat_analysis.in')        
      READ(1,'(/)')
      READ (1,*) NMAT 
      READ(1,'(//)')

      DO J = 1,NMAT
        READ(1,*) E1(J)
        READ(1,*) E2(J) 
        READ(1,*) NU12(J) 
        READ(1,*) G12(J) 
        READ(1,*) THICK(J)
        READ(1,*) EP1MAX(J)
        READ(1,*) EP2MAX(J)
        READ(1,*) G12MAX(J)
        READ(1,*) RHO(J)
        READ(1,*) COST(J)
        READ(1,'(//)')

      END DO

      READ(1,*) ALPHA
      READ(1,'(/)')
      READ(1,*) UNBAL
      READ(1,*) LAMBDA
      READ(1,'(///)')

      READ(1,*) NX
      READ(1,*) NY
      READ(1,*) NXY
      READ(1,*) M
      READ(1,*) N
      READ(1,'(///)')

      DO J = 1,NMAT

        WRITE(2,500)'PROPERTIES FOR MATERIAL',J,' (psi):'
        WRITE(2,510)'E1','E2','G12','NU12'        
        WRITE(2,520)E1(J),E2(J),G12(J),NU12(J) 
        WRITE(2,*)
        WRITE(2,530)'PLY THICK (in.)','DENSITY (lb/in^3)',
     $  '  COST FACTOR'
        WRITE(2,540)THICK(J),RHO(J),COST(J)
        WRITE(2,*)
        WRITE(2,530)'EPSILON_1-MAX   ','EPSILON_2-MAX     ', 
     $ 'GAMMA_12-MAX'
        WRITE(2,540)EP1MAX(J),EP2MAX(J),G12MAX(J)

      END DO

      WRITE(2,*)
      WRITE(2,550)'PLATE LOADINGS (lbs/in, lbs):' 
      WRITE(2,560)'NX','NY','NXY'
      WRITE(2,570)NX,NY,NXY 
      WRITE(2,*)
      WRITE(2,*)
      WRITE(2,580)'PLATE DIMENSIONS (in.)     MODE SHAPES', 
     $ 'WEIGHTING FACTOR '
      WRITE(2,*)
      WRITE(2,590)'X        Y','M        N              ALPHA'
      WRITE(2,600)X(1),Y(1),M,N,ALPHA 
      WRITE(2,*)

!
! check value of weighting factor, alpha
! 

      IF((ALPHA .GT. 1) .OR. (ALPHA .LT. 0)) THEN
       WRITE(*,*)'ALPHA = ',ALPHA
       WRITE(*,*)'WEIGHTING FUNCTION MUST BE 0<ALPHA<1'
       WRITE(*,*)'PROGRAM TERMINATED.'
       STOP
      END IF


500   FORMAT(/,A,1X,I2,A,/)
510   FORMAT(5X,A,11X,A,10X,A,7X,A,/,50('-'))
520   FORMAT(3(E10.4,3X),F6.4)
530   FORMAT(3(A,5X),/,70('-')) 
540   FORMAT(E9.3,12X,E9.3,15X,E9.3,/)
550   FORMAT(A,/)
560   FORMAT(3X,A,11X,A,10X,A,/,45('-'))
570   FORMAT(E9.3,3X,E9.3,3X,E9.3)
580   FORMAT(A,9X,A)
590   FORMAT(4X,A,14X,A,/,70('-'))
600   FORMAT(F8.4,2X,F8.4,9X,I2,7X,I2,13X,F6.4,/,/,/,/)
610   FORMAT(30X,A,/,/,/,/)

C
C calculate the stiffness invariants given the material properties
C
      CALL UCALC(U,E1,E2,NU12,G12,NMAT)
      END IF

C
C  need to read in the laminate number as 'I' (whatever the loop variable is)
C                  the total number of plies for a laminate, NPLY (.....%laminate_size) 
C     

C
C  decode the laminate string from integers to angle values
C

      CALL DECODE(IS,ISD,NPLY,I,NUM,MATL,MATLD,NMAT,
     $              N15S,N30S,N45S,N60S,N75S)

C
C  calculate the Q-Bar, A, B, and D Matrices for laminate using classical lamination theory
C

      CALL CLT(I,A,D,U,ISD,THICK,NUM,MATLD,LAMTH)

C
C calculate the A-inverse, and D-inverse matrices
C

      CALL MATINV(A,AI)
      CALL MATINV(D,DI)

C
C calculate no. of contiguous plies in excess of 4 for a laminate 
C this is not used as a constrain currently but may be
C implemented if desired
C
      MN = M*N
      CALL CRSHR(I,D,NX,NY,NXY,X,Y,M,N,CRSH,MN,BLCHECK,
     $           LAMTH)

C calculate the m. of s. for the critical shear load

      CALL MSCRSH(I,CRSH,NX,NY,NXY,CRSHMS,BLCHECK)

C calculate epsilon-1, epsilon-2, and gamma-12 for each ply in
C the laminate

      CALL STRAIN(I,EP1,EP2,GAM12,AI,DI,NX,NY,NXY,
     $               ISD,NUM)

C calculate the margin of safety for the ply that undergoes
C the largest strain in a laminate
C (CONSTRAINT #3)

      CALL MSSTR(I,EP1,EP2,GAM12,NUM,EP1MAX,EP2MAX,
     $               G12MAX,STRMS,DIRCR,PLYCR,MATLD)

C calculate the fitness of each laminate in the population

      CALL FITSUB(I,STRMS,FIT,CRSHMS,VIOL,THICK,LAMBDA,
     $                  N15S,N30S,N45S,N60S,N75S,UNBAL,NMAT,
     $                  LAYCOST,MATCOST,MATLD,NUM,WTMAT1,WTMAT2,
     $                  WTTOT,RHO,X,Y,COST,M1COST,M2COST,ISD,
     $                  ALPHA)

      if ( WRITE)  then
      CALL OUT1(ISD,I,NUM,CRSH,CRSHMS,EP1,EP2,GAM12,FIT,
     $      STRMS,DIRCR,PLYCR,E1,E2,G12,NU12,
     $      NX,NY,NXY,X,Y,EP1MAX,EP2MAX,G12MAX,THICK,
     $      VIOL,MATLD,NMAT,BLCHECK,MATCOST,M1COST,M2COST,
     $      LAYCOST,WTMAT1,WTMAT2,WTTOT,ALPHA)
      end if
      

      END


C----------------------- END OF MAIN PROGRAM --------------------------

C----------------------------------------------------------------------

      SUBROUTINE UCALC(U,E1,E2,NU12,G12,NMAT)

*************************************************************
C        INPUTS:
C          E1(J),E2(J),NU12(J),G12(J) - properties for each
C                                     - of the two materials

C        OUTPUTS:  
C          U(I,J) - stiffness invariants, J - represents 
C                 - material type
*************************************************************

      DOUBLE PRECISION U(5,10),E1(10),E2(10),NU12(10),
     $                 NU21,G12(10),Q11,Q12,Q22,Q66     
      INTEGER NMAT,I


C calculate poisson's ratio

      DO 5, I = 1,NMAT

       NU21 = (E2(I)/E1(I))*NU12(I)

C calculate the reduced stiffness matrix  

       Q11 = E1(I)/(1.D+00-NU12(I)*NU21)
       Q12 = (NU12(I)*E2(I))/(1.D+00-NU12(I)*NU21)
       Q22 = E2(I)/(1.D+00-NU12(I)*NU21)
       Q66 = G12(I)

C calculate the stiffness invariants for each material

       U(1,I) = (1.D+00/8.D+00)*(3.D+00*Q11 + 3.D+00*Q22 +
     $            2.D+00*Q12 + 4.D+00*Q66)
       U(2,I) = (1.D+00/2.D+00)*(Q11-Q22)
       U(3,I) = (1.D+00/8.D+00)*(Q11 + Q22 - 2.D+00*Q12 - 4.D+00*Q66)
       U(4,I) = (1.D+00/8.D+00)*(Q11 + Q22 + 6.D+00*Q12 - 4.D+00*Q66)
       U(5,I) = (1.D+00/8.D+00)*(Q11 + Q22 - 2.D+00*Q12 + 4.D+00*Q66)

5     CONTINUE

      RETURN
      END

C----------------------------------------------------------------------

C----------------------------------------------------------------------
      SUBROUTINE DECODE(IS,ISD,NPLY,I,NUM,MATL,MATLD,NMAT,
     $                  N15S,N30S,N45S,N60S,N75S)

***********************************************************************
C     INPUTS:
C             NPLY - the maximum number of single plies in half
C                   - the laminate
C                 I - the Ith laminate in the population 
C           IS(I,J) - the array which stores laminate (I's)
C                   - coded ply orientations 
C         MATL(I,J) - the array that stores laminate (I's)
C                   - material type   
C              NMAT - the number of different types of materials

C     OUTPUTS:
C             ISD(J) - the array which stores a laminate (I's)
C                    - decoded ply orientations (ply orientations  
C                    - are in degrees)
C           MATLD(J) - the array that stores laminate (I's) material
C                    - type 
C            NUM(I)  - stores the no. of filled single plies in 
C                    - laminate (I)
C            NXYS(J) - stores the number of +/- XY degree plies
C                    - for each different material
***********************************************************************

      INTEGER POPMAX,PLYMAX,TOTMAX
      PARAMETER (POPMAX = 500,PLYMAX = 200,TOTMAX = 400)
      INTEGER I,NPLY,ISD(TOTMAX),NUM(TOTMAX),NMAT,IS(POPMAX,PLYMAX),
     $        L,LL,J,MATL(POPMAX,PLYMAX),MATLD(TOTMAX),
     $        N15S(10),N30S(10),N45S(10),N60S(10),N75S(10),
     $        PLMIN15(10),PLMIN30(10),PLMIN45(10),PLMIN60(10),
     $        PLMIN75(10)

C initially set the number of filled plies to zero

      NUM(I) = 0
      L = 0
      LL = NPLY

      DO 5, J = 1,NMAT
         PLMIN15(J) = 1
         N15S(J) = 0 
         PLMIN30(J) = 1
         N30S(J) = 0 
         PLMIN45(J) = 1
         N45S(J) = 0 
         PLMIN60(J) = 1
         N60S(J) = 0 
         PLMIN75(J) = 1
         N75S(J) = 0 
5     CONTINUE

      DO 10, J = NPLY,1,-1

C if stack (J) is empty, set the max no. of available stacks to one less

            IF(IS(I,J) .EQ. 0) THEN
               LL = LL -1
               GOTO 10

C if stack (J) in laminate (I) is coded as 1, set the appropriate two 
C plies (assuming symmetry) to 0 degrees

            ELSE IF(IS(I,J) .EQ. 1) THEN
               L = L + 1  
               ISD(L) = 0
               ISD(2*LL+1-L) = 0
               NUM(I) = NUM(I) + 2 

C if stack (J) in laminate (I) is coded as 2, set the appropriate two 
C plies (assuming symmetry) to 15 or -15 degrees

           ELSE IF(IS(I,J) .EQ. 2) THEN
               L = L + 1  
               ISD(L) = 15*PLMIN15(MATL(I,J))
               ISD(2*LL+1-L) = 15*PLMIN15(MATL(I,J))
               NUM(I) = NUM(I) + 2
               IF(PLMIN15(MATL(I,J)) .EQ. 1) THEN
                 PLMIN15(MATL(I,J)) = -1
               ELSE
                 PLMIN15(MATL(I,J)) = 1
               END IF

C count the number of +/- 15 degree plies for each material 
C in the laminate

               N15S(MATL(I,J)) = N15S(MATL(I,J)) + 1

C if stack (J) in laminate (I) is coded as 3, set the appropriate two 
C plies (assuming symmetry) to 30 or -30 degrees

           ELSE IF(IS(I,J) .EQ. 3) THEN
               L = L + 1  
               ISD(L) = 30*PLMIN30(MATL(I,J))
               ISD(2*LL+1-L) = 30*PLMIN30(MATL(I,J))
               NUM(I) = NUM(I) + 2
               IF(PLMIN30(MATL(I,J)) .EQ. 1) THEN
                 PLMIN30(MATL(I,J)) = -1
               ELSE
                 PLMIN30(MATL(I,J)) = 1
               END IF

C count the number of +/- 30 degree plies for each material 
C in the laminate

               N30S(MATL(I,J)) = N30S(MATL(I,J)) + 1

C if stack (J) in laminate (I) is coded as 4, set the appropriate two 
C plies (assuming symmetry) to 45 or -45 degrees

           ELSE IF(IS(I,J) .EQ. 4) THEN
               L = L + 1  
               ISD(L) = 45*PLMIN45(MATL(I,J))
               ISD(2*LL+1-L) = 45*PLMIN45(MATL(I,J))
               NUM(I) = NUM(I) + 2
               IF(PLMIN45(MATL(I,J)) .EQ. 1) THEN
                 PLMIN45(MATL(I,J)) = -1
               ELSE
                 PLMIN45(MATL(I,J)) = 1
               END IF

C count the number of +/- 45 degree plies for each material 
C in the laminate

               N45S(MATL(I,J)) = N45S(MATL(I,J)) + 1

C if stack (J) in laminate (I) is coded as 5, set the appropriate two 
C plies (assuming symmetry) to 60 or -60 degrees

           ELSE IF(IS(I,J) .EQ. 5) THEN
               L = L + 1  
               ISD(L) = 60*PLMIN60(MATL(I,J))
               ISD(2*LL+1-L) = 60*PLMIN60(MATL(I,J))
               NUM(I) = NUM(I) + 2
               IF(PLMIN60(MATL(I,J)) .EQ. 1) THEN
                 PLMIN60(MATL(I,J)) = -1
               ELSE
                 PLMIN60(MATL(I,J)) = 1
               END IF

C count the number of +/- 60 degree plies for each material 
C in the laminate

               N60S(MATL(I,J)) = N60S(MATL(I,J)) + 1

C if stack (J) in laminate (I) is coded as 6, set the appropriate two 
C plies (assuming symmetry) to 75 or -75 degrees

           ELSE IF(IS(I,J) .EQ. 6) THEN
               L = L + 1  
               ISD(L) = 75*PLMIN75(MATL(I,J))
               ISD(2*LL+1-L) = 75*PLMIN75(MATL(I,J))
               NUM(I) = NUM(I) + 2
               IF(PLMIN75(MATL(I,J)) .EQ. 1) THEN
                 PLMIN75(MATL(I,J)) = -1
               ELSE
                 PLMIN75(MATL(I,J)) = 1
               END IF

C count the number of +/- 75 degree plies for each material 
C in the laminate

               N75S(MATL(I,J)) = N75S(MATL(I,J)) + 1


C if stack (J) in laminate (I) is coded as 7, set the appropriate two 
C plies (assuming symmetry) to 90 degrees

            ELSE IF(IS(I,J) .EQ. 7) THEN           
               L = L + 1  
               ISD(L) = 90
               ISD(2*LL+1-L) = 90
               NUM(I) = NUM(I) + 2
            END IF

C fill array that stores the 2 different material types

            IF(MATL(I,J) .EQ. 1) THEN
               MATLD(L) = 1
               MATLD(2*LL+1-L) = 1
            ELSE IF(MATL(I,J) .EQ. 2) THEN
               MATLD(L) = 2
               MATLD(2*LL+1-L) = 2
            END IF             
         
10       CONTINUE

      RETURN
      END 

C----------------------------------------------------------------------


C--------------------------------------------------------------------

      SUBROUTINE CLT(I,A,D,U,ISD,THICK,NUM,MATLD,LAMTH)

***********************************************************************
C     INPUTS:
C             I - the Ith laminate in the population 
C        ISD(J) - the array which stores laminate (I's)
C               - decoded ply orientations (J's)
C      MATLD(J) - the array that stores laminate (I's) material
C               - type 
C         THICK - the thickness of each individual ply
C       NUM(I)  - stores the no. of filled single plies in laminate (I)
C          U(5) - stiffness invariants

C     OUTPUTS:
C            A(3,3) - stores the A-matrix for laminate(I)
C            D(3,3) - stores the D-matrix for laminate(I)
C          LAMTH(I) - stores the total laminate thickess
***********************************************************************

      INTEGER TOTMAX,POPMAX
      PARAMETER(TOTMAX = 400, POPMAX = 500)      
      DOUBLE PRECISION A(3,3),D(3,3),QB(3,3),U(5,10),PI,C,
     $ THICK(10),LL,LLL,LAMTH(POPMAX)          
      INTEGER NUM(TOTMAX),ISD(TOTMAX),I,J,K,L,MATLD(TOTMAX)
      
      PI = 3.141592654D+00
      C = 180.0D+00
      LL = 0.D+00

C calculate thickness of the laminate

      DO 1, J = 1,NUM(I)/2
        LL = LL + THICK(MATLD(J))
1     CONTINUE


      LAMTH(I) = 2.D+00*LL
      LLL = LL

      DO 25, J = 1,NUM(I)

        LLL = LLL - THICK(MATLD(J))

C     calculate Q-BAR matrix terms (classical lamination theory)

        QB(1,1) = U(1,MATLD(J)) + U(2,MATLD(J))*COS((2*ISD(J))*PI/C)
     $            + U(3,MATLD(J))*COS(DBLE(4*ISD(J))*PI/C)

        QB(1,2) = U(4,MATLD(J)) -  
     $            U(3,MATLD(J))*COS(DBLE(4*ISD(J))*PI/C)

        QB(1,3) = 0.5D+00*U(2,MATLD(J))*SIN(DBLE(2*ISD(J))*PI/C) +
     $                    U(3,MATLD(J))*SIN(DBLE(4*ISD(J))*PI/C)

        QB(2,1) = QB(1,2)

        QB(2,2) = U(1,MATLD(J)) -   
     $            U(2,MATLD(J))*COS(DBLE(2*ISD(J))*PI/C) +
     $            U(3,MATLD(J))*COS(DBLE(4*ISD(J))*PI/C)
  
        QB(2,3) = 0.5D+00*U(2,MATLD(J))*SIN(DBLE(2*ISD(J))*PI/C) - 
     $                    U(3,MATLD(J))*SIN(DBLE(4*ISD(J))*PI/C) 

        QB(3,1) = QB(1,3)

        QB(3,2) = QB(2,3)

        QB(3,3) = U(5,MATLD(J)) - 
     $            U(3,MATLD(J))*COS(DBLE(4*ISD(J))*PI/C) 

C set all A and D matrix terms initially to zero

         IF(J .EQ. 1) THEN

           DO 10, K = 1,3
            DO 5, L = 1,3
               A(K,L) = 0.D+00
               D(K,L) = 0.D+00
5            CONTINUE
10         CONTINUE   

         END IF

         DO 20, K = 1,3
          DO 15, L = 1,3

C calculate A matrix (classical lamination theory)

            A(K,L) = A(K,L) + THICK(MATLD(J))*QB(K,L)

C calculate D matrix (classical lamination theory)

            D(K,L) = D(K,L) + (LL**3 - LLL**3)*QB(K,L)

15       CONTINUE
20     CONTINUE

       LL = LLL

25    CONTINUE

C make sure to multiply the D matrix by 1/3

           DO 35,K = 1,3
             DO 30, L = 1,3
               D(K,L) = (1.D+00/3.D+00)*D(K,L)
30         CONTINUE 
35       CONTINUE 

        RETURN
        END
        
C----------------------------------------------------------------------

C---------------------------------------------------------------------- 
      SUBROUTINE MATINV(A,AI)

***********************************************************************
C     INPUTS:
C       A(3,3) - stores the A-matrix for laminate(I)
C       D(3,3) - stores the D-matrix for laminate(I)

C     OUTPUTS:
C      AI(3,3) - stores the inverse of the A-matrix for laminate(I)
C      DI(3,3) - stores the inverse of the D-matrix for laminate(I)
***********************************************************************

*** WARNING: THIS SUBROUTINE INVERTS SYMMETRIC MATRICES ONLY ***

      INTEGER L
      PARAMETER(L = 3)      
      DOUBLE PRECISION A(3,3),AI(3,3),MAT(L,L),DD
      INTEGER I,II,J,K,LL

      DO 10, II = 1,L
        DO 5, LL = 1,L
          MAT(II,LL) = A(II,LL)
5       CONTINUE
10    CONTINUE
         
      DO 35, K=1,L

         DD=MAT(K,K)
         DO 15, J=1,L
          MAT(K,J)=-MAT(K,J)/DD
15       CONTINUE

          DO 30, I=1,L

            IF (K.EQ.I) GOTO 25

            DO 20, J=1,L
              IF (K.NE.J) MAT(I,J)=MAT(I,J)+MAT(I,K)*MAT(K,J)
20          CONTINUE

C array MAT also stores the inverse values  

25          MAT(I,K)=MAT(I,K)/DD

30        CONTINUE 

            MAT(K,K)=1.0D0/DD

35    CONTINUE

C translate inverse values stored in MAT to the respective values
C in the AI and DI arrays

      DO 45, I = 1,3
        DO 40, J = 1,3

          AI(I,J) = MAT(I,J)

40      CONTINUE
45    CONTINUE

      RETURN
      END
 
C----------------------------------------------------------------------

C----------------------------------------------------------------------      

      SUBROUTINE CRSHR(AA,D,NX,NY,NXY,X,Y,M,N,CRSH,MN,
     $                 BLCHECK,LAMTH)
     
C  this subroutine calculates the critical shear buckling load of 
C  specially orthotropic plates (D16 and D26 not equal to zero)
C  the analysis uses Classical Lamination Theory and the Galerkin
C  energy method to calculate the critical buckling loads: NXcr, NYcr, 
C  and NXYcr.

C  the Galerkin energy method was taken from 
C  Whitney: "Structural Analysis of Anisotropic Plates"
C            Chapter 6.5 

C  assumptions: i) plate is symmetric about the mid-plane ([B] = 0)
C              ii) simply supported edges     

      

      INTEGER POPMAX
      PARAMETER(POPMAX = 500)      

      DOUBLE PRECISION D(3,3),NX,NY,NXY,R,CC(100),DD(100),EE(100),
     $                 MM(100,100),MAT1(100,100),MAT2(100,100),
     $                 NUMER,DENOM,PI,ALFR(100),ALFI(100),BETA(100),
     $                 Z(100,100),LAMRL(100),LAMIM(100),NXCR,NYCR,
     $                 NXYCR,CRLAMDA,MINVAL,CRSH(POPMAX),FF(100,100),
     $                 LAMTH(POPMAX),X(POPMAX),Y(POPMAX)
     
      INTEGER I,J,M,N,MN,MATZ,II,JJ,P1,P2,COND1,
     $        COND2,IERR,CRM,CRN,AA,P1B,BLCHECK(POPMAX),CRP




      IF(M .GT. 10) THEN
       WRITE(3,*)'M TOO LARGE. IT HAS BEEN SET TO THE MAX VALUE OF 10'
       M = 10
      END IF

      IF(N .GT. 10) THEN
       WRITE(3,*)'N TOO LARGE. IT HAS BEEN SET TO THE MAX VALUE OF 10'
       N = 10
      END IF

C
C initialize all matrices
C

      DO 35, I = 1,MN
       DO 30, J = 1,MN
        FF(I,J) = 0.D+00
        MM(I,J) = 0.D+00
        MAT1(I,J) = 0.D+00
        MAT2(I,J) = 0.D+00
30      CONTINUE
35    CONTINUE

      MATZ = 1    
      PI = DACOS(-1.D+00) 
      R = X(AA)/Y(AA)                
      P1 = 0                       
                                                   
      DO 65, I = 1,M
        DO 60, J = 1,N
          P1 = P1+1
          CC(P1) = (PI**4)*(D(1,1)*I**4 + 
     $             2*(D(1,2) + 2*D(3,3))*(I**2)*(J**2)*(R**2) + 
     $             D(2,2)*(J**4)*(R**4))
          DD(P1) = -1*(PI**2)*(NX*(I**2)*(X(AA)**2) + 
     $             NY*(J**2)*(X(AA)**2)*(R**2)) 
          EE(P1) = 32.D+00*I*J*(R**3)*(Y(AA)**2)*NXY
C
C changes for anisotrpic plates made here - addition of FF(I,J) array
C
          P1B = 0
          DO 45, II = 1,M
            DO 40, JJ = 1,N
              P1B = P1B + 1
              FF(P1,P1B) = -32.D+00*I*J*R*(PI**2)*
     $                      ((II**2+I**2)*D(1,3) +
     $                      (R**2)*(JJ**2+J**2)*D(2,3))
40          CONTINUE
45        CONTINUE

C
C end changes
C 

          P2 = 0
          DO 55, II = 1,M
            COND1 = I+II
            DO 50, JJ = 1,N
              COND2 = J+JJ
              P2 = P2+1
              IF((MOD(COND1,2) .EQ. 0) .OR. (MOD(COND2,2) .EQ. 0)) THEN
                MM(P1,P2) = 0.D+00
              ELSE 
                NUMER = II*JJ
                DENOM = ( (I**2)-(II**2) )*( (J**2)-(JJ**2) )
                MM(P1,P2) = DBLE(NUMER)/DBLE(DENOM)
              END IF
50          CONTINUE
55        CONTINUE

60      CONTINUE
65    CONTINUE 

C
C      IF(FMT3 .EQ. 1) THEN
C        WRITE(4,*)'[M]:'
C        WRITE(4,*)
C        DO 70, I = 1,MN
C          WRITE(4,500)(MM(I,J), J=1,MN)
C70      CONTINUE       
C        WRITE(4,*)
C      END IF
C

C
C changes made for anisotropic plates - FF(I,J)
C
      
      DO 85, I = 1,MN
        DO 80, J = 1,MN
          IF(I .EQ. J) THEN
            MAT1(I,J) = CC(I) + FF(I,J)*MM(I,J)
          ELSE
            MAT1(I,J) = FF(I,J)*MM(I,J)
          END IF
80      CONTINUE
85    CONTINUE


C
C end changes
C

      DO 95, I = 1,MN
        DO 90, J = 1,MN
          IF(I .EQ. J) THEN
            MAT2(I,J) = DD(I) + EE(I)*MM(I,J)
          ELSE
            MAT2(I,J) = EE(I)*MM(I,J)
          END IF
90      CONTINUE
95    CONTINUE  

C
C output arrays
C

C
C      IF(FMT3 .EQ. 1) THEN
C        WRITE(4,*)'[MAT1]:'
C        WRITE(4,*)
C        DO 90, I = 1,MN
C          WRITE(4,500)(MAT1(I,J), J=1,MN)
C90      CONTINUE       
C        WRITE(4,*)
C
C        WRITE(4,*)'[MAT2]:'
C        WRITE(4,*)
C        DO 95, I = 1,MN
C          WRITE(4,500)(MAT2(I,J), J=1,MN)
C95      CONTINUE       
C        WRITE(4,*)
C      END IF
C

      CALL rgg(100,MN,MAT1,MAT2,ALFR,ALFI,BETA,MATZ,Z,IERR)


      IF(IERR .NE. 0) THEN
        WRITE(3,*)'***** RGG ERROR:',IERR,'*****'
      END IF

      P1 = 0 
      MINVAL = 1.D+50

      DO 105, I = 1,M
        DO 100, J = 1,N
          P1 = P1 + 1
          IF(DABS(BETA(P1)) .LT. 1.D-06) THEN
            WRITE(3,*)'ERROR IN EIGEN SOLVER. CHECK SIGNS ON NX,NY,NXY'
            WRITE(3,*)'PROGAM TERMINATED'
            STOP
          ELSE 
            LAMRL(P1) = ALFR(P1)/BETA(P1)
            LAMIM(P1) = ALFI(P1)/BETA(P1)
C
C            IF(FMT3 .EQ. 1) THEN
C              WRITE(4,510)'M =',I,'N =',J,LAMRL(P1),LAMIM(P1),'i' 
C            END IF
C
          END IF

          IF(LAMRL(P1) .LT. 0.D+00) GOTO 100
          IF(LAMRL(P1) .LT. MINVAL) THEN
            CRLAMDA = LAMRL(P1)         
            MINVAL = LAMRL(P1)
            CRP = P1 
            CRM = I
            CRN = J
          END IF  

100     CONTINUE
105   CONTINUE        

      NXCR = CRLAMDA*NX
      NYCR = CRLAMDA*NY
      NXYCR = CRLAMDA*NXY

C
C store critical shear value in array
C

      CRSH(AA) = NXCR
      BLCHECK(AA) = 1

      IF((ABS(CRSH(AA)) .GT. ABS(NYCR)) .AND.
     $   (NYCR .NE. 0)) THEN
        CRSH(AA) = NYCR
        BLCHECK(AA) = 2
      END IF
 
      IF((ABS(CRSH(AA)) .GT. ABS(NXYCR)).AND.
     $   (NXYCR .NE. 0)) THEN
        CRSH(AA) = NXYCR
        BLCHECK(AA) = 3
      END IF 

C
C output eigen value analysis
C

       WRITE(4,*)
       WRITE(4,*)'EIGEN VECTORS:'
       WRITE(4,*)'CRP',CRP
       DO 110, I = 1,MN
         WRITE(4,*) I
         WRITE(4,500)Z(CRP,I)
110    CONTINUE   

       WRITE(4,*)
       WRITE(4,*)
       WRITE(4,*)'CRITICAL VALUES FOR M, N, AND LAMDA'
       WRITE(4,*)
       WRITE(4,520)'M =',CRM,'N =',CRN,'LAMDA =',CRLAMDA 
       WRITE(4,*)
       WRITE(4,*)'CRITICAL BUCKLING LOADS: NX, NY, AND NXY.'
       WRITE(4,*)
       WRITE(4,530)'NXCR =',NXCR
       WRITE(4,530)'NYCR =',NYCR
       WRITE(4,530)'NXYCR =',NXYCR


500   FORMAT(100(F8.4,1X,:))
510   FORMAT(A,1X,I3,2X,A,1X,I3,3X,E13.6,1X,E13.6,A)                
520   FORMAT(A,1X,I3,2X,A,1X,I3,3X,A,1X,E11.6)                
530   FORMAT(A,1X,E13.6)


      RETURN               
      END
C----------------------------------------------------------------------

C---------------------------------------------------------------------

      SUBROUTINE MSCRSH(I,CRSH,NX,NY,NXY,CRSHMS,BLCHECK)
***********************************************************************
C     INPUTS:
C              I - the Ith laminate in the population
C        CRSH(I) - the critical shear load of the plate (laminate (I))
C            NXY - the shear load applied to the plate

C     OUTPUTS:
C      CRSHMS(I) - the critical shear load m. of s.
***********************************************************************

      INTEGER POPMAX
      PARAMETER (POPMAX = 500)      
      DOUBLE PRECISION CRSH(POPMAX),CRSHMS(POPMAX),NXY,
     $                 NX,NY
      INTEGER I,BLCHECK(POPMAX)      

C calculate the margin of safety for the critical shear load
 
      IF(BLCHECK(I) .EQ. 1) THEN
       CRSHMS(I) = 1.D+00 - ABS(NX/CRSH(I))
      ELSE IF(BLCHECK(I) .EQ. 2) THEN
       CRSHMS(I) = 1.D+00 - ABS(NY/CRSH(I))
      ELSE IF(BLCHECK(I) .EQ. 3) THEN
       CRSHMS(I) = 1.D+00 - ABS(NXY/CRSH(I))
      END IF

      RETURN
      END     

C----------------------------------------------------------------------

C----------------------------------------------------------------------
      SUBROUTINE STRAIN(I,EP1,EP2,GAM12,AI,DI,NX,NY,NXY,
     $                  ISD,NUM)

***********************************************************************
C     INPUTS:
C           I - the Ith laminate in the population
C   NX,NY,NXY - load applied to the laminate
C   MX,MY,MXY - moments applied to the laminate
C     AI(3,3) - the inverse of the A-matrix
C     DI(3,3) - the inverse of the D-matrix
C     NUM(I)  - stores the no. of filled single plies in laminate (I)
C    ISD(I,J) - the array which stores laminate (I's)
C             - decoded ply orientations 

C     OUTPUTS:
C     EP1(I,J) - stores epsilon-1 strain for ply (J) in laminate (I)
C     EP2(I,J) - stores epsilon-2 strain for ply (J) in laminate (I)
C   GAM12(I,J) - stores gamma-12 strain for ply (J) in laminate (I)
***********************************************************************

      INTEGER POPMAX,TOTMAX
      PARAMETER (POPMAX = 500,TOTMAX = 400)
      DOUBLE PRECISION EP1(POPMAX,TOTMAX),
     $ EP2(POPMAX,TOTMAX),GAM12(POPMAX,TOTMAX),
     $ AI(3,3),DI(3,3),T(3,3),NX,NY,NXY,PI,C,
     $ EPX,EPY,GAMXY
       INTEGER NUM(TOTMAX),ISD(TOTMAX),I,J


      PI = 3.141592654D+00
      C = 180.0D+00

C calculate the off-axes strains using classical lamination theory

      EPX = AI(1,1)*NX + AI(1,2)*NY + AI(1,3)*NXY 
      EPY = AI(2,1)*NX + AI(2,2)*NY + AI(2,3)*NXY 
      GAMXY = AI(3,1)*NX + AI(3,2)*NY + AI(3,3)*NXY 

C
C initially set all prinicipal strains to zero
C

      DO 5, J = 1,(NUM(I)/2)
         EP1(I,J) = 0.D+00
         EP2(I,J) = 0.D+00
         GAM12(I,J) = 0.D+00
5     CONTINUE

C              
C using the transformation matrix - calcuate the principal
C direction strains 
C

      DO 10, J = 1,(NUM(I)/2)
         T(1,1) = ( COS(DBLE(ISD(J))*PI/C) )**2
         T(1,2) = ( SIN(DBLE(ISD(J))*PI/C) )**2
         T(1,3) = ( (COS(DBLE(ISD(J))*PI/C))*
     $           (SIN(DBLE(ISD(J))*PI/C)) )
         T(2,1) = ( SIN(DBLE(ISD(J))*PI/C) )**2
         T(2,2) = ( COS(DBLE(ISD(J))*PI/C) )**2
         T(2,3) = -1.0D+00*( (COS(DBLE(ISD(J))*PI/C))*
     $           (SIN(DBLE(ISD(J))*PI/C)) )
         T(3,1) = -2.D+00*( (COS(DBLE(ISD(J))*PI/C))*
     $           (SIN(DBLE(ISD(J))*PI/C)) )
         T(3,2) = 2.D+00*( (COS(DBLE(ISD(J))*PI/C))*
     $           (SIN(DBLE(ISD(J))*PI/C)) )
         T(3,3) = ( COS(DBLE(ISD(J))*PI/C) )**2 - 
     $           ( SIN(DBLE(ISD(J))*PI/C) )**2

         EP1(I,J) = T(1,1)*EPX + T(1,2)*EPY + T(1,3)*GAMXY 
         EP2(I,J) = T(2,1)*EPX + T(2,2)*EPY + T(2,3)*GAMXY 
         GAM12(I,J) = T(3,1)*EPX + T(3,2)*EPY + T(3,3)*GAMXY 

10    CONTINUE

      RETURN
      END

C----------------------------------------------------------------------

C----------------------------------------------------------------------      

      SUBROUTINE MSSTR(I,EP1,EP2,GAM12,NUM,EP1MAX,EP2MAX,
     $                  G12MAX,STRMS,DIRCR,PLYCR,MATLD)  

***********************************************************************
C  INPUTS:
C            I - the Ith laminate in the population
C     EP1(I,J) - stores epsilon-1 strain for ply (J) in laminate (I)
C     EP2(I,J) - stores epsilon-2 strain for ply (J) in laminate (I)
C   GAM12(I,J) - stores gamma-12 strain for ply (J) in laminate (I)
C      NUM(I)  - stores the no. of filled single plies in laminate (I)
C    EP1MAX(I) - the maximum allowable epsilon-1 for mat'l (I) 
C    EP2MAX(I) - the maximum allowable epsilon-2 for mat'l (I)
C    G12MAX(I) - the maximum allowable gamma-12 for mat'l (I)
C     MATLD(J) - the array that stores laminate (I's) material
C              - type 

C  OUTPUTS:
C     STRMS(I) - min among plies: (all. strain/ ply strain)
C     DIRCR(I) - the principal direction that STRMS(I) 
C              - occurs in (1,2,1-2) 
C     PLYCR(I) - the position in laminate (I) of the ply that 
C              - undergoes the largest strain
***********************************************************************

      INTEGER POPMAX,TOTMAX
      PARAMETER (POPMAX = 500,TOTMAX = 400)
      DOUBLE PRECISION EP1(POPMAX,TOTMAX),EP2(POPMAX,TOTMAX),
     $ GAM12(POPMAX,TOTMAX),STRMS(POPMAX),
     $ EP1MAX(10),EP2MAX(10),G12MAX(10),CALC 
      INTEGER NUM(TOTMAX),DIRCR(POPMAX),PLYCR(POPMAX),I,J,
     $ MATLD(TOTMAX)

C set STRMS(I) to 0

      STRMS(I) = 0.D+00

      DO 10, J = 1,(NUM(I)/2) 

C calculate the ratio of maximum allowable strain in the 
C 1-direction to the strain in the 1-direction for
C ply (J) in laminate (I)

        CALC =ABS(EP1(I,J))/EP1MAX(MATLD(J))

C check to see whether CALC is less than the value stored in 
C STRMS(I) 
C if it is, then replace STRMS(I) with CALC
C check to see what strain direction is critical
C check to see what ply is being analysed

        IF(CALC .GT. STRMS(I)) THEN 
           STRMS(I) = CALC
           DIRCR(I) = 1
           PLYCR(I) = J
        END IF

C calculate the ratio of maximum allowable strain in the 
C 2-direction to the strain in the 2-direction for
C ply (J) in laminate (I)

        CALC =ABS(EP2(I,J))/EP2MAX(MATLD(J))

C check to see whether CALC is less than the value stored in 
C STRMS(I) 
C if it is, then replace STRMS(I) with CALC
C check to see what strain direction is critical
C check to see what ply is being analysed

        IF(CALC .GT. STRMS(I)) THEN 
           STRMS(I) = CALC
           DIRCR(I) = 2
           PLYCR(I) = J
        END IF 

C calculate the ratio of maximum allowable shear strain  
C to the shear strain for ply (J) in laminate (I)

        CALC =ABS(GAM12(I,J))/G12MAX(MATLD(J))
 
C check to see whether CALC is less than the value stored in 
C STRMS(I) 
C if it is, then replace STRMS(I) with CALC 
C check to see what strain direction is critical
C check to see what ply is being analysed

       IF(CALC .GT. STRMS(I)) THEN
           STRMS(I) = CALC
           DIRCR(I) = 3
           PLYCR(I) = J
       END IF
     
10    CONTINUE

C calculate the the margin of safety by subtracting 1 from STRMS(I)

      STRMS(I) = 1.D+00 - STRMS(I)  

      RETURN
      END

C----------------------------------------------------------------------

C---------------------------------------------------------------------

      SUBROUTINE FITSUB(I,STRMS,FIT,CRSHMS,VIOL,THICK,LAMBDA,
     $                  N15S,N30S,N45S,N60S,N75S,UNBAL,NMAT,
     $                  LAYCOST,MATCOST,MATLD,NUM,WTMAT1,WTMAT2,
     $                  WTTOT,RHO,X,Y,COST,M1COST,M2COST,ISD,
     $                  ALPHA)
**********************************************************************
C     INPUTS:
C            I - the Ith laminate in the population 
C     STRMS(I) - min among plies: (all. strain/ ply strain) 
C        THICK - the thickness of an individual ply 
C     LAMTH(I) - stores the thickness for laminate (I)
C       LAMBDA - exponential factor to help ga code converge 
C      N45S(I) - the number of +/- 45 plies for each material
C              - type in the laminate
C        UNBAL - a constant used to penalize unbalanced laminates
C         NMAT - the number of different types of materials

C    OUTPUTS:
C       FIT(I) - stores the fitness for laminate (I)
C         VIOL - stores design feasibility
**********************************************************************

C  *************************** NOTE *******************************
C the fitness function works by scaling all constraints and 
C then determining which constraint violates its boundaries the
C most (designated - CSTVMX).
C if no constraints are violated, the value for the fitness is just the 
C value of the thickness of the laminate
C *******************************************************************

      INTEGER POPMAX,PLYMAX,TOTMAX
      PARAMETER (POPMAX = 500, PLYMAX = 200, TOTMAX = 400)

      DOUBLE PRECISION FIT(POPMAX),CRSHMS(POPMAX),
     $ STRMS(POPMAX),CSTVMX,AVGWT,LAYCOST(POPMAX),
     $ THICK(10),LAMBDA,PENLTY(10),UNBAL,AVGLAY,
     $ MATCOST(POPMAX),WTMAT1(POPMAX),WTMAT2(POPMAX),
     $ WTTOT(POPMAX),RHO(10),COST(10),M1COST(POPMAX),
     $ M2COST(POPMAX),C1(7),C2(7),C3,RFACT,PLYS,AVGMAT,
     $ CSTPLY(10),WTPLY(10),OBJCTV,ALPHA,PLYWT(10),TC,
     $ X(POPMAX),Y(POPMAX)

      INTEGER I,VIOL(POPMAX),N15S(10),N30S(10),N45S(10),
     $ N60S(10),N75S(10),NMAT,J,MATLD(TOTMAX),NUM(TOTMAX),
     $ ISD(TOTMAX),TEMP,IANG
    

C
C  Calculate the layup cost for a laminate
C

      C1(1) = 3.9358D-04
      C1(2) = 3.9035D-04
      C1(3) = 3.5032D-04
      C1(4) = 2.9377D-04
      C1(5) = 2.0773D-04
      C1(6) = 1.0469D-04
      C1(7) = 0.D+00

      DO 5, J = 1,7
       TEMP = 8-J
       C2(J) = C1(TEMP)
5     CONTINUE

      C3 = 1.9241D-05
      LAYCOST(I) = 0.D+00

      DO 10, J = 1,NUM(I)

       IF(ABS(ISD(J)) .EQ. 0) THEN
        IANG = 1
       ELSE IF(ABS(ISD(J)) .EQ. 15) THEN
        IANG = 2 
       ELSE IF(ABS(ISD(J)) .EQ. 30) THEN
        IANG = 3 
       ELSE IF(ABS(ISD(J)) .EQ. 45) THEN
        IANG = 4 
       ELSE IF(ABS(ISD(J)) .EQ. 60) THEN
        IANG = 5 
       ELSE IF(ABS(ISD(J)) .EQ. 75) THEN
        IANG = 6 
       ELSE IF(ABS(ISD(J)) .EQ. 90) THEN
        IANG = 7 
       END IF

C  
C       RFACT = THICK(MATLD(J))/PLYS) 
C
C
C       RFACT is the ratio of the ply thickness to the thickness
C       of the material (PLYS) that layup machine can dispense.
C       Since a value for PLYS is not available we will assume a
C       value of 1 for RFACT 
C
C
       RFACT = 1.D+00

       LAYCOST(I) = LAYCOST(I) + 
     $ RFACT*(C1(IANG)*X(I) + C2(IANG)*Y(I) + C3*X(I)*Y(I))


10    CONTINUE

C
C calculate the average layup cost of a ply
C

      AVGLAY = LAYCOST(I)/NUM(I)

C
C  Calculate the cost and weight of a single ply 
C

      DO 15, J=1,NMAT
       WTPLY(J) = RHO(J)*X(I)*Y(I)*THICK(J)
       PLYWT(J) = WTPLY(J)
       CSTPLY(J) = COST(J)*WTPLY(J)
15    CONTINUE
 
C
C  Calculate the average weight and material cost of a ply
C


      AVGWT = 0.D+00
      AVGMAT = 0.D+00
      DO 20, J = 1,NMAT    
       AVGWT = AVGWT + WTPLY(J) 
       AVGMAT = AVGMAT + CSTPLY(J) 
20    CONTINUE

      AVGWT = AVGWT/DBLE(NMAT)
      AVGMAT = AVGMAT/DBLE(NMAT)

C
C      Calculate the weight of the laminate
C
       
      WTMAT1(I) = 0.D+00
      WTMAT2(I) = 0.D+00

      DO 25, J = 1,NUM(I)
        IF (MATLD(J) .EQ. 1) THEN
         WTMAT1(I) = WTMAT1(I) + WTPLY(1) 
        ELSE IF (MATLD(J) .EQ. 2) THEN
         WTMAT2(I) = WTMAT2(I) + WTPLY(2)
       END IF
25    CONTINUE

      WTTOT(I) = WTMAT1(I) + WTMAT2(I)

C
C  Calculate the material cost of the laminate
C
      M1COST(I) = COST(1)*WTMAT1(I)
      M2COST(I) = COST(2)*WTMAT2(I)
      MATCOST(I) = M1COST(I) + M2COST(I)


C
C calculate the modified alpha value for cost
C and weight combination
C
       TC = MATCOST(I) + LAYCOST(I)

C
C scale appropriate quantities
C calculate the objective function and scale with alpha
C


      CSTPLY(1)= (1.D+00-ALPHA)*CSTPLY(1)
      CSTPLY(2)= (1.D+00-ALPHA)*CSTPLY(2)
      AVGMAT = (1.D+00-ALPHA)*AVGMAT
      AVGLAY = (1.D+00-ALPHA)*AVGLAY
      PLYWT(1) = ALPHA*PLYWT(1)
      PLYWT(2) = ALPHA*PLYWT(2)
      AVGWT = ALPHA*AVGWT
      OBJCTV = ALPHA*WTTOT(I) +
     $        (1.D+00-ALPHA)*TC


C
C is the number of +/- theta plies an even number?
C
      DO 30, J = 1,NMAT
        PENLTY(J) = 0 
        IF(MOD(N15S(J),2) .NE. 0) THEN
          PENLTY(J) = PENLTY(J) + UNBAL
        END IF
        IF(MOD(N30S(J),2) .NE. 0) THEN
          PENLTY(J) = PENLTY(J) + UNBAL
        END IF
        IF(MOD(N45S(J),2) .NE. 0) THEN
          PENLTY(J) = PENLTY(J) + UNBAL
        END IF
        IF(MOD(N60S(J),2) .NE. 0) THEN
          PENLTY(J) = PENLTY(J) + UNBAL
        END IF
        IF(MOD(N75S(J),2) .NE. 0) THEN
          PENLTY(J) = PENLTY(J) + UNBAL
        END IF   
30    CONTINUE


C CSTVMX is set equal to a large number

      CSTVMX = 1000000.D+00

      IF(CRSHMS(I) .LT. CSTVMX) CSTVMX = CRSHMS(I)
      IF(STRMS(I) .LT. CSTVMX) CSTVMX = STRMS(I)

C if no constraint is violated, then the fitness is set equal to the
C negative thickness plus a small bonus function, and a penalty
C function if the number of +/- theta plies are not equal. 

      IF(CSTVMX .GE. 0.D+00) THEN
        FIT(I) =  -1.D+00*OBJCTV
     $            - (PLYWT(1) + CSTPLY(1) + AVGLAY)*PENLTY(1)
     $            - (PLYWT(2) + CSTPLY(2) + AVGLAY)*PENLTY(2) 
     $            + CSTVMX*(AVGWT + AVGLAY + AVGMAT)
        VIOL(I) = 2

C if a constraint is violated, then the fitness is set equal to the  
C negative thickness and multiplied by the largest
C constraint violation. A small step function is also added to
C make sure that very thin laminates get penalized appropriately, 
C and a penalty function is implemented if the number of +/- 45
C plies is not equal

      ELSE

        FIT(I) = -1.D+00*OBJCTV*((1.D+00-CSTVMX)**LAMBDA)
     $            - (PLYWT(1) + CSTPLY(1) + AVGLAY)*PENLTY(1)
     $            - (PLYWT(2) + CSTPLY(2) + AVGLAY)*PENLTY(2) 
        VIOL(I) = 1

      END IF

      RETURN
      END

C----------------------------------------------------------------------

C----------------------------------------------------------------------      

      SUBROUTINE OUT1(ISD,I,NUM,CRSH,CRSHMS,EP1,EP2,GAM12,FIT,
     $      STRMS,DIRCR,PLYCR,E1,E2,G12,NU12,
     $      NX,NY,NXY,X,Y,EP1MAX,EP2MAX,G12MAX,THICK,
     $      VIOL,MATLD,NMAT,BLCHECK,MATCOST,M1COST,M2COST,
     $      LAYCOST,WTMAT1,WTMAT2,WTTOT,ALPHA)

***********************************************************************
C     INPUTS:
C         NPLY - the maximum number of single plies in half
C              - the laminate
C            I - the Ith laminate in the population 
C     ISD(I,J) - the array which stores laminate (I's)
C              - decoded ply orientations 
C     MATLD(J) - the array that stores laminate (I's) material
C              - type 
C      NUM(I)  - stores the no. of filled single plies in laminate (I)
C  E1(J),E2(J) - Young's moduli for each material
C       G12(J) - shear modulus for each material
C      NU12(J) - poisson's ratio for each material  
C    NX,NY,NXY - applied in-plane load
C    MX,MY,MXY - applied moment 
C       FIT(I) - stores the fitness for laminate (I)
C      CRSH(I) - the critical shear load of the plate (laminate (I))
C    CRSHMS(I) - the critical shear load m. of s.
C     EP1(I,J) - stores epsilon-1 strain for ply (J) in laminate (I)
C     EP2(I,J) - stores epsilon-2 strain for ply (J) in laminate (I)
C   GAM12(I,J) - stores gamma-12 strain for ply (J) in laminate (I)
C     STRMS(I) - min among plies: (all. strain/ ply strain) 
C     DIRCR(I) - the principal dir. that STRMS(I) occurs in (1,2,1-2) 
C     PLYCR(I) - the position in the laminate of the ply with 
C              - max quantity (all. strain/ply strain) 
C        GENER - the number of generations completed
C         VIOL - stores design feasibility
C         NAMT - the number of different materials
C           PC - probability of repair
C           PM - probability of mutation
C         PADD - probability of adding a ply
C         PDEL - probability of deleting a ply
C       PMMATL - probability of mutating material array
C         PREP - probability of repair
C           PP - probability of permutation
***********************************************************************
  
      INTEGER POPMAX,TOTMAX
      PARAMETER (POPMAX = 500,TOTMAX = 400)


      DOUBLE PRECISION CRSH(POPMAX),CRSHMS(POPMAX),
     $ EP1(POPMAX,TOTMAX),EP2(POPMAX,TOTMAX),         
     $ GAM12(POPMAX,TOTMAX),FIT(POPMAX),STRMS(POPMAX),         
     $ E1(10),E2(10),G12(10),NU12(10),NX,NY,NXY,
     $ EP1MAX(10),EP2MAX(10),G12MAX(10),THICK(10),
     $ MATCOST(POPMAX),M1COST(POPMAX),M2COST(POPMAX),
     $ LAYCOST(POPMAX),TOTCOST(POPMAX),WTMAT1(POPMAX),WTMAT2(POPMAX),
     $ WTTOT(POPMAX),LAMTH,DENOM,ALPHA,X(POPMAX),Y(POPMAX)
  

      INTEGER NUM(TOTMAX),ISD(TOTMAX),
     $ DIRCR(POPMAX),PLYCR(POPMAX),I,J,
     $ VIOL(POPMAX),NMAT,MATLD(TOTMAX),K,
     $ BLCHECK(POPMAX)

     
       TOTCOST(I) = MATCOST(I) + LAYCOST(I)  

C
C calculate laminate thickness
C

      LAMTH = 0.D+00
      DO 5, K = 1,NUM(I)
       LAMTH = LAMTH + THICK(MATLD(K))
5     CONTINUE

C
C output results on all laminates
C

C
C write stacking sequence for laminate (I)
C

      WRITE(2,500)'LAMINATE:',I,'(HALF-LAMINATE STACKING',
     $             'SEQ. - SYMMETRIC)'
      WRITE(2,510)'ORIENT.:',(ISD(J),J = NUM(I)/2+1,NUM(I))
      WRITE(2,510)'MATL.:  ',(MATLD(J), J = NUM(I)/2+1,NUM(I))

       WRITE(2,*)
       WRITE(2,*)
       WRITE(2,520)'         - PLATE DIMENSIONS -'
       WRITE(2,*)'X DIMENSION         Y DIMENSION'
       WRITE(2,*)X(I),Y(I)
       WRITE(2,*)
       WRITE(2,520)'         - MANUFACTURING COST -'
       WRITE(2,*)
       WRITE(2,530)'MATERIAL 1    MATERIAL 2   TOTAL MAT COST',
     $     'LAYUP COST    LAMINATE COST'
       WRITE(2,540)M1COST(I),M2COST(I),MATCOST(I),
     $             LAYCOST(I),TOTCOST(I)
       WRITE(2,*)

       WRITE(2,520)'         - MATERIAL WEIGHT -'
       WRITE(2,*)
       WRITE(2,530)'MATERIAL 1    MATERIAL 2',
     $             'LAMINATE WEIGHT    LAMINATE THICKNESS'
       WRITE(2,550)WTMAT1(I),WTMAT2(I),WTTOT(I),LAMTH


       WRITE(2,520)'    - OBJECTIVE FUNCTION RATIOS -'
       WRITE(2,*)
       DENOM = (1.D+00-ALPHA)*TOTCOST(I) + ALPHA*WTTOT(I)

       WRITE(2,560)'ALPHA (CST)','ALPHA (WT)',
     $   'MCOST (% OF)','LCOST (% OF)','WT (% OF)'
       WRITE(2,570)1.D00-ALPHA,ALPHA,
     $ ( (1.D+00-ALPHA)*MATCOST(I)/DENOM )*100.D+00,
     $ ( (1.D+00-ALPHA)*LAYCOST(I)/DENOM )*100.D+00,
     $ ( ALPHA*WTTOT(I)/DENOM )*100.D+00
       WRITE(2,*)

       WRITE(2,520)'    - FITNESS AND CONSTRAINTS -'
       WRITE(2,*)
       IF(BLCHECK(I) .EQ. 1) THEN
        WRITE(2,580)'FITNESS','CRITICAL LOAD (X) & M.S.'
       ELSE IF(BLCHECK(I) .EQ. 2) THEN
        WRITE(2,580)'FITNESS','CRITICAL LOAD (Y) & M.S.'
       ELSE IF(BLCHECK(I) .EQ. 3) THEN
        WRITE(2,580)'FITNESS','CRITICAL LOAD (XY) & M.S.'
       END IF

       WRITE(2,590)FIT(I),CRSH(I),CRSHMS(I)
       WRITE(2,*)
       WRITE(2,600)'PLY-BY-PLY STRAINS:','EP1','EP2','GAM12'

        DO 10, K = (NUM(I)/2)+1,NUM(I)
          WRITE(2,610)'PLY:',K,EP1(I,K),EP2(I,K),GAM12(I,K)
10      CONTINUE  

       WRITE(2,*)
       WRITE(2,620)'PLY STRAIN MARGIN OF SAFETY:',STRMS(I)

      IF(DIRCR(I) .EQ. 1) THEN
          WRITE(2,630)'PLY',PLYCR(I),'IS CRITICAL IN EPSILON-1 STRAIN'
          WRITE(2,640)'PLY STRAIN: ',EP1(I,PLYCR(I))
       ELSE IF(DIRCR(I) .EQ. 2) THEN
          WRITE(2,630)'PLY',PLYCR(I),'IS CRITICAL IN EPSILON-2 STRAIN'
          WRITE(2,640)'PLY STRAIN: ',EP2(I,PLYCR(I))
       ELSE IF(DIRCR(I) .EQ. 3) THEN
          WRITE(2,630)'PLY',PLYCR(I),'IS CRITICAL IN GAMMA-12 STRAIN'
          WRITE(2,640)'PLY STRAIN: ',GAM12(I,PLYCR(I))
       END IF

       IF(VIOL(I) .EQ. 1) THEN
         WRITE(2,*)
         WRITE(2,*)'CONSTRAINT(S) VIOLATED!!'
       END IF  

      WRITE(2,*) 
 
C+++++++++++++++++++++++++ FORMAT STATEMENTS +++++++++++++++++++++++

500   FORMAT(3X,A,1X,I4,2X,A,1X,A,/,72('='))
510   FORMAT(A,100(I3,1X,:))
520   FORMAT(16X,A)
530   FORMAT(A,4X,A,/,75('='))
540   FORMAT(F10.6,4X,F10.6,4X,F10.6,6X,F10.6,6X,F10.6/)
550   FORMAT(F10.6,4X,F10.6,4X,F10.6,13X,F10.6,/)
560   FORMAT(A,3X,A,3X,A,3X,A,3X,A,/,72('='))
570   FORMAT(F7.5,8X,F7.5,6X,F7.3,9X,F7.3,6X,F7.3)
580   FORMAT(1X,A,7X,A,/,72('='))
590   FORMAT(F14.6,8X,F8.3,3X,F8.3)

600   FORMAT(A,8X,A,10X,A,8X,A,/,70('='))
610   FORMAT(A,I2,15X,F11.8,3X,F11.8,3X,F11.8)
620   FORMAT(A,2X,F7.3,/)
630   FORMAT(A,1X,I3,1X,A)
640   FORMAT(A,1X,F16.8)



      RETURN
      END 
C----------------------------------------------------------------------

C--------------------------------------------------------------------



      subroutine rgg(nm,n,a,b,alfr,alfi,beta,matz,z,ierr)
c
      integer n,nm,ierr,matz
      double precision a(nm,100),b(nm,100),alfr(100),alfi(100),
     $   beta(100),z(nm,100)
      logical tf
c
c     this subroutine calls the recommended sequence of
c     subroutines from the eigensystem subroutine package (eispack)
c     to find the eigenvalues and eigenvectors (if desired)
c     for the real general generalized eigenproblem  ax = (lambda)bx.
c
c     on input
c
c        nm  must be set to the row dimension of the two-dimensional
c        array parameters as declared in the calling program
c        dimension statement.
c
c        n  is the order of the matrices  a  and  b.
c
c        a  contains a real general matrix.
c
c        b  contains a real general matrix.
c
c        matz  is an integer variable set equal to zero if
c        only eigenvalues are desired.  otherwise it is set to
c        any non-zero integer for both eigenvalues and eigenvectors.
c
c     on output
c
c        alfr  and  alfi  contain the real and imaginary parts,
c        respectively, of the numerators of the eigenvalues.
c
c        beta  contains the denominators of the eigenvalues,
c        which are thus given by the ratios  (alfr+i*alfi)/beta.
c        complex conjugate pairs of eigenvalues appear consecutively
c        with the eigenvalue having the positive imaginary part first.
c
c        z  contains the real and imaginary parts of the eigenvectors
c        if matz is not zero.  if the j-th eigenvalue is real, the
c        j-th column of  z  contains its eigenvector.  if the j-th
c        eigenvalue is complex with positive imaginary part, the
c        j-th and (j+1)-th columns of  z  contain the real and
c        imaginary parts of its eigenvector.  the conjugate of this
c        vector is the eigenvector for the conjugate eigenvalue.
c
c        ierr  is an integer output variable set equal to an error
c           completion code described in the documentation for qzit.
c           the normal completion code is zero.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      if (n .le. nm) go to 10
      ierr = 10 * n
      go to 50
c
   10 if (matz .ne. 0) go to 20
c     .......... find eigenvalues only ..........
      tf = .false.
      call  qzhes(nm,n,a,b,tf,z)
      call  qzit(nm,n,a,b,0.0d0,tf,z,ierr)
      call  qzval(nm,n,a,b,alfr,alfi,beta,tf,z)
      go to 50
c     .......... find both eigenvalues and eigenvectors ..........
   20 tf = .true.
      call  qzhes(nm,n,a,b,tf,z)
      call  qzit(nm,n,a,b,0.0d0,tf,z,ierr)
      call  qzval(nm,n,a,b,alfr,alfi,beta,tf,z)
      if (ierr .ne. 0) go to 50
      call  qzvec(nm,n,a,b,alfr,alfi,beta,z)
   50 return
      end
      double precision function epslon (x)
      double precision x
c
c     estimate unit roundoff in quantities of size x.
c
      double precision a,b,c,eps
c
c     this program should function properly on all systems
c     satisfying the following two assumptions,
c        1.  the base used in representing floating point
c            numbers is not a power of three.
c        2.  the quantity  a  in statement 10 is represented to 
c            the accuracy used in floating point variables
c            that are stored in memory.
c     the statement number 10 and the go to 10 are intended to
c     force optimizing compilers to generate code satisfying 
c     assumption 2.
c     under these assumptions, it should be true that,
c            a  is not exactly equal to four-thirds,
c            b  has a zero for its last bit or digit,
c            c  is not exactly equal to one,
c            eps  measures the separation of 1.0 from
c                 the next larger floating point number.
c     the developers of eispack would appreciate being informed
c     about any systems where these assumptions do not hold.
c
c     this version dated 4/6/83.
c
      a = 4.0d0/3.0d0
   10 b = a - 1.0d0
      c = b + b + b
      eps = dabs(c-1.0d0)
      if (eps .eq. 0.0d0) go to 10
      epslon = eps*dabs(x)
      return
      end
      subroutine qzhes(nm,n,a,b,matz,z)
c
      integer i,j,k,l,n,lb,l1,nm,nk1,nm1,nm2
      double precision a(nm,n),b(nm,n),z(nm,n)
      double precision r,s,t,u1,u2,v1,v2,rho
      logical matz
c
c     this subroutine is the first step of the qz algorithm
c     for solving generalized matrix eigenvalue problems,
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
c
c     this subroutine accepts a pair of real general matrices and
c     reduces one of them to upper hessenberg form and the other
c     to upper triangular form using orthogonal transformations.
c     it is usually followed by  qzit,  qzval  and, possibly,  qzvec.
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 matrices.
c
c        a contains a real general matrix.
c
c        b contains a real general matrix.
c
c        matz should be set to .true. if the right hand transformations
c          are to be accumulated for later use in computing
c          eigenvectors, and to .false. otherwise.
c
c     on output
c
c        a has been reduced to upper hessenberg form.  the elements
c          below the first subdiagonal have been set to zero.
c
c        b has been reduced to upper triangular form.  the elements
c          below the main diagonal have been set to zero.
c
c        z contains the product of the right hand transformations if
c          matz has been set to .true.  otherwise, z is not referenced.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
c     .......... initialize z ..........
      if (.not. matz) go to 10
c
      do 3 j = 1, n
c
         do 2 i = 1, n
            z(i,j) = 0.0d0
    2    continue
c
         z(j,j) = 1.0d0
    3 continue
c     .......... reduce b to upper triangular form ..........
   10 if (n .le. 1) go to 170
      nm1 = n - 1
c
      do 100 l = 1, nm1
         l1 = l + 1
         s = 0.0d0
c
         do 20 i = l1, n
            s = s + dabs(b(i,l))
   20    continue
c
         if (s .eq. 0.0d0) go to 100
         s = s + dabs(b(l,l))
         r = 0.0d0
c
         do 25 i = l, n
            b(i,l) = b(i,l) / s
            r = r + b(i,l)**2
   25    continue
c
         r = dsign(dsqrt(r),b(l,l))
         b(l,l) = b(l,l) + r
         rho = r * b(l,l)
c
         do 50 j = l1, n
            t = 0.0d0
c
            do 30 i = l, n
               t = t + b(i,l) * b(i,j)
   30       continue
c
            t = -t / rho
c
            do 40 i = l, n
               b(i,j) = b(i,j) + t * b(i,l)
   40       continue
c
   50    continue
c
         do 80 j = 1, n
            t = 0.0d0
c
            do 60 i = l, n
               t = t + b(i,l) * a(i,j)
   60       continue
c
            t = -t / rho
c
            do 70 i = l, n
               a(i,j) = a(i,j) + t * b(i,l)
   70       continue
c
   80    continue
c
         b(l,l) = -s * r
c
         do 90 i = l1, n
            b(i,l) = 0.0d0
   90    continue
c
  100 continue
c     .......... reduce a to upper hessenberg form, while
c                keeping b triangular ..........
      if (n .eq. 2) go to 170
      nm2 = n - 2
c
      do 160 k = 1, nm2
         nk1 = nm1 - k
c     .......... for l=n-1 step -1 until k+1 do -- ..........
         do 150 lb = 1, nk1
            l = n - lb
            l1 = l + 1
c     .......... zero a(l+1,k) ..........
            s = dabs(a(l,k)) + dabs(a(l1,k))
            if (s .eq. 0.0d0) go to 150
            u1 = a(l,k) / s
            u2 = a(l1,k) / s
            r = dsign(dsqrt(u1*u1+u2*u2),u1)
            v1 =  -(u1 + r) / r
            v2 = -u2 / r
            u2 = v2 / v1
c
            do 110 j = k, n
               t = a(l,j) + u2 * a(l1,j)
               a(l,j) = a(l,j) + t * v1
               a(l1,j) = a(l1,j) + t * v2
  110       continue
c
            a(l1,k) = 0.0d0
c
            do 120 j = l, n
               t = b(l,j) + u2 * b(l1,j)
               b(l,j) = b(l,j) + t * v1
               b(l1,j) = b(l1,j) + t * v2
  120       continue
c     .......... zero b(l+1,l) ..........
            s = dabs(b(l1,l1)) + dabs(b(l1,l))
            if (s .eq. 0.0d0) go to 150
            u1 = b(l1,l1) / s
            u2 = b(l1,l) / s
            r = dsign(dsqrt(u1*u1+u2*u2),u1)
            v1 =  -(u1 + r) / r
            v2 = -u2 / r
            u2 = v2 / v1
c
            do 130 i = 1, l1
               t = b(i,l1) + u2 * b(i,l)
               b(i,l1) = b(i,l1) + t * v1
               b(i,l) = b(i,l) + t * v2
  130       continue
c
            b(l1,l) = 0.0d0
c
            do 140 i = 1, n
               t = a(i,l1) + u2 * a(i,l)
               a(i,l1) = a(i,l1) + t * v1
               a(i,l) = a(i,l) + t * v2
  140       continue
c
            if (.not. matz) go to 150
c
            do 145 i = 1, n
               t = z(i,l1) + u2 * z(i,l)
               z(i,l1) = z(i,l1) + t * v1
               z(i,l) = z(i,l) + t * v2
  145       continue
c
  150    continue
c
  160 continue
c
  170 return
      end
      subroutine qzit(nm,n,a,b,eps1,matz,z,ierr)
c
      integer i,j,k,l,n,en,k1,k2,ld,ll,l1,na,nm,ish,itn,its,km1,lm1,
     x        enm2,ierr,lor1,enorn
      double precision a(nm,n),b(nm,n),z(nm,n)
      double precision r,s,t,a1,a2,a3,ep,sh,u1,u2,u3,v1,v2,v3,ani,a11,
     x       a12,a21,a22,a33,a34,a43,a44,bni,b11,b12,b22,b33,b34,
     x       b44,epsa,epsb,eps1,anorm,bnorm,epslon
      logical matz,notlas
c
c     this subroutine is the second step of the qz algorithm
c     for solving generalized matrix eigenvalue problems,
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart,
c     as modified in technical note nasa tn d-7305(1973) by ward.
c
c     this subroutine accepts a pair of real matrices, one of them
c     in upper hessenberg form and the other in upper triangular form.
c     it reduces the hessenberg matrix to quasi-triangular form using
c     orthogonal transformations while maintaining the triangular form
c     of the other matrix.  it is usually preceded by  qzhes  and
c     followed by  qzval  and, possibly,  qzvec.
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 matrices.
c
c        a contains a real upper hessenberg matrix.
c
c        b contains a real upper triangular matrix.
c
c        eps1 is a tolerance used to determine negligible elements.
c          eps1 = 0.0 (or negative) may be input, in which case an
c          element will be neglected only if it is less than roundoff
c          error times the norm of its matrix.  if the input eps1 is
c          positive, then an element will be considered negligible
c          if it is less than eps1 times the norm of its matrix.  a
c          positive value of eps1 may result in faster execution,
c          but less accurate results.
c
c        matz should be set to .true. if the right hand transformations
c          are to be accumulated for later use in computing
c          eigenvectors, and to .false. otherwise.
c
c        z contains, if matz has been set to .true., the
c          transformation matrix produced in the reduction
c          by  qzhes, if performed, or else the identity matrix.
c          if matz has been set to .false., z is not referenced.
c
c     on output
c
c        a has been reduced to quasi-triangular form.  the elements
c          below the first subdiagonal are still zero and no two
c          consecutive subdiagonal elements are nonzero.
c
c        b is still in upper triangular form, although its elements
c          have been altered.  the location b(n,1) is used to store
c          eps1 times the norm of b for later use by  qzval  and  qzvec.
c
c        z contains the product of the right hand transformations
c          (for both steps) if matz has been set to .true..
c
c        ierr is set to
c          zero       for normal return,
c          j          if the limit of 30*n iterations is exhausted
c                     while the j-th eigenvalue is being sought.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      ierr = 0
c     .......... compute epsa,epsb ..........
      anorm = 0.0d0
      bnorm = 0.0d0
c
      do 30 i = 1, n
         ani = 0.0d0
         if (i .ne. 1) ani = dabs(a(i,i-1))
         bni = 0.0d0
c
         do 20 j = i, n
            ani = ani + dabs(a(i,j))
            bni = bni + dabs(b(i,j))
   20    continue
c
         if (ani .gt. anorm) anorm = ani
         if (bni .gt. bnorm) bnorm = bni
   30 continue
c
      if (anorm .eq. 0.0d0) anorm = 1.0d0
      if (bnorm .eq. 0.0d0) bnorm = 1.0d0
      ep = eps1
      if (ep .gt. 0.0d0) go to 50
c     .......... use roundoff level if eps1 is zero ..........
      ep = epslon(1.0d0)
   50 epsa = ep * anorm
      epsb = ep * bnorm
c     .......... reduce a to quasi-triangular form, while
c                keeping b triangular ..........
      lor1 = 1
      enorn = n
      en = n
      itn = 30*n
c     .......... begin qz step ..........
   60 if (en .le. 2) go to 1001
      if (.not. matz) enorn = en
      its = 0
      na = en - 1
      enm2 = na - 1
   70 ish = 2
c     .......... check for convergence or reducibility.
c                for l=en step -1 until 1 do -- ..........
      do 80 ll = 1, en
         lm1 = en - ll
         l = lm1 + 1
         if (l .eq. 1) go to 95
         if (dabs(a(l,lm1)) .le. epsa) go to 90
   80 continue
c
   90 a(l,lm1) = 0.0d0
      if (l .lt. na) go to 95
c     .......... 1-by-1 or 2-by-2 block isolated ..........
      en = lm1
      go to 60
c     .......... check for small top of b ..........
   95 ld = l
  100 l1 = l + 1
      b11 = b(l,l)
      if (dabs(b11) .gt. epsb) go to 120
      b(l,l) = 0.0d0
      s = dabs(a(l,l)) + dabs(a(l1,l))
      u1 = a(l,l) / s
      u2 = a(l1,l) / s
      r = dsign(dsqrt(u1*u1+u2*u2),u1)
      v1 = -(u1 + r) / r
      v2 = -u2 / r
      u2 = v2 / v1
c
      do 110 j = l, enorn
         t = a(l,j) + u2 * a(l1,j)
         a(l,j) = a(l,j) + t * v1
         a(l1,j) = a(l1,j) + t * v2
         t = b(l,j) + u2 * b(l1,j)
         b(l,j) = b(l,j) + t * v1
         b(l1,j) = b(l1,j) + t * v2
  110 continue
c
      if (l .ne. 1) a(l,lm1) = -a(l,lm1)
      lm1 = l
      l = l1
      go to 90
  120 a11 = a(l,l) / b11
      a21 = a(l1,l) / b11
      if (ish .eq. 1) go to 140
c     .......... iteration strategy ..........
      if (itn .eq. 0) go to 1000
      if (its .eq. 10) go to 155
c     .......... determine type of shift ..........
      b22 = b(l1,l1)
      if (dabs(b22) .lt. epsb) b22 = epsb
      b33 = b(na,na)
      if (dabs(b33) .lt. epsb) b33 = epsb
      b44 = b(en,en)
      if (dabs(b44) .lt. epsb) b44 = epsb
      a33 = a(na,na) / b33
      a34 = a(na,en) / b44
      a43 = a(en,na) / b33
      a44 = a(en,en) / b44
      b34 = b(na,en) / b44
      t = 0.5d0 * (a43 * b34 - a33 - a44)
      r = t * t + a34 * a43 - a33 * a44
      if (r .lt. 0.0d0) go to 150
c     .......... determine single shift zeroth column of a ..........
      ish = 1
      r = dsqrt(r)
      sh = -t + r
      s = -t - r
      if (dabs(s-a44) .lt. dabs(sh-a44)) sh = s
c     .......... look for two consecutive small
c                sub-diagonal elements of a.
c                for l=en-2 step -1 until ld do -- ..........
      do 130 ll = ld, enm2
         l = enm2 + ld - ll
         if (l .eq. ld) go to 140
         lm1 = l - 1
         l1 = l + 1
         t = a(l,l)
         if (dabs(b(l,l)) .gt. epsb) t = t - sh * b(l,l)
         if (dabs(a(l,lm1)) .le. dabs(t/a(l1,l)) * epsa) go to 100
  130 continue
c
  140 a1 = a11 - sh
      a2 = a21
      if (l .ne. ld) a(l,lm1) = -a(l,lm1)
      go to 160
c     .......... determine double shift zeroth column of a ..........
  150 a12 = a(l,l1) / b22
      a22 = a(l1,l1) / b22
      b12 = b(l,l1) / b22
      a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11)
     x     / a21 + a12 - a11 * b12
      a2 = (a22 - a11) - a21 * b12 - (a33 - a11) - (a44 - a11)
     x     + a43 * b34
      a3 = a(l1+1,l1) / b22
      go to 160
c     .......... ad hoc shift ..........
  155 a1 = 0.0d0
      a2 = 1.0d0
      a3 = 1.1605d0
  160 its = its + 1
      itn = itn - 1
      if (.not. matz) lor1 = ld
c     .......... main loop ..........
      do 260 k = l, na
         notlas = k .ne. na .and. ish .eq. 2
         k1 = k + 1
         k2 = k + 2
         km1 = max0(k-1,l)
         ll = min0(en,k1+ish)
         if (notlas) go to 190
c     .......... zero a(k+1,k-1) ..........
         if (k .eq. l) go to 170
         a1 = a(k,km1)
         a2 = a(k1,km1)
  170    s = dabs(a1) + dabs(a2)
         if (s .eq. 0.0d0) go to 70
         u1 = a1 / s
         u2 = a2 / s
         r = dsign(dsqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 180 j = km1, enorn
            t = a(k,j) + u2 * a(k1,j)
            a(k,j) = a(k,j) + t * v1
            a(k1,j) = a(k1,j) + t * v2
            t = b(k,j) + u2 * b(k1,j)
            b(k,j) = b(k,j) + t * v1
            b(k1,j) = b(k1,j) + t * v2
  180    continue
c
         if (k .ne. l) a(k1,km1) = 0.0d0
         go to 240
c     .......... zero a(k+1,k-1) and a(k+2,k-1) ..........
  190    if (k .eq. l) go to 200
         a1 = a(k,km1)
         a2 = a(k1,km1)
         a3 = a(k2,km1)
  200    s = dabs(a1) + dabs(a2) + dabs(a3)
         if (s .eq. 0.0d0) go to 260
         u1 = a1 / s
         u2 = a2 / s
         u3 = a3 / s
         r = dsign(dsqrt(u1*u1+u2*u2+u3*u3),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         v3 = -u3 / r
         u2 = v2 / v1
         u3 = v3 / v1
c
         do 210 j = km1, enorn
            t = a(k,j) + u2 * a(k1,j) + u3 * a(k2,j)
            a(k,j) = a(k,j) + t * v1
            a(k1,j) = a(k1,j) + t * v2
            a(k2,j) = a(k2,j) + t * v3
            t = b(k,j) + u2 * b(k1,j) + u3 * b(k2,j)
            b(k,j) = b(k,j) + t * v1
            b(k1,j) = b(k1,j) + t * v2
            b(k2,j) = b(k2,j) + t * v3
  210    continue
c
         if (k .eq. l) go to 220
         a(k1,km1) = 0.0d0
         a(k2,km1) = 0.0d0
c     .......... zero b(k+2,k+1) and b(k+2,k) ..........
  220    s = dabs(b(k2,k2)) + dabs(b(k2,k1)) + dabs(b(k2,k))
         if (s .eq. 0.0d0) go to 240
         u1 = b(k2,k2) / s
         u2 = b(k2,k1) / s
         u3 = b(k2,k) / s
         r = dsign(dsqrt(u1*u1+u2*u2+u3*u3),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         v3 = -u3 / r
         u2 = v2 / v1
         u3 = v3 / v1
c
         do 230 i = lor1, ll
            t = a(i,k2) + u2 * a(i,k1) + u3 * a(i,k)
            a(i,k2) = a(i,k2) + t * v1
            a(i,k1) = a(i,k1) + t * v2
            a(i,k) = a(i,k) + t * v3
            t = b(i,k2) + u2 * b(i,k1) + u3 * b(i,k)
            b(i,k2) = b(i,k2) + t * v1
            b(i,k1) = b(i,k1) + t * v2
            b(i,k) = b(i,k) + t * v3
  230    continue
c
         b(k2,k) = 0.0d0
         b(k2,k1) = 0.0d0
         if (.not. matz) go to 240
c
         do 235 i = 1, n
            t = z(i,k2) + u2 * z(i,k1) + u3 * z(i,k)
            z(i,k2) = z(i,k2) + t * v1
            z(i,k1) = z(i,k1) + t * v2
            z(i,k) = z(i,k) + t * v3
  235    continue
c     .......... zero b(k+1,k) ..........
  240    s = dabs(b(k1,k1)) + dabs(b(k1,k))
         if (s .eq. 0.0d0) go to 260
         u1 = b(k1,k1) / s
         u2 = b(k1,k) / s
         r = dsign(dsqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 250 i = lor1, ll
            t = a(i,k1) + u2 * a(i,k)
            a(i,k1) = a(i,k1) + t * v1
            a(i,k) = a(i,k) + t * v2
            t = b(i,k1) + u2 * b(i,k)
            b(i,k1) = b(i,k1) + t * v1
            b(i,k) = b(i,k) + t * v2
  250    continue
c
         b(k1,k) = 0.0d0
         if (.not. matz) go to 260
c
         do 255 i = 1, n
            t = z(i,k1) + u2 * z(i,k)
            z(i,k1) = z(i,k1) + t * v1
            z(i,k) = z(i,k) + t * v2
  255    continue
c
  260 continue
c     .......... end qz step ..........
      go to 70
c     .......... set error -- all eigenvalues have not
c                converged after 30*n iterations ..........
 1000 ierr = en
c     .......... save epsb for use by qzval and qzvec ..........
 1001 if (n .gt. 1) b(n,1) = epsb
      return
      end
      subroutine qzval(nm,n,a,b,alfr,alfi,beta,matz,z)
c
      integer i,j,n,en,na,nm,nn,isw
      double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n)
      double precision c,d,e,r,s,t,an,a1,a2,bn,cq,cz,di,dr,ei,ti,tr,u1,
     x       u2,v1,v2,a1i,a11,a12,a2i,a21,a22,b11,b12,b22,sqi,sqr,
     x       ssi,ssr,szi,szr,a11i,a11r,a12i,a12r,a22i,a22r,epsb
      logical matz
c
c     this subroutine is the third step of the qz algorithm
c     for solving generalized matrix eigenvalue problems,
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
c
c     this subroutine accepts a pair of real matrices, one of them
c     in quasi-triangular form and the other in upper triangular form.
c     it reduces the quasi-triangular matrix further, so that any
c     remaining 2-by-2 blocks correspond to pairs of complex
c     eigenvalues, and returns quantities whose ratios give the
c     generalized eigenvalues.  it is usually preceded by  qzhes
c     and  qzit  and may be followed by  qzvec.
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 matrices.
c
c        a contains a real upper quasi-triangular matrix.
c
c        b contains a real upper triangular matrix.  in addition,
c          location b(n,1) contains the tolerance quantity (epsb)
c          computed and saved in  qzit.
c
c        matz should be set to .true. if the right hand transformations
c          are to be accumulated for later use in computing
c          eigenvectors, and to .false. otherwise.
c
c        z contains, if matz has been set to .true., the
c          transformation matrix produced in the reductions by qzhes
c          and qzit, if performed, or else the identity matrix.
c          if matz has been set to .false., z is not referenced.
c
c     on output
c
c        a has been reduced further to a quasi-triangular matrix
c          in which all nonzero subdiagonal elements correspond to
c          pairs of complex eigenvalues.
c
c        b is still in upper triangular form, although its elements
c          have been altered.  b(n,1) is unaltered.
c
c        alfr and alfi contain the real and imaginary parts of the
c          diagonal elements of the triangular matrix that would be
c          obtained if a were reduced completely to triangular form
c          by unitary transformations.  non-zero values of alfi occur
c          in pairs, the first member positive and the second negative.
c
c        beta contains the diagonal elements of the corresponding b,
c          normalized to be real and non-negative.  the generalized
c          eigenvalues are then the ratios ((alfr+i*alfi)/beta).
c
c        z contains the product of the right hand transformations
c          (for all three steps) if matz has been set to .true.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      epsb = b(n,1)
      isw = 1
c     .......... find eigenvalues of quasi-triangular matrices.
c                for en=n step -1 until 1 do -- ..........
      do 510 nn = 1, n
         en = n + 1 - nn
         na = en - 1
         if (isw .eq. 2) go to 505
         if (en .eq. 1) go to 410
         if (a(en,na) .ne. 0.0d0) go to 420
c     .......... 1-by-1 block, one real root ..........
  410    alfr(en) = a(en,en)
         if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en)
         beta(en) = dabs(b(en,en))
         alfi(en) = 0.0d0
         go to 510
c     .......... 2-by-2 block ..........
  420    if (dabs(b(na,na)) .le. epsb) go to 455
         if (dabs(b(en,en)) .gt. epsb) go to 430
         a1 = a(en,en)
         a2 = a(en,na)
         bn = 0.0d0
         go to 435
  430    an = dabs(a(na,na)) + dabs(a(na,en)) + dabs(a(en,na))
     x      + dabs(a(en,en))
         bn = dabs(b(na,na)) + dabs(b(na,en)) + dabs(b(en,en))
         a11 = a(na,na) / an
         a12 = a(na,en) / an
         a21 = a(en,na) / an
         a22 = a(en,en) / an
         b11 = b(na,na) / bn
         b12 = b(na,en) / bn
         b22 = b(en,en) / bn
         e = a11 / b11
         ei = a22 / b22
         s = a21 / (b11 * b22)
         t = (a22 - e * b22) / b22
         if (dabs(e) .le. dabs(ei)) go to 431
         e = ei
         t = (a11 - e * b11) / b11
  431    c = 0.5d0 * (t - s * b12)
         d = c * c + s * (a12 - e * b12)
         if (d .lt. 0.0d0) go to 480
c     .......... two real roots.
c                zero both a(en,na) and b(en,na) ..........
         e = e + (c + dsign(dsqrt(d),c))
         a11 = a11 - e * b11
         a12 = a12 - e * b12
         a22 = a22 - e * b22
         if (dabs(a11) + dabs(a12) .lt.
     x       dabs(a21) + dabs(a22)) go to 432
         a1 = a12
         a2 = a11
         go to 435
  432    a1 = a22
         a2 = a21
c     .......... choose and apply real z ..........
  435    s = dabs(a1) + dabs(a2)
         u1 = a1 / s
         u2 = a2 / s
         r = dsign(dsqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 440 i = 1, en
            t = a(i,en) + u2 * a(i,na)
            a(i,en) = a(i,en) + t * v1
            a(i,na) = a(i,na) + t * v2
            t = b(i,en) + u2 * b(i,na)
            b(i,en) = b(i,en) + t * v1
            b(i,na) = b(i,na) + t * v2
  440    continue
c
         if (.not. matz) go to 450
c
         do 445 i = 1, n
            t = z(i,en) + u2 * z(i,na)
            z(i,en) = z(i,en) + t * v1
            z(i,na) = z(i,na) + t * v2
  445    continue
c
  450    if (bn .eq. 0.0d0) go to 475
         if (an .lt. dabs(e) * bn) go to 455
         a1 = b(na,na)
         a2 = b(en,na)
         go to 460
  455    a1 = a(na,na)
         a2 = a(en,na)
c     .......... choose and apply real q ..........
  460    s = dabs(a1) + dabs(a2)
         if (s .eq. 0.0d0) go to 475
         u1 = a1 / s
         u2 = a2 / s
         r = dsign(dsqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 470 j = na, n
            t = a(na,j) + u2 * a(en,j)
            a(na,j) = a(na,j) + t * v1
            a(en,j) = a(en,j) + t * v2
            t = b(na,j) + u2 * b(en,j)
            b(na,j) = b(na,j) + t * v1
            b(en,j) = b(en,j) + t * v2
  470    continue
c
  475    a(en,na) = 0.0d0
         b(en,na) = 0.0d0
         alfr(na) = a(na,na)
         alfr(en) = a(en,en)
         if (b(na,na) .lt. 0.0d0) alfr(na) = -alfr(na)
         if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en)
         beta(na) = dabs(b(na,na))
         beta(en) = dabs(b(en,en))
         alfi(en) = 0.0d0
         alfi(na) = 0.0d0
         go to 505
c     .......... two complex roots ..........
  480    e = e + c
         ei = dsqrt(-d)
         a11r = a11 - e * b11
         a11i = ei * b11
         a12r = a12 - e * b12
         a12i = ei * b12
         a22r = a22 - e * b22
         a22i = ei * b22
         if (dabs(a11r) + dabs(a11i) + dabs(a12r) + dabs(a12i) .lt.
     x       dabs(a21) + dabs(a22r) + dabs(a22i)) go to 482
         a1 = a12r
         a1i = a12i
         a2 = -a11r
         a2i = -a11i
         go to 485
  482    a1 = a22r
         a1i = a22i
         a2 = -a21
         a2i = 0.0d0
c     .......... choose complex z ..........
  485    cz = dsqrt(a1*a1+a1i*a1i)
         if (cz .eq. 0.0d0) go to 487
         szr = (a1 * a2 + a1i * a2i) / cz
         szi = (a1 * a2i - a1i * a2) / cz
         r = dsqrt(cz*cz+szr*szr+szi*szi)
         cz = cz / r
         szr = szr / r
         szi = szi / r
         go to 490
  487    szr = 1.0d0
         szi = 0.0d0
  490    if (an .lt. (dabs(e) + ei) * bn) go to 492
         a1 = cz * b11 + szr * b12
         a1i = szi * b12
         a2 = szr * b22
         a2i = szi * b22
         go to 495
  492    a1 = cz * a11 + szr * a12
         a1i = szi * a12
         a2 = cz * a21 + szr * a22
         a2i = szi * a22
c     .......... choose complex q ..........
  495    cq = dsqrt(a1*a1+a1i*a1i)
         if (cq .eq. 0.0d0) go to 497
         sqr = (a1 * a2 + a1i * a2i) / cq
         sqi = (a1 * a2i - a1i * a2) / cq
         r = dsqrt(cq*cq+sqr*sqr+sqi*sqi)
         cq = cq / r
         sqr = sqr / r
         sqi = sqi / r
         go to 500
  497    sqr = 1.0d0
         sqi = 0.0d0
c     .......... compute diagonal elements that would result
c                if transformations were applied ..........
  500    ssr = sqr * szr + sqi * szi
         ssi = sqr * szi - sqi * szr
         i = 1
         tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21
     x      + ssr * a22
         ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22
         dr = cq * cz * b11 + cq * szr * b12 + ssr * b22
         di = cq * szi * b12 + ssi * b22
         go to 503
  502    i = 2
         tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21
     x      + cq * cz * a22
         ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21
         dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22
         di = -ssi * b11 - sqi * cz * b12
  503    t = ti * dr - tr * di
         j = na
         if (t .lt. 0.0d0) j = en
         r = dsqrt(dr*dr+di*di)
         beta(j) = bn * r
         alfr(j) = an * (tr * dr + ti * di) / r
         alfi(j) = an * t / r
         if (i .eq. 1) go to 502
  505    isw = 3 - isw
  510 continue
      b(n,1) = epsb
c
      return
      end
      subroutine qzvec(nm,n,a,b,alfr,alfi,beta,z)
c
      integer i,j,k,m,n,en,ii,jj,na,nm,nn,isw,enm2
      double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n)
      double precision d,q,r,s,t,w,x,y,di,dr,ra,rr,sa,ti,tr,t1,t2,w1,x1,
     x       zz,z1,alfm,almi,almr,betm,epsb
c
c     this subroutine is the optional fourth step of the qz algorithm
c     for solving generalized matrix eigenvalue problems,
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
c
c     this subroutine accepts a pair of real matrices, one of them in
c     quasi-triangular form (in which each 2-by-2 block corresponds to
c     a pair of complex eigenvalues) and the other in upper triangular
c     form.  it computes the eigenvectors of the triangular problem and
c     transforms the results back to the original coordinate system.
c     it is usually preceded by  qzhes,  qzit, and  qzval.
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 matrices.
c
c        a contains a real upper quasi-triangular matrix.
c
c        b contains a real upper triangular matrix.  in addition,
c          location b(n,1) contains the tolerance quantity (epsb)
c          computed and saved in  qzit.
c
c        alfr, alfi, and beta  are vectors with components whose
c          ratios ((alfr+i*alfi)/beta) are the generalized
c          eigenvalues.  they are usually obtained from  qzval.
c
c        z contains the transformation matrix produced in the
c          reductions by  qzhes,  qzit, and  qzval, if performed.
c          if the eigenvectors of the triangular problem are
c          desired, z must contain the identity matrix.
c
c     on output
c
c        a is unaltered.  its subdiagonal elements provide information
c           about the storage of the complex eigenvectors.
c
c        b has been destroyed.
c
c        alfr, alfi, and beta are unaltered.
c
c        z contains the real and imaginary parts of the eigenvectors.
c          if alfi(i) .eq. 0.0, the i-th eigenvalue is real and
c            the i-th column of z contains its eigenvector.
c          if alfi(i) .ne. 0.0, the i-th eigenvalue is complex.
c            if alfi(i) .gt. 0.0, the eigenvalue is the first of
c              a complex pair and the i-th and (i+1)-th columns
c              of z contain its eigenvector.
c            if alfi(i) .lt. 0.0, the eigenvalue is the second of
c              a complex pair and the (i-1)-th and i-th columns
c              of z contain the conjugate of its eigenvector.
c          each eigenvector is normalized so that the modulus
c          of its largest component is 1.0 .
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      epsb = b(n,1)
      isw = 1
c     .......... for en=n step -1 until 1 do -- ..........
      do 800 nn = 1, n
         en = n + 1 - nn
         na = en - 1
         if (isw .eq. 2) go to 795
         if (alfi(en) .ne. 0.0d0) go to 710
c     .......... real vector ..........
         m = en
         b(en,en) = 1.0d0
         if (na .eq. 0) go to 800
         alfm = alfr(m)
         betm = beta(m)
c     .......... for i=en-1 step -1 until 1 do -- ..........
         do 700 ii = 1, na
            i = en - ii
            w = betm * a(i,i) - alfm * b(i,i)
            r = 0.0d0
c
            do 610 j = m, en
  610       r = r + (betm * a(i,j) - alfm * b(i,j)) * b(j,en)
c
            if (i .eq. 1 .or. isw .eq. 2) go to 630
            if (betm * a(i,i-1) .eq. 0.0d0) go to 630
            zz = w
            s = r
            go to 690
  630       m = i
            if (isw .eq. 2) go to 640
c     .......... real 1-by-1 block ..........
            t = w
            if (w .eq. 0.0d0) t = epsb
            b(i,en) = -r / t
            go to 700
c     .......... real 2-by-2 block ..........
  640       x = betm * a(i,i+1) - alfm * b(i,i+1)
            y = betm * a(i+1,i)
            q = w * zz - x * y
            t = (x * s - zz * r) / q
            b(i,en) = t
            if (dabs(x) .le. dabs(zz)) go to 650
            b(i+1,en) = (-r - w * t) / x
            go to 690
  650       b(i+1,en) = (-s - y * t) / zz
  690       isw = 3 - isw
  700    continue
c     .......... end real vector ..........
         go to 800
c     .......... complex vector ..........
  710    m = na
         almr = alfr(m)
         almi = alfi(m)
         betm = beta(m)
c     .......... last vector component chosen imaginary so that
c                eigenvector matrix is triangular ..........
         y = betm * a(en,na)
         b(na,na) = -almi * b(en,en) / y
         b(na,en) = (almr * b(en,en) - betm * a(en,en)) / y
         b(en,na) = 0.0d0
         b(en,en) = 1.0d0
         enm2 = na - 1
         if (enm2 .eq. 0) go to 795
c     .......... for i=en-2 step -1 until 1 do -- ..........
         do 790 ii = 1, enm2
            i = na - ii
            w = betm * a(i,i) - almr * b(i,i)
            w1 = -almi * b(i,i)
            ra = 0.0d0
            sa = 0.0d0
c
            do 760 j = m, en
               x = betm * a(i,j) - almr * b(i,j)
               x1 = -almi * b(i,j)
               ra = ra + x * b(j,na) - x1 * b(j,en)
               sa = sa + x * b(j,en) + x1 * b(j,na)
  760       continue
c
            if (i .eq. 1 .or. isw .eq. 2) go to 770
            if (betm * a(i,i-1) .eq. 0.0d0) go to 770
            zz = w
            z1 = w1
            r = ra
            s = sa
            isw = 2
            go to 790
  770       m = i
            if (isw .eq. 2) go to 780
c     .......... complex 1-by-1 block ..........
            tr = -ra
            ti = -sa
  773       dr = w
            di = w1
c     .......... complex divide (t1,t2) = (tr,ti) / (dr,di) ..........
  775       if (dabs(di) .gt. dabs(dr)) go to 777
            rr = di / dr
            d = dr + di * rr
            t1 = (tr + ti * rr) / d
            t2 = (ti - tr * rr) / d
            go to (787,782), isw
  777       rr = dr / di
            d = dr * rr + di
            t1 = (tr * rr + ti) / d
            t2 = (ti * rr - tr) / d
            go to (787,782), isw
c     .......... complex 2-by-2 block ..........
  780       x = betm * a(i,i+1) - almr * b(i,i+1)
            x1 = -almi * b(i,i+1)
            y = betm * a(i+1,i)
            tr = y * ra - w * r + w1 * s
            ti = y * sa - w * s - w1 * r
            dr = w * zz - w1 * z1 - x * y
            di = w * z1 + w1 * zz - x1 * y
            if (dr .eq. 0.0d0 .and. di .eq. 0.0d0) dr = epsb
            go to 775
  782       b(i+1,na) = t1
            b(i+1,en) = t2
            isw = 1
            if (dabs(y) .gt. dabs(w) + dabs(w1)) go to 785
            tr = -ra - x * b(i+1,na) + x1 * b(i+1,en)
            ti = -sa - x * b(i+1,en) - x1 * b(i+1,na)
            go to 773
  785       t1 = (-r - zz * b(i+1,na) + z1 * b(i+1,en)) / y
            t2 = (-s - zz * b(i+1,en) - z1 * b(i+1,na)) / y
  787       b(i,na) = t1
            b(i,en) = t2
  790    continue
c     .......... end complex vector ..........
  795    isw = 3 - isw
  800 continue
c     .......... end back substitution.
c                transform to original coordinate system.
c                for j=n step -1 until 1 do -- ..........
      do 880 jj = 1, n
         j = n + 1 - jj
c
         do 880 i = 1, n
            zz = 0.0d0
c
            do 860 k = 1, j
  860       zz = zz + z(i,k) * b(k,j)
c
            z(i,j) = zz
  880 continue
c     .......... normalize so that modulus of largest
c                component of each vector is 1.
c                (isw is 1 initially from before) ..........
      do 950 j = 1, n
         d = 0.0d0
         if (isw .eq. 2) go to 920
         if (alfi(j) .ne. 0.0d0) go to 945
c
         do 890 i = 1, n
            if (dabs(z(i,j)) .gt. d) d = dabs(z(i,j))
  890    continue
c
         do 900 i = 1, n
  900    z(i,j) = z(i,j) / d
c
         go to 950
c
  920    do 930 i = 1, n
            r = dabs(z(i,j-1)) + dabs(z(i,j))
            if (r .ne. 0.0d0) r = r * dsqrt((z(i,j-1)/r)**2
     x                                     +(z(i,j)/r)**2)
            if (r .gt. d) d = r
  930    continue
c
         do 940 i = 1, n
            z(i,j-1) = z(i,j-1) / d
            z(i,j) = z(i,j) / d
  940    continue
c
  945    isw = 3 - isw
  950 continue
c
      return
      end