subroutine cchdc ( a, lda, p, work, ipvt, job, info ) !*****************************************************************************80 ! !! CCHDC: Cholesky decomposition of a Hermitian positive definite matrix. ! ! Discussion: ! ! A pivoting option allows the user to estimate the condition of a ! Hermitian positive definite matrix or determine the rank of a ! Hermitian positive semidefinite matrix. ! ! For Hermitian positive definite matrices, INFO = P is the normal return. ! ! For pivoting with Hermitian positive semidefinite matrices, INFO will ! in general be less than P. However, INFO may be greater than ! the rank of A, since rounding error can cause an otherwise zero ! element to be positive. Indefinite systems will always cause ! INFO to be less than P. ! ! Modified: ! ! 02 May 2007 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,P). On input, A contains the matrix ! whose decomposition is to be computed. Only the upper half of A ! need be stored. The lower part of the array A is not referenced. ! On output, A contains in its upper half the Cholesky factor ! of the matrix A as it has been permuted by pivoting. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer P, the order of the matrix. ! ! Workspace, complex WORK(P). ! ! Input/output, integer IPVT(P). IPVT is not referenced if JOB == 0. ! On input, IPVT contains integers that control the selection of the ! pivot elements, if pivoting has been requested. Each diagonal element ! A(K,K) is placed in one of three classes according to the input ! value of IPVT(K): ! IPVT(K) > 0, X(K) is an initial element. ! IPVT(K) == 0, X(K) is a free element. ! IPVT(K) < 0, X(K) is a final element. ! Before the decomposition is computed, initial elements are moved by ! symmetric row and column interchanges to the beginning of the array A ! and final elements to the end. Both initial and final elements ! are frozen in place during the computation and only free elements ! are moved. At the K-th stage of the reduction, if A(K,K) is occupied ! by a free element, it is interchanged with the largest free element ! A(L,L) with K <= L. ! On output, IPVT(K) contains the index of the diagonal element ! of A that was moved into the J-th position, if pivoting was requested. ! ! Input, integer JOB, specifies whether column pivoting is to be done. ! 0, no pivoting is done. ! nonzero, pivoting is done. ! ! Output, integer INFO, contains the index of the last positive ! diagonal element of the Cholesky factor. ! implicit none integer lda integer p complex a(lda,p) integer i integer info integer j integer job integer ipvt(p) integer k integer kb integer l real maxdia integer maxl logical negk integer pl integer plp1 integer pu logical swapk complex t complex temp complex work(p) pl = 1 pu = 0 info = p if ( job /= 0 ) then ! ! Pivoting has been requested. Rearrange the elements according to IPVT. ! do k = 1, p swapk = ( 0 < ipvt(k) ) negk = ( ipvt(k) < 0 ) if ( negk ) then ipvt(k) = -k else ipvt(k) = k end if if ( swapk ) then if ( k /= pl ) then call cswap ( pl-1, a(1:pl-1,k), 1, a(1:pl-1,pl), 1 ) t = a(k,k) a(k,k) = a(pl,pl) a(pl,pl) = t a(pl,k) = conjg ( a(pl,k) ) plp1 = pl + 1 do j = plp1, p if ( j < k ) then t = conjg ( a(pl,j) ) a(pl,j) = conjg ( a(j,k) ) a(j,k) = t else if ( j /= k ) then t = a(k,j) a(k,j) = a(pl,j) a(pl,j) = t end if end do ipvt(k) = ipvt(pl) ipvt(pl) = k end if pl = pl + 1 end if end do pu = p do kb = pl, p k = p - kb + pl if ( ipvt(k) < 0 ) then ipvt(k) = -ipvt(k) if ( pu /= k ) then call cswap ( k-1, a(1:k-1,k), 1, a(1:k-1,pu), 1 ) t = a(k,k) a(k,k) = a(pu,pu) a(pu,pu) = t a(k,pu) = conjg ( a(k,pu) ) do j = k + 1, p if ( j < pu ) then t = conjg ( a(k,j) ) a(k,j) = conjg ( a(j,pu) ) a(j,pu) = t else if ( j /= pu ) then t = a(k,j) a(k,j) = a(pu,j) a(pu,j) = t end if end do i = ipvt(k) ipvt(k) = ipvt(pu) ipvt(pu) = i end if pu = pu - 1 end if end do end if do k = 1, p ! ! Reduction loop. ! maxdia = real ( a(k,k) ) maxl = k ! ! Determine the pivot element. ! if ( pl <= k .and. k < pu ) then do l = k+1, pu if ( maxdia < real ( a(l,l) ) ) then maxdia = real ( a(l,l) ) maxl = l end if end do end if ! ! Quit if the pivot element is not positive. ! if ( maxdia <= 0.0E+00 ) then info = k - 1 return end if ! ! Start the pivoting and update IPVT. ! if ( k /= maxl ) then call cswap ( k-1, a(1:k-1,k), 1, a(1:k-1,maxl), 1 ) a(maxl,maxl) = a(k,k) a(k,k) = cmplx ( maxdia, 0.0E+00 ) i = ipvt(maxl) ipvt(maxl) = ipvt(k) ipvt(k) = i a(k,maxl) = conjg ( a(k,maxl) ) end if ! ! Reduction step. Pivoting is contained across the rows. ! work(k) = cmplx ( sqrt ( real ( a(k,k) ) ), 0.0E+00 ) a(k,k) = work(k) do j = k+1, p if ( k /= maxl ) then if ( j < maxl ) then t = conjg ( a(k,j) ) a(k,j) = conjg ( a(j,maxl) ) a(j,maxl) = t else if ( j /= maxl ) then t = a(k,j) a(k,j) = a(maxl,j) a(maxl,j) = t end if end if a(k,j) = a(k,j) / work(k) work(j) = conjg ( a(k,j) ) temp = -a(k,j) a(k+1:j,j) = a(k+1:j,j) + temp * work(k+1:j) end do end do return end subroutine cchdd ( r, ldr, p, x, z, ldz, nz, y, rho, c, s, info ) !*****************************************************************************80 ! !! CCHDD downdates an augmented Cholesky decomposition. ! ! Discussion: ! ! CCHDD downdates an augmented Cholesky decomposition or the ! triangular factor of an augmented QR decomposition. ! Specifically, given an upper triangular matrix R of order P, a ! row vector X, a column vector Z, and a scalar Y, CCHDD ! determines a unitary matrix U and a scalar ZETA such that ! ! ( R Z ) ( RR ZZ ) ! U * ( ) = ( ), ! ( 0 ZETA ) ( X Y ) ! ! where RR is upper triangular. If R and Z have been obtained ! from the factorization of a least squares problem, then ! RR and ZZ are the factors corresponding to the problem ! with the observation (X,Y) removed. In this case, if RHO ! is the norm of the residual vector, then the norm of ! the residual vector of the downdated problem is ! sqrt ( RHO**2 - ZETA**2 ). ! CCHDD will simultaneously downdate several triplets (Z,Y,RHO) ! along with R. ! ! For a less terse description of what CCHDD does and how ! it may be applied, see the LINPACK guide. ! ! The matrix U is determined as the product U(1)*...*U(P) ! where U(I) is a rotation in the (P+1,I)-plane of the ! form ! ! ( C(I) -conjg ( S(I) ) ) ! ( ). ! ( S(I) C(I) ) ! ! The rotations are chosen so that C(I) is real. ! ! The user is warned that a given downdating problem may ! be impossible to accomplish or may produce ! inaccurate results. For example, this can happen ! if X is near a vector whose removal will reduce the ! rank of R. Beware. ! ! Modified: ! ! 15 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex R(LDR,P); on input, the upper triangular matrix ! that is to be downdated. On output, the downdated matrix. The ! part of R below the diagonal is not referenced. ! ! Input, integer LDR, the leading dimension of R. P <= LDR. ! ! Input, integer P, the order of the matrix. ! ! Input, complex X(P), the row vector that is to ! be removed from R. ! ! Input/output, complex Z(LDZ,NZ); on input, an array of NZ ! P-vectors which are to be downdated along with R. On output, ! the downdated vectors. ! ! Input, integer LDZ, the leading dimension of Z. P <= LDZ. ! ! Input, integer NZ, the number of vectors to be downdated. ! NZ may be zero, in which case Z, Y, and R are not referenced. ! ! Input, complex Y(NZ), the scalars for the downdating ! of the vectors Z. ! ! Input/output, real RHO(NZ). On input, the norms of the residual ! vectors that are to be downdated. On output, the downdated norms. ! ! Output, real C(P), the cosines of the transforming rotations. ! ! Output, complex S(P), the sines of the transforming rotations. ! ! Output, integer INFO: ! 0, if the entire downdating was successful. ! -1, if R could not be downdated. In this case, all quantities ! are left unaltered. ! 1, if some RHO could not be downdated. The offending RHO's are ! set to -1. ! implicit none integer ldr integer ldz integer nz integer p real a real alpha real azeta complex b real c(p) complex cdotc integer i integer ii integer info integer j real norm complex r(ldr,p) real rho(nz) complex s(p) real scale real scnrm2 complex t complex x(p) complex xx complex y(nz) complex z(ldz,nz) complex zeta ! ! Solve the system hermitian(R) * A = X, placing the result in S. ! info = 0 s(1) = conjg ( x(1) ) / conjg ( r(1,1) ) ! ! The FORTRAN90 DOT_PRODUCT automatically conjugates the first vector. ! do j = 2, p s(j) = conjg ( x(j) ) - dot_product ( r(1:j-1,j), s(1:j-1) ) s(j) = s(j) / conjg ( r(j,j) ) end do norm = scnrm2 ( p, s, 1 ) if ( 1.0E+00 <= norm ) then info = -1 return end if alpha = sqrt ( 1.0E+00 - norm**2 ) ! ! Determine the transformations. ! do ii = 1, p i = p - ii + 1 scale = alpha + abs ( s(i) ) a = alpha / scale b = s(i) / scale norm = sqrt ( a**2 + ( real ( b ) )**2 + ( aimag ( b ) )**2 ) c(i) = a / norm s(i) = conjg ( b ) / norm alpha = scale * norm end do ! ! Apply the transformations to R. ! do j = 1, p xx = cmplx ( 0.0E+00, 0.0E+00 ) do ii = 1, j i = j - ii + 1 t = c(i) * xx + s(i) * r(i,j) r(i,j) = c(i) * r(i,j) - conjg ( s(i) ) * xx xx = t end do end do ! ! If required, downdate Z and RHO. ! do j = 1, nz zeta = y(j) do i = 1, p z(i,j) = ( z(i,j) - conjg ( s(i) ) * zeta ) / c(i) zeta = c(i) * zeta - s(i) * z(i,j) end do azeta = abs ( zeta ) if ( rho(j) < azeta ) then info = 1 rho(j) = -1.0E+00 else rho(j) = rho(j) * sqrt ( 1.0E+00 - ( azeta / rho(j) )**2 ) end if end do return end subroutine cchex ( r, ldr, p, k, l, z, ldz, nz, c, s, job ) !*****************************************************************************80 ! !! CCHEX updates a Cholesky factorization. ! ! Discussion: ! ! CCHEX updates a Cholesky factorization ! ! A = hermitian(R) * R ! ! of a positive definite matrix A of order P under diagonal ! permutations of the form ! ! E' * A * E ! ! where E is a permutation matrix. Specifically, given ! an upper triangular matrix R and a permutation matrix ! E (which is specified by K, L, and JOB), CCHEX determines ! a unitary matrix U such that ! ! U * R * E = RR, ! ! where RR is upper triangular. At the user's option, the ! transformation U will be multiplied into the array Z. ! ! If A = hermitian(X)*X, so that R is the triangular part of the ! QR factorization of X, then RR is the triangular part of the ! QR factorization of X * E, that is, X with its columns permuted. ! ! For a less terse description of what CCHEX does and how ! it may be applied, see the LINPACK guide. ! ! The matrix Q is determined as the product U(L-K)*...*U(1) ! of plane rotations of the form ! ! ( C(I) S(I) ) ! ( ) , ! ( -conjg(S(i)) C(I) ) ! ! where C(I) is real, the rows these rotations operate on ! are described below. ! ! There are two types of permutations, which are determined ! by the value of job. ! ! JOB = 1, right circular shift: ! The columns are rearranged in the following order. ! ! 1, ..., K-1, L, K, K+1, ..., L-1, L+1, ..., P. ! ! U is the product of L-K rotations U(I), where U(I) ! acts in the (L-I,L-I+1)-plane. ! ! JOB = 2, left circular shift: ! The columns are rearranged in the following order ! ! 1, ..., K-1, K+1, K+2, ..., L, L, L+1, ..., P. ! ! U is the product of L-K rotations U(I), where U(I) ! acts in the (K+I-1,K+I)-plane. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex R(LDR,P); On input, the upper triangular factor ! that is to be updated. On output, the updated factor. Elements ! below the diagonal are not referenced. ! ! Input, integer LDR, the leading dimension of R, which is at least P. ! ! Input, integer P, the order of the matrix. ! ! Input, integer K, the first column to be permuted. ! ! Input, integer L, the last column to be permuted. ! L must be strictly greater than K. ! ! Input/output, complex Z(LDZ,NZ); on input, an array of NZ P-vectors into ! which the transformation U is multiplied. On output, the updated ! matrix. Z is not referenced if NZ = 0. ! ! Input, integer LDZ, the leading dimension of Z, which must ! be at least P. ! ! Input, integer NZ, the number of columns of the matrix Z. ! ! Output, real C(P), the cosines of the transforming rotations. ! ! Output, complex S(P), the sines of the transforming rotations. ! ! Input, integer JOB, determines the type of permutation. ! 1, right circular shift. ! 2, left circular shift. ! implicit none integer ldr integer ldz integer nz integer p real c(p) integer i integer ii integer il integer iu integer j integer jj integer job integer k integer l complex r(ldr,p) complex s(p) complex t complex z(ldz,nz) if ( job == 1 ) then ! ! Right circular shift. ! ! Reorder the columns. ! do i = 1, l ii = l - i + 1 s(i) = r(ii,l) end do do jj = k, l - 1 j = l - 1 - jj + k r(1:j,j+1) = r(1:j,j) r(j+1,j+1) = cmplx ( 0.0E+00, 0.0E+00 ) end do do i = 1, k - 1 ii = l - i + 1 r(i,k) = s(ii) end do ! ! Calculate the rotations. ! t = s(1) do i = 1, l - k call crotg ( s(i+1), t, c(i), s(i) ) t = s(i+1) end do r(k,k) = t do j = k+1, p il = max ( 1, l - j + 1 ) do ii = il, l - k i = l - ii t = c(ii) * r(i,j) + s(ii) * r(i+1,j) r(i+1,j) = c(ii) * r(i+1,j) - conjg ( s(ii) ) * r(i,j) r(i,j) = t end do end do ! ! If required, apply the transformations to Z. ! do j = 1, nz do ii = 1, l - k i = l - ii t = c(ii) * z(i,j) + s(ii) * z(i+1,j) z(i+1,j) = c(ii) * z(i+1,j) - conjg ( s(ii) ) * z(i,j) z(i,j) = t end do end do else ! ! Left circular shift. ! ! Reorder the columns. ! do i = 1, k ii = l - k + i s(ii) = r(i,k) end do do j = k, l - 1 r(1:j,j) = r(1:j,j+1) jj = j - k + 1 s(jj) = r(j+1,j+1) end do do i = 1, k ii = l - k + i r(i,l) = s(ii) end do r(k+1:l,l) = cmplx ( 0.0E+00, 0.0E+00 ) ! ! Reduction loop. ! do j = k, p ! ! Apply the rotations. ! if ( j /= k ) then iu = min ( j - 1, l - 1 ) do i = k, iu ii = i - k + 1 t = c(ii) * r(i,j) + s(ii) * r(i+1,j) r(i+1,j) = c(ii) * r(i+1,j) - conjg ( s(ii) ) * r(i,j) r(i,j) = t end do end if if ( j < l ) then jj = j - k + 1 t = s(jj) call crotg ( r(j,j), t, c(jj), s(jj) ) end if end do ! ! Apply the rotations to Z. ! do j = 1, nz do i = k, l - 1 ii = i - k + 1 t = c(ii) * z(i,j) + s(ii) * z(i+1,j) z(i+1,j) = c(ii) * z(i+1,j) - conjg ( s(ii) ) * z(i,j) z(i,j) = t end do end do end if return end subroutine cchud ( r, ldr, p, x, z, ldz, nz, y, rho, c, s ) !*****************************************************************************80 ! !! CCHUD updates an augmented Cholesky decomposition. ! ! Discussion: ! ! CCHUD updates an augmented Cholesky decomposition of the ! triangular part of an augmented QR decomposition. Specifically, ! given an upper triangular matrix R of order P, a row vector ! X, a column vector Z, and a scalar Y, CCHUD determines a ! unitary matrix U and a scalar ZETA such that ! ! ( R Z ) ( RR ZZ ) ! U * ( ) = ( ), ! ( X Y ) ( 0 ZETA ) ! ! where RR is upper triangular. If R and Z have been ! obtained from the factorization of a least squares ! problem, then RR and ZZ are the factors corresponding to ! the problem with the observation (X,Y) appended. In this ! case, if RHO is the norm of the residual vector, then the ! norm of the residual vector of the updated problem is ! sqrt ( RHO**2 + ZETA**2 ). CCHUD will simultaneously update ! several triplets (Z,Y,RHO). ! ! For a less terse description of what CCHUD does and how ! it may be applied see the LINPACK guide. ! ! The matrix U is determined as the product U(P)*...*U(1), ! where U(I) is a rotation in the (I,P+1) plane of the ! form ! ! ( C(I) S(I) ) ! ( ). ! ( -conjg ( S(I) ) C(I) ) ! ! The rotations are chosen so that C(I) is real. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex R(LDR,P), the upper triangular matrix ! that is to be updated. The part of R below the diagonal is ! not referenced. ! ! Input, integer LDR, the leading dimension of R. ! P <= LDR. ! ! Input, integer P, the order of the matrix. ! ! Input, complex X(P), the row to be added to R. ! ! Input/output, complex Z(LDZ,NZ), NZ P-vectors to ! be updated with R. ! ! Input, integer LDZ, the leading dimension of Z. ! P <= LDZ. ! ! Integer, integer NZ, the number of vectors to be updated. ! NZ may be zero, in which case Z, Y, and RHO are not referenced. ! ! Input, complex Y(NZ), the scalars for updating the vectors Z. ! ! Input/output, real RHO(NZ); on input, the norms of the residual ! vectors that are to be updated. If RHO(J) is negative, it is ! left unaltered. On output, the updated values. ! ! Output, real C(P). the cosines of the transforming rotations. ! ! Output, complex S(P), the sines of the transforming rotations. ! implicit none integer ldr integer ldz integer nz integer p real azeta real c(p) integer i integer j complex r(ldr,p) real rho(nz) complex s(p) real scale complex t complex x(p) complex xj complex y(nz) complex z(ldz,nz) complex zeta ! ! Update R. ! do j = 1, p xj = x(j) ! ! Apply the previous rotations. ! do i = 1, j - 1 t = c(i) * r(i,j) + s(i) * xj xj = c(i) * xj - conjg ( s(i) ) * r(i,j) r(i,j) = t end do ! ! Compute the next rotation. ! call crotg ( r(j,j), xj, c(j), s(j) ) end do ! ! If required, update Z and RHO. ! do j = 1, nz zeta = y(j) do i = 1, p t = c(i) * z(i,j) + s(i) * zeta zeta = c(i) * zeta - conjg ( s(i) ) * z(i,j) z(i,j) = t end do azeta = abs ( zeta ) if ( azeta /= 0.0E+00 .and. 0.0E+00 <= rho(j) ) then scale = azeta + rho(j) rho(j) = scale * sqrt ( ( azeta / scale )**2 + ( rho(j) / scale )**2 ) end if end do return end subroutine cgbco ( abd, lda, n, ml, mu, ipvt, rcond, z ) !*****************************************************************************80 ! !! CGBCO factors a complex band matrix and estimates its condition. ! ! Discussion: ! ! If RCOND is not needed, CGBFA is slightly faster. ! ! To solve A*X = B, follow CGBCO by CGBSL. ! ! To compute inverse(A)*C, follow CGBCO by CGBSL. ! ! To compute determinant(A), follow CGBCO by CGBDI. ! ! Band storage: ! ! If A is a band matrix, the following program segment ! will set up the input. ! ! ml = (band width below the diagonal) ! mu = (band width above the diagonal) ! m = ml + mu + 1 ! do j = 1, n ! i1 = max ( 1, j - mu ) ! i2 = min ( n, j + ml ) ! do i = i1, i2 ! k = i - j + m ! abd(k,j) = a(i,j) ! end do ! end do ! ! This uses rows ML+1 through 2*ML+MU+1 of ABD. ! In addition, the first ML rows in ABD are used for ! elements generated during the triangularization. ! The total number of rows needed in ABD is 2*ML+MU+1. ! The ML+MU by ML+MU upper left triangle and the ! ML by ML lower right triangle are not referenced. ! ! Example: ! ! If the original matrix A is ! ! 11 12 13 0 0 0 ! 21 22 23 24 0 0 ! 0 32 33 34 35 0 ! 0 0 43 44 45 46 ! 0 0 0 54 55 56 ! 0 0 0 0 65 66 ! ! Then N = 6, ML = 1, MU = 2, 5 <= LDA and ABD should contain ! ! * * * + + + ! * * 13 24 35 46 ! * 12 23 34 45 56 ! 11 22 33 44 55 66 ! 21 32 43 54 65 * ! ! * = not used, ! + = used for pivoting. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex ABD(LDA,N), on input, contains the matrix in ! band storage. The columns of the matrix are stored in the columns ! of ABD and the diagonals of the matrix are stored in rows ML+1 ! through 2*ML+MU+1 of ABD. On output, an upper triangular matrix ! in band storage and the multipliers which were used to obtain it. ! The factorization can be written A = L*U where L is a product of ! permutation and unit lower triangular matrices and U is upper triangular. ! ! Input, integer LDA, the leading dimension of ABD. ! LDA must be at least 2*ML+MU+1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, the number of diagonals below the main diagonal. ! 0 <= ML < N. ! ! Input, integer MU, the number of diagonals above the main diagonal. ! 0 <= MU < N. ! More efficient if ML <= MU. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! epsilon may cause relative perturbations in X of size (EPSILON/RCOND). ! If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate ! underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is close to a singular matrix, then Z is an ! approximate null vector in the sense that ! norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ). ! implicit none integer lda integer n complex abd(lda,n) real anorm real cabs1 complex cdotc complex csign1 complex ek integer info integer ipvt(n) integer is integer j integer ju integer k integer l integer la integer lm integer lz integer m integer ml integer mm integer mu real rcond real s real scasum real sm complex t complex wk complex wkm real ynorm complex z(n) ! ! Compute 1-norm of A. ! anorm = 0.0E+00 l = ml + 1 is = l + mu do j = 1, n anorm = max ( anorm, scasum ( l, abd(is:is+l-1,j), 1 ) ) if ( ml + 1 < is ) then is = is - 1 end if if ( j <= mu ) then l = l + 1 end if if ( n - ml <= j ) then l = l - 1 end if end do ! ! Factor ! call cgbfa ( abd, lda, n, ml, mu, ipvt, info ) ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and hermitian(A)*Y = E. ! ! Hermitian(A) is the conjugate transpose of A. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where hermitian(U)*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve hermitian(U) * W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) m = ml + mu + 1 ju = 0 do k = 1, n if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csign1 ( ek, -z(k) ) end if if ( cabs1 ( abd(m,k) ) < cabs1 ( ek - z(k) ) ) then s = cabs1 ( abd(m,k) ) / cabs1 ( ek - z(k) ) z(1:n) = z(1:n) * s ek = cmplx ( s, 0.0E+00 ) * ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1 ( wk ) sm = cabs1 ( wkm ) if ( cabs1 ( abd(m,k) ) /= 0.0E+00 ) then wk = wk / conjg ( abd(m,k) ) wkm = wkm / conjg ( abd(m,k) ) else wk = cmplx ( 1.0E+00, 0.0E+00 ) wkm = cmplx ( 1.0E+00, 0.0E+00 ) end if ju = min ( max ( ju, mu + ipvt(k) ), n ) mm = m if ( k+1 <= ju ) then do j = k+1, ju mm = mm - 1 sm = sm + cabs1 ( z(j) + wkm * conjg ( abd(mm,j) ) ) z(j) = z(j) + wk * conjg ( abd(mm,j) ) s = s + cabs1 ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm mm = m do j = k+1, ju mm = mm - 1 z(j) = z(j) + t * conjg ( abd(mm,j) ) end do end if end if z(k) = wk end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ! ! Solve hermitian(L) * Y = W. ! do k = n, 1, -1 lm = min ( ml, n - k ) if ( k < n ) then ! z(k) = z(k) + cdotc ( lm, abd(m+1:m+lm,k), 1, z(k+1:k+lm), 1 ) z(k) = z(k) + dot_product ( abd(m+1:m+lm,k), z(k+1:k+lm) ) end if if ( 1.0E+00 < cabs1 ( z(k) ) ) then s = 1.0E+00 / cabs1 ( z(k) ) z(1:n) = z(1:n) * s end if l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = 1.0E+00 ! ! Solve L * V = Y. ! do k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t lm = min ( ml, n - k ) if ( k < n ) then z(k+1:k+lm) = z(k+1:k+lm) + t * abd(m+1:m+lm,k) end if if ( 1.0E+00 < cabs1 ( z(k) ) ) then s = 1.0E+00 / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm ! ! Solve U * Z = W. ! do k = n, 1, -1 if ( cabs1 ( abd(m,k) ) < cabs1 ( z(k) ) ) then s = cabs1 ( abd(m,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if if ( cabs1 ( abd(m,k) ) /= 0.0E+00 ) then z(k) = z(k) / abd(m,k) else z(k) = cmplx ( 1.0E+00, 0.0E+00 ) end if lm = min ( k, m ) - 1 la = m - lm lz = k - lm t = -z(k) z(lz:lz+lm-1) = z(lz:lz+lm-1) + t * abd(la:la+lm-1,k) end do ! ! Make ZNORM = 1. ! s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine cgbdi ( abd, lda, n, ml, mu, ipvt, det ) !*****************************************************************************80 ! !! CGBDI computes the determinant of a band matrix factored by CGBCO or CGBFA. ! ! Discussion: ! ! If the inverse is needed, use CGBSL N times. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex ABD(LDA,N), the output from CGBCO or CGBFA. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, the number of diagonals below the main diagonal. ! ! Input, integer MU, the number of diagonals above the main diagonal. ! ! Input, integer IPVT(N), the pivot vector from CGBCO or CGBFA. ! ! Output, complex DET(2), determinant of original matrix. ! Determinant = DET(1) * 10.0**DET(2) with 1.0 <= cabs1 ( DET(1) ) < 10.0 ! or DET(1) = 0.0. Also, DET(2) is strictly real. ! implicit none integer lda integer n complex abd(lda,n) real cabs1 complex det(2) integer i integer ipvt(n) integer m integer ml integer mu m = ml + mu + 1 det(1) = cmplx ( 1.0E+00, 0.0E+00 ) det(2) = cmplx ( 0.0E+00, 0.0E+00 ) do i = 1, n if ( ipvt(i) /= i ) then det(1) = -det(1) end if det(1) = det(1) * abd(m,i) if ( cabs1 ( det(1) ) == 0.0E+00 ) then exit end if do while ( cabs1 ( det(1) ) < 1.0E+00 ) det(1) = det(1) * cmplx ( 10.0E+00, 0.0E+00 ) det(2) = det(2) - cmplx ( 1.0E+00, 0.0E+00 ) end do do while ( 10.0E+00 <= cabs1 ( det(1) ) ) det(1) = det(1) / cmplx ( 10.0E+00, 0.0E+00 ) det(2) = det(2) + cmplx ( 1.0E+00, 0.0E+00 ) end do end do return end subroutine cgbfa ( abd, lda, n, ml, mu, ipvt, info ) !*****************************************************************************80 ! !! CGBFA factors a complex band matrix by elimination. ! ! Discussion: ! ! CGBFA is usually called by CGBCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Band storage: ! ! If A is a band matrix, the following program segment ! will set up the input. ! ! ml = (band width below the diagonal) ! mu = (band width above the diagonal) ! m = ml + mu + 1 ! do j = 1, n ! i1 = max ( 1, j - mu ) ! i2 = min ( n, j + ml ) ! do i = i1, i2 ! k = i - j + m ! abd(k,j) = a(i,j) ! end do ! end do ! ! This uses rows ML+1 through 2*ML+MU+1 of ABD. ! In addition, the first ML rows in ABD are used for ! elements generated during the triangularization. ! The total number of rows needed in ABD is 2*ML+MU+1. ! The ML+MU by ML+MU upper left triangle and the ! ML by ML lower right triangle are not referenced. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex ABD(LDA,N), on input, contains the matrix in ! band storage. The columns of the matrix are stored in the columns ! of ABD and the diagonals of the matrix are stored in rows ML+1 ! through 2*ML+MU+1 of ABD. On output, an upper triangular matrix ! in band storage and the multipliers which were used to obtain it. ! The factorization can be written A = L*U where L is a product of ! permutation and unit lower triangular matrices and U is upper triangular. ! ! Input, integer LDA, the leading dimension of ABD. ! LDA must be at least 2*ML+MU+1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, the number of diagonals below the main diagonal. ! 0 <= ML < N. ! ! Input, integer MU, the number of diagonals above the main diagonal. ! 0 <= MU < N. More efficient if ML <= MU. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO. ! 0, normal value. ! K, if U(K,K) == 0.0. This is not an error condition for this ! subroutine, but it does indicate that CGBSL will divide by zero if ! called. Use RCOND in CGBCO for a reliable indication of singularity. ! implicit none integer lda integer n complex abd(lda,n) real cabs1 integer i integer i0 integer icamax integer info integer ipvt(n) integer j integer j0 integer j1 integer ju integer jz integer k integer l integer lm integer m integer ml integer mm integer mu complex t m = ml + mu + 1 info = 0 ! ! Zero initial fill-in columns. ! j0 = mu + 2 j1 = min ( n, m ) - 1 do jz = j0, j1 i0 = m + 1 - jz do i = i0, ml abd(i,jz) = cmplx ( 0.0E+00, 0.0E+00 ) end do end do jz = j1 ju = 0 ! ! Gaussian elimination with partial pivoting. ! do k = 1, n-1 ! ! Zero next fill-in column ! jz = jz + 1 if ( jz <= n ) then abd(1:ml,jz) = cmplx ( 0.0E+00, 0.0E+00 ) end if ! ! Find L = pivot index. ! lm = min ( ml, n - k ) l = icamax ( lm+1, abd(m:m+lm,k), 1 ) + m - 1 ipvt(k) = l + k - m ! ! Zero pivot implies this column already triangularized. ! if ( cabs1 ( abd(l,k) ) == 0.0E+00 ) then info = k cycle end if ! ! Interchange if necessary. ! if ( l /= m ) then t = abd(l,k) abd(l,k) = abd(m,k) abd(m,k) = t end if ! ! Compute multipliers. ! t = - cmplx ( 1.0E+00, 0.0E+00 ) / abd(m,k) abd(m+1:m+lm,k) = abd(m+1:m+lm,k) * t ! ! Row elimination with column indexing. ! ju = min ( max ( ju, mu + ipvt(k) ), n ) mm = m do j = k+1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if ( l /= mm ) then abd(l,j) = abd(mm,j) abd(mm,j) = t end if abd(mm+1:mm+lm,j) = abd(mm+1:mm+lm,j) + t * abd(m+1:m+lm,k) end do end do ipvt(n) = n if ( cabs1 ( abd(m,n) ) == 0.0E+00 ) then info = n end if return end subroutine cgbsl ( abd, lda, n, ml, mu, ipvt, b, job ) !*****************************************************************************80 ! !! CGBSL solves a complex band system factored by CGBCO or CGBFA. ! ! Discussion: ! ! CGBSL can solve A * X = B or hermitan ( A ) * X = B. ! ! A division by zero will occur if the input factor contains a ! zero on the diagonal. Technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of LDA. It will not occur if the subroutines are ! called correctly and if CGBCO has set 0.0 < RCOND ! or CGBFA has set INFO = 0. ! ! To compute inverse ( A ) * C where C is a matrix with P columns: ! ! call cgbco(abd,lda,n,ml,mu,ipvt,rcond,z) ! ! if ( rcond is not too small ) then ! do j = 1, p ! call cgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) ! end do ! end if ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex ABD(LDA,N), the output from CGBCO or CGBFA. ! ! Input, integer LDA, the leading dimension of ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, the number of diagonals below the main diagonal. ! ! Input, integer MU, the number of diagonals above the main diagonal. ! ! Input, integer IPVT(N), the pivot vector from CGBCO or CGBFA. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! ! Input, integer JOB. ! 0, to solve A*x = b, ! nonzero, to solve hermitian(A)*x = b, where hermitian(A) is the ! conjugate transpose. ! implicit none integer lda integer n complex abd(lda,n) complex b(n) complex cdotc integer ipvt(n) integer job integer k integer l integer la integer lb integer lm integer m integer ml integer mu complex t m = mu + ml + 1 if ( job == 0 ) then ! ! JOB = 0, solve A * X = B. ! ! First solve L * Y = B. ! if ( ml /= 0 ) then do k = 1, n - 1 lm = min ( ml, n - k ) l = ipvt(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if b(k+1:k+lm) = b(k+1:k+lm) + t * abd(m+1:m+lm,k) end do end if ! ! Now solve U * X = Y. ! do k = n, 1, -1 b(k) = b(k) / abd(m,k) lm = min ( k, m ) - 1 la = m - lm lb = k - lm t = -b(k) b(lb:lb+lm-1) = b(lb:lb+lm-1) + t * abd(la:la+lm-1,k) end do else ! ! JOB = nonzero, solve hermitian(A) * X = B. ! ! First solve hermitian ( U ) * Y = B. ! do k = 1, n lm = min ( k, m ) - 1 la = m - lm lb = k - lm ! t = cdotc ( lm, abd(la:la+lm-1,k), 1, b(lb:lb+lm-1), 1 ) t = dot_product ( abd(la:la+lm-1,k), b(lb:lb+lm-1) ) b(k) = ( b(k) - t ) / conjg ( abd(m,k) ) end do ! ! Now solve hermitian ( L ) * X = Y. ! if ( ml /= 0 ) then do k = n-1, 1, -1 lm = min ( ml, n - k ) ! b(k) = b(k) + cdotc ( lm, abd(m+1:m+lm,k), 1, b(k+1:k+lm), 1 ) b(k) = b(k) + dot_product ( abd(m+1:m+lm,k), b(k+1:k+lm) ) l = ipvt(k) if ( l /= k ) then t = b(l) b(l) = b(k) b(k) = t end if end do end if end if return end subroutine cgeco ( a, lda, n, ipvt, rcond, z ) !*****************************************************************************80 ! !! CGECO factors a complex matrix and estimates its condition. ! ! Discussion: ! ! If RCOND is not needed, CGEFA is slightly faster. ! ! To solve A*X = B, follow CGECO by CGESL. ! ! To compute inverse(A)*C, follow CGECO by CGESL. ! ! To compute determinant(A), follow CGECO by CGEDI. ! ! To compute inverse(A), follow CGECO by CGEDI. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complexA(LDA,N), on input, the matrix to be factored. ! On output, an upper triangular matrix and the multipliers which were ! used to obtain it. The factorization can be written A = L*U where ! L is a product of permutation and unit lower triangular matrices ! and U is upper triangular. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of A. ! For the system A*X = B, relative perturbations in A and B of size ! EPSILON may cause relative perturbations in X of size (EPSILON/RCOND). ! If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate ! underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is close to a singular matrix, then Z is ! an approximate null vector in the sense that ! norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ). ! implicit none integer lda integer n complex a(lda,n) real anorm real cabs1 complex cdotc complex csign1 complex ek integer info integer ipvt(n) integer j integer k integer l real rcond real s real scasum real sm complex t complex wk complex wkm real ynorm complex z(n) ! ! Compute the 1-norm of A. ! anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, scasum ( n, a(1:n,j), 1 ) ) end do ! ! Factor. ! call cgefa ( a, lda, n, ipvt, info ) ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and hermitian(A)*Y = E. ! ! Hermitian(A) is the conjugate transpose of A. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where hermitian(U)*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve hermitian(U)*W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) do k = 1, n if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csign1 ( ek, -z(k) ) end if if ( cabs1 ( a(k,k) ) < cabs1 ( ek - z(k) ) ) then s = cabs1 ( a(k,k) ) / cabs1 ( ek - z(k) ) z(1:n) = z(1:n) * s ek = cmplx ( s, 0.0E+00 ) * ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1 ( wk ) sm = cabs1 ( wkm ) if ( cabs1 ( a(k,k) ) /= 0.0E+00 ) then wk = wk / conjg ( a(k,k) ) wkm = wkm / conjg ( a(k,k) ) else wk = cmplx ( 1.0E+00, 0.0E+00 ) wkm = cmplx ( 1.0E+00, 0.0E+00 ) end if do j = k+1, n sm = sm + cabs1 ( z(j) + wkm * conjg ( a(k,j) ) ) z(j) = z(j) + wk * conjg ( a(k,j) ) s = s + cabs1 ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm z(k+1:n) = z(k+1:n) + t * conjg ( a(k,k+1:n) ) end if z(k) = wk end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ! ! Solve hermitian(L) * Y = W. ! do k = n, 1, -1 if ( k < n ) then ! z(k) = z(k) + cdotc ( n-k, a(k+1:n,k), 1, z(k+1:n), 1 ) z(k) = z(k) + dot_product ( a(k+1:n,k), z(k+1:n) ) end if if ( 1.0E+00 < cabs1 ( z(k) ) ) then s = 1.0E+00 / cabs1 ( z(k) ) z(1:n) = z(1:n) * s end if l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = 1.0E+00 ! ! Solve L * V = Y. ! do k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t if ( k < n ) then z(k+1:n) = z(k+1:n) + t * a(k+1:n,k) end if if ( 1.0E+00 < cabs1 ( z(k) ) ) then s = 1.0E+00 / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm ! ! Solve U * Z = V. ! do k = n, 1, -1 if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) ) then s = cabs1 ( a(k,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if if ( cabs1 ( a(k,k) ) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = cmplx ( 1.0E+00, 0.0E+00 ) end if t = -z(k) z(1:k-1) = z(1:k-1) + t * a(1:k-1,k) end do ! ! Make ZNORM = 1. ! s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine cgedi ( a, lda, n, ipvt, det, work, job ) !*****************************************************************************80 ! !! CGEDI computes the determinant and inverse of a matrix. ! ! Discussion: ! ! The matrix must have been factored by CGECO or CGEFA. ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal and the inverse is requested. ! It will not occur if the subroutines are called correctly ! and if CGECO has set 0.0 < RCOND or CGEFA has set ! INFO == 0. ! ! Modified: ! ! 29 April 2007 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the factor information ! from CGECO or CGEFA. On output, the inverse matrix, if it ! was requested, ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CGECO or CGEFA. ! ! Output, complex DET(2), the determinant of the original matrix, ! if requested. Otherwise not referenced. ! Determinant = DET(1) * 10.0**DET(2) with ! 1.0 <= cabs1 ( DET(1) ) < 10.0 or DET(1) == 0.0. ! Also, DET(2) is strictly real. ! ! Workspace, complex WORK(N). ! ! Input, integer JOB. ! 11, both determinant and inverse. ! 01, inverse only. ! 10, determinant only. ! implicit none integer lda integer n complex a(lda,n) real cabs1 complex det(2) integer i integer ipvt(n) integer j integer job integer k integer l complex t complex work(n) ! ! Compute the determinant. ! if ( job / 10 /= 0 ) then det(1) = cmplx ( 1.0E+00, 0.0E+00 ) det(2) = cmplx ( 0.0E+00, 0.0E+00 ) do i = 1, n if ( ipvt(i) /= i ) then det(1) = -det(1) end if det(1) = det(1) * a(i,i) if ( cabs1 ( det(1) ) == 0.0E+00 ) then exit end if do while ( cabs1 ( det(1) ) < 1.0E+00 ) det(1) = det(1) * cmplx ( 10.0E+00, 0.0E+00 ) det(2) = det(2) - cmplx ( 1.0E+00, 0.0E+00 ) end do do while ( 10.0E+00 <= cabs1 ( det(1) ) ) det(1) = det(1) / cmplx ( 10.0E+00, 0.0E+00 ) det(2) = det(2) + cmplx ( 1.0E+00, 0.0E+00 ) end do end do end if ! ! Compute inverse(U). ! if ( mod ( job, 10 ) /= 0 ) then do k = 1, n a(k,k) = cmplx ( 1.0E+00, 0.0E+00 ) / a(k,k) t = -a(k,k) a(1:k-1,k) = a(1:k-1,k) * t do j = k+1, n t = a(k,j) a(k,j) = cmplx ( 0.0E+00, 0.0E+00 ) a(1:k,j) = a(1:k,j) + t * a(1:k,k) end do end do ! ! Form inverse(U) * inverse(L). ! do k = n-1, 1, -1 work(k+1:n) = a(k+1:n,k) a(k+1:n,k) = cmplx ( 0.0E+00, 0.0E+00 ) do j = k+1, n t = work(j) a(1:n,k) = a(1:n,k) + t * a(1:n,j) end do l = ipvt(k) if ( l /= k ) then work(1:n) = a(1:n,k) a(1:n,k) = a(1:n,l) a(1:n,l) = work(1:n) end if end do end if return end subroutine cgefa ( a, lda, n, ipvt, info ) !*****************************************************************************80 ! !! CGEFA factors a complex matrix by Gaussian elimination. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the matrix to be factored. ! On output, an upper triangular matrix and the multipliers which were ! used to obtain it. The factorization can be written A = L*U where ! L is a product of permutation and unit lower triangular matrices and ! U is upper triangular. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO, ! 0, normal value. ! K, if U(K,K) == 0.0. This is not an error condition for this ! subroutine, but it does indicate that CGESL or CGEDI will divide by zero ! if called. Use RCOND in CGECO for a reliable indication of singularity. ! implicit none integer lda integer n complex a(lda,n) real cabs1 integer icamax integer info integer ipvt(n) integer j integer k integer l complex t ! ! Gaussian elimination with partial pivoting. ! info = 0 do k = 1, n - 1 ! ! Find L = pivot index. ! l = icamax ( n-k+1, a(k:n,k), 1 ) + k - 1 ipvt(k) = l ! ! Zero pivot implies this column already triangularized. ! if ( cabs1 ( a(l,k) ) == 0.0E+00 ) then info = k cycle end if ! ! Interchange if necessary. ! if ( l /= k ) then t = a(l,k) a(l,k) = a(k,k) a(k,k) = t end if ! ! Compute multipliers ! t = - cmplx ( 1.0E+00, 0.0E+00 ) / a(k,k) a(k+1:n,k) = a(k+1:n,k) * t ! ! Row elimination with column indexing ! do j = k+1, n t = a(l,j) if ( l /= k ) then a(l,j) = a(k,j) a(k,j) = t end if a(k+1:n,j) = a(k+1:n,j) + t * a(k+1:n,k) end do end do ipvt(n) = n if ( cabs1 ( a(n,n) ) == 0.0E+00 ) then info = n end if return end subroutine cgesl ( a, lda, n, ipvt, b, job ) !*****************************************************************************80 ! !! CGESL solves a complex system factored by CGECO or CGEFA. ! ! Discussion: ! ! A division by zero will occur if the input factor contains a ! zero on the diagonal. Technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of LDA. It will not occur if the subroutines are ! called correctly and if CGECO has set 0.0 < RCOND ! or CGEFA has set INFO == 0. ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call cgeco(a,lda,n,ipvt,rcond,z) ! ! if (rcond is not too small) then ! do j = 1, p ! call cgesl ( a, lda, n, ipvt, c(1,j), 0 ) ! end do ! end if ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex A(LDA,N), the factored matrix information, ! as output from CGECO or CGEFA. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CGECO or CGEFA. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! ! Input, integer JOB. ! 0, to solve A*X = B, ! nonzero, to solve hermitian(A)*X = B where hermitian(A) is the ! conjugate transpose. ! implicit none integer lda integer n complex a(lda,n) complex b(n) complex cdotc integer ipvt(n) integer job integer k integer l complex t if ( job == 0 ) then ! ! JOB = 0, solve A * X = B. ! ! First solve L * Y = B. ! do k = 1, n-1 l = ipvt(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if b(k+1:n) = b(k+1:n) + t * a(k+1:n,k) end do ! ! Now solve U * X = Y. ! do k = n, 1, -1 b(k) = b(k) / a(k,k) t = -b(k) b(1:k-1) = b(1:k-1) + t * a(1:k-1,k) end do else ! ! JOB nonzero, solve hermitian(A) * X = B. ! ! First solve hermitian(U) * Y = B. ! do k = 1, n ! t = cdotc ( k-1, a(1:k-1,k), 1, b(1:k-1), 1 ) t = dot_product ( a(1:k-1,k), b(1:k-1) ) b(k) = ( b(k) - t ) / conjg ( a(k,k) ) end do ! ! Now solve hermitian(L) * X = Y. ! do k = n-1, 1, -1 ! b(k) = b(k) + cdotc ( n-k, a(k+1:n,k), 1, b(k+1:n), 1 ) b(k) = b(k) + dot_product ( a(k+1:n,k), b(k+1:n) ) l = ipvt(k) if ( l /= k ) then t = b(l) b(l) = b(k) b(k) = t end if end do end if return end subroutine cgtsl ( n, c, d, e, b, info ) !*****************************************************************************80 ! !! CGTSL solves a complex general tridiagonal system. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input/output, complex C(N); on input, the subdiagonal of the tridiagonal ! matrix in entries C(2:N). On output, C has been overwritten. ! ! Input/output, complex D(N); on input, the diagonal of the tridiagonal ! matrix. On output, D has been overwritten. ! ! Input/output, complex E(N); on input, the superdiagonal of the tridiagonal ! matrix in entries E(1:N-1). On output, E has been overwritten. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! ! Output, integer INFO. ! 0, normal value. ! K, if the K-th element of the diagonal becomes exactly zero. The ! subroutine returns when this is detected. ! implicit none integer n complex b(n) complex c(n) real cabs1 complex d(n) complex e(n) integer info integer k complex t info = 0 c(1) = d(1) if ( 1 <= n-1 ) then d(1) = e(1) e(1) = cmplx ( 0.0E+00, 0.0E+00 ) e(n) = cmplx ( 0.0E+00, 0.0E+00 ) do k = 1, n-1 if ( cabs1 ( c(k) ) <= cabs1 ( c(k+1) ) ) then t = b(k) b(k) = b(k+1) b(k+1) = t t = c(k) c(k) = c(k+1) c(k+1) = t t = d(k) d(k) = d(k+1) d(k+1) = t t = e(k) e(k) = e(k+1) e(k+1) = t end if if ( cabs1 ( c(k) ) == 0.0E+00 ) then info = k return end if t = -c(k+1) / c(k) c(k+1) = d(k+1) + t * d(k) d(k+1) = e(k+1) + t * e(k) e(k+1) = cmplx ( 0.0E+00, 0.0E+00 ) b(k+1) = b(k+1) + t * b(k) end do end if if ( cabs1 ( c(n) ) == 0.0E+00 ) then info = n return end if ! ! Back solve. ! b(n) = b(n) / c(n) if ( 1 < n ) then b(n-1) = ( b(n-1) - d(n-1) * b(n) ) / c(n-1) do k = n-2, 1, -1 b(k) = ( b(k) - d(k) * b(k+1) - e(k) * b(k+2) ) / c(k) end do end if return end subroutine chico ( a, lda, n, ipvt, rcond, z ) !*****************************************************************************80 ! !! CHICO factors a complex hermitian matrix and estimates its condition. ! ! Discussion: ! ! If RCOND is not needed, CHIFA is slightly faster. ! ! To solve A*X = B, follow CHICO by CHISL. ! ! To compute inverse(A)*C, follow CHICO by CHISL. ! ! To compute inverse(A), follow CHICO by CHIDI. ! ! To compute determinant(A), follow CHICO by CHIDI. ! ! To compute inertia(A), follow CHICO by CHIDI. ! ! Modified: ! ! 01 May 2007 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the hermitian matrix to be ! factored. On output, a block diagonal matrix and the multipliers ! which were used to obtain it. The factorization can be written ! A = U*D*hermitian(U) where U is a product of permutation and unit ! upper triangular matrices, hermitian(U) is the conjugate transpose ! of U, and D is block diagonal with 1 by 1 and 2 by 2 blocks. ! Only the diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of ! the matrix. For the system A*X = B, relative perturbations in A and B ! of size EPSILON may cause relative perturbations in X of size ! (EPSILON/RCOND). If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is close to a singular matrix, then Z is an ! approximate null vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! implicit none integer lda integer n complex a(lda,n) complex ak complex akm1 real anorm complex bk complex bkm1 real cabs1 complex cdotc complex csign1 complex denom complex ek integer i integer info integer ipvt(n) integer j integer k integer kp integer kps integer ks real rcond real s real scasum complex t real ynorm complex z(n) ! ! Find norm of A using only upper half. ! do j = 1, n z(j) = cmplx ( scasum ( j, a(1:j,j), 1 ), 0.0E+00 ) do i = 1, j - 1 z(i) = cmplx ( real ( z(i) ) + cabs1 ( a(i,j) ), 0.0E+00 ) end do end do anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, real ( z(j) ) ) end do ! ! Factor. ! call chifa ( a, lda, n, ipvt, info ) ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where U*D*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve U*D*W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) k = n do while ( 0 < k ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if kp = abs ( ipvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then t = z(kps) z(kps) = z(kp) z(kp) = t end if if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csign1 ( ek, z(k) ) end if z(k) = z(k) + ek z(1:k-ks) = z(1:k-ks) + z(k) * a(1:k-ks,k) if ( ks /= 1 ) then if ( cabs1 ( z(k-1) ) /= 0.0E+00 ) then ek = csign1 ( ek, z(k-1) ) end if z(k-1) = z(k-1) + ek z(1:k-ks) = z(1:k-ks) + z(k-1) * a(1:k-ks,k-1) end if if ( ks /= 2 ) then if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) ) then s = cabs1 ( a(k,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ek = cmplx ( s, 0.0E+00 ) * ek end if if ( cabs1 ( a(k,k) ) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = cmplx ( 1.0E+00, 0.0E+00 ) end if else ak = a(k,k) / conjg ( a(k-1,k) ) akm1 = a(k-1,k-1) / a(k-1,k) bk = z(k) / conjg ( a(k-1,k) ) bkm1 = z(k-1) / a(k-1,k) denom = ak * akm1 - cmplx ( 1.0E+00, 0.0E+00 ) z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ! ! Solve hermitian(U) * Y = W. ! k = 1 do while ( k <= n ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= 1 ) then ! z(k) = z(k) + cdotc ( k-1, a(1:k-1,k), 1, z(1:k-1), 1 ) z(k) = z(k) + dot_product ( a(1:k-1,k), z(1:k-1) ) if ( ks == 2 ) then ! z(k+1) = z(k+1) + cdotc ( k-1, a(1:k-1,k+1), 1, z(1:k-1), 1 ) z(k+1) = z(k+1) + dot_product ( a(1:k-1,k+1), z(1:k-1) ) end if kp = abs ( ipvt(k) ) if ( kp /= k ) then t = z(k) z(k) = z(kp) z(kp) = t end if end if k = k + ks end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = 1.0E+00 ! ! Solve U*D*V = Y. ! k = n do while ( 0 < k ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= ks ) then kp = abs ( ipvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then t = z(kps) z(kps) = z(kp) z(kp) = t end if z(1:k-ks) = z(1:k-ks) + z(k) * a(1:k-ks,k) if ( ks == 2 ) then z(1:k-ks) = z(1:k-ks) + z(k-1) * a(1:k-ks,k-1) end if end if if ( ks /= 2 ) then if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) ) then s = cabs1 ( a(k,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if if ( cabs1 ( a(k,k) ) /= 0.0E+00 ) then z(k) = z(k) / a(k,k) else z(k) = cmplx ( 1.0E+00, 0.0E+00 ) end if else ak = a(k,k) / conjg ( a(k-1,k) ) akm1 = a(k-1,k-1) / a(k-1,k) bk = z(k) / conjg ( a(k-1,k) ) bkm1 = z(k-1) / a(k-1,k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm ! ! Solve hermitian(U) * Z = V. ! k = 1 do while ( k <= n ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= 1 ) then ! z(k) = z(k) + cdotc ( k-1, a(1:k-1,k), 1, z(1:k-1), 1 ) z(k) = z(k) + dot_product ( a(1:k-1,k), z(1:k-1) ) if ( ks == 2 ) then ! z(k+1) = z(k+1) + cdotc ( k-1, a(1:k-1,k+1), 1, z(1:k-1), 1 ) z(k+1) = z(k+1) + dot_product ( a(1:k-1,k+1), z(1:k-1) ) end if kp = abs ( ipvt(k) ) if ( kp /= k ) then t = z(k) z(k) = z(kp) z(kp) = t end if end if k = k + ks end do ! ! Make ZNORM = 1. ! s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine chidi ( a, lda, n, ipvt, det, inert, work, job ) !*****************************************************************************80 ! !! CHIDI computes the determinant and inverse of a matrix factored by CHIFA. ! ! Discussion: ! ! CHIDI computes the determinant, inertia (number of positive, zero, ! and negative eigenvalues) and inverse of a complex hermitian matrix ! using the factors from CHIFA. ! ! A division by zero may occur if the inverse is requested ! and CHICO has set RCOND == 0.0 or CHIFA has set INFO /= 0. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the factored matrix ! from CHIFA. On output, if the inverse was requested, A contains ! the inverse matrix. The strict lower triangle of A is never ! referenced. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CHIFA. ! ! Output, real DET(2), the determinant of the original matrix. ! Determinant = DET(1) * 10.0**DET(2) with 1.0 <= abs ( DET(1) ) < 10.0 ! or DET(1) = 0.0. ! ! Output, integer INERT(3), the inertia of the original matrix. ! INERT(1) = number of positive eigenvalues. ! INERT(2) = number of negative eigenvalues. ! INERT(3) = number of zero eigenvalues. ! ! Workspace, complex WORK(N), ! ! Input, integer JOB, has the decimal expansion ABC where: ! if C /= 0, the inverse is computed, ! if B /= 0, the determinant is computed, ! if A /= 0, the inertia is computed. ! For example, JOB = 111 gives all three. ! implicit none integer lda integer n complex a(lda,n) real ak complex akkp1 real akp1 complex cdotc real d real det(2) integer job integer inert(3) integer ipvt(n) integer j integer jb integer k integer km1 integer ks integer kstep logical nodet logical noert logical noinv real t complex temp complex work(n) noinv = mod ( job, 10 ) == 0 nodet = mod ( job, 100 ) / 10 == 0 noert = mod ( job, 1000 ) / 100 == 0 if ( .not. nodet .or. .not. noert ) then if ( .not. noert ) then inert(1:3) = 0 end if if ( .not. nodet ) then det(1) = 1.0E+00 det(2) = 0.0E+00 end if t = 0.0E+00 do k = 1, n d = real ( a(k,k) ) ! ! Check if 1 by 1. ! if ( ipvt(k) <= 0 ) then ! ! 2 by 2 block ! Use DET = ( D / T * C - T ) * T, T = abs ( S ) ! to avoid underflow/overflow troubles. ! Take two passes through scaling. Use T for flag. ! if ( t == 0.0E+00 ) then t = abs ( a(k,k+1) ) d = ( d / t ) * real ( a(k+1,k+1) ) - t else d = t t = 0.0E+00 end if end if if ( .not. noert ) then if ( 0.0E+00 < d ) then inert(1) = inert(1) + 1 else if ( d < 0.0E+00 ) then inert(2) = inert(2) + 1 else if ( d == 0.0E+00 ) then inert(3) = inert(3) + 1 end if end if if ( .not. nodet ) then det(1) = det(1) * d if ( det(1) /= 0.0E+00 ) then do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = det(1) * 10.0E+00 det(2) = det(2) - 1.0E+00 end do do while ( 10.0E+00 <= abs ( det(1) ) ) det(1) = det(1) / 10.0E+00 det(2) = det(2) + 1.0E+00 end do end if end if end do end if ! ! Compute inverse(A). ! if ( .not. noinv ) then k = 1 do while ( k <= n ) km1 = k - 1 if ( 0 <= ipvt(k) ) then ! ! 1 by 1 ! a(k,k) = cmplx ( 1.0E+00 / real ( a(k,k) ), 0.0E+00 ) if ( 1 <= km1 ) then work(1:km1) = a(1:km1,k) do j = 1, km1 a(j,k) = cdotc ( j, a(1:j,j), 1, work, 1 ) a(1:j-1,k) = a(1:j-1,k) + work(j) * a(1:j-1,j) end do a(k,k) = a(k,k) + cmplx ( & real ( cdotc ( km1, work, 1, a(1:km1,k), 1 ) ), 0.0E+00 ) end if kstep = 1 else ! ! 2 by 2 ! t = abs ( a(k,k+1) ) ak = real ( a(k,k) ) / t akp1 = real ( a(k+1,k+1) ) / t akkp1 = a(k,k+1) / t d = t * ( ak * akp1 - 1.0E+00 ) a(k,k) = cmplx ( akp1 / d, 0.0E+00 ) a(k+1,k+1) = cmplx ( ak / d, 0.0E+00 ) a(k,k+1) = -akkp1 / d if ( 1 <= km1 ) then work(1:km1) = a(1:km1,k+1) do j = 1, km1 a(j,k+1) = cdotc ( j, a(1:j,j), 1, work, 1 ) a(1:j-1,k+1) = a(1:j-1,k+1) + work(j) * a(1:j-1,j) end do a(k+1,k+1) = a(k+1,k+1) + cmplx ( & real ( cdotc ( km1, work, 1, a(1:km1,k+1), 1 ) ), 0.0E+00 ) a(k,k+1) = a(k,k+1) + cdotc ( km1, a(1:km1,k), 1, a(1:km1,k+1), 1 ) work(1:km1) = a(1:km1,k) do j = 1, km1 a(j,k) = cdotc ( j, a(1:j,j), 1, work, 1 ) a(1:j-1,k) = a(1:j-1,k) + work(j) * a(1:j-1,j) end do a(k,k) = a(k,k) + cmplx ( & real ( cdotc ( km1, work, 1, a(1:km1,k), 1 ) ), 0.0E+00 ) end if kstep = 2 end if ! ! Swap ! ks = abs ( ipvt(k) ) if ( ks /= k ) then call cswap ( ks, a(1:ks,ks), 1, a(1:ks,k), 1 ) do j = k, ks, -1 temp = conjg ( a(j,k) ) a(j,k) = conjg ( a(ks,j) ) a(ks,j) = temp end do if ( kstep /= 1 ) then temp = a(ks,k+1) a(ks,k+1) = a(k,k+1) a(k,k+1) = temp end if end if k = k + kstep end do end if return end subroutine chifa ( a, lda, n, ipvt, info ) !*****************************************************************************80 ! !! CHIFA factors a complex hermitian matrix. ! ! Discussion: ! ! CHIFA performs the factoring by elimination with symmetric pivoting. ! ! To solve A*X = B, follow CHIFA by CHISL. ! ! To compute inverse(A)*C, follow CHIFA by CHISL. ! ! To compute determinant(A), follow CHIFA by CHIDI. ! ! To compute inertia(A), follow CHIFA by CHIDI. ! ! To compute inverse(A), follow CHIFA by CHIDI. ! ! Modified: ! ! 01 May 2007 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the hermitian matrix to be ! factored. On output, a block diagonal matrix and the multipliers which ! were used to obtain it. The factorization can be written ! A = U*D*hermitian(U) where U is a product of permutation and unit upper ! triangular matrices, hermitian(U) is the conjugate transpose of U, and ! D is block diagonal with 1 by 1 and 2 by 2 blocks. Only the diagonal ! and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO. ! 0, normal value. ! K, if the K-th pivot block is singular. This is not an error condition ! for this subroutine, but it does indicate that CHISL or CHIDI may ! divide by zero if called. ! implicit none integer lda integer n complex a(lda,n) real absakk complex ak complex akm1 real alpha complex bk complex bkm1 real cabs1 real colmax complex denom integer icamax integer imax integer info integer ipvt(n) integer j integer jj integer jmax integer k integer km1 integer km2 integer kstep complex mulk complex mulkm1 real rowmax logical swap complex t ! ! Initialize. ! ! ALPHA is used in choosing pivot block size. ! alpha = ( 1.0E+00 + sqrt ( 17.0E+00 ) ) / 8.0E+00 info = 0 ! ! Main loop on K, which goes from N to 1. ! k = n do ! ! Leave the loop if K = 0 or K = 1. ! if ( k == 0 ) then exit end if if ( k == 1 ) then ipvt(1) = 1 if ( cabs1 ( a(1,1) ) == 0.0E+00 ) then info = 1 end if exit end if ! ! This section of code determines the kind of ! elimination to be performed. When it is completed, ! KSTEP will be set to the size of the pivot block, and ! SWAP will be set to TRUE if an interchange is ! required. ! km1 = k - 1 absakk = cabs1 ( a(k,k) ) ! ! Determine the largest off-diagonal element in column K. ! imax = icamax ( k-1, a(1:k-1,k), 1 ) colmax = cabs1 ( a(imax,k) ) if ( alpha * colmax <= absakk ) then kstep = 1 swap = .false. else ! ! Determine the largest off-diagonal element in row IMAX. ! rowmax = 0.0E+00 do j = imax + 1, k rowmax = max ( rowmax, cabs1 ( a(imax,j) ) ) end do if ( imax /= 1 ) then jmax = icamax ( imax-1, a(1:imax-1,imax), 1 ) rowmax = max ( rowmax, cabs1 ( a(jmax,imax) ) ) end if if ( alpha * rowmax <= cabs1 ( a(imax,imax) ) ) then kstep = 1 swap = .true. else if ( alpha * colmax * ( colmax / rowmax ) <= absakk ) then kstep = 1 swap = .false. else kstep = 2 swap = ( imax /= km1 ) end if end if ! ! Column K is zero. Set INFO and iterate the loop. ! if ( max ( absakk, colmax ) == 0.0E+00 ) then ipvt(k) = k info = k k = k - kstep cycle end if if ( kstep /= 2 ) then ! ! 1 x 1 pivot block. ! if ( swap ) then call cswap ( imax, a(1:imax,imax), 1, a(1:imax,k), 1 ) do jj = imax, k j = k + imax - jj t = conjg ( a(j,k) ) a(j,k) = conjg ( a(imax,j) ) a(imax,j) = t end do end if ! ! Perform the elimination. ! do jj = 1, km1 j = k - jj mulk = -a(j,k) / a(k,k) t = conjg ( mulk ) a(1:j,j) = a(1:j,j) + t * a(1:j,k) a(j,j) = cmplx ( real ( a(j,j) ), 0.0E+00 ) a(j,k) = mulk end do ! ! Set the pivot array. ! ipvt(k) = k if ( swap ) then ipvt(k) = imax end if else ! ! 2 x 2 pivot block. ! if ( swap ) then call cswap ( imax, a(1:imax,imax), 1, a(1:imax,k-1), 1 ) do jj = imax, km1 j = km1 + imax - jj t = conjg ( a(j,k-1) ) a(j,k-1) = conjg ( a(imax,j) ) a(imax,j) = t end do t = a(k-1,k) a(k-1,k) = a(imax,k) a(imax,k) = t end if ! ! Perform the elimination. ! km2 = k - 2 if ( 0 < k - 2 ) then ak = a(k,k) / a(k-1,k) akm1 = a(k-1,k-1) / conjg ( a(k-1,k) ) denom = 1.0E+00 - ak * akm1 do jj = 1, k - 2 j = km1 - jj bk = a(j,k) / a(k-1,k) bkm1 = a(j,k-1) / conjg ( a(k-1,k) ) mulk = ( akm1 * bk - bkm1 ) / denom mulkm1 = ( ak * bkm1 - bk ) / denom t = conjg ( mulk ) a(1:j,j) = a(1:j,j) + t * a(1:j,k) t = conjg ( mulkm1 ) a(1:j,j) = a(1:j,j) + t * a(1:j,k-1) a(j,k) = mulk a(j,k-1) = mulkm1 a(j,j) = cmplx ( real ( a(j,j) ), 0.0E+00 ) end do end if ! ! Set the pivot array. ! if ( swap ) then ipvt(k) = -imax else ipvt(k) = 1 - k end if ipvt(k-1) = ipvt(k) end if k = k - kstep end do return end subroutine chisl ( a, lda, n, ipvt, b ) !*****************************************************************************80 ! !! CHISL solves a complex hermitian system factored by CHIFA. ! ! Discussion: ! ! A division by zero may occur if CHICO has set RCOND == 0.0 ! or CHIFA has set INFO /= 0. ! ! To compute inverse(A) * C where C is a matrix with P columns ! ! call chifa(a,lda,n,ipvt,info) ! if ( info == 0 ) then ! do j = 1, p ! call chisl(a,lda,n,ipvt,c(1,j)) ! end do ! end if ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex A(LDA,N), the output from CHIFA. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CHIFA. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! implicit none integer lda integer n complex a(lda,n) complex ak complex akm1 complex b(n) complex bk complex bkm1 complex cdotc complex denom integer ipvt(n) integer k integer kp complex t ! ! Loop backward applying the transformations and D inverse to B. ! k = n do while ( 0 < k ) ! ! 1 x 1 pivot block. ! if ( 0 <= ipvt(k) ) then if ( k /= 1 ) then kp = ipvt(k) if ( kp /= k ) then t = b(k) b(k) = b(kp) b(kp) = t end if b(1:k-1) = b(1:k-1) + a(1:k-1,k) * b(k) end if ! ! Apply D inverse. ! b(k) = b(k) / a(k,k) k = k - 1 ! ! 2 x 2 pivot block. ! else if ( k /= 2 ) then kp = abs ( ipvt(k) ) if ( kp /= k - 1 ) then t = b(k-1) b(k-1) = b(kp) b(kp) = t end if b(1:k-2) = b(1:k-2) + b(k) * a(1:k-2,k) b(1:k-2) = b(1:k-2) + b(k-1) * a(1:k-2,k-1) end if ! ! Apply D inverse. ! ak = a(k,k) / conjg ( a(k-1,k) ) akm1 = a(k-1,k-1) / a(k-1,k) bk = b(k) / conjg ( a(k-1,k) ) bkm1 = b(k-1) / a(k-1,k) denom = ak * akm1 - 1.0E+00 b(k) = ( akm1 * bk - bkm1 ) / denom b(k-1) = ( ak * bkm1 - bk ) / denom k = k - 2 end if end do ! ! Loop forward applying the transformations. ! k = 1 do while ( k <= n ) ! ! 1 x 1 pivot block. ! if ( 0 <= ipvt(k) ) then if ( k /= 1 ) then b(k) = b(k) + cdotc ( k-1, a(1:k-1,k), 1, b(1:k-1), 1 ) kp = ipvt(k) if ( kp /= k ) then t = b(k) b(k) = b(kp) b(kp) = t end if end if k = k + 1 else ! ! 2 x 2 pivot block. ! if ( k /= 1 ) then b(k) = b(k) + cdotc ( k-1, a(1:k-1,k), 1, b(1:k-1), 1 ) b(k+1) = b(k+1) + cdotc ( k-1, a(1:k-1,k+1), 1, b(1:k-1), 1 ) kp = abs ( ipvt(k) ) if ( kp /= k ) then t = b(k) b(k) = b(kp) b(kp) = t end if end if k = k + 2 end if end do return end subroutine chpco ( ap, n, ipvt, rcond, z ) !*****************************************************************************80 ! !! CHPCO factors a complex hermitian packed matrix and estimates its condition. ! ! Discussion: ! ! If RCOND is not needed, CHPFA is slightly faster. ! ! To solve A*X = B, follow CHPCO by CHPSL. ! ! To compute inverse(A)*C, follow CHPCO by CHPSL. ! ! To compute inverse(A), follow CHPCO by CHPDI. ! ! To compute determinant(A), follow CHPCO by CHPDI. ! ! To compute inertia(A), follow CHPCO by CHPDI. ! ! Packed storage: ! ! The following program segment will pack the upper ! triangle of a hermitian matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex AP(N*(N+1)/2); on input, the packed form of a ! hermitian matrix A. The columns of the upper triangle are stored ! sequentially in a one-dimensional array of length N*(N+1)/2. On ! output, a block diagonal matrix and the multipliers which were used ! to obtain it stored in packed form. The factorization can be written ! A = U*D*hermitian(U) where U is a product of permutation and unit ! upper triangular matrices, hermitian(U) is the conjugate transpose ! of U, and D is block diagonal with 1 by 1 and 2 by 2 blocks. ! ! Input, integer N, the order of the matrix. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, real RCOND, an estimate of the reciprocal condition of ! the matrix. For the system A*X = B, relative perturbations in A and B ! of size EPSILON may cause relative perturbations in X of size ! (EPSILON/RCOND). If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is close to a singular matrix, then Z is an ! approximate null vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! implicit none integer n complex ak complex akm1 real anorm complex ap((n*(n+1))/2) complex bk complex bkm1 real cabs1 complex cdotc complex csign1 complex denom complex ek integer i integer ij integer ik integer ikm1 integer ikp1 integer info integer ipvt(n) integer j integer j1 integer k integer kk integer km1k integer km1km1 integer kp integer kps integer ks real rcond real s real scasum complex t real ynorm complex z(n) ! ! Find norm of A using only upper half. ! j1 = 1 do j = 1, n z(j) = cmplx ( scasum ( j, ap(j1:j1+j-1), 1 ), 0.0E+00 ) ij = j1 j1 = j1 + j do i = 1, j-1 z(i) = cmplx ( real ( z(i) ) + cabs1 ( ap(ij) ), 0.0E+00 ) ij = ij + 1 end do end do anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, real ( z(j) ) ) end do ! ! Factor. ! call chpfa ( ap, n, ipvt, info ) ! ! RCOND = 1/(norm(A) * (estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where U*D*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve U*D*W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) k = n ik = ( n * ( n - 1 ) ) / 2 do while ( 0 < k ) kk = ik + k ikm1 = ik - ( k - 1 ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if kp = abs ( ipvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then t = z(kps) z(kps) = z(kp) z(kp) = t end if if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csign1 ( ek, z(k) ) end if z(k) = z(k) + ek z(1:k-ks) = z(1:k-ks) + z(k) * ap(ik+1:ik+k-ks) if ( ks /= 1 ) then if ( cabs1 ( z(k-1) ) /= 0.0E+00 ) then ek = csign1 ( ek, z(k-1) ) end if z(k-1) = z(k-1) + ek z(1:k-ks) = z(1:k-ks) + z(k-1) * ap(ikm1+1:ikm1+k-ks) end if if ( ks /= 2 ) then if ( cabs1 ( ap(kk) ) < cabs1 ( z(k) ) ) then s = cabs1 ( ap(kk) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ek = cmplx ( s, 0.0E+00 ) * ek end if if ( cabs1 ( ap(kk) ) /= 0.0E+00 ) then z(k) = z(k) / ap(kk) else z(k) = cmplx ( 1.0E+00, 0.0E+00 ) end if else km1k = ik + k - 1 km1km1 = ikm1 + k - 1 ak = ap(kk) / conjg ( ap(km1k) ) akm1 = ap(km1km1) / ap(km1k) bk = z(k) / conjg ( ap(km1k) ) bkm1 = z(k-1) / ap(km1k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks ik = ik - k if ( ks == 2 ) then ik = ik - ( k + 1 ) end if end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ! ! Solve hermitian(U) * Y = W. ! k = 1 ik = 0 do while ( k <= n ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= 1 ) then z(k) = z(k) + cdotc ( k-1, ap(ik+1:ik+k-1), 1, z(1:k-1), 1 ) ikp1 = ik + k if ( ks == 2 ) then z(k+1) = z(k+1) + cdotc ( k-1, ap(ikp1+1:ikp1+k-1), 1, z(1:k-1), 1 ) end if kp = abs ( ipvt(k) ) if ( kp /= k ) then t = z(k) z(k) = z(kp) z(kp) = t end if end if ik = ik + k if ( ks == 2 ) then ik = ik + ( k + 1 ) end if k = k + ks end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = 1.0E+00 ! ! Solve U*D*V = Y. ! k = n ik = ( n * ( n - 1 ) ) / 2 do while ( 0 < k ) kk = ik + k ikm1 = ik - ( k - 1 ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= ks ) then kp = abs ( ipvt(k) ) kps = k + 1 - ks if ( kp /= kps ) then t = z(kps) z(kps) = z(kp) z(kp) = t end if z(1:k-ks) = z(1:k-ks) + z(k) * ap(ik+1:ik+k-ks) if ( ks == 2 ) then z(1:k-ks) = z(1:k-ks) + z(k-1) * ap(ikm1+1:ikm1+k-ks) end if end if if ( ks /= 2 ) then if ( cabs1 ( ap(kk) ) < cabs1 ( z(k) ) ) then s = cabs1 ( ap(kk) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if if ( cabs1 ( ap(kk) ) /= 0.0E+00 ) then z(k) = z(k) / ap(kk) else z(k) = cmplx ( 1.0E+00, 0.0E+00 ) end if else km1k = ik + k - 1 km1km1 = ikm1 + k - 1 ak = ap(kk) / conjg ( ap(km1k) ) akm1 = ap(km1km1) / ap(km1k) bk = z(k) / conjg ( ap(km1k) ) bkm1 = z(k-1) / ap(km1k) denom = ak * akm1 - 1.0E+00 z(k) = ( akm1 * bk - bkm1 ) / denom z(k-1) = ( ak * bkm1 - bk ) / denom end if k = k - ks ik = ik - k if ( ks == 2 ) then ik = ik - ( k + 1 ) end if end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm ! ! Solve hermitian(U) * Z = V. ! k = 1 ik = 0 do while ( k <= n ) if ( ipvt(k) < 0 ) then ks = 2 else ks = 1 end if if ( k /= 1 ) then z(k) = z(k) + cdotc ( k-1, ap(ik+1:ik+k-1), 1, z(1:k-1), 1 ) ikp1 = ik + k if ( ks == 2 ) then z(k+1) = z(k+1) + cdotc ( k-1, ap(ikp1+1:ikp1+k-1), 1, z(1:k-1), 1 ) end if kp = abs ( ipvt(k) ) if ( kp /= k ) then t = z(k) z(k) = z(kp) z(kp) = t end if end if ik = ik + k if ( ks == 2 ) then ik = ik + ( k + 1 ) end if k = k + ks end do ! ! Make ZNORM = 1. ! s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine chpdi ( ap, n, ipvt, det, inert, work, job ) !*****************************************************************************80 ! !! CHPDI: determinant, inertia and inverse of a complex hermitian matrix. ! ! Discussion: ! ! The routine uses the factors from CHPFA. ! ! The matrix is stored in packed form. ! ! A division by zero will occur if the inverse is requested and CHPCO has ! set RCOND == 0.0 or CHPFA has set INFO /= 0. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex AP(N*(N+1)/2); on input, the factored matrix ! from CHPFA. If the inverse was requested, then on output, AP contains ! the upper triangle of the inverse of the original matrix, stored in packed ! form. The columns of the upper triangle are stored sequentially in a ! one-dimensional array. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CHPFA. ! ! Output, real DET(2), if requested, the determinant of the original matrix. ! Determinant = DET(1) * 10.0**DET(2) with 1.0 <= abs ( DET(1) ) < 10.0 ! or DET(1) = 0.0. ! ! Output, integer INERT(3), if requested, the inertia of the original matrix. ! INERT(1) = number of positive eigenvalues. ! INERT(2) = number of negative eigenvalues. ! INERT(3) = number of zero eigenvalues. ! ! Workspace, complex WORK(N). ! ! Input, integer JOB, has the decimal expansion ABC where: ! if C /= 0, the inverse is computed, ! if B /= 0, the determinant is computed, ! if A /= 0, the inertia is computed. ! For example, JOB = 111 gives all three. ! implicit none integer n real ak complex akkp1 real akp1 complex ap((n*(n+1))/2) complex cdotc real d real det(2) integer ij integer ik integer ikp1 integer iks integer inert(3) integer ipvt(n) integer j integer jb integer jk integer jkp1 integer job integer k integer kk integer kkp1 integer km1 integer ks integer ksj integer kskp1 integer kstep logical nodet logical noert logical noinv real t complex temp complex work(n) noinv = mod ( job, 10 ) == 0 nodet = mod ( job, 100 ) / 10 == 0 noert = mod ( job, 1000 ) / 100 == 0 if ( .not. nodet .or. .not. noert ) then if ( .not. noert ) then inert(1:3) = 0 end if if ( .not. nodet ) then det(1) = 1.0E+00 det(2) = 0.0E+00 end if t = 0.0E+00 ik = 0 do k = 1, n kk = ik + k d = real ( ap(kk) ) ! ! Check if 1 by 1 ! if ( ipvt(k) <= 0 ) then ! ! 2 by 2 block ! Use DET (D S; S C) = ( D / T * C - T ) * T, T = abs ( S ) ! to avoid underflow/overflow troubles. ! Take two passes through scaling. Use T for flag. ! if ( t == 0.0E+00 ) then ikp1 = ik + k kkp1 = ikp1 + k t = abs ( ap(kkp1) ) d = ( d / t ) * real ( ap(kkp1+1) ) - t else d = t t = 0.0E+00 end if end if if ( .not. noert ) then if ( 0.0E+00 < d ) then inert(1) = inert(1) + 1 else if ( d < 0.0E+00 ) then inert(2) = inert(2) + 1 else if ( d == 0.0E+00 ) then inert(3) = inert(3) + 1 end if end if if ( .not. nodet ) then det(1) = det(1) * d if ( det(1) /= 0.0E+00 ) then do while ( abs ( det(1) ) < 1.0E+00 ) det(1) = det(1) * 10.0E+00 det(2) = det(2) - 1.0E+00 end do do while ( 10.0E+00 <= abs ( det(1) ) ) det(1) = det(1) / 10.0E+00 det(2) = det(2) + 1.0E+00 end do end if end if ik = ik + k end do end if ! ! Compute inverse(A). ! if ( .not. noinv ) then k = 1 ik = 0 do while ( k <= n ) km1 = k - 1 kk = ik + k ikp1 = ik + k kkp1 = ikp1 + k ! ! 1 by 1 ! if ( 0 <= ipvt(k) ) then ap(kk) = cmplx ( 1.0E+00 / real ( ap(kk) ), 0.0E+00 ) if ( 1 <= km1 ) then work(1:km1) = ap(ik+1:ik+km1) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = cdotc ( j, ap(ij+1:ij+j), 1, work, 1 ) ap(ik+1:ik+j-1) = ap(ik+1:ik+j-1) + work(j) * ap(ij+1:ij+j-1) ij = ij + j end do ap(kk) = ap(kk) + cmplx ( real ( & cdotc ( km1, work, 1, ap(ik+1:ik+km1), 1 ) ), 0.0E+00 ) end if kstep = 1 ! ! 2 by 2 ! else t = abs ( ap(kkp1) ) ak = real ( ap(kk) ) / t akp1 = real ( ap(kkp1+1) ) / t akkp1 = ap(kkp1) / t d = t * ( ak * akp1 - 1.0E+00 ) ap(kk) = cmplx ( akp1 / d, 0.0E+00 ) ap(kkp1+1) = cmplx ( ak / d, 0.0E+00 ) ap(kkp1) = -akkp1 / d if ( 1 <= km1 ) then work(1:km1) = ap(ikp1+1:ikp1+km1) ij = 0 do j = 1, km1 jkp1 = ikp1 + j ap(jkp1) = cdotc ( j, ap(ij+1:ij+j), 1, work, 1 ) ap(ikp1+1:ikp1+j-1) = ap(ikp1+1:ikp1+j-1) & + work(j) * ap(ij+1:ij+j-1) ij = ij + j end do ap(kkp1+1) = ap(kkp1+1) + cmplx ( real ( cdotc ( km1, work, 1, & ap(ikp1+1), 1 ) ), 0.0E+00 ) ap(kkp1) = ap(kkp1) & + cdotc ( km1, ap(ik+1:ik+km1), 1, ap(ikp1+1:ikp1+km1), 1 ) work(1:km1) = ap(ik+1:ik+km1) ij = 0 do j = 1, km1 jk = ik + j ap(jk) = cdotc ( j, ap(ij+1:ij+j), 1, work, 1 ) ap(ik+1:ik+j-1) = ap(ik+1:ik+j-1) + work(j) * ap(ij+1:ij+j-1) ij = ij + j end do ap(kk) = ap(kk) & + cmplx ( real ( cdotc ( km1, work, 1, ap(ik+1:ik+km1), 1 ) ), & 0.0E+00 ) end if kstep = 2 end if ! ! Swap ! ks = abs ( ipvt(k) ) if ( ks /= k ) then iks = ( ks * ( ks - 1 ) ) / 2 call cswap ( ks, ap(iks+1:iks+ks), 1, ap(ik+1:ik+ks), 1 ) ksj = ik + ks do jb = ks, k j = k + ks - jb jk = ik + j temp = conjg ( ap(jk) ) ap(jk) = conjg ( ap(ksj) ) ap(ksj) = temp ksj = ksj - ( j - 1 ) end do if ( kstep /= 1 ) then kskp1 = ikp1 + ks temp = ap(kskp1) ap(kskp1) = ap(kkp1) ap(kkp1) = temp end if end if ik = ik + k if ( kstep == 2 ) then ik = ik + k + 1 end if k = k + kstep end do end if return end subroutine chpfa ( ap, n, ipvt, info ) !*****************************************************************************80 ! !! CHPFA factors a complex hermitian packed matrix. ! ! Discussion: ! ! To solve A*X = B, follow CHPFA by CHPSL. ! ! To compute inverse(A)*C, follow CHPFA by CHPSL. ! ! To compute determinant(A), follow CHPFA by CHPDI. ! ! To compute inertia(A), follow CHPFA by CHPDI. ! ! To compute inverse(A), follow CHPFA by CHPDI. ! ! Packed storage: ! ! The following program segment will pack the upper ! triangle of a hermitian matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex AP(N*(N+1)/2); on input, the packed form ! of a hermitian matrix. The columns of the upper triangle are ! stored sequentially in a one-dimensional array. On output, a ! block diagonal matrix and the multipliers which were used to ! obtain it stored in packed form. The factorization can be ! written A = U*D*hermitian(U) where U is a product of permutation ! and unit upper triangular matrices , hermitian(U) is the ! conjugate transpose of U, and D is block diagonal with 1 by 1 ! and 2 by 2 blocks. ! ! Input, integer N, the order of the matrix. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO. ! 0, normal value. ! K, if the K-th pivot block is singular. This is not an error condition ! for this subroutine, but it does indicate that CHPSL or CHPDI may divide ! by zero if called. ! implicit none integer n real absakk complex ak complex akm1 real alpha complex ap((n*(n+1))/2) complex bk complex bkm1 real cabs1 real colmax complex denom integer icamax integer ij integer ijj integer ik integer ikm1 integer im integer imax integer imim integer imj integer imk integer info integer ipvt(n) integer j integer jj integer jk integer jkm1 integer jmax integer jmim integer k integer kk integer km1 integer km1k integer km1km1 integer km2 integer kstep complex mulk complex mulkm1 real rowmax logical swap complex t ! ! Initialize. ! ! ALPHA is used in choosing pivot block size. ! alpha = ( 1.0E+00 + sqrt ( 17.0E+00 ) ) / 8.0E+00 info = 0 ! ! Main loop on K, which goes from N to 1. ! k = n ik = ( n * ( n - 1 ) ) / 2 do ! ! Leave the loop if K = 0 or K = 1. ! if ( k == 0 ) then exit end if if ( k == 1 ) then ipvt(1) = 1 if ( cabs1 ( ap(1) ) == 0.0E+00 ) then info = 1 end if exit end if ! ! This section of code determines the kind of ! elimination to be performed. When it is completed, ! KSTEP will be set to the size of the pivot block, and ! SWAP will be set to TRUE if an interchange is ! required. ! km1 = k - 1 kk = ik + k absakk = cabs1 ( ap(kk) ) ! ! Determine the largest off-diagonal element in column K. ! imax = icamax ( k-1, ap(ik+1:ik+k-1), 1 ) imk = ik + imax colmax = cabs1 ( ap(imk) ) if ( alpha * colmax <= absakk ) then kstep = 1 swap = .false. ! ! Determine the largest off-diagonal element in row IMAX. ! else rowmax = 0.0E+00 im = ( imax * ( imax - 1 ) ) / 2 imj = im + 2 * imax do j = imax + 1, k rowmax = max ( rowmax, cabs1 ( ap(imj) ) ) imj = imj + j end do if ( imax /= 1 ) then jmax = icamax ( imax-1, ap(im+1:im+imax-1), 1 ) jmim = jmax + im rowmax = max ( rowmax, cabs1 ( ap(jmim) ) ) end if imim = imax + im if ( alpha * rowmax <= cabs1 ( ap(imim) ) ) then kstep = 1 swap = .true. else if ( alpha * colmax * ( colmax / rowmax ) <= absakk ) then kstep = 1 swap = .false. else kstep = 2 swap = ( imax /= km1 ) end if end if ! ! Column K is zero. Set INFO and iterate the loop. ! if ( max ( absakk, colmax ) == 0.0E+00 ) then ipvt(k) = k info = k ik = ik - ( k - 1 ) if ( kstep == 2 ) then ik = ik - ( k - 2 ) end if k = k - kstep cycle end if if ( kstep /= 2 ) then ! ! 1 x 1 pivot block. ! if ( swap ) then call cswap ( imax, ap(im+1:im+imax), 1, ap(ik+1:ik+imax), 1 ) imj = ik + imax do jj = imax, k j = k + imax - jj jk = ik + j t = conjg ( ap(jk) ) ap(jk) = conjg ( ap(imj) ) ap(imj) = t imj = imj - ( j - 1 ) end do end if ! ! Perform the elimination. ! ij = ik - ( k - 1 ) do jj = 1, km1 j = k - jj jk = ik + j mulk = -ap(jk) / ap(kk) t = conjg ( mulk ) ap(ij+1:ij+j) = ap(ij+1:ij+j) + t * ap(ik+1:ik+j) ijj = ij + j ap(ijj) = cmplx ( real ( ap(ijj) ), 0.0E+00 ) ap(jk) = mulk ij = ij - ( j - 1 ) end do ! ! Set the pivot array. ! if ( swap ) then ipvt(k) = imax else ipvt(k) = k end if ! ! 2 x 2 pivot block. ! else km1k = ik + k - 1 ikm1 = ik - ( k - 1 ) if ( swap ) then call cswap ( imax, ap(im+1:im+imax), 1, ap(ikm1+1:ikm1+imax), 1 ) imj = ikm1 + imax do jj = imax, km1 j = km1 + imax - jj jkm1 = ikm1 + j t = conjg ( ap(jkm1) ) ap(jkm1) = conjg ( ap(imj) ) ap(imj) = t imj = imj - ( j - 1 ) end do t = ap(km1k) ap(km1k) = ap(imk) ap(imk) = t end if ! ! Perform the elimination. ! km2 = k - 2 if ( km2 /= 0 ) then ak = ap(kk) / ap(km1k) km1km1 = ikm1 + k - 1 akm1 = ap(km1km1) / conjg ( ap(km1k) ) denom = 1.0E+00 - ak * akm1 ij = ik - ( k - 1 ) - ( k - 2 ) do jj = 1, km2 j = km1 - jj jk = ik + j bk = ap(jk) / ap(km1k) jkm1 = ikm1 + j bkm1 = ap(jkm1) / conjg ( ap(km1k) ) mulk = ( akm1 * bk - bkm1 ) / denom mulkm1 = ( ak * bkm1 - bk ) / denom t = conjg ( mulk ) ap(ij+1:ij+j) = ap(ij+1:ij+j) + t * ap(ik+1:ik+j) t = conjg ( mulkm1 ) ap(ij+1:ij+j) = ap(ij+1:ij+j) + t * ap(ikm1+1:ikm1+j) ap(jk) = mulk ap(jkm1) = mulkm1 ijj = ij + j ap(ijj) = cmplx ( real ( ap(ijj) ), 0.0E+00 ) ij = ij - ( j - 1 ) end do end if ! ! Set the pivot array. ! if ( swap ) then ipvt(k) = -imax else ipvt(k) = 1 - k end if ipvt(k-1) = ipvt(k) end if ik = ik - ( k - 1 ) if ( kstep == 2 ) then ik = ik - ( k - 2 ) end if k = k - kstep end do return end subroutine chpsl ( ap, n, ipvt, b ) !*****************************************************************************80 ! !! CHPSL solves a complex hermitian system factored by CHPFA. ! ! Discussion: ! ! A division by zero may occur if CHPCO set RCOND to 0.0 ! or CHPFA set INFO nonzero. ! ! To compute ! ! inverse ( A ) * C ! ! where C is a matrix with P columns ! ! call chpfa(ap,n,ipvt,info) ! ! if ( info == 0 ) then ! do j = 1, p ! call chpsl(ap,n,ipvt,c(1,j)) ! end do ! end if ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex AP(N*(N+1)/2), the output from CHPFA. ! ! Input, integer N, the order of the matrix. ! ! Input, integer IPVT(N), the pivot vector from CHPFA. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! implicit none integer n complex ak complex akm1 complex ap((n*(n+1))/2) complex b(n) complex bk complex bkm1 complex cdotc complex denom integer ik integer ikm1 integer ikp1 integer ipvt(n) integer k integer kk integer km1k integer km1km1 integer kp complex t ! ! Loop backward applying the transformations and inverse ( D ) to B. ! k = n ik = ( n * ( n - 1 ) ) / 2 do while ( 0 < k ) kk = ik + k ! ! 1 x 1 pivot block. ! if ( 0 <= ipvt(k) ) then if ( k /= 1 ) then kp = ipvt(k) if ( kp /= k ) then t = b(k) b(k) = b(kp) b(kp) = t end if b(1:k-1) = b(1:k-1) + b(k) * ap(ik+1:ik+k-1) end if ! ! Apply D inverse. ! b(k) = b(k) / ap(kk) k = k - 1 ik = ik - k else ! ! 2 x 2 pivot block. ! ikm1 = ik - ( k - 1 ) if ( k /= 2 ) then kp = abs ( ipvt(k) ) if ( kp /= k - 1 ) then t = b(k-1) b(k-1) = b(kp) b(kp) = t end if b(1:k-2) = b(1:k-2) + b(k) * ap(ik+1:ik+k-2) b(1:k-2) = b(1:k-2) + b(k-1) * ap(ikm1+1:ikm1+k-2) end if ! ! Apply D inverse. ! km1k = ik + k - 1 kk = ik + k ak = ap(kk) / conjg ( ap(km1k) ) km1km1 = ikm1 + k - 1 akm1 = ap(km1km1) / ap(km1k) bk = b(k) / conjg ( ap(km1k) ) bkm1 = b(k-1) / ap(km1k) denom = ak * akm1 - 1.0E+00 b(k) = ( akm1 * bk - bkm1 ) / denom b(k-1) = ( ak * bkm1 - bk ) / denom k = k - 2 ik = ik - ( k + 1 ) - k end if end do ! ! Loop forward applying the transformations. ! k = 1 ik = 0 do while ( k <= n ) ! ! 1 x 1 pivot block. ! if ( 0 <= ipvt(k) ) then if ( k /= 1 ) then b(k) = b(k) + cdotc ( k-1, ap(ik+1:ik+k-1), 1, b(1:k-1), 1 ) kp = ipvt(k) if ( kp /= k ) then t = b(k) b(k) = b(kp) b(kp) = t end if end if ik = ik + k k = k + 1 ! ! 2 x 2 pivot block. ! else if ( k /= 1 ) then b(k) = b(k) + cdotc ( k-1, ap(ik+1:ik+k-1), 1, b(1:k-1), 1 ) ikp1 = ik + k b(k+1) = b(k+1) + cdotc ( k-1, ap(ikp1+1:ikp1+k-1), 1, b(1:k-1), 1 ) kp = abs ( ipvt(k) ) if ( kp /= k ) then t = b(k) b(k) = b(kp) b(kp) = t end if end if ik = ik + k + k + 1 k = k + 2 end if end do return end subroutine cpbco ( abd, lda, n, m, rcond, z, info ) !*****************************************************************************80 ! !! CPBCO factors a complex hermitian positive definite band matrix. ! ! Discussion: ! ! The routine also estimates the condition number of the matrix. ! ! If RCOND is not needed, CPBFA is slightly faster. ! ! To solve A*X = B, follow CPBCO by CPBSL. ! ! To compute inverse(A)*C, follow CPBCO by CPBSL. ! ! To compute determinant(A), follow CPBCO by CPBDI. ! ! Band storage: ! ! If A is a hermitian positive definite band matrix, ! the following program segment will set up the input. ! ! m = (band width above diagonal) ! do j = 1, n ! i1 = max ( 1, j-m ) ! do i = i1, j ! k = i-j+m+1 ! abd(k,j) = a(i,j) ! end do ! end do ! ! This uses M+1 rows of A, except for the M by M ! upper left triangle, which is ignored. ! ! Example: ! ! If the original matrix is ! ! 11 12 13 0 0 0 ! 12 22 23 24 0 0 ! 13 23 33 34 35 0 ! 0 24 34 44 45 46 ! 0 0 35 45 55 56 ! 0 0 0 46 56 66 ! ! then N = 6, M = 2 and ABD should contain ! ! * * 13 24 35 46 ! * 12 23 34 45 56 ! 11 22 33 44 55 66 ! ! Modified: ! ! 13 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex ABD(LDA,N); on input, the matrix to be factored. ! The columns of the upper triangle are stored in the columns of ABD, ! and the diagonals of the upper triangle are stored in the rows of ABD. ! On output, an upper triangular matrix R, stored in band form, so that ! A = hermitian(R) * R. If INFO /= 0, the factorization is not complete. ! ! Input, integer LDA, the leading dimension of ABD. ! LDA must be at least M+1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! 0 <= M < N. ! ! Output, real RCOND, an estimate of the reciprocal condition of ! the matrix. For the system A*X = B, relative perturbations in A and B ! of size EPSILON may cause relative perturbations in X of size ! (EPSILON/RCOND). If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is singular to working precision, then Z is ! an approximate null vector in the sense that ! norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ). ! If INFO /= 0, Z is unchanged. ! ! Output, integer INFO. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is not ! positive definite. ! implicit none integer lda integer n complex abd(lda,n) real anorm real cabs1 complex cdotc complex csign1 complex ek integer i integer info integer j integer j2 integer k integer kp1 integer l integer la integer lb integer lm integer m integer mu real rcond real s real scasum real sm complex t complex wk complex wkm real ynorm complex z(n) rcond = 0.0E+00 ! ! Find the norm of A. ! do j = 1, n l = min ( j, m + 1 ) mu = max ( m + 2 - j, 1 ) z(j) = cmplx ( scasum ( l, abd(mu:mu+l-1,j), 1 ), 0.0E+00 ) k = j - l do i = mu, m k = k + 1 z(k) = cmplx ( real ( z(k) ) + cabs1 ( abd(i,j) ), 0.0E+00 ) end do end do anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, real ( z(j) ) ) end do ! ! Factor. ! call cpbfa ( abd, lda, n, m, info ) if ( info /= 0 ) then return end if ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where hermitian(R)*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve hermitian(R)*W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) do k = 1, n if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csign1 ( ek, -z(k) ) end if if ( real ( abd(m+1,k) ) < cabs1 ( ek - z(k) ) ) then s = real ( abd(m+1,k) ) / cabs1 ( ek - z(k) ) z(1:n) = z(1:n) * s ek = cmplx ( s, 0.0E+00 ) * ek end if wk = ek - z(k) wkm = - ek - z(k) s = cabs1 ( wk ) sm = cabs1 ( wkm ) wk = wk / abd(m+1,k) wkm = wkm / abd(m+1,k) j2 = min ( k + m, n ) i = m + 1 if ( k+1 <= j2 ) then do j = k+1, j2 i = i - 1 sm = sm + cabs1 ( z(j) + wkm * conjg ( abd(i,j) ) ) z(j) = z(j) + wk * conjg ( abd(i,j) ) s = s + cabs1 ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm i = m + 1 do j = k+1, j2 i = i - 1 z(j) = z(j) + t * conjg ( abd(i,j) ) end do end if end if z(k) = wk end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ! ! Solve R * Y = W. ! do k = n, 1, -1 if ( real ( abd(m+1,k) ) < cabs1 ( z(k) ) ) then s = real ( abd(m+1,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s end if z(k) = z(k) / abd(m+1,k) lm = min ( k - 1, m ) la = m + 1 - lm lb = k - lm t = -z(k) z(lb:lb+lm-1) = z(lb:lb+lm-1) + t * abd(la:la+lm-1,k) end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = 1.0E+00 ! ! Solve hermitian(R)*V = Y. ! do k = 1, n lm = min ( k - 1, m ) la = m + 1 - lm lb = k - lm z(k) = z(k) - cdotc ( lm, abd(la:la+lm-1,k), 1, z(lb:lb+lm-1), 1 ) if ( real ( abd(m+1,k) ) < cabs1 ( z(k) ) ) then s = real ( abd(m+1,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if z(k) = z(k) / abd(m+1,k) end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm ! ! Solve R * Z = W. ! do k = n, 1, -1 if ( real ( abd(m+1,k) ) < cabs1 ( z(k) ) ) then s = real ( abd(m+1,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if z(k) = z(k) / abd(m+1,k) lm = min ( k - 1, m ) la = m + 1 - lm lb = k - lm t = -z(k) z(lb:lb+lm-1) = z(lb:lb+lm-1) + t * abd(la:la+lm-1,k) end do ! ! Make ZNORM = 1. ! s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine cpbdi ( abd, lda, n, m, det ) !*****************************************************************************80 ! !! CPBDI gets the determinant of a hermitian positive definite band matrix. ! ! Discussion: ! ! CPBDI uses the factors computed by CPBCO or CPBFA. ! ! If the inverse is needed, use CPBSL N times. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex ABD(LDA,N), the output from CPBCO or CPBFA. ! ! Input, integer LDA, the leading dimension of the array ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! ! Output, real DET(2), the determinant of the original matrix in the ! form determinant = DET(1) * 10.0**DET(2) with 1.0 <= DET(1) < 10.0 ! or DET(1) == 0.0. ! implicit none integer lda integer n complex abd(lda,n) real det(2) integer i integer m det(1) = 1.0E+00 det(2) = 0.0E+00 do i = 1, n det(1) = det(1) * real ( abd(m+1,i) )**2 if ( det(1) == 0.0E+00 ) then exit end if do while ( det(1) < 1.0E+00 ) det(1) = det(1) * 10.0E+00 det(2) = det(2) - 1.0E+00 end do do while ( 10.0E+00 <= det(1) ) det(1) = det(1) / 10.0E+00 det(2) = det(2) + 1.0E+00 end do end do return end subroutine cpbfa ( abd, lda, n, m, info ) !*****************************************************************************80 ! !! CPBFA factors a complex hermitian positive definite band matrix. ! ! Discussion: ! ! CPBFA is usually called by CPBCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Band storage: ! ! If A is a hermitian positive definite band matrix, ! the following program segment will set up the input. ! ! m = (band width above diagonal) ! do j = 1, n ! i1 = max ( 1, j-m ) ! do i = i1, j ! k = i-j+m+1 ! abd(k,j) = a(i,j) ! end do ! end do ! ! Modified: ! ! 23 March 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex ABD(LDA,N); on input, the matrix to be factored. ! The columns of the upper triangle are stored in the columns of ABD ! and the diagonals of the upper triangle are stored in the rows of ABD. ! On output, an upper triangular matrix R, stored in band form, so that ! A = hermitian(R)*R. ! ! Input, integer LDA, the leading dimension of ABD. ! LDA must be at least M+1. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! 0 <= M < N. ! ! Output, integer INFO. ! 0, for normal return. ! K, if the leading minor of order K is not positive definite. ! implicit none integer lda integer n complex abd(lda,n) complex cdotc integer m integer ik integer info integer j integer jk integer k integer mu real s complex t info = 0 do j = 1, n s = 0.0E+00 ik = m + 1 jk = max ( j - m, 1 ) mu = max ( m + 2 - j, 1 ) do k = mu, m t = abd(k,j) - cdotc ( k-mu, abd(ik:ik+k-mu-1,jk), 1, abd(mu:mu+k-mu-1,j), 1 ) t = t / abd(m+1,jk) abd(k,j) = t s = s + real ( t * conjg ( t ) ) ik = ik - 1 jk = jk + 1 end do s = real ( abd(m+1,j) ) - s if ( s <= 0.0E+00 .or. aimag ( abd(m+1,j) ) /= 0.0E+00 ) then info = j exit end if abd(m+1,j) = cmplx ( sqrt ( s ), 0.0E+00 ) end do return end subroutine cpbsl ( abd, lda, n, m, b ) !*****************************************************************************80 ! !! CPBSL solves a complex hermitian positive definite band system. ! ! Discussion: ! ! The system matrix must have been factored by CPBCO or CPBFA. ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal. Technically this indicates ! singularity but it is usually caused by improper subroutine ! arguments. It will not occur if the subroutines are called ! correctly and INFO == 0. ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call cpbco(abd,lda,n,rcond,z,info) ! ! if (rcond is too small .or. info /= 0) then ! error ! end if ! ! do j = 1, p ! call cpbsl(abd,lda,n,c(1,j)) ! end do ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex ABD(LDA,N), the output from CPBCO or CPBFA. ! ! Input, integer LDA, the leading dimension of ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer M, the number of diagonals above the main diagonal. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! implicit none integer lda integer n complex abd(lda,n) complex b(n) complex cdotc integer k integer la integer lb integer lm integer m complex t ! ! Solve hermitian(R) * Y = B. ! do k = 1, n lm = min ( k - 1, m ) la = m + 1 - lm lb = k - lm t = cdotc ( lm, abd(la:la+lm-1,k), 1, b(lb:lb+lm-1), 1 ) b(k) = ( b(k) - t ) / abd(m+1,k) end do ! ! Solve R * X = Y. ! do k = n, 1, -1 lm = min ( k - 1, m ) la = m + 1 - lm lb = k - lm b(k) = b(k) / abd(m+1,k) t = -b(k) b(lb:lb+lm-1) = b(lb:lb+lm-1) + t * abd(la:la+lm-1,k) end do return end subroutine cpoco ( a, lda, n, rcond, z, info ) !*****************************************************************************80 ! !! CPOCO factors a complex hermitian positive definite matrix. ! ! Discussion: ! ! The routine also estimates the condition of the matrix. ! ! If RCOND is not needed, CPOFA is slightly faster. ! ! To solve A*X = B, follow CPOCO by CPOSL. ! ! To compute inverse(A)*C, follow CPOCO by CPOSL. ! ! To compute determinant(A), follow CPOCO by CPODI. ! ! To compute inverse(A), follow CPOCO by CPODI. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the hermitian matrix to be ! factored. On output, an upper triangular matrix R so that ! A = hermitian(R)*R ! where hermitian(R) is the conjugate transpose. The strict lower ! triangle is unaltered. If INFO /= 0, the factorization is not complete. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, real RCOND, an estimate of the reciprocal condition of ! the matrix. For the system A*X = B, relative perturbations in A and B ! of size EPSILON may cause relative perturbations in X of size ! (EPSILON/RCOND). If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is close to a singular matrix, then Z is an ! approximate null vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! ! Output, integer INFO. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is not ! positive definite. ! implicit none integer lda integer n complex a(lda,n) real anorm real cabs1 complex cdotc complex csign1 complex ek integer i integer info integer j integer k integer kp1 real rcond real s real scasum real sm complex t complex wk complex wkm real ynorm complex z(n) rcond = 0.0E+00 ! ! Find norm of A using only upper half. ! do j = 1, n z(j) = cmplx ( scasum ( j, a(1:j,j), 1 ), 0.0E+00 ) do i = 1, j - 1 z(i) = cmplx ( real ( z(i) ) + cabs1 ( a(i,j) ), 0.0E+00 ) end do end do anorm = 0.0E+00 do j = 1, n anorm = max ( anorm, real ( z(j) ) ) end do ! ! Factor. ! call cpofa ( a, lda, n, info ) if ( info /= 0 ) then return end if ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where hermitian(R)*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve hermitian(R)*W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) do k = 1, n if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csign1 ( ek, -z(k) ) end if if ( real ( a(k,k) ) < cabs1 ( ek - z(k) ) ) then s = real ( a(k,k) ) / cabs1 ( ek - z(k) ) z(1:n) = z(1:n) * s ek = cmplx ( s, 0.0E+00 ) * ek end if wk = ek - z(k) wkm = -ek - z(k) s = cabs1 ( wk ) sm = cabs1 ( wkm ) wk = wk / a(k,k) wkm = wkm / a(k,k) kp1 = k + 1 if ( kp1 <= n ) then do j = kp1, n sm = sm + cabs1 ( z(j) + wkm * conjg ( a(k,j) ) ) z(j) = z(j) + wk * conjg ( a(k,j) ) s = s + cabs1 ( z(j) ) end do if ( s < sm ) then t = wkm - wk wk = wkm z(kp1:n) = z(kp1:n) + t * conjg ( a(k,kp1:n) ) end if end if z(k) = wk end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ! ! Solve R * Y = W. ! do k = n, 1, -1 if ( real ( a(k,k) ) < cabs1 ( z(k) ) ) then s = real ( a(k,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s end if z(k) = z(k) / a(k,k) t = -z(k) z(1:k-1) = z(1:k-1) + t * a(1:k-1,k) end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = 1.0E+00 ! ! Solve hermitian(R) * V = Y. ! do k = 1, n z(k) = z(k) - cdotc ( k-1, a(1:k-1,k), 1, z(1:k-1), 1 ) if ( real ( a(k,k) ) < cabs1 ( z(k) ) ) then s = real ( a(k,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if z(k) = z(k) / a(k,k) end do s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm ! ! Solve R * Z = V. ! do k = n, 1, -1 if ( real ( a(k,k) ) < cabs1 ( z(k) ) ) then s = real ( a(k,k) ) / cabs1 ( z(k) ) z(1:n) = z(1:n) * s ynorm = s * ynorm end if z(k) = z(k) / a(k,k) t = -z(k) z(1:k-1) = z(1:k-1) + t * a(1:k-1,k) end do ! ! Make ZNORM = 1. ! s = 1.0E+00 / scasum ( n, z, 1 ) z(1:n) = z(1:n) * s ynorm = s * ynorm if ( anorm /= 0.0E+00 ) then rcond = ynorm / anorm else rcond = 0.0E+00 end if return end subroutine cpodi ( a, lda, n, det, job ) !*****************************************************************************80 ! !! CPODI: determinant, inverse of a complex hermitian positive definite matrix. ! ! Discussion: ! ! The matrix is assumed to have been factored by CPOCO, CPOFA or CQRDC. ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal and the inverse is requested. ! It will not occur if the subroutines are called correctly ! and if CPOCO or CPOFA has set INFO == 0. ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); on input, the output A from CPOCO or ! CPOFA, or the output X from CQRDC. On output, if CPOCO or CPOFA was ! used to factor A, then CPODI produces the upper half of inverse(A). ! If CQRDC was used to decompose X, then CPODI produces the upper half ! of inverse(hermitian(X)*X) where hermitian(X) is the conjugate transpose. ! Elements of A below the diagonal are unchanged. ! If the units digit of JOB is zero, A is unchanged. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, real DET(2), if requested, the determinant of A or of ! hermitian(X)*X. Determinant = DET(1) * 10.0**DET(2) with ! 1.0 <= abs ( DET(1) ) < 10.0 or DET(1) = 0.0. ! ! Input, integer JOB. ! 11, both determinant and inverse. ! 01, inverse only. ! 10, determinant only. ! implicit none integer lda integer n complex a(lda,n) real det(2) integer i integer j integer job integer k complex t ! ! Compute determinant ! if ( ( job / 10 ) /= 0 ) then det(1) = 1.0E+00 det(2) = 0.0E+00 do i = 1, n det(1) = det(1) * real ( a(i,i) )**2 if ( det(1) == 0.0E+00 ) then exit end if do while ( det(1) < 1.0E+00 ) det(1) = det(1) * 10.0E+00 det(2) = det(2) - 1.0E+00 end do do while ( 10.0E+00 <= det(1) ) det(1) = det(1) / 10.0E+00 det(2) = det(2) + 1.0E+00 end do end do end if ! ! Compute inverse(R). ! if ( mod ( job, 10 ) /= 0 ) then do k = 1, n a(k,k) = cmplx ( 1.0E+00, 0.0E+00 ) / a(k,k) t = -a(k,k) a(1:k-1,k) = a(1:k-1,k) * t do j = k+1, n t = a(k,j) a(k,j) = cmplx ( 0.0E+00, 0.0E+00 ) a(1:k,j) = a(1:k,j) + t * a(1:k,k) end do end do ! ! Form inverse(R) * hermitian(inverse(R)). ! do j = 1, n do k = 1, j-1 t = conjg ( a(k,j) ) a(1:k,k) = a(1:k,k) + t * a(1:k,j) end do t = conjg ( a(j,j) ) a(1:j,j) = a(1:j,j) * t end do end if return end subroutine cpofa ( a, lda, n, info ) !*****************************************************************************80 ! !! CPOFA factors a complex hermitian positive definite matrix. ! ! Discussion: ! ! CPOFA is usually called by CPOCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! (time for CPOCO) = (1 + 18/N) * (time for CPOFA). ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex A(LDA,N); On input, the hermitian matrix to be ! factored. On output, an upper triangular matrix R so that ! A = hermitian(R)*R ! where hermitian(R) is the conjugate transpose. The strict lower ! triangle is unaltered. If INFO /= 0, the factorization is not ! complete. Only the diagonal and upper triangle are used. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Output, integer INFO. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is ! not positive definite. ! implicit none integer lda integer n complex a(lda,n) complex cdotc integer info integer j integer k real s complex t info = 0 do j = 1, n s = 0.0E+00 do k = 1, j-1 t = a(k,j) - cdotc ( k-1, a(1:k-1,k), 1, a(1:k-1,j), 1 ) t = t / a(k,k) a(k,j) = t s = s + real ( t * conjg ( t ) ) end do s = real ( a(j,j) ) - s if ( s <= 0.0E+00 .or. aimag ( a(j,j) ) /= 0.0E+00 ) then info = j exit end if a(j,j) = cmplx ( sqrt ( s ), 0.0E+00 ) end do return end subroutine cposl ( a, lda, n, b ) !*****************************************************************************80 ! !! CPOSL solves a complex hermitian positive definite system. ! ! Discussion: ! ! CPOSL uses the factors computed by CPOCO or CPOFA. ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal. Technically this indicates ! singularity but it is usually caused by improper subroutine ! arguments. It will not occur if the subroutines are called ! correctly and INFO == 0. ! ! To compute inverse(A) * C where C is a matrix with p columns ! ! call cpoco(a,lda,n,rcond,z,info) ! ! if (rcond is too small or info is not zero ) then ! error ! end if ! ! do j = 1, p ! call cposl(a,lda,n,c(1,j)) ! end do ! ! Modified: ! ! 07 May 2006 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input, complex A(LDA,N), the output from CPOCO or CPOFA. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer N, the order of the matrix. ! ! Input/output, complex B(N). On input, the right hand side. ! On output, the solution. ! implicit none integer lda integer n complex a(lda,n) complex b(n) complex cdotc integer k complex t ! ! Solve hermitian(R) * Y = B. ! do k = 1, n t = cdotc ( k-1, a(1:k-1,k), 1, b(1:k-1), 1 ) b(k) = ( b(k) - t ) / a(k,k) end do ! ! Solve R * X = Y. ! do k = n, 1, -1 b(k) = b(k) / a(k,k) b(1:k-1) = b(1:k-1) - a(1:k-1,k) * b(k) end do return end subroutine cppco ( ap, n, rcond, z, info ) !*****************************************************************************80 ! !! CPPCO factors a complex hermitian positive definite matrix. ! ! Discussion: ! ! The routine also estimates the condition of the matrix. ! ! The matrix is stored in packed form. ! ! If RCOND is not needed, CPPFA is slightly faster. ! ! To solve A*X = B, follow CPPCO by CPPSL. ! ! To compute inverse(A)*C, follow CPPCO by CPPSL. ! ! To compute determinant(A), follow CPPCO by CPPDI. ! ! To compute inverse(A), follow CPPCO by CPPDI. ! ! Packed storage: ! ! The following program segment will pack the upper ! triangle of a hermitian matrix. ! ! k = 0 ! do j = 1, n ! do i = 1, j ! k = k + 1 ! ap(k) = a(i,j) ! end do ! end do ! ! Modified: ! ! 13 May 2007 ! ! Author: ! ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, complex AP(N*(N+1)/2); on input, the packed form of a ! hermitian matrix. The columns of the upper triangle are stored ! sequentially in a one-dimensional array. On output, an upper ! triangular matrix R, stored in packed form, so that A = hermitian(R) * R. ! If INFO /= 0 , the factorization is not complete. ! ! Input, integer N, the order of the matrix. ! ! Output, real RCOND, an estimate of the reciprocal condition of ! the matrix. For the system A*X = B, relative perturbations in A and B ! of size EPSILON may cause relative perturbations in X of size ! (EPSILON/RCOND). If RCOND is so small that the logical expression ! 1.0 + RCOND == 1.0 ! is true, then A may be singular to working precision. In particular, ! RCOND is zero if exact singularity is detected or the estimate underflows. ! ! Workspace, complex Z(N), a work vector whose contents are usually ! unimportant. If A is close to a singular matrix, then Z is an ! approximate null vector in the sense that ! norm(A*Z) = RCOND * norm(A) * norm(Z). ! ! Output, integer INFO. ! 0, for normal return. ! K, signals an error condition. The leading minor of order K is not ! positive definite. ! implicit none integer n real anorm complex ap((n*(n+1))/2) real cabs1 complex cdotc complex csign1 complex ek integer i integer ij integer info integer j integer j1 integer k integer kj integer kk real rcond real s real scasum real sm complex t complex wk complex wkm real ynorm complex z(n) rcond = 0.0E+00 ! ! Find norm of A. ! j1 = 1 do j = 1, n z(j) = cmplx ( scasum ( j, ap(j1:j1+j-1), 1 ), 0.0E+00 ) ij = j1 j1 = j1 + j do i = 1, j-1 z(i) = cmplx ( real ( z(i) ) + cabs1 ( ap(ij) ), 0.0E+00 ) ij = ij + 1 end do end do anorm = maxval ( real ( z(1:n) ) ) ! ! Factor. ! call cppfa ( ap, n, info ) if ( info /= 0 ) then return end if ! ! RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))). ! ! Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E. ! ! The components of E are chosen to cause maximum local ! growth in the elements of W where hermitian(R)*W = E. ! ! The vectors are frequently rescaled to avoid overflow. ! ! Solve hermitian(R)*W = E. ! ek = cmplx ( 1.0E+00, 0.0E+00 ) z(1:n) = cmplx ( 0.0E+00, 0.0E+00 ) kk = 0 do k = 1, n kk = kk + k if ( cabs1 ( z(k) ) /= 0.0E+00 ) then ek = csig