home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / mateigen / hsbg.for < prev    next >
Text File  |  1985-11-29  |  3KB  |  125 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE HSBG
  5. C
  6. C        PURPOSE
  7. C           TO REDUCE A REAL MATRIX INTO UPPER ALMOST TRIANGULAR FORM
  8. C
  9. C        USAGE
  10. C           CALL HSBG(N,A,IA)
  11. C
  12. C        DESCRIPTION OF THE PARAMETERS
  13. C           N      ORDER OF THE MATRIX
  14. C           A      THE INPUT MATRIX, N BY N
  15. C           IA     SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY
  16. C                  A IN THE CALLING PROGRAM WHEN THE MATRIX IS IN
  17. C                  DOUBLE SUBSCRIPTED DATA STORAGE MODE.  IA=N WHEN
  18. C                  THE MATRIX IS IN SSP VECTOR STORAGE MODE.
  19. C
  20. C        REMARKS
  21. C           THE HESSENBERG FORM REPLACES THE ORIGINAL MATRIX IN THE
  22. C           ARRAY A.
  23. C
  24. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  25. C           NONE
  26. C
  27. C        METHOD
  28. C           SIMILARITY TRANSFORMATIONS USING ELEMENTARY ELIMINATION
  29. C           MATRICES, WITH PARTIAL PIVOTING.
  30. C
  31. C        REFERENCES
  32. C           J.H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM -
  33. C           CLARENDON PRESS, OXFORD, 1965.
  34. C
  35. C     ..................................................................
  36. C
  37.       SUBROUTINE HSBG(N,A,IA)
  38.       DIMENSION A(1)
  39.       DOUBLE PRECISION S
  40.       L=N
  41.       NIA=L*IA
  42.       LIA=NIA-IA
  43. C
  44. C        L IS THE ROW INDEX OF THE ELIMINATION
  45. C
  46.    20 IF(L-3) 360,40,40
  47.    40 LIA=LIA-IA
  48.       L1=L-1
  49.       L2=L1-1
  50. C
  51. C        SEARCH FOR THE PIVOTAL ELEMENT IN THE LTH ROW
  52. C
  53.       ISUB=LIA+L
  54.       IPIV=ISUB-IA
  55.       PIV=ABS(A(IPIV))
  56.       IF(L-3) 90,90,50
  57.    50 M=IPIV-IA
  58.       DO 80 I=L,M,IA
  59.       T=ABS(A(I))
  60.       IF(T-PIV) 80,80,60
  61.    60 IPIV=I
  62.       PIV=T
  63.    80 CONTINUE
  64.    90 IF(PIV) 100,320,100
  65.   100 IF(PIV-ABS(A(ISUB))) 180,180,120
  66. C
  67. C        INTERCHANGE THE COLUMNS
  68. C
  69.   120 M=IPIV-L
  70.       DO 140 I=1,L
  71.       J=M+I
  72.       T=A(J)
  73.       K=LIA+I
  74.       A(J)=A(K)
  75.   140 A(K)=T
  76. C
  77. C        INTERCHANGE THE ROWS
  78. C
  79.       M=L2-M/IA
  80.       DO 160 I=L1,NIA,IA
  81.       T=A(I)
  82.       J=I-M
  83.       A(I)=A(J)
  84.   160 A(J)=T
  85. C
  86. C        TERMS OF THE ELEMENTARY TRANSFORMATION
  87. C
  88.   180 DO 200 I=L,LIA,IA
  89.   200 A(I)=A(I)/A(ISUB)
  90. C
  91. C        RIGHT TRANSFORMATION
  92. C
  93.       J=-IA
  94.       DO 240 I=1,L2
  95.       J=J+IA
  96.       LJ=L+J
  97.       DO 220 K=1,L1
  98.       KJ=K+J
  99.       KL=K+LIA
  100.   220 A(KJ)=A(KJ)-A(LJ)*A(KL)
  101.   240 CONTINUE
  102. C
  103. C        LEFT TRANSFORMATION
  104. C
  105.       K=-IA
  106.       DO 300 I=1,N
  107.       K=K+IA
  108.       LK=K+L1
  109.       S=A(LK)
  110.       LJ=L-IA
  111.       DO 280 J=1,L2
  112.       JK=K+J
  113.       LJ=LJ+IA
  114.   280 S=S+A(LJ)*A(JK)*1.0D0
  115.   300 A(LK)=S
  116. C
  117. C        SET THE LOWER PART OF THE MATRIX TO ZERO
  118. C
  119.       DO 310 I=L,LIA,IA
  120.   310 A(I)=0.0
  121.   320 L=L1
  122.       GO TO 20
  123.   360 RETURN
  124.       END
  125.