cppfa.f 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  1. *DECK CPPFA
  2. SUBROUTINE CPPFA (AP, N, INFO)
  3. C***BEGIN PROLOGUE CPPFA
  4. C***PURPOSE Factor a complex Hermitian positive definite matrix stored
  5. C in packed form.
  6. C***LIBRARY SLATEC (LINPACK)
  7. C***CATEGORY D2D1B
  8. C***TYPE COMPLEX (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 CPPFA factors a complex Hermitian positive definite matrix
  15. C stored in packed form.
  16. C
  17. C CPPFA is usually called by CPPCO, but it can be called
  18. C directly with a saving in time if RCOND is not needed.
  19. C (Time for CPPCO) = (1 + 18/N)*(Time for CPPFA) .
  20. C
  21. C On Entry
  22. C
  23. C AP COMPLEX (N*(N+1)/2)
  24. C the packed form of a Hermitian 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 = CTRANS(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 Hermitian 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 CDOTC
  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 CPPFA
  68. INTEGER N,INFO
  69. COMPLEX AP(*)
  70. C
  71. COMPLEX CDOTC,T
  72. REAL S
  73. INTEGER J,JJ,JM1,K,KJ,KK
  74. C***FIRST EXECUTABLE STATEMENT CPPFA
  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) - CDOTC(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 + REAL(T*CONJG(T))
  90. 10 CONTINUE
  91. 20 CONTINUE
  92. JJ = JJ + J
  93. S = REAL(AP(JJ)) - S
  94. IF (S .LE. 0.0E0 .OR. AIMAG(AP(JJ)) .NE. 0.0E0) GO TO 40
  95. AP(JJ) = CMPLX(SQRT(S),0.0E0)
  96. 30 CONTINUE
  97. INFO = 0
  98. 40 CONTINUE
  99. RETURN
  100. END