ctrsl.f 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150
  1. *DECK CTRSL
  2. SUBROUTINE CTRSL (T, LDT, N, B, JOB, INFO)
  3. C***BEGIN PROLOGUE CTRSL
  4. C***PURPOSE Solve a system of the form T*X=B or CTRANS(T)*X=B, where
  5. C T is a triangular matrix. Here CTRANS(T) is the conjugate
  6. C transpose.
  7. C***LIBRARY SLATEC (LINPACK)
  8. C***CATEGORY D2C3
  9. C***TYPE COMPLEX (STRSL-S, DTRSL-D, CTRSL-C)
  10. C***KEYWORDS LINEAR ALGEBRA, LINPACK, TRIANGULAR LINEAR SYSTEM,
  11. C TRIANGULAR MATRIX
  12. C***AUTHOR Stewart, G. W., (U. of Maryland)
  13. C***DESCRIPTION
  14. C
  15. C CTRSL solves systems of the form
  16. C
  17. C T * X = B
  18. C or
  19. C CTRANS(T) * X = B
  20. C
  21. C where T is a triangular matrix of order N. Here CTRANS(T)
  22. C denotes the conjugate transpose of the matrix T.
  23. C
  24. C On Entry
  25. C
  26. C T COMPLEX(LDT,N)
  27. C T contains the matrix of the system. The zero
  28. C elements of the matrix are not referenced, and
  29. C the corresponding elements of the array can be
  30. C used to store other information.
  31. C
  32. C LDT INTEGER
  33. C LDT is the leading dimension of the array T.
  34. C
  35. C N INTEGER
  36. C N is the order of the system.
  37. C
  38. C B COMPLEX(N).
  39. C B contains the right hand side of the system.
  40. C
  41. C JOB INTEGER
  42. C JOB specifies what kind of system is to be solved.
  43. C If JOB is
  44. C
  45. C 00 solve T*X = B, T lower triangular,
  46. C 01 solve T*X = B, T upper triangular,
  47. C 10 solve CTRANS(T)*X = B, T lower triangular,
  48. C 11 solve CTRANS(T)*X = B, T upper triangular.
  49. C
  50. C On Return
  51. C
  52. C B B contains the solution, if INFO .EQ. 0.
  53. C Otherwise B is unaltered.
  54. C
  55. C INFO INTEGER
  56. C INFO contains zero if the system is nonsingular.
  57. C Otherwise INFO contains the index of
  58. C the first zero diagonal element of T.
  59. C
  60. C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
  61. C Stewart, LINPACK Users' Guide, SIAM, 1979.
  62. C***ROUTINES CALLED CAXPY, CDOTC
  63. C***REVISION HISTORY (YYMMDD)
  64. C 780814 DATE WRITTEN
  65. C 890831 Modified array declarations. (WRB)
  66. C 890831 REVISION DATE from Version 3.2
  67. C 891214 Prologue converted to Version 4.0 format. (BAB)
  68. C 900326 Removed duplicate information from DESCRIPTION section.
  69. C (WRB)
  70. C 920501 Reformatted the REFERENCES section. (WRB)
  71. C***END PROLOGUE CTRSL
  72. INTEGER LDT,N,JOB,INFO
  73. COMPLEX T(LDT,*),B(*)
  74. C
  75. C
  76. COMPLEX CDOTC,TEMP
  77. INTEGER CASE,J,JJ
  78. COMPLEX ZDUM
  79. REAL CABS1
  80. CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
  81. C***FIRST EXECUTABLE STATEMENT CTRSL
  82. C
  83. C CHECK FOR ZERO DIAGONAL ELEMENTS.
  84. C
  85. DO 10 INFO = 1, N
  86. IF (CABS1(T(INFO,INFO)) .EQ. 0.0E0) GO TO 150
  87. 10 CONTINUE
  88. INFO = 0
  89. C
  90. C DETERMINE THE TASK AND GO TO IT.
  91. C
  92. CASE = 1
  93. IF (MOD(JOB,10) .NE. 0) CASE = 2
  94. IF (MOD(JOB,100)/10 .NE. 0) CASE = CASE + 2
  95. GO TO (20,50,80,110), CASE
  96. C
  97. C SOLVE T*X=B FOR T LOWER TRIANGULAR
  98. C
  99. 20 CONTINUE
  100. B(1) = B(1)/T(1,1)
  101. IF (N .LT. 2) GO TO 40
  102. DO 30 J = 2, N
  103. TEMP = -B(J-1)
  104. CALL CAXPY(N-J+1,TEMP,T(J,J-1),1,B(J),1)
  105. B(J) = B(J)/T(J,J)
  106. 30 CONTINUE
  107. 40 CONTINUE
  108. GO TO 140
  109. C
  110. C SOLVE T*X=B FOR T UPPER TRIANGULAR.
  111. C
  112. 50 CONTINUE
  113. B(N) = B(N)/T(N,N)
  114. IF (N .LT. 2) GO TO 70
  115. DO 60 JJ = 2, N
  116. J = N - JJ + 1
  117. TEMP = -B(J+1)
  118. CALL CAXPY(J,TEMP,T(1,J+1),1,B(1),1)
  119. B(J) = B(J)/T(J,J)
  120. 60 CONTINUE
  121. 70 CONTINUE
  122. GO TO 140
  123. C
  124. C SOLVE CTRANS(T)*X=B FOR T LOWER TRIANGULAR.
  125. C
  126. 80 CONTINUE
  127. B(N) = B(N)/CONJG(T(N,N))
  128. IF (N .LT. 2) GO TO 100
  129. DO 90 JJ = 2, N
  130. J = N - JJ + 1
  131. B(J) = B(J) - CDOTC(JJ-1,T(J+1,J),1,B(J+1),1)
  132. B(J) = B(J)/CONJG(T(J,J))
  133. 90 CONTINUE
  134. 100 CONTINUE
  135. GO TO 140
  136. C
  137. C SOLVE CTRANS(T)*X=B FOR T UPPER TRIANGULAR.
  138. C
  139. 110 CONTINUE
  140. B(1) = B(1)/CONJG(T(1,1))
  141. IF (N .LT. 2) GO TO 130
  142. DO 120 J = 2, N
  143. B(J) = B(J) - CDOTC(J-1,T(1,J),1,B(1),1)
  144. B(J) = B(J)/CONJG(T(J,J))
  145. 120 CONTINUE
  146. 130 CONTINUE
  147. 140 CONTINUE
  148. 150 CONTINUE
  149. RETURN
  150. END