Google
 

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