bnslv.f 2.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
  1. *DECK BNSLV
  2. SUBROUTINE BNSLV (W, NROWW, NROW, NBANDL, NBANDU, B)
  3. C***BEGIN PROLOGUE BNSLV
  4. C***SUBSIDIARY
  5. C***PURPOSE Subsidiary to BINT4 and BINTK
  6. C***LIBRARY SLATEC
  7. C***TYPE SINGLE PRECISION (BNSLV-S, DBNSLV-D)
  8. C***AUTHOR (UNKNOWN)
  9. C***DESCRIPTION
  10. C
  11. C BNSLV is the BANSLV routine from
  12. C * A Practical Guide to Splines * by C. de Boor
  13. C
  14. C Companion routine to BNFAC . It returns the solution X of the
  15. C linear system A*X = B in place of B , given the LU-factorization
  16. C for A in the work array W from BNFAC.
  17. C
  18. C ***** I N P U T ******
  19. C W, NROWW,NROW,NBANDL,NBANDU.....Describe the LU-factorization of a
  20. C banded matrix A of order NROW as constructed in BNFAC .
  21. C For details, see BNFAC .
  22. C B.....Right side of the system to be solved .
  23. C
  24. C ***** O U T P U T ******
  25. C B.....Contains the solution X , of order NROW .
  26. C
  27. C ***** M E T H O D ******
  28. C (With A = L*U, as stored in W,) the unit lower triangular system
  29. C L(U*X) = B is solved for Y = U*X, and Y stored in B . Then the
  30. C upper triangular system U*X = Y is solved for X . The calcul-
  31. C ations are so arranged that the innermost loops stay within columns.
  32. C
  33. C***SEE ALSO BINT4, BINTK
  34. C***ROUTINES CALLED (NONE)
  35. C***REVISION HISTORY (YYMMDD)
  36. C 800901 DATE WRITTEN
  37. C 890531 Changed all specific intrinsics to generic. (WRB)
  38. C 890831 Modified array declarations. (WRB)
  39. C 891214 Prologue converted to Version 4.0 format. (BAB)
  40. C 900328 Added TYPE section. (WRB)
  41. C***END PROLOGUE BNSLV
  42. C
  43. INTEGER NBANDL, NBANDU, NROW, NROWW, I, J, JMAX, MIDDLE, NROWM1
  44. REAL W(NROWW,*), B(*)
  45. C***FIRST EXECUTABLE STATEMENT BNSLV
  46. MIDDLE = NBANDU + 1
  47. IF (NROW.EQ.1) GO TO 80
  48. NROWM1 = NROW - 1
  49. IF (NBANDL.EQ.0) GO TO 30
  50. C FORWARD PASS
  51. C FOR I=1,2,...,NROW-1, SUBTRACT RIGHT SIDE(I)*(I-TH COLUMN
  52. C OF L ) FROM RIGHT SIDE (BELOW I-TH ROW) .
  53. DO 20 I=1,NROWM1
  54. JMAX = MIN(NBANDL,NROW-I)
  55. DO 10 J=1,JMAX
  56. B(I+J) = B(I+J) - B(I)*W(MIDDLE+J,I)
  57. 10 CONTINUE
  58. 20 CONTINUE
  59. C BACKWARD PASS
  60. C FOR I=NROW,NROW-1,...,1, DIVIDE RIGHT SIDE(I) BY I-TH DIAG-
  61. C ONAL ENTRY OF U, THEN SUBTRACT RIGHT SIDE(I)*(I-TH COLUMN
  62. C OF U) FROM RIGHT SIDE (ABOVE I-TH ROW).
  63. 30 IF (NBANDU.GT.0) GO TO 50
  64. C A IS LOWER TRIANGULAR .
  65. DO 40 I=1,NROW
  66. B(I) = B(I)/W(1,I)
  67. 40 CONTINUE
  68. RETURN
  69. 50 I = NROW
  70. 60 B(I) = B(I)/W(MIDDLE,I)
  71. JMAX = MIN(NBANDU,I-1)
  72. DO 70 J=1,JMAX
  73. B(I-J) = B(I-J) - B(I)*W(MIDDLE-J,I)
  74. 70 CONTINUE
  75. I = I - 1
  76. IF (I.GT.1) GO TO 60
  77. 80 B(1) = B(1)/W(MIDDLE,1)
  78. RETURN
  79. END