1      subroutine chifa(a,lda,n,kpvt,info)
2      integer lda,n,kpvt(1),info
3      complex a(lda,1)
4c
5c     chifa factors a complex hermitian matrix by elimination
6c     with symmetric pivoting.
7c
8c     to solve  a*x = b , follow chifa by chisl.
9c     to compute  inverse(a)*c , follow chifa by chisl.
10c     to compute  determinant(a) , follow chifa by chidi.
11c     to compute  inertia(a) , follow chifa by chidi.
12c     to compute  inverse(a) , follow chifa by chidi.
13c
14c     on entry
15c
16c        a       complex(lda,n)
17c                the hermitian matrix to be factored.
18c                only the diagonal and upper triangle are used.
19c
20c        lda     integer
21c                the leading dimension of the array  a .
22c
23c        n       integer
24c                the order of the matrix  a .
25c
26c     on return
27c
28c        a       a block diagonal matrix and the multipliers which
29c                were used to obtain it.
30c                the factorization can be written  a = u*d*ctrans(u)
31c                where  u  is a product of permutation and unit
32c                upper triangular matrices , ctrans(u) is the
33c                conjugate transpose of  u , and  d  is block diagonal
34c                with 1 by 1 and 2 by 2 blocks.
35c
36c        kpvt    integer(n)
37c                an integer vector of pivot indices.
38c
39c        info    integer
40c                = 0  normal value.
41c                = k  if the k-th pivot block is singular. this is
42c                     not an error condition for this subroutine,
43c                     but it does indicate that chisl or chidi may
44c                     divide by zero if called.
45c
46c     linpack. this version dated 08/14/78 .
47c     james bunch, univ. calif. san diego, argonne nat. lab.
48c
49c     subroutines and functions
50c
51c     blas caxpy,cswap,icamax
52c     fortran abs,aimag,amax1,cmplx,conjg,real,sqrt
53c
54c     internal variables
55c
56      complex ak,akm1,bk,bkm1,denom,mulk,mulkm1,t
57      real absakk,alpha,colmax,rowmax
58      integer imax,imaxp1,j,jj,jmax,k,km1,km2,kstep,icamax
59      logical swap
60c
61      complex zdum
62      real cabs1
63      cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum))
64c
65c     initialize
66c
67c     alpha is used in choosing pivot block size.
68      alpha = (1.0e0 + sqrt(17.0e0))/8.0e0
69c
70      info = 0
71c
72c     main loop on k, which goes from n to 1.
73c
74      k = n
75   10 continue
76c
77c        leave the loop if k=0 or k=1.
78c
79c     ...exit
80         if (k .eq. 0) go to 200
81         if (k .gt. 1) go to 20
82            kpvt(1) = 1
83            if (cabs1(a(1,1)) .eq. 0.0e0) info = 1
84c     ......exit
85            go to 200
86   20    continue
87c
88c        this section of code determines the kind of
89c        elimination to be performed.  when it is completed,
90c        kstep will be set to the size of the pivot block, and
91c        swap will be set to .true. if an interchange is
92c        required.
93c
94         km1 = k - 1
95         absakk = cabs1(a(k,k))
96c
97c        determine the largest off-diagonal element in
98c        column k.
99c
100         imax = icamax(k-1,a(1,k),1)
101         colmax = cabs1(a(imax,k))
102         if (absakk .lt. alpha*colmax) go to 30
103            kstep = 1
104            swap = .false.
105         go to 90
106   30    continue
107c
108c           determine the largest off-diagonal element in
109c           row imax.
110c
111            rowmax = 0.0e0
112            imaxp1 = imax + 1
113            do 40 j = imaxp1, k
114               rowmax = amax1(rowmax,cabs1(a(imax,j)))
115   40       continue
116            if (imax .eq. 1) go to 50
117               jmax = icamax(imax-1,a(1,imax),1)
118               rowmax = amax1(rowmax,cabs1(a(jmax,imax)))
119   50       continue
120            if (cabs1(a(imax,imax)) .lt. alpha*rowmax) go to 60
121               kstep = 1
122               swap = .true.
123            go to 80
124   60       continue
125            if (absakk .lt. alpha*colmax*(colmax/rowmax)) go to 70
126               kstep = 1
127               swap = .false.
128            go to 80
129   70       continue
130               kstep = 2
131               swap = imax .ne. km1
132   80       continue
133   90    continue
134         if (amax1(absakk,colmax) .ne. 0.0e0) go to 100
135c
136c           column k is zero.  set info and iterate the loop.
137c
138            kpvt(k) = k
139            info = k
140         go to 190
141  100    continue
142         if (kstep .eq. 2) go to 140
143c
144c           1 x 1 pivot block.
145c
146            if (.not.swap) go to 120
147c
148c              perform an interchange.
149c
150               call cswap(imax,a(1,imax),1,a(1,k),1)
151               do 110 jj = imax, k
152                  j = k + imax - jj
153                  t = conjg(a(j,k))
154                  a(j,k) = conjg(a(imax,j))
155                  a(imax,j) = t
156  110          continue
157  120       continue
158c
159c           perform the elimination.
160c
161            do 130 jj = 1, km1
162               j = k - jj
163               mulk = -a(j,k)/a(k,k)
164               t = conjg(mulk)
165               call caxpy(j,t,a(1,k),1,a(1,j),1)
166               a(j,j) = cmplx(real(a(j,j)),0.0e0)
167               a(j,k) = mulk
168  130       continue
169c
170c           set the pivot array.
171c
172            kpvt(k) = k
173            if (swap) kpvt(k) = imax
174         go to 190
175  140    continue
176c
177c           2 x 2 pivot block.
178c
179            if (.not.swap) go to 160
180c
181c              perform an interchange.
182c
183               call cswap(imax,a(1,imax),1,a(1,k-1),1)
184               do 150 jj = imax, km1
185                  j = km1 + imax - jj
186                  t = conjg(a(j,k-1))
187                  a(j,k-1) = conjg(a(imax,j))
188                  a(imax,j) = t
189  150          continue
190               t = a(k-1,k)
191               a(k-1,k) = a(imax,k)
192               a(imax,k) = t
193  160       continue
194c
195c           perform the elimination.
196c
197            km2 = k - 2
198            if (km2 .eq. 0) go to 180
199               ak = a(k,k)/a(k-1,k)
200               akm1 = a(k-1,k-1)/conjg(a(k-1,k))
201               denom = 1.0e0 - ak*akm1
202               do 170 jj = 1, km2
203                  j = km1 - jj
204                  bk = a(j,k)/a(k-1,k)
205                  bkm1 = a(j,k-1)/conjg(a(k-1,k))
206                  mulk = (akm1*bk - bkm1)/denom
207                  mulkm1 = (ak*bkm1 - bk)/denom
208                  t = conjg(mulk)
209                  call caxpy(j,t,a(1,k),1,a(1,j),1)
210                  t = conjg(mulkm1)
211                  call caxpy(j,t,a(1,k-1),1,a(1,j),1)
212                  a(j,k) = mulk
213                  a(j,k-1) = mulkm1
214                  a(j,j) = cmplx(real(a(j,j)),0.0e0)
215  170          continue
216  180       continue
217c
218c           set the pivot array.
219c
220            kpvt(k) = 1 - k
221            if (swap) kpvt(k) = -imax
222            kpvt(k-1) = kpvt(k)
223  190    continue
224         k = k - kstep
225      go to 10
226  200 continue
227      return
228      end
229