subroutine zpbfa(abd,lda,n,m,info) integer lda,n,m,info complex*16 abd(lda,1) c c zpbfa factors a complex*16 hermitian positive definite matrix c stored in band form. c c zpbfa is usually called by zpbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry c c abd complex*16(lda, n) c the matrix to be factored. the columns of the upper c triangle are stored in the columns of abd and the c diagonals of the upper triangle are stored in the c rows of abd . see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. m + 1 . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c 0 .le. m .lt. n . c c on return c c abd an upper triangular matrix r , stored in band c form, so that a = ctrans(r)*r . c c info integer c = 0 for normal return. c = k if the leading minor of order k is not c positive definite. c c band storage c c if a is a hermitian positive definite band matrix, c the following program segment will set up the input. c c m = (band width above diagonal) c do 20 j = 1, n c i1 = max0(1, j-m) c do 10 i = i1, j c k = i-j+m+1 c abd(k,j) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas zdotc c fortran dcmplx,dconjg,max0,dsqrt c c internal variables c complex*16 zdotc,t double precision s integer ik,j,jk,k,mu double precision dreal,dimag complex*16 zdumr,zdumi dreal(zdumr) = zdumr dimag(zdumi) = (0.0d0,-1.0d0)*zdumi c begin block with ...exits to 40 c c do 30 j = 1, n info = j s = 0.0d0 ik = m + 1 jk = max0(j-m,1) mu = max0(m+2-j,1) if (m .lt. mu) go to 20 do 10 k = mu, m t = abd(k,j) - zdotc(k-mu,abd(ik,jk),1,abd(mu,j),1) t = t/abd(m+1,jk) abd(k,j) = t s = s + dreal(t*dconjg(t)) ik = ik - 1 jk = jk + 1 10 continue 20 continue s = dreal(abd(m+1,j)) - s c ......exit if (s .le. 0.0d0 .or. dimag(abd(m+1,j)) .ne. 0.0d0) * go to 40 abd(m+1,j) = dcmplx(dsqrt(s),0.0d0) 30 continue info = 0 40 continue return end