      SUBROUTINE  MtxCnv( A, B, N, N1, W, IER )                          EAB00100
C***********************************************************************EAB00200
C*                                                                     *EAB00300
C*  PURPOSE                                                            *EAB00400
C*      SOLVES GENERALIZED EIGENVALUE PROBLEM                          *EAB00500
C*          (A) X = LAMBDA (B) X                                       *EAB00600
C*      FOR REAL SYMMETRIC MATRICES (A) AND (B),                       *EAB00700
C*      THE LATTER BEING POSITIVE DEFINITE.                            *EAB00800
C*      DOUBLE PRECISION.                                              *EAB00900
C*                                                                     *EAB01000
C*  USAGE                                                              *EAB01100
C*      CALL  MtxCnv( A, B, N, N1, W, IER )                            *EAB01200
C*                                                                     *EAB01300
C*  DESCRIPTION OF PARAMETERS                                          *EAB01400
C*    (INPUT)                                                          *EAB01500
C*    A(N1,N)  - 2-DIM ARRAY CONTAINING THE REAL SYMMETRIC MATRIX.     *EAB01600
C*    B(N1,N)  - 2-DIM ARRAY CONTAINING THE REAL SYMMETRIC             *EAB01700
C*               POSITIVE-DEFINITE MATRIX.                             *EAB01800
C*                                                                     *EAB01900
C*      N      - ORDER OF THE MATRIX.                                  *EAB02000
C*      N1     - INTEGER WHICH SPECIFIES THE SIZE OF THE ARRAYS        *EAB02100
C*               A, B, AND V.                                          *EAB02200
C*                                                                     *EAB02600
C*    (WORK AREA)                                                      *EAB03200
C*     W(N,7)  - WORK AREA USED IN COMPUTATION OF EIGENVALUES AND      *EAB03300
C*               EIGENVECTORS                                          *EAB03400
C*     IER     - RETURN CODE                                           *EAB03410
C*               0: NORMAL RETURN                                      *EAB03420
C*               NONZERO: ERROR RETURN                                 *EAB03430
C*                                                                     *EAB03500
C*  REMARK                                                             *EAB03600
C*      1.  THE ARRAYS  A  AND  B  WILL NOT BE CONSERVED.              *EAB03700
C*      2.  IT SUFFICES TO FILL A(I,J) AND B(I,J)  FOR  I .GE. J.      *EAB03800
C*                                                                     *EAB03900
C*  SUBROUTINE REQUIRED                                                *EAB04000
C*      EIGRS                                                          *EAB04100
C*                                                                     *EAB04200
C*  COPYRIGHT                                                          *EAB04300
C*      Y. ONODERA,    AUGUST, 1989.                                   *EAB04400
C*                                                                     *EAB04500
C***********************************************************************EAB04600
      IMPLICIT REAL*8  (A-H,O-Z)                                        EAB04700
      IMPLICIT INTEGER (I-N)                                            EAB04800
      real*8 A(N1,N), B(N1,N), W(N,7)                                     EAB04900

      IF ( N1 .LE. 0 )                    GOTO 8                        EAB05100
      IF ( N  .LE. 0  .OR. N .GT. N1 )    GOTO 8                        EAB05200
        GOTO 10                                                         EAB05500
    8 WRITE(6,9) N, N1                                                  EAB05600
    9 FORMAT(' (SUBR EIGAB) INVALID ARGUMENT.  N, N1 =',                EAB05700
     & 4I5)                                                             EAB05800
      IER = 2                                                           EAB05850
      RETURN                                                            EAB05900
   10 CONTINUE                                                          EAB06000
      IF ( N .NE. 1 ) GOTO 20                                           EAB06100
      T = B(1,1)                                                        EAB06200
      IF ( T .LE. 0 ) GOTO 15                                           EAB06300
