dhels.f 3.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798
  1. *DECK DHELS
  2. SUBROUTINE DHELS (A, LDA, N, Q, B)
  3. C***BEGIN PROLOGUE DHELS
  4. C***SUBSIDIARY
  5. C***PURPOSE Internal routine for DGMRES.
  6. C***LIBRARY SLATEC (SLAP)
  7. C***CATEGORY D2A4, D2B4
  8. C***TYPE DOUBLE PRECISION (SHELS-S, DHELS-D)
  9. C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
  10. C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
  11. C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
  12. C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
  13. C Seager, Mark K., (LLNL), seager@llnl.gov
  14. C Lawrence Livermore National Laboratory
  15. C PO Box 808, L-60
  16. C Livermore, CA 94550 (510) 423-3141
  17. C***DESCRIPTION
  18. C This routine is extracted from the LINPACK routine SGESL with
  19. C changes due to the fact that A is an upper Hessenberg matrix.
  20. C
  21. C DHELS solves the least squares problem:
  22. C
  23. C MIN(B-A*X,B-A*X)
  24. C
  25. C using the factors computed by DHEQR.
  26. C
  27. C *Usage:
  28. C INTEGER LDA, N
  29. C DOUBLE PRECISION A(LDA,N), Q(2*N), B(N+1)
  30. C
  31. C CALL DHELS(A, LDA, N, Q, B)
  32. C
  33. C *Arguments:
  34. C A :IN Double Precision A(LDA,N)
  35. C The output from DHEQR which contains the upper
  36. C triangular factor R in the QR decomposition of A.
  37. C LDA :IN Integer
  38. C The leading dimension of the array A.
  39. C N :IN Integer
  40. C A is originally an (N+1) by N matrix.
  41. C Q :IN Double Precision Q(2*N)
  42. C The coefficients of the N Givens rotations
  43. C used in the QR factorization of A.
  44. C B :INOUT Double Precision B(N+1)
  45. C On input, B is the right hand side vector.
  46. C On output, B is the solution vector X.
  47. C
  48. C***SEE ALSO DGMRES
  49. C***ROUTINES CALLED DAXPY
  50. C***REVISION HISTORY (YYMMDD)
  51. C 890404 DATE WRITTEN
  52. C 890404 Previous REVISION DATE
  53. C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
  54. C 890922 Numerous changes to prologue to make closer to SLATEC
  55. C standard. (FNF)
  56. C 890929 Numerous changes to reduce SP/DP differences. (FNF)
  57. C 910411 Prologue converted to Version 4.0 format. (BAB)
  58. C 910502 Added C***FIRST EXECUTABLE STATEMENT line. (FNF)
  59. C 910506 Made subsidiary to DGMRES. (FNF)
  60. C 920511 Added complete declaration section. (WRB)
  61. C***END PROLOGUE DHELS
  62. C The following is for optimized compilation on LLNL/LTSS Crays.
  63. CLLL. OPTIMIZE
  64. C .. Scalar Arguments ..
  65. INTEGER LDA, N
  66. C .. Array Arguments ..
  67. DOUBLE PRECISION A(LDA,*), B(*), Q(*)
  68. C .. Local Scalars ..
  69. DOUBLE PRECISION C, S, T, T1, T2
  70. INTEGER IQ, K, KB, KP1
  71. C .. External Subroutines ..
  72. EXTERNAL DAXPY
  73. C***FIRST EXECUTABLE STATEMENT DHELS
  74. C
  75. C Minimize(B-A*X,B-A*X). First form Q*B.
  76. C
  77. DO 20 K = 1, N
  78. KP1 = K + 1
  79. IQ = 2*(K-1) + 1
  80. C = Q(IQ)
  81. S = Q(IQ+1)
  82. T1 = B(K)
  83. T2 = B(KP1)
  84. B(K) = C*T1 - S*T2
  85. B(KP1) = S*T1 + C*T2
  86. 20 CONTINUE
  87. C
  88. C Now solve R*X = Q*B.
  89. C
  90. DO 40 KB = 1, N
  91. K = N + 1 - KB
  92. B(K) = B(K)/A(K,K)
  93. T = -B(K)
  94. CALL DAXPY(K-1, T, A(1,K), 1, B(1), 1)
  95. 40 CONTINUE
  96. RETURN
  97. C------------- LAST LINE OF DHELS FOLLOWS ----------------------------
  98. END