Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dsinv.ssp
There are 2 other files named dsinv.ssp in the archive. Click here to see a list.
C                                                                       DSIN  10
C     ..................................................................DSIN  20
C                                                                       DSIN  30
C        SUBROUTINE DSINV                                               DSIN  40
C                                                                       DSIN  50
C        PURPOSE                                                        DSIN  60
C           INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX           DSIN  70
C                                                                       DSIN  80
C        USAGE                                                          DSIN  90
C           CALL DSINV(A,N,EPS,IER)                                     DSIN 100
C                                                                       DSIN 110
C        DESCRIPTION OF PARAMETERS                                      DSIN 120
C           A      - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN    DSIN 130
C                    SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT     DSIN 140
C                    MATRIX.                                            DSIN 150
C                    ON RETURN A CONTAINS THE RESULTANT UPPER           DSIN 160
C                    TRIANGULAR MATRIX IN DOUBLE PRECISION.             DSIN 170
C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.      DSIN 180
C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED      DSIN 190
C                    AS RELATIVE TOLERANCE FOR TEST ON LOSS OF          DSIN 200
C                    SIGNIFICANCE.                                      DSIN 210
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS         DSIN 220
C                    IER=0  - NO ERROR                                  DSIN 230
C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-  DSIN 240
C                             TER N OR BECAUSE SOME RADICAND IS NON-    DSIN 250
C                             POSITIVE (MATRIX A IS NOT POSITIVE        DSIN 260
C                             DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-  DSIN 270
C                             FICANCE)                                  DSIN 280
C                    IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFI-  DSIN 290
C                             CANCE. THE RADICAND FORMED AT FACTORIZA-  DSIN 300
C                             TION STEP K+1 WAS STILL POSITIVE BUT NO   DSIN 310
C                             LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).  DSIN 320
C                                                                       DSIN 330
C        REMARKS                                                        DSIN 340
C           THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE  DSIN 350
C           STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.DSIN 360
C           IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-  DSIN 370
C           LAR MATRIX IS STORED COLUMNWISE TOO.                        DSIN 380
C           THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL  DSIN 390
C           CALCULATED RADICANDS ARE POSITIVE.                          DSIN 400
C                                                                       DSIN 410
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DSIN 420
C           DMFSD                                                       DSIN 430
C                                                                       DSIN 440
C        METHOD                                                         DSIN 450
C           SOLUTION IS DONE USING FACTORIZATION BY SUBROUTINE DMFSD.   DSIN 460
C                                                                       DSIN 470
C     ..................................................................DSIN 480
C                                                                       DSIN 490
      SUBROUTINE DSINV(A,N,EPS,IER)                                     DSIN 500
C                                                                       DSIN 510
C                                                                       DSIN 520
      DIMENSION A(1)                                                    DSIN 530
      DOUBLE PRECISION A,DIN,WORK                                       DSIN 540
C                                                                       DSIN 550
C        FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMFSD            DSIN 560
C        A = TRANSPOSE(T) * T                                           DSIN 570
      CALL DMFSD(A,N,EPS,IER)                                           DSIN 580
      IF(IER) 9,1,1                                                     DSIN 590
C                                                                       DSIN 600
C        INVERT UPPER TRIANGULAR MATRIX T                               DSIN 610
C        PREPARE INVERSION-LOOP                                         DSIN 620
    1 IPIV=N*(N+1)/2                                                    DSIN 630
      IND=IPIV                                                          DSIN 640
C                                                                       DSIN 650
C        INITIALIZE INVERSION-LOOP                                      DSIN 660
      DO 6 I=1,N                                                        DSIN 670
      DIN=1.D0/A(IPIV)                                                  DSIN 680
      A(IPIV)=DIN                                                       DSIN 690
      MIN=N                                                             DSIN 700
      KEND=I-1                                                          DSIN 710
      LANF=N-KEND                                                       DSIN 720
      IF(KEND) 5,5,2                                                    DSIN 730
    2 J=IND                                                             DSIN 740
C                                                                       DSIN 750
C        INITIALIZE ROW-LOOP                                            DSIN 760
      DO 4 K=1,KEND                                                     DSIN 770
      WORK=0.D0                                                         DSIN 780
      MIN=MIN-1                                                         DSIN 790
      LHOR=IPIV                                                         DSIN 800
      LVER=J                                                            DSIN 810
C                                                                       DSIN 820
C        START INNER LOOP                                               DSIN 830
      DO 3 L=LANF,MIN                                                   DSIN 840
      LVER=LVER+1                                                       DSIN 850
      LHOR=LHOR+L                                                       DSIN 860
    3 WORK=WORK+A(LVER)*A(LHOR)                                         DSIN 870
C        END OF INNER LOOP                                              DSIN 880
C                                                                       DSIN 890
      A(J)=-WORK*DIN                                                    DSIN 900
    4 J=J-MIN                                                           DSIN 910
C        END OF ROW-LOOP                                                DSIN 920
C                                                                       DSIN 930
    5 IPIV=IPIV-MIN                                                     DSIN 940
    6 IND=IND-1                                                         DSIN 950
C        END OF INVERSION-LOOP                                          DSIN 960
C                                                                       DSIN 970
C        CALCULATE INVERSE(A) BY MEANS OF INVERSE(T)                    DSIN 980
C        INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T))                DSIN 990
C        INITIALIZE MULTIPLICATION-LOOP                                 DSIN1000
      DO 8 I=1,N                                                        DSIN1010
      IPIV=IPIV+I                                                       DSIN1020
      J=IPIV                                                            DSIN1030
C                                                                       DSIN1040
C        INITIALIZE ROW-LOOP                                            DSIN1050
      DO 8 K=I,N                                                        DSIN1060
      WORK=0.D0                                                         DSIN1070
      LHOR=J                                                            DSIN1080
C                                                                       DSIN1090
C        START INNER LOOP                                               DSIN1100
      DO 7 L=K,N                                                        DSIN1110
      LVER=LHOR+K-I                                                     DSIN1120
      WORK=WORK+A(LHOR)*A(LVER)                                         DSIN1130
    7 LHOR=LHOR+L                                                       DSIN1140
C        END OF INNER LOOP                                              DSIN1150
C                                                                       DSIN1160
      A(J)=WORK                                                         DSIN1170
    8 J=J+K                                                             DSIN1180
C        END OF ROW- AND MULTIPLICATION-LOOP                            DSIN1190
C                                                                       DSIN1200
    9 RETURN                                                            DSIN1210
      END                                                               DSIN1220