SUBROUTINE CDINVS (
*
*  inputs
*
     :                  MTXIN, N, NMAX,
*
*  outputs
*
     :                  MTXOUT, STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  Invert a symmetric matrix of double precision.
*
*  Description:
*  ------------
*  Maximum order of the input matrix is 30.
*
*  FORTRAN name: CDINVS.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       None
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       None
*  SDAS:
*       UMSPUT
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       02-15-87     J.-C. HSU        design and coding
*     2       08-20-87     J.-C. HSU        F77 standard
*-------------------------------------------------------------------------------
*
*== input:
*                                       --order of input matrix,
        INTEGER                 N,
*                                       --(maximum) order of the matrix array
*                                       --as declared in the calling routine
     :                          NMAX
*                                       --input matrix
        DOUBLE PRECISION        MTXIN(NMAX, NMAX)
*
*== output:
*                                       --output matrix (inverse of MTXIN)
        DOUBLE PRECISION        MTXOUT(NMAX, NMAX)
*                                       --error status
        INTEGER                 STATUS
*
*== local:
*                                       --loop indices
        INTEGER                 J, K, L, M, KM1, KP1,
*                                       --error status
     :                          STATOK
*                                       --error message
        CHARACTER*130           CONTXT, MESS
*                                       --largest (absolute) diagonal element
        DOUBLE PRECISION        BIG,
*                                       --working arrays
     :                          P(30), Q(30), R(30),
*                                       --absolute value of diagonal elements
     :                          AB
*=========================begin hsp.inc=========================================
*                                       --status return code
        INTEGER         OK, ERRNUM(20)
        INTEGER         DEST, PRIO
        DATA OK /0/
        DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
     :               711, 712, 713, 714, 715, 716, 717, 718, 719, 720/
*                                       --message destination and priority
        DATA DEST, PRIO /1, 0/
*=========================end hsp.inc===========================================
*------------------------------------------------------------------------------
*
*  check if the order is out of range
*
        IF (N. GT. 30) THEN
            STATUS = ERRNUM(1)
            CONTXT = 'matrix size is larger than 30 by 30'
            GO TO 999
        ELSE IF (N .GT. NMAX) THEN
            STATUS = ERRNUM(2)
            CONTXT = 'matrix size is larger than what is specified in'
     :                  // ' the calling routine'
            GO TO 999
        END IF
*
*  initilization
*
        DO 20 M = 1, N
            R(M) = 1.D0
            DO 10 K = 1, N
                MTXOUT(K, M) = MTXIN(K, M)
   10       CONTINUE
   20   CONTINUE
*
*  look for the largest (absolute) diagonal element
*
        DO 80 M = 1, N
            BIG = 0.
            DO 30 L = 1, N
                AB = ABS (MTXOUT(L, L))
                IF (AB .GT. BIG .AND. R(L) .NE. 0.) THEN
                    BIG = AB
                    K = L
                END IF
   30       CONTINUE
*
*  check if the largest diagonal element is zero
*
            IF (BIG .EQ. 0.) THEN
                STATUS = ERRNUM(3)
                CONTXT = 'matrix to be inverted is ill-defined'
                GO TO 999
            END IF
*
            R(K) = 0
            Q(K) = 1.D0 / MTXOUT(K, K)
            P(K) = 1.D0
            MTXOUT(K, K) = 0.
*
            KM1 = K - 1
            IF (KM1 .NE. 0.) THEN
                DO 40 L = 1, KM1
                    P(L) = MTXOUT(L, K)
                    IF (R(L) .NE. 0.) THEN
                        Q(L) = -MTXOUT(L, K) * Q(K)
                    ELSE
                        Q(L) = MTXOUT(L, K) * Q(K)
                    END IF
*
                    MTXOUT(L, K) = 0.
   40           CONTINUE
            END IF
*
            KP1 = K + 1
            IF (KP1 .LE. N) THEN
                DO 50 L = KP1, N
                    IF (R(L) .NE. 0) THEN
                        P(L) = MTXOUT(K, L)
                    ELSE
                        P(L) = -MTXOUT(K, L)
                    END IF
*
                    Q(L) = -MTXOUT(K, L) * Q(K)
                    MTXOUT(K, L) = 0.
   50           CONTINUE
            END IF
*
            DO 70 L = 1, N
                DO 60 K = L, N
                    MTXOUT(L, K) = MTXOUT(L, K) + P(L) * Q(K)
   60           CONTINUE
   70       CONTINUE
   80   CONTINUE
*
*  equate the symmetric terms
*
        M = N + 1
        L = N
        DO 100 K = 2, N
            M = M -1
            L = L -1
            DO 90 J = 1, L
                MTXOUT(M, J) = MTXOUT(J, M)
   90       CONTINUE
  100   CONTINUE
*
        STATUS = OK
        GO TO 1000
*
*  write error message
*
  999   MESS = 'CDINVS: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,    STATOK)
*
 1000   RETURN
        END