1 subroutine cppfa(ap,n,info) 2 integer n,info 3 complex ap(1) 4c 5c cppfa factors a complex hermitian positive definite matrix 6c stored in packed form. 7c 8c cppfa is usually called by cppco, but it can be called 9c directly with a saving in time if rcond is not needed. 10c (time for cppco) = (1 + 18/n)*(time for cppfa) . 11c 12c on entry 13c 14c ap complex (n*(n+1)/2) 15c the packed form of a hermitian matrix a . the 16c columns of the upper triangle are stored sequentially 17c in a one-dimensional array of length n*(n+1)/2 . 18c see comments below for details. 19c 20c n integer 21c the order of the matrix a . 22c 23c on return 24c 25c ap an upper triangular matrix r , stored in packed 26c form, so that a = ctrans(r)*r . 27c 28c info integer 29c = 0 for normal return. 30c = k if the leading minor of order k is not 31c positive definite. 32c 33c 34c packed storage 35c 36c the following program segment will pack the upper 37c triangle of a hermitian matrix. 38c 39c k = 0 40c do 20 j = 1, n 41c do 10 i = 1, j 42c k = k + 1 43c ap(k) = a(i,j) 44c 10 continue 45c 20 continue 46c 47c linpack. this version dated 08/14/78 . 48c cleve moler, university of new mexico, argonne national lab. 49c 50c subroutines and functions 51c 52c blas cdotc 53c fortran aimag,cmplx,conjg,real,sqrt 54c 55c internal variables 56c 57 complex cdotc,t 58 real s 59 integer j,jj,jm1,k,kj,kk 60c begin block with ...exits to 40 61c 62c 63 jj = 0 64 do 30 j = 1, n 65 info = j 66 s = 0.0e0 67 jm1 = j - 1 68 kj = jj 69 kk = 0 70 if (jm1 .lt. 1) go to 20 71 do 10 k = 1, jm1 72 kj = kj + 1 73 t = ap(kj) - cdotc(k-1,ap(kk+1),1,ap(jj+1),1) 74 kk = kk + k 75 t = t/ap(kk) 76 ap(kj) = t 77 s = s + real(t*conjg(t)) 78 10 continue 79 20 continue 80 jj = jj + j 81 s = real(ap(jj)) - s 82c ......exit 83 if (s .le. 0.0e0 .or. aimag(ap(jj)) .ne. 0.0e0) go to 40 84 ap(jj) = cmplx(sqrt(s),0.0e0) 85 30 continue 86 info = 0 87 40 continue 88 return 89 end 90