sppfa.f 3.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  1. *DECK SPPFA
  2. SUBROUTINE SPPFA (AP, N, INFO)
  3. C***BEGIN PROLOGUE SPPFA
  4. C***PURPOSE Factor a real symmetric positive definite matrix stored in
  5. C packed form.
  6. C***LIBRARY SLATEC (LINPACK)
  7. C***CATEGORY D2B1B
  8. C***TYPE SINGLE PRECISION (SPPFA-S, DPPFA-D, CPPFA-C)
  9. C***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, PACKED,
  10. C POSITIVE DEFINITE
  11. C***AUTHOR Moler, C. B., (U. of New Mexico)
  12. C***DESCRIPTION
  13. C
  14. C SPPFA factors a real symmetric positive definite matrix
  15. C stored in packed form.
  16. C
  17. C SPPFA is usually called by SPPCO, but it can be called
  18. C directly with a saving in time if RCOND is not needed.
  19. C (Time for SPPCO) = (1 + 18/N)*(Time for SPPFA) .
  20. C
  21. C On Entry
  22. C
  23. C AP REAL (N*(N+1)/2)
  24. C the packed form of a symmetric matrix A . The
  25. C columns of the upper triangle are stored sequentially
  26. C in a one-dimensional array of length N*(N+1)/2 .
  27. C See comments below for details.
  28. C
  29. C N INTEGER
  30. C the order of the matrix A .
  31. C
  32. C On Return
  33. C
  34. C AP an upper triangular matrix R , stored in packed
  35. C form, so that A = TRANS(R)*R .
  36. C
  37. C INFO INTEGER
  38. C = 0 for normal return.
  39. C = K if the leading minor of order K is not
  40. C positive definite.
  41. C
  42. C
  43. C Packed Storage
  44. C
  45. C The following program segment will pack the upper
  46. C triangle of a symmetric matrix.
  47. C
  48. C K = 0
  49. C DO 20 J = 1, N
  50. C DO 10 I = 1, J
  51. C K = K + 1
  52. C AP(K) = A(I,J)
  53. C 10 CONTINUE
  54. C 20 CONTINUE
  55. C
  56. C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
  57. C Stewart, LINPACK Users' Guide, SIAM, 1979.
  58. C***ROUTINES CALLED SDOT
  59. C***REVISION HISTORY (YYMMDD)
  60. C 780814 DATE WRITTEN
  61. C 890831 Modified array declarations. (WRB)
  62. C 890831 REVISION DATE from Version 3.2
  63. C 891214 Prologue converted to Version 4.0 format. (BAB)
  64. C 900326 Removed duplicate information from DESCRIPTION section.
  65. C (WRB)
  66. C 920501 Reformatted the REFERENCES section. (WRB)
  67. C***END PROLOGUE SPPFA
  68. INTEGER N,INFO
  69. REAL AP(*)
  70. C
  71. REAL SDOT,T
  72. REAL S
  73. INTEGER J,JJ,JM1,K,KJ,KK
  74. C***FIRST EXECUTABLE STATEMENT SPPFA
  75. JJ = 0
  76. DO 30 J = 1, N
  77. INFO = J
  78. S = 0.0E0
  79. JM1 = J - 1
  80. KJ = JJ
  81. KK = 0
  82. IF (JM1 .LT. 1) GO TO 20
  83. DO 10 K = 1, JM1
  84. KJ = KJ + 1
  85. T = AP(KJ) - SDOT(K-1,AP(KK+1),1,AP(JJ+1),1)
  86. KK = KK + K
  87. T = T/AP(KK)
  88. AP(KJ) = T
  89. S = S + T*T
  90. 10 CONTINUE
  91. 20 CONTINUE
  92. JJ = JJ + J
  93. S = AP(JJ) - S
  94. IF (S .LE. 0.0E0) GO TO 40
  95. AP(JJ) = SQRT(S)
  96. 30 CONTINUE
  97. INFO = 0
  98. 40 CONTINUE
  99. RETURN
  100. END