c      E(1) = A(1,1) / T                                                 EAB06400
c      IF ( NV .EQ. 1 ) V(1,1) = 1                                       EAB06500
      RETURN                                                            EAB06600
   15 WRITE(6,16)                                                       EAB06700
   16 FORMAT( ' (SUBR EIGAB) MATRIX B IS NOT POSITIVE DEFINITE.')       EAB06800
      IER = 1                                                           EAB06850
      RETURN                                                            EAB06900
C ***************************************************************       EAB07000
C             MODIFIED CHOLESKY FACTORIZATION STARTS:                   EAB07100
C               (B) = (L TRANSPOSED) (D) (L),                           EAB07200
C             WHERE (D) IS DIAGONAL MATRIX,                             EAB07300
C                   (L) IS LOWER TRIANGULAR MATRIX.                     EAB07400
C ***************************************************************       EAB07500
   20 DO 24 K=N,1,-1                                                    EAB07600
       SUM = 0                                                          EAB07700
       DO 21 I=K+1,N                                                    EAB07800
        SUM = SUM + W(I,1) * B(I,K)**2                                  EAB07900
   21  CONTINUE                                                         EAB08000
       D = B(K,K) - SUM                                                 EAB08100
       IF ( D .LE. 0 ) GOTO 15                                          EAB08200
       W(K,1) = D                                                       EAB08300
       D = 1 / D                                                        EAB08400
       B(K,K) = SQRT(D)                                                 EAB08500
       DO 23 J=1,K-1                                                    EAB08600
        SUM = 0                                                         EAB08700
        DO 22 I=K+1,N                                                   EAB08800
         SUM = SUM + W(I,1) * B(I,K) * B(I,J)                           EAB08900
   22   CONTINUE                                                        EAB09000
        B(K,J) = ( B(K,J) - SUM ) * D                                   EAB09100
   23  CONTINUE                                                         EAB09200
   24 CONTINUE                                                          EAB09300
C ***************************************************************       EAB09400
C             TRANSFORM THE MATRIX                                      EAB09500
C             (A NEW) = (L TRANSPOSED)**(-1) (A) (L)**(-1)              EAB09600
C ***************************************************************       EAB09700
      DO 36 K=N-1,1,-1                                                  EAB09800
       DO 32 J=1,K-1                                                    EAB09900
        SUM = 0                                                         EAB10000
        DO 31 I=K+1,N                                                   EAB10100
         SUM = SUM + B(I,K) * A(I,J)                                    EAB10200
   31   CONTINUE                                                        EAB10300
        A(K,J) = A(K,J) - SUM                                           EAB10400
   32  CONTINUE                                                         EAB10500
       SUMK = 0                                                         EAB10600
       DO 35 J=K+1,N                                                    EAB10700
        SUM = 0                                                         EAB10800
        DO 33 I=K+1,J-1                                                 EAB10900
         SUM = SUM + B(I,K) * A(J,I)                                    EAB11000
   33   CONTINUE                                                        EAB11100
        DO 34 I=J,N                                                     EAB11200
         SUM = SUM + B(I,K) * A(I,J)                                    EAB11300
   34   CONTINUE                                                        EAB11400
        T = A(J,K) - SUM                                                EAB11500
        SUMK = SUMK + B(J,K) * ( A(J,K) + T )                           EAB11600
        A(J,K) = T                                                      EAB11700
   35  CONTINUE                                                         EAB11800
       A(K,K) = A(K,K) - SUMK                                           EAB11900
   36 CONTINUE                                                          EAB12000
C ***************************************************************       EAB12100
C             PRE-MULTIPLY AND POST-MULTIPLY (A) WITH (D)**(-1/2)       EAB12200
C ***************************************************************       EAB12300
      DO 42 I=1,N                                                       EAB12400
       A(I,I) = A(I,I) / W(I,1)                                         EAB12500
       T=B(I,I)                                                         EAB12600
       DO 41 J=1,I-1                                                    EAB12700
        A(I,J) = A(I,J) * B(J,J) * T                                    EAB12800
        a(j,i) = a(i,j)
   41  CONTINUE                                                         EAB12900
   42 CONTINUE                                                          EAB13000
      RETURN                                                            EAB15300
      END                                                               EAB15400